%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter11;
PROC IMPORT OUT= WORK.pressure 
            DATAFILE= "&location\chapter11.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH11TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Blood Pressure Example. Refer to Table 11.1, page 427";run;
symbol i=rlcli95 v=dot;
proc gplot;
plot diastole*age;
title "Plot of Blood Pressure by Age. Refer to Figure 11.1(a), page 428";run;
proc reg;
model diastole=age;
output out=rpressur r=resid p=predict;
proc gplot;
plot resid*age;
title "Plot of Residual by Age. Refer to Figure 11.1(b), page 428";run;
data rpressur;set rpressur;
absresid=abs(resid);
proc gplot;
plot absresid*age;
title "Plot of Absolute Value of Resid by Age. Refer to Figure 11.1(c), page 428";run;
proc reg;
model diastole=age;
title "Initial Regression weight, Refer to Equation 11.18, page 428";run;
proc reg data=rpressur;
model absresid=age;
output out=rpressur p=weight;
title "Regression of Age on Absolute value of Residual. Refer to Eq. 11.19, page 428";run;
proc print;title "initial weights";run;
data rpressur;set rpressur;
weight=1/(weight*weight);
proc print;
title "individual estimated weights. See Discussion, page 428";run;
proc glm data=rpressur;
weight weight;
model diastole=age/clparm alpha=.05;
title "WLS estimates. Refer to Equation 11.20, page 428";run;

PROC IMPORT OUT= WORK.bodyfat 
            DATAFILE= "&location\chapter11.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH07TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Bodyfat data. Refer to page 434";run;
proc standard mean=0 std=1 data=bodyfat out=bodyfatz;var bodyfat triceps thigh midarm;
symbol1 i=spline c=blue line=2 v=dot;
symbol2 i=spline c=red line=3 v=dot;
symbol3 i=spline c=green line=1 v=dot;
proc reg data = bodyfatz outest = temp outstb noprint;
  model bodyfat = triceps thigh midarm/ ridge = (0.001 to 0.1 by .001) outvif  ;
  plot / ridgeplot vref=0;
title "Ridge estimated standardized regression weights. Refer to Figure 11.3, page 435";run;

/*Or, alternately, you can create the estimates and plot them via proc plot like this:
  (Which means you can make the tables indicated in Tables 11.2 and 11.3 on page 4343)*/
proc reg data=bodyfatz outvif outest=ridgev 
ridge=0 .002 .004 .006 .008 .01 .02 .03 .04 .05 .1 .5 1;
model bodyfat=triceps thigh midarm;

data betas;set ridgev;
file print;
if _type_="RIDGE";
put "ridge=" _ridge_ " triceps=" triceps 8.3 " thigh=" thigh 8.3 " midarm=" midarm 8.3;
title "Ridge estimated standardized regression weights. Refer to Table 11.2, page 434";run;
data VIFs;set ridgev;
file print;
if _type_="RIDGEVIF";
put "ridge=" _ridge_ " VIF(triceps)=" triceps 8.3 " VIF(thigh)=" thigh 8.3 " VIF(midarm)=" midarm 8.3;
title "Ridge estimated Variance Inflation Factors. Refer to Table 11.3, page 434";run;
symbol1 i=spline c=blue line=2 v=dot;
symbol2 i=spline c=red line=3 v=dot;
symbol3 i=spline c=green line=1 v=dot;
data betas;set betas;
if _ridge_=0 then logval=log10(_ridge_+.001);else
logval=log10(_ridge_);
axis1 value=('.001' '.01' '.1' '1') ;
proc gplot data=betas;
plot (triceps thigh midarm)*logval/overlay vref=0 haxis=axis1;
title "Plot of standardized regression weights by ridge coefficients. See Figure 11.3, page 435";run;

proc reg data=bodyfatz outest=ridgefinal outvif
ridge=.02;
model bodyfat=triceps thigh midarm/stb;
title "Final Ridge Solution. Refer to page 435";run;
proc print data=ridgefinal;run;
proc means data=bodyfat;run;
data ridgefinal;set ridgefinal;
if _type_="RIDGE";
sy=5.1061856;
striceps=5.0232591;
sthigh=5.2346115;
smidarm=3.6471474;
mtriceps=25.3050000;
mthigh=51.1700000;
mmidarm=27.6200000;
my=20.1950000;
triceps=triceps*sy/striceps;
thigh=thigh*sy/sthigh;
midarm=midarm*sy/smidarm;
int=my-triceps*mtriceps-thigh*mthigh-midarm*mmidarm;
proc print;
title "unstandardized coefficients for ridge solution. Refer to page 435 (Differences due to rounding error in book)";run;
data huberbi;
do u=-5 to 5 by .02;
if abs(u)<=1.345 then huber=1;else huber=1.345/abs(u);
if abs(u)<=4.685 then bisquare=(1-(u/4.685)**2)**2;else bisquare=0;
output;
end;
symbol1 i=spline c=red v=none;
symbol2 i=spline c=blue v=none;
proc gplot;
plot (bisquare huber)*u/overlay;
title "huber and bisquare weight functions. Refer to Figures 11.4a and 11.4b, page 440";run;
PROC IMPORT OUT= WORK.ets 
            DATAFILE= "&location\chapter11.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH11TA04$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Math Proficiency Example. Refer to Table 11.4, page 443";run;
symbol1 i=rl v=dot color=black;
proc gplot data=ets;
plot mathprof*homelib;
title "Linear Regression Plot. Refer to Figure 11.5(a) page 442";run;
data etsshort;set ets;if state="Guam" or state="D.C." or state="Virgin Islands" then delete;
proc sort data=etsshort out=etsshort;by homelib;*Data should be sorted by predictor first;
proc loess data=etsshort;
model mathprof=homelib/direct smooth=.01 .03 .05 .1 .2 .3 .4 .5 .6 .7 .8 .9 1
dfmethod=exact;
ods output outputstatistics=loessout;
proc print data=loessout;run;
proc sort data=loessout out=loessout;by smoothingparameter homelib;
symbol1 c=black v=dot i=rlclm95;
symbol2 c=red v=point i=join l=1;
symbol3 c=blue v=none i=rlcli95;
proc gplot data=loessout;by smoothingparameter;
plot (depvar pred depvar)*homelib/overlay;run;
title 'Compare with Figure 11.5 page 442' ;run;
proc loess data=etsshort;
model mathprof=homelib/degree=2 smooth=.2
;
ods output outputstatistics=loessout;
proc sort data=loessout out=loessout;by smoothingparameter homelib;
proc print data=loessout;run;

symbol1 c=black v=dot i=rlclm95;
symbol2 c=red v=point i=join l=1;
symbol3 c=blue v=none i=rlcli95;
proc gplot data=loessout;
plot (depvar pred depvar)*homelib/overlay;run;
title "Refer Figure 11.5 A, page 442 (Abscissa it text should be X2 instead of X3)";
title2 "Outliers were removed from Loess regressions";run;
proc reg data=ets;
model mathprof=homelib/stb;
plot residual.*homelib;
title "Residuals from Linear Regression. Refer to Figure 11.5(b), page 442";run;
symbol1 i=rq v=dot color=black;
proc gplot data=ets;
plot mathprof*homelib;
title "Plot of OSL quadratic Fit. Refer to Figure 11.5(c), page 442";run;
proc reg data=ets;
model mathprof=homelib;
plot cookd.*obs.;
title "Plot of Cook's D by observation- Refer to Figure 11.5 (d), page 442";run;
data ets;set ets;
chome=homelib;
proc standard data=ets out=ets mean=0;var chome;
data ets;set ets;
chome2=chome*chome;
symbol1 i=rl c=black v=dot;
symbol1 i=rl c=red v=dot;
%include "&location\mad.sas";
%include "&location\robust_hubert.sas";
%include "&location\robust_bisquare.sas";
title "Robust quadratic fit and robust weights, Refer to Figure 11.5e and 11.5f, page 442";
%robust_hubert(ets, mathprof, chome chome2, 0.000005, 8);run;
title "Same analysis, but using Bisquare (not shown in text) ";
%robust_bisquare(ets, mathprof, chome chome2, 0.000005, 8);
/*Note for the Future (well, this spring, really)
  This code will not work unless you have SAS version 9.2, but I provide it just so you can see what
  it looks like. This is much easier than using Michael Friendly's macro & allows you to draw loess
  regression lines in the graph as well, as shown in Figure 11.6, page 446*/
/*
proc sgscatter data=ets
;
plot (mathprof parents homelib reading tvwatch absences)*
(mathprof parents homelib reading tvwatch absences) / loess;
run;
*/
proc corr data=ets;var mathprof parents homelib reading tvwatch absences;
title "Initial correlation matrix. Refer to Figure 11.6b, page 446";run;
%include "C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter11\scatter.sas";
%scatter(data = ets, var = mathprof parents homelib reading tvwatch absences);run;
proc reg;
model mathprof=parents homelib reading tvwatch absences/stb influence;
id state;
plot predicted.*(parents homelib reading tvwatch absences);
title "Initial regression model. Refer to Equation 11.52, page 445";
title2 "and discussion on pages 446-447";run;
proc reg;
model mathprof=parents homelib reading tvwatch absences/selection=rsquare cp adjrsq rmse;
title "all possible subsets regression. Refer to Figure 11.7, p. 448";run;
title "Robust quadratic fit and robust weights, Refer to page 447";
%robust_hubert(ets, mathprof, homelib reading tvwatch, 0.000005, 8);run;
proc reg;
model mathprof=homelib reading tvwatch;
title "Ordinary Least Squares regression for comparison. Refer to page 447";run;

proc reg data=_tempout_;
weight _w2_;
model mathprof=homelib reading tvwatch;
output out=etsw p=pmathprof ;
proc corr spearman;var mathprof pmathprof;
title "Spearman correlation between observed and predicted math proficiency. See page 448";run;
PROC IMPORT OUT= WORK.insure 
            DATAFILE= "&location\chapter11.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH11TA07$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Life Insurance Example. Refer to Table 11.7, p. 451";run;

proc loess data=insure;
model insurance=risk income/direct smooth=.5
dfmethod=exact interp=cubic;
ods output outputstatistics=loessout;

data insure;merge insure loessout;
proc g3grid data=insure out=gridvals;
grid income*risk=pred;
run;
proc g3d data=gridvals;
   plot income*risk=pred;
title "3d grid of predicted insurance based on risk and income";
run;

proc gcontour data=gridvals;
   plot income*risk=pred;
title "Contour plot of Loess regression. Compare with Figure 11.8, page 452";run;
   
proc means data=gridvals;class income;var risk pred;
title "note the interpolated values are what you can use for slices- these are not exact slices";run;
data smallg;set gridvals;
if income=2.8 or (income>6 and income<7) or income=9.1;
proc sort data=smallg out=smallg;by income;
symbol1 i=spline v=dot c=blue;

proc gplot;by income;
plot pred*risk;
title "slices of Loess regressions. Compare (loosely) with Figure 11.8(b), page 452";run;

PROC IMPORT OUT= WORK.Steroid 
            DATAFILE= "&location\chapter11.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH11TA08$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Steroid Level example. Compare with Table 11.8, page 453";run;
data steroid;set steroid;
if age>=8 and age<9 then region=1;else
if age>=9 and age<10 then region=2;else
if age>=10 and age<13 then region=3;else
if age>=13 and age<14 then region=4;else
if age>=14 and age<25 then region=5;
proc means;class region;var steroid;
title "Final predicted values. Compare with Table 11.8, page 453 columns 4-6";
run;