/*
Matrix Activity #3

Obtain the LS regression coefficients, and
the associated covariance matrix, for the bodyfat 
percentage example  from class notes (p. 33) 
Use all three predictors.
Use the regression to estimate the mean bodyfat
among a population with x1=20, x2=43, x3=29.
Report a standard error.  Do this using PROC REG, 
then IML.
*/

options ls=75 nodate;
data bodyfats;
  input 
  x1    x2    x3    y;
  x0=1;
  cards;
  19.5  43.1  29.1  11.9
  24.7  49.8  28.2  22.8
  30.7  51.9  37.0  18.7
  29.8  54.3  31.1  20.1
  19.1  42.2  30.9  12.9
  25.6  53.9  23.7  21.7
  31.4  58.5  27.6  27.1
  27.9  52.1  30.6  25.4
  22.1  49.9  23.2  21.3
  25.5  53.5  24.8  19.3
  31.1  56.6  30.0  25.4
  30.4  56.7  28.3  27.2
  18.7  46.5  23.0  11.7
  19.7  44.2  28.6  17.8
  14.6  42.7  21.3  12.8
  29.5  54.4  30.1  23.9
  27.7  55.3  25.7  22.6
  30.2  58.6  24.6  25.4
  22.7  48.2  27.1  14.8
  25.2  51.0  27.5  21.1
  20     43    29    .
;

proc reg ;
   model y=x1-x3/p cli clm; /* full model */
   output out=two p=p lclm=lcb uclm=ucl;
   model y=x2/p cli clm; /* reduced model */
   output out=three p=p lclm=lcb uclm=ucl;
run;
proc g3d data=two;
   title "not easy to estimate dependence on both x1 and x2";
   scatter x1*x2=p;
run;  
data two; set two; if y=.; proc print;  title "full model w/  x1,x2,x3"; run;
data three; set three; if y=.; proc print ; title "reduced model w/  x2"; run;

data bodyfats;  set bodyfats;  if y>.; run;

PROC IML;

USE bodyfats;

READ ALL VAR{y} INTO y;
READ ALL VAR{x0 x1 x2 x3} INTO X;  /* (full model) */

n=nrow(X);
p=ncol(X);
XpX=t(X) * X;
XpXi=inv(XpX);
betahat = XpXi * (t(X)*y);
yhat = X * betahat;
residuals = y - yhat;
SSE = t(residuals)*residuals;
MSE = 1/(n-p) * SSE;
Sigmahat = MSE * XpXi;

a={1,20,43,29};
/*(population with x1=20, x2=43, x3=29)*/
muhat=t(a) * betahat;
semuhat=sqrt(t(a)*Sigmahat * a);  
print muhat;
print semuhat;

/* Reduced Model */
READ ALL VAR{x0 x2 } INTO Xr;

p=ncol(Xr);
XpX=t(Xr) * Xr;
XpXi=inv(XpX);
betahat = XpXi * (t(Xr)*y);
yhat = X * betahat;
residuals = y - yhat;
SSE = t(residuals)*residuals;
MSE = 1/(n-p) * SSE;
Sigmahat = MSE * XpXi;

a={1,43};
/*(population with x1=20, x2=43, x3=29)*/
muhatr=t(a) * betahat;
semuhatr=sqrt(t(a)*Sigmahat * a);  
print muhatr;
print semuhatr;
end;
quit;