%let location=C:\Documents and Settings\woodph\My Documents\glm2\kutnersolutions\chapter24;
PROC IMPORT OUT= WORK.learn 
            DATAFILE= "&location\chapter24.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH24TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value agefmt
1='1=Young '
2='2=Middle'
3='3=Old';
value sexfmt
1='1=Male '
2='2=Female';

data learn;set learn;
format age agefmt. sex sexfmt.;
proc means;
title 'data listing. Refer to Table 24.1, page 994';run;
proc means;class sex age iq;var time;
title 'data listing via proc means. Refer to Table 24.1, page 994';run;
proc sort data=learn out=learn;by iq;
symbol1 i=join v=dot;
symbol2 i=join v=circle;
symbol3 i=join v=square;
proc gplot;by iq;
plot time*sex=age/vaxis=5 to 25 by 5;
title "sex by age plots by iq. Refer to Figure 24.1, page 998";run;
proc means;class sex age iq;var time;run;
proc glm;class sex age iq;
model time=sex age iq sex*age sex*iq age*iq/solution;
title 'initial GLM (not really realistic because of the pattern of interaction)';run;
PROC IMPORT OUT= WORK.TABLE24_2 
            DATAFILE= "&location\chapter24.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH24TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;class sex age iq;var time;
title 'data listing via proc means. Refer to Table 24.2, page 999';run;
data table24_2;set table24_2;
format age agefmt. sex sexfmt.;
proc sort data=table24_2 out=table24_2;by iq;
symbol1 i=join v=circle;
symbol2 i=join v=dot;
symbol3 i=join v=star;
proc gplot;by iq;
plot time*sex=age/vaxis=5 to 25 by 5 vminor=0 haxis=1 to 2 by 1 hminor=0;
title 'Refer to Figure 24.1 figs a & b, page 1000';
*1000 pages of statistics text- is there no end to delight?;run;
proc sort data=table24_2 out=table24age_2;by age;
proc gplot;by age;
plot time*sex=iq/vaxis=5 to 25 by 5 vminor=0 haxis=1 to 2 by 1 hminor=0;
title "Refer to Figure 24.1, Figs 24.2c, 24.2d, & 24.2e";run;
PROC IMPORT OUT= WORK.exercise 
            DATAFILE= "&location\chapter24.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH24TA04$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title 'data. refer to Table 24.2, page 1004';run;
proc sort data=exercise out=exercise;by sex smoking;
proc univariate;var stress;class bmi smoking;by sex;
histogram stress/normal;inset std mean median;
title 'histograms of stress levels';run;
proc boxplot;by sex smoking;
plot stress*bmi;
title 'box plots by sex and smoking';run;

proc glm data=exercise;class sex smoking bmi;
model stress=sex|smoking|bmi;
output out=outglm r=resid p=predict;
title 'reference plot on page 1007';run;

proc capability data=outglm noprint;
  qqplot resid;
title 'qq plot for overall model. Refer to Figure 24.6, page 1007';run;
proc format;
value fatfmt
1="1=Low Fat "
2="2=High Fat";
value smokefmt
1="1=Light Smoking"
2="2=Heavy Smoking";
data outglm;set outglm;
format sex sexfmt. bmi fatfmt. smoking smokefmt.;
proc sort data=outglm out=outglm;by sex;
proc gplot;by sex;
plot predict*bmi=smoking/vaxis=0 to 30 by 10 haxis=1 to 2 by 1 hminor=0 vminor=0;
title "Mean plots. Refer to Figure 24.5, page 1007";run;
ods output estimates=estims;
ods output OverallANOVA=anova;
proc glm data=exercise;class sex bmi smoking;
model stress=sex|smoking|bmi;
estimate 'l1' smoking 1 -1 bmi*smoking 1 -1 0 0/e;
estimate 'l2' smoking 1 -1 bmi*smoking 0 0 1 -1/e;
estimate 'l3' sex 1 -1/e;
title "Refer to Figure 24.7, page 1011 for the anova summary table.";
title2 "Refer to page 1018 for the contrasts of treatment means";run;
run;
data _null_;
set anova;
if source='Error' then call symput('Df',DF);
run;
data estims;set estims;
drop dependent tvalue probt;
bonf=tinv(1-(.05/6),&df);
lower=estimate-bonf*stderr;
upper=estimate+bonf*stderr;
proc print data=estims;
title "bonferroni contrasts of comparisons. Refer to page 1018";
run; 

proc glm data=exercise;class sex bmi smoking;
model stress=sex|smoking|bmi;
lsmeans smoking*bmi/pdiff cl alpha=0.0083333333333333333333333333333333;
title 'assuming common error';run;
proc glm data=exercise;class bmi smoking;
model stress=sex smoking|bmi;
lsmeans smoking*bmi/out=lsout;
data lsout;set lsout;
format sex sexfmt. bmi fatfmt. smoking smokefmt.;
proc gplot;
plot lsmean*smoking=bmi/vminor=0 hminor=0 vaxis=0 to 30 by 10 haxis=1 to 2 by 1;
title "refer to Figure 24.8b, page 1019";run;
PROC IMPORT OUT= WORK.exercise 
            DATAFILE= "&location\chapter24.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH24TA07$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc glm data=exercise;class bmi smoking;
model stress=sex|smoking|bmi;
title "GLM for unbalanced data. REfer to page 1020";
proc reg data=exercise;
model stress=sm1--triple;
title "same analysis but using coded regression vectors Compare with Above and refer to page 1020";run;
/*For the power analyses discussed at the end, you'd use GLMPower in SAS instead*/