%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\chapter23;
PROC IMPORT OUT= WORK.bone 
            DATAFILE= "&location\chapter23.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH23TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value bonefmt
1='1=Severely Depressed  '
2='2=Moderately Depressed'
3='3=Mildly Depressed';
value sexfmt
1='1=Male'
2='2=Female';
data bone;set bone;
format sex sexfmt. bone bonefmt.;
proc sort data=bone out=bone;by sex bone;var hormone;
title "Bone Development data. Refer to Table 23.1, page 954";run;
proc means;
title 'Summary data output from bone data. ';run;
proc sort data=bone out=bone;by sex;
proc univariate;var hormone;class bone;by sex;
histogram hormone/normal;
title 'histograms of hormone levels (Not shown, but standard diagnostics)';run;
proc boxplot;by sex;
plot hormone*bone;
title "box plots by sex (Not shown, but supplemental to Figure 23.1, page 955)";run;
proc means;class sex bone;var hormone;
output out=outmeans mean=mhormone;
title "computing raw cell means";
data outmeans;set outmeans;
if sex=. or bone=. then delete;
proc print;run;
axis1 order=0 to 2.5 by .5;
symbol1 i=join;
proc gplot;
plot mhormone*bone=sex/hminor=0 vminor=0 vaxis=axis1;
title "Means of Bone Development on Hormone by Sex. Refer to Figure 23.1, page 955";run;
/* Now, it's important to understand the educational thrust of the material in Chapter 23.
The first discussion talks about how regular coded vectors, such as found in Table 23.2,
can be used to test hypotheses with unbalanced N just the same as for the balanced case,
it's just that the coefficient's aren't interpretable to the observed means in the design,
which is the rationale for the weighting strategies used. Table 23.3, for example,
shows these regressions and the SS associated with the full and various reduced models.
In practice, we wouldn't "do" the analyses this way, given that these questions are 
easily addressed within Proc GLM. I therefore present the "easy" way to do the analyses
and then, for completeness calculate the vectors and regressions which produce the
materials described on pages 956 and 957. I personally hope that this ia a bit of
overkill and, if the interpretation of the regression coefficients isn't obvious,
that should serve as an indication that further study of earlier materials on
coding is needed.*/
proc glm data=bone;class bone sex;
model hormone=bone|sex;
output out=outglm r=resid p=predict;
title 'Unbalanced two factor GLM "the easy way" reference material on page 957 & 959';run;
symbol1 v=dot;
proc capability data=outglm noprint;
  qqplot resid;
title 'qq plot for first model';run;
proc sort data=outglm out=outglm;by sex;
proc capability data=outglm noprint;
  qqplot resid;by sex;
title 'qq plot for first (full) model by sex';run;

proc glm data=bone;class bone sex;
model hormone=bone|sex/solution;
lsmeans bone/adjust=tukey pdiff cl alpha=.1;
lsmeans bone/alpha=.05 cl stderr adjust=t;
title "First LSMEANS is tukey pairwise cf page 963";
title2 "Second LSMEANS is Single DF test on page 964.";
title3 "You'll have to divide the .9 by 0.18399502 values in printout to replicate the t statistic in the book";
title4 "(You can't expect SAS to do *everything for you after all!)";run; 
data bone;set bone;
if sex=1 then x1=1;else if sex=2 then x1=-1;
label x1='Effect Coded Sex';
if bone=1 then do;x2=1;x3=0;end;else
if bone=2 then do;x2=0;x3=1;end;else
if bone=3 then do;x2=-1;x3=-1;end;
label x2='Severe vs. mild'
 x3='moderate vs. mild';
x1x2=x1*x2;
x1x3=x1*x3;
label x1x2='sex by (sev. vs. mild)'
 x1x3='sex by (mod. vs. mild)';
proc print;var sex bone hormone x1 x2 x3 x1x2 x1x3;
title "Data for Regression fits- Refer to Table 23.2, page 956";run;
proc reg;
model hormone= x1 x2 x3 x1x2 x1x3;
title "Full Model: Refer to Table 23.3a, page 957";
model hormone= x1 x2 x3 ;
title "Growth and Sex Main Effects Only: Refer to Table 23.3b, page 957";
model hormone= x2 x3 x1x2 x1x3;
title "Growth and Sex interactions Only: Refer to Table 23.33, page 957";
model hormone= x1 x1x2 x1x3;
title "Growth and Sex interactions Only: Refer to Table 23.33, page 957";
run;
/*From there, you'd simply cut and paste the incremental SS to calculate the F-statistics shown in the
  book.
For illustrative purposes, it's interesting to compare these results with what would happen
if you cast the hypotheses are relative to a group considered a control.
Suppose that you wanted to consider the mildly depressed group a control- Here is what would result.*/
data dbone;set bone;
/*Dummy coding on bone development and tests of interaction vectors*/
if bone=1 then do;bone1=1;bone2=0;end;
if bone=2 then do;bone1=0;bone2=1;end;
if bone=3 then do;bone1=0;bone2=0;end;
if sex=2 then sex=0;
sbone1=bone1*sex;
sbone2=bone2*sex;
proc reg;
model hormone=bone1 bone2/partialr2 ss1;
boneme: test bone1=bone2=0;
model hormone=bone1 bone2 sex/partialr2 ss1;
t1sex: test sex=0;
model hormone=bone1 bone2 sex sbone1 sbone2/partialr2 ss1 seqb;
inttest: test sbone1=sbone2=0;
run;
proc reg;
model hormone=bone1 bone2 sex sbone1 sbone2;
interactions: test sbone1=sbone2=0;
mains: test sex=0;
mainb: test sbone1=sbone2=0;
title 'type III ss via proc reg';run;
/*It's frequently helpful to use proc standard to center or even standardize your data in cases where
the design is unbalanced, as shown here:*/
proc standard mean=0 std=1 data=bone out=stanbone;var bone1 bone2 sex;
title 'centering predictors in an unbalanced design';
data stanbone;set stanbone;
sbone1=bone1*sex;
sbone2=bone2*sex;
run;
/*The rest of your analysis would follow as before- Centering helps to remove nonessential collinearity
  but should interpreted carefully*/
/*Missing observations in randomized designs. Page 967*/
PROC IMPORT OUT= WORK.table23_6 
            DATAFILE= "&location\chapter23.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH23TA06$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "Data. Refer to Table 23.6a, page 968";run;
data table23_6;set table23_6;
if block=1 then do;x1=1;x2=0;end;
if block=2 then do;x1=0;x2=1;end;
if block=3 then do;x1=-1;x2=-1;end;
if trt=1 then do;x3=1;x4=0;end;
if trt=2 then do;x3=0;x4=1;end;
if trt=3 then do;x3=-1;x4=-1;end;
label
x1='block 1 vs. 3'
x2='block 2 vs. 3'
x3='treatment 1 vs. 3'
x4='treatment 2 vs. 3';

proc print;var trt block y x1 x2 x3 x4;
proc glm;class trt block;
model y=trt block/solution;
title "Easy way using GLM Refer to Table 23.7a and pages 968-969";run;
proc glm;
model y=x1 x2 x3 x4/solution clparm;
title "regression approach- refer to Table 23.7b, page 969";run;
proc reg;
model y=x1 x2 x3 x4/covb corrb;
title "regression approach- refer to Table 23.7c, page 969";run;
proc glm;class block trt;
model y=block trt;
lsmeans trt/pdiff tdiff cl alpha=.05;
title "LS means comparison. Refer to page 970";run;
/* Now, for the weighting strategies given, you'd use the estimate statement in GLM.
 The data discussed aren't available in raw form, just summarized in Table 19.11
(if you can find it, let me know and I'll add it here)
but you should be able to simply take the coefficients they describe and plug them into estimate*/
PROC IMPORT OUT= WORK.TABLE23_8 
            DATAFILE= "&location\chapter23.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH23TA08$"; 
     GETNAMES=YES;
     MIXED=YES;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "Unequal weighting of means via regression. Refer to Table 23.8, page 973";run;
proc reg;
model y=z1-z4/stb noint;
title "regression weights for weighted contrasts. refer to Page 974";run;
proc reg;
model y=x5-x6 z1-z4/stb noint;
test X5=X6=0;
title "regression weights for weighted contrasts. refer to Page 974";run;
data bone;set bone;
if sex=1 and bone=1 then weight=3/14;
if sex=1 and bone=2 then weight=2/14;
if sex=1 and bone=3 then weight=2/14;
if sex=2 and bone=1 then weight=1/14;
if sex=2 and bone=2 then weight=3/14;
if sex=2 and bone=3 then weight=3/14;

proc glm data=bone;class bone sex;weight weight;
model hormone=sex bone sex*bone/solution;
means sex/cldiff t;
lsmeans sex/obsmargins om;
title "Followup tests of means assuming proportional weights. Refer to Page 979-980";run;
proc glm data=bone;class bone sex;
model hormone=sex bone /solution;
means sex/cldiff t;
lsmeans sex/obsmargins om adjust=t pdiff tdiff;
lsmeans bone/obsmargins om adjust=t pdiff tdiff;
title "Followup tests of means assuming proportional representation. Refer to Page 979-980";run;
/*There is, of course, a way to do these comparisons via an appropriate contrast statement which 
  would follow the discussion in the book more closely. I've yet to figure that out yet, though.
  What is presented above, though, probably gives you more appropriate analyses for many applied
  situations*/