options ps=54 ls=80 nodate nonumber;  

proc iml;
XY1 = {43.7 19 177,  43.7 43 279, 43.7 56 346, 54.6 13 160,  54.6 19 193,  54.6 43 280, 54.6 56 335,                                    
55.7 13 169, 55.7 26 212, 55.7 34.5 244, 55.7 43 285, 58.8 13 181, 58.8 43 298, 60.5 19 212, 60.5 43 317,                               
60.5 56 347, 61.9 13 186, 61.9 19 216, 61.9 34.5 265, 61.9 43 306, 61.9 56 348, 66.7 13 209,                                            
66.7 43 324,  66.7 56 352};

/*  XY1 is a matrix containing n=24 rows and 3 columns (X1,X2,Y) we must form the X and Y matrices from it
This is similar to when we input a raw data frame (dataset) and create a matrix  */   

X0 = j(nrow(XY1),1,1);   /* Create a Column of 1s for the intercept (same number of rows as XY1) */
X012 = X0||XY1(|,1|)||XY1(|,2|);   
Y = XY1(|,3|);            /* Create Y from the third column of XY1   */
X01 = X012(|,1:2|);
X02=X012(|,1|)||X012(|,3|);

P0 = X0 * inv(X0`*X0) * X0`;
P01 = X01 * inv(X01`*X01) * X01`;
P02 = X02 * inv(X02`*X02) * X02`;
P012 = X012 * inv(X012`*X012) * X012`;

betahat = inv(X012`*X012)*X012`*Y;
s2 = Y`*(i(nrow(Y))-P012)*Y/(nrow(Y)-nrow(betahat));
V_betahat = s2 * inv(X012`*X012);

s_betahat = sqrt(diag(V_betahat))*j(3,1,1);
t_betahat = betahat/s_betahat;
p_betahat = 2*(1-probt(abs(t_betahat),nrow(Y)-nrow(betahat)));

print betahat s_betahat t_betahat p_betahat;

R0 = Y` * P0 * Y;
R01 = Y` * P01 * Y;
R02 = Y` * P02 * Y;
R012 = Y` * P012 * Y;

R_1_0 = R01-R0;
R_2_01 = R012-R01;

R_2_0 = R02-R0;
R_1_02 = R012-R02;

print R_1_0 R_2_01;
print R_2_0 R_1_02;




run;                                                                                                                                    
stop;