%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter10;
PROC IMPORT OUT= WORK.Life 
            DATAFILE= "&location\chapter10.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH10TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Life Insurance Example- Refer to Table 10.1, page 387";run;
proc corr;
title "Basic correlations- Refer to page 386";run;
proc reg;
model insure=income risk;
title "Basic regression equation- Refer to Equation 10.3, page 386";run;
proc reg;
model insure=income risk;
plot residual.*income;
title "residual plot- Refer to Figure 10.3a, page 387";run;
proc reg;
model insure=risk;
output out=resids r=Y_risk;
proc reg data=resids;
model income=risk;
output out=resids r=X1_risk;
proc gplot;
plot y_risk*X1_risk;
title "Added variable plot for X1- Refer to Figure 10.3b, page 387";run;
PROC IMPORT OUT= WORK.bodyfat 
            DATAFILE= "&location\chapter10.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH07TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Bodyfat example from Chapter 7. Refer to Table 7.1, page 257";run;
proc reg;
model bodyfat=triceps thigh;
title "Regression function. Refer to Table 7.2, page 258 and page 388";run;
proc reg;
model bodyfat=triceps thigh;
plot residual.*(triceps thigh)/vref=0;
title "Residual plot against triceps Refer to Figure 10.4a and 10.4c, page 389";run;
proc reg;
model bodyfat triceps=thigh;
output out=r10_4 r=rbf_th rtr_th;
symbol i=rl v=dot line=2;
proc gplot;
plot rbf_th*rtr_th/vref=0;
title "Added Variable plot for triceps, Refer to Figure 10.4b, page 389";run;
run;
proc reg;
model bodyfat thigh=triceps;
output out=r10_4 r=rbf_tr rth_tr;
symbol i=rl v=dot line=2;
proc gplot;
plot rbf_tr*rth_tr/vref=0;
title "Added Variable plot for thigh, Refer to Figure 10.4d, page 389";run;
PROC IMPORT OUT= WORK.HatMat 
            DATAFILE= "&location\chapter10.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH10TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "Illustration of Hat Matrix data. Refer to Table 10.2, page 393";run;
proc reg;
model y=x1 x2;
output out=hatmat p=Ypred r=ei h=hii ;
data hatmat;set hatmat;
s2e=ei*ei;
proc print;
title "Data and Basic Results. Refer to Table 102, page 393";run;
proc reg;
model y=x1 x2/all;run;
proc iml;
use hatmat var {x1 x2};
read all into raw;
use hatmat var{y};
read all into y;
nrow=nrow(raw);
ones=j(nrow,1,1);
X=ones||raw;
parms=ncol(X);
print "Adjoined X matrix, Refer to Page 393" X;
hat=X*inv(X`*X)*X`;
print "Hat matrix. Refer to formula 10.10, page 392 and Table 10.2(b), page 383" hat;
yhat=hat*y;
print "Yhat via hat matrix. Refer to formula 10.11, page 392 and page 393" yhat;
identity=i(nrow);
e=(Identity-hat)*y;
print "error via Hat matrix and Y. Refer to formula 10.12, page 392" e;
mse=e`*e/(nrow-parms);
print "mse calculation- Refer to formula 10.13 and 6.32" mse;
vare_i=mse#(identity-hat);
print "variance of individual errors via MSE, see Equations 10.16a and 10.16b page 392 and 10.2c, page 393" vare_i;
proc reg data=bodyfat;
model bodyfat=thigh triceps/influence;
output out=bodyfatr r=residual h=hat student=student rstudent=rstudent dffits=dffits cookd=cookd p=pbodyfat;
title "residuals, studentized residuals, deleted studentized residuals dffits, Cook's D & dfbetas";
title2 "Refer to pp. 394-405";run;
data bodyfatr;set bodyfatr;
obs+1;
prob_t=1-probt(abs(rstudent),20-2-1);*two tailed test, therefore uninterested in sign of t;
critt=tinv(1-(.10/2/20),20-3-1);*two tailed test, family alpha of .10, 20 observations;
proc print data=bodyfatr;var residual hat rstudent prob_t critt;
title "residuals, hat, and studentized deleted residuals for bodyfat. Refer to Table 10.3, page 397";
title2 "see also critical t under Bonferroni correction with alpha =.10 page 397";run;
symbol1 v=dot i=none c=red;
symbol2 v=circle i=none c=green;
symbol3 v=square i=none c=blue;
symbol4 v=plus i=none c=black;
proc means data=hatmat;run;
proc gplot data=hatmat;
plot x2*x1=hii/vref=23.5 href=14;
title "Refer to Figure 10.6, page 398 (x marks the centroid)";run;
data bodyfat;set bodyfat;
n+1;
proc gplot data=bodyfat;
plot thigh*triceps=n;
title "observational plot for bodyfat. Refer to Table 10.7, page 399";run;
proc insight data=bodyfat;
title "Often it's easier to look at these annotations via Insight Refer to page 399";run;
proc print data=bodyfatr;var dffits;
title "dffits- see page 401";run;
proc reg data=bodyfat;
model bodyfat=thigh triceps/influence;
title "Dffits, and DFBETA's: Refer to table 10.4, page 402, 404, & 405";run;
proc print data=bodyfatr;var cookd;
title "Cook's d, Refer to Table 10.4, page 402";run;
axis1 order=-4.5 to 4.5 by .5;
axis2 order=10 to 30 by 5;
proc gplot data=bodyfatr;
bubble residual*pbodyfat=cookd/vaxis=axis1 haxis=axis2;
title "Bubble plot of Cook's D. Refer to Figure 10.8(a) page 404";run;
symbol1 i=join;
axis1 order=0 to .5 by .1;
axis2 order=0 to 25 by 5;
proc gplot data=bodyfatr;
plot cookd*obs/vaxis=axis1 haxis=axis2;
title "Index Influence Plot of Cook's D. Refer Figure 10.8(b) page 404";run;
proc reg data=bodyfat;
model bodyfat=thigh triceps/influence;
output out=bodyfat1 p=predy;
reweight obs.=3;
model bodyfat=thigh triceps/influence;
output out=bodyfat2 p=predy3;
data bodyfat1;merge bodyfat1 bodyfat2;
abspdiff=abs(predy3-predy)*100;
obs+1;
proc means;var abspdiff;
title "Average Absolute percent different in predicted Y scores, refer to page 406";
proc sort data=bodyfat1 out=bodyfat1 ;by descending abspdiff;
proc print noobs;var obs abspdiff;
title "Absolute percent difference in predicted Y scores by observation, refer to page 406";run;
proc standard mean=0 std=1 data=bodyfat out=bodyfatz;var bodyfat thigh triceps midarm;
proc reg data=bodyfatz outest=vifs outvif ridge=0;
model bodyfat=thigh triceps midarm;
title "Quick way to get VIF's, refer to page 408";run;
proc print data=vifs;
title "Variance Influation Factors Refer to Table 10.5, page 409";run;
data _null_;set vifs;
max=max(of thigh triceps midarm);
mean=mean(of thigh triceps midarm);
file print;
if _type_="PARMS" then put "thigh beta=" thigh 8.2 " triceps beta=" triceps 8.2 " midarm beta=" midarm 8.2;
if _type_="RIDGEVIF" then do;put "thigh VIF=" thigh 8.2 " triceps VIF=" triceps 8.2 " midarm VIF=" midarm 8.2;
 put "Maximum VIF is: " max 8.2;
 put "Average VIF is: " mean 8.2;
 end;
run;
PROC IMPORT OUT= WORK.surgery 
            DATAFILE= "&location\chapter10.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH09TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Surgery Unit Example- Data Summary- Refer to Page 410";run;
data surgery;set surgery;
xclot=clotting;
xprognos=prognos;
xenzyme=enzyme;
xalcheav=alcheav;
proc standard mean=0 data=surgery out=surgery;var xclot xprognos xenzyme xalcheav;
data surgery;set surgery;
clot_pro=xclot*xprognos;
clot_enz=xclot*xenzyme;
clot_alc=xclot*xalcheav;
pro_enz=xprognos*xenzyme;
pro_alc=xprognos*xalcheav;
enz_alc=xenzyme*xalcheav;
proc reg;var _numeric_;
model LnSurviv=clotting prognos enzyme alcheav;
output out=rsurg r=residual;
symbol1 i=rlcli95;
proc gplot;
plot residual*(clot_pro clot_enz clot_alc pro_enz pro_alc enz_alc);
title "added variable plots for interactions (not shown in text). Refer to Page 410";run;
proc reg;
model LnSurviv=clotting prognos enzyme alcheav clot_pro--enz_alc;
noints: test clot_pro--enz_alc=0;
title "test of no interactions in surgery data. Refer to page 410.";run;
proc reg;var _numeric_;
model LnSurviv=clotting prognos enzyme alcheav;
plot residual.*predicted. residual.*age residual.*nqq.;
title "Residual and Added Variable Plots for Surgical Unit example. Refer to Figures 10.9 a, b & d, page 411";run;
proc reg data=rsurg;
model age=clotting prognos enzyme alcheav;
output out=rsurg r=ageresid;
proc gplot;
plot residual*ageresid;
title "Added Variable Plot for Age. Refer to Figure 10.9c, page 411";run;
proc reg data=rsurg;
model LnSurviv=clotting prognos enzyme alcheav age;
title "Formal Test of Dropping Age from the Model. Refer to Page 411";run;
data rsurg;set rsurg;
nresid=residual;
proc rank normal=blom data=rsurg out=rsurg;var nresid;
proc corr;var residual nresid;
title "Correlation between Normal theory residual and observed residual. Refer to page 411";run;
proc reg data=surgery outest=vifs outvif ridge=0;
model LnSurviv=clotting prognos enzyme alcheav;
title "Quick way to get VIF's, refer to page 412";run;
proc print data=vifs;
title "Variance Influation Factors Refer to page 412";run;
data vifs;set vifs;keep clotting prognos enzyme alcheav;
if _type_="RIDGEVIF";
proc transpose data=vifs out=tvifs prefix=vif;
proc print;var _Name_ vif1;
title "VIF's in Table format- Refer to Page 412";run;
proc reg data=surgery;
model LnSurviv=clotting prognos enzyme alcheav;
output out=surgeryr r=residual h=hat rstudent=rstudent dffits=dffits cookd=cookd p=pbodyfat;
title "residuals,  deleted studentized residuals dffits, & Cook's D Refer to Page 412";
data surgeryr;set surgeryr;
obs+1;
symbol1 i=join c=black;
proc gplot;
plot (rstudent hat cookd dffits)*obs;
title "Studentized Deleted Residuals, lever, Cook's D, & DFFITS. Refer to Figure 10.10a-d, page 413";run;
data surgeryr;set surgeryr;
probf=probf(cookd,5,49);
proc print;var cookd probf;
title "probability levels of Cook's D. Refer to Page 412";run;
proc reg data=surgery;
model LnSurviv=clotting prognos enzyme alcheav;
output out=surgery1 p=predy;
reweight obs.=17;
model LnSurviv=clotting prognos enzyme alcheav;
output out=surgery2 p=predy17;
data bodyfat1;merge surgery1 surgery2;
abspdiff=abs(predy17-predy)*100;
obs+1;
proc means;var abspdiff;
title "Average Absolute percent difference in predicted Y scores, refer to page 412";
proc sort data=bodyfat1 out=bodyfat1;by descending abspdiff;
proc print;var obs abspdiff predy17 predy Lnsurviv clotting prognos enzyme alcheav;
title "Absolute Absolute percent differences in descending order, refer to page 412";run;