%let location=C:\Documents and Settings\woodph\My Documents\glm2\kutnersolutions\chapter28;
PROC IMPORT OUT= WORK.taste 
            DATAFILE= "&location\chapter28.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH28TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title 'grand mean';run;
proc means;class  product;var rating;
title 'data listing by product';run;
proc means;class block;var rating;
title 'data listing by block';run;
proc means mean;class block product;var rating;
title "Balanced incomplete block design for five treatments. Refer to Figure 28.1, p. 1174";run;
proc glm;class block;
model rating=block/solution;
title "Note, SAS assumes an effect coding relative to the last group (block ten in this example)";run;
proc glm data=taste;class block product;
model rating=block product;
lsmeans product/cl pdiff adjust=tukey out=lsmeans;
output out=resids p=prating r=rrating;
title 'least square means & glm Refer Figure 28.2, p. 1176 & page 1181 for LSmeans';run;
proc capability noprint;
qqplot rrating;
title 'normal probability plot- reference page 1178';
proc gplot;
plot rrating*prating;
title 'residuals by predicted plot. reference page 1128';
proc print data=lsmeans;run;
symbol1 i=join v=dot;
proc gplot data=lsmeans;
plot lsmean*product/vaxis=1 to 10 by 1 vminor=0 haxis=1 to 5 by 1 hminor=0;
title "least squares means by product. Refer to Figure 28.4, page 1181";run;
/*Or, if you want, you can do these analyses via the regression approach with coded effect vectors*/
proc print;
title "coding scheme for regression model. Refer to pp. 1177-1178";run;
proc reg data=taste;
model rating=b1-b9 p1-p4/clm;
product: test p1=p2=p3=p4=0;
block: test b1=b2=b3=b4=b5=b6=b7=b8=b9=0;
title 'regression approach to design. reference page 1176';run;
/*Or, if you want to do it the old fashioned way:*/
proc reg data=taste;
model rating=b1-b9;
title "Refer to Figure 28.2(b), page 1176";run;
proc reg data=taste;
model rating=p1-p4;
title "Refer to Figure 28.2(c), page 1176";run;
/*Background Music example. p. 1187*/
PROC IMPORT OUT= WORK.music 
            DATAFILE= "&location\chapter28.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH28TA04$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value musicfmt
1="1=A"
2="2=B"
3="3=C"
4="4=D"
5="5=E";
data music;set music;format music musicfmt.;
proc means;class music;var perform;
title 'marginal means by type of music. Refer to bottom of Table 28.4, p. 1187';
proc means;class day;var perform;
title 'marginal means by day of week. Refer to Last Column of Table 28.4, p. 1187';
proc means;class week;var perform;
title 'marginal means by week. Refer to Last Row of Table 28.4, p. 1187';run;
proc means mean;class week day;
title "Data listing. Refer to Table 28.4, p. 1187";run;
proc glm data=music;class music week day;
model perform=week day music;
lsmeans music/adjust=tukey pdiff cl alpha=.1;
output out=resids r=rmusic p=pmusic rstudent=stud;
title 'lsmeans for groups- reference page 1192';run;
proc capability data=resids;
qqplot rmusic;
title 'probability plot. Refer to Figure 28.5(b), p. 1191';run;
symbol i=none v=dot;
proc gplot;
plot rmusic*pmusic/vaxis=-10 to 10 by 5 vminor=0 haxis=0 to 40 by 10 hminor=0;
title "Residual Plot. Refer to Figure 28.5(a), p. 1191";run;
proc univariate;
histogram perform stud;inset mean std;
title 'residuals and leverage';run;

/*tukey additivity*/
proc glm data=music;class music week day;
model perform=week day music;
  ods output overallanova=overall modelanova=model;
run;
quit;
ods listing;
ods output close;
data _null_;
  set overall;
  if source='Corrected Total' then call symput('overall', ss);
run;
data _null_;
  set model ;
  if hypothesistype=1 and source='week' then call symput('ssa1', ss);
  if hypothesistype=1 and source='day' then call symput('ssa2', ss);
  if hypothesistype=1 and source='music' then call symput('ssb', ss);
  if hypothesistype=1 and source='week' then call symput('dfa1', df);
  if hypothesistype=1 and source='day' then call symput('dfa2', df);
  if hypothesistype=1 and source='music' then call symput('dfb', df);
run;
proc means data=music;var perform;
output out=meant mean=meant;
title 'calculate grand mean';run;
data _null;set meant;
call symput('meant',meant);run;
proc means data=music;class week;var perform;
output out=meanr mean=meanr;
title 'calculate means by week';run;
proc means data=music;class day;var perform;
output out=meanr1 mean=meanr1;
title 'calculate means by day';run;

proc means data=music;class music;var perform;
output out=means mean=means;
title 'calculate means by music';run;
proc sort data=means out=means;by music;
proc sort data=music out=music;by music;
data allmeans;merge music means;by music;
proc sort data=meanr out=meanr;by week;
proc sort data=allmeans out=allmeans;by week;
data allmeant;merge allmeans meanr;by week;
meant=&meant;
if week>. and music>.;
drop _type_ _freq_;
proc sort data=allmeant out=allmeant;by day;
proc sort data=meanr1 out=meanr1;by day;
data allmeant;merge allmeant meanr1;by day;
meant=&meant;
if day>. and music>.;
drop _type_ _freq_;

ss=(meanr1-meant)*(meanr - meant)*(means - meant)*perform;
proc means data=allmeant;var ss;
output out=totss sum=sumss;
title 'calculate sstotal';
data _null_;set totss;
call symput('total',sumss);run;

data final;
  msa1 = &ssa1/(&dfa1+1);
  msa2 = &ssa1/(&dfa2+1);
  msb = &ssb/(&dfb+1);
  ssab = (&total*&total) / ( msa1*msa2*msb );
  ssrem = &overall - &ssa1 - &ssa2 - &ssb - ssab;
  f = ssab/( ssrem/((&dfa1+1)*(&dfa2+1)*(&dfb+1) - (&dfa1+1) - (&dfa2+1) - (&dfb+1)) );
  p_value = 1- cdf('F',f,  1, (&dfa1+1)*(&dfa2+1)*(&dfb+1) - (&dfa1+1) - (&dfa2+1) - (&dfb+1) );
run;
proc print data=final;
title 'Tukey additivity test. Talked about on p. 1191 (and here it is for you to look at too!)';run;
/*Retraining program experiment. p. 1196*/
PROC IMPORT OUT= WORK.Train 
            DATAFILE= "&location\chapter28.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH28TA08$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value methfmt
1="1=A"
2="2=B"
3="3=C";
data train;set train;
format method methfmt.;
proc means;class iq;var score;
title 'means by iq';run;
proc means;class age;var score;
title 'means by age';run;
proc means;class method;var score;
title 'means by method';run;
proc univariate;class method;
histogram score;
title 'first exploration of mean differences';
proc univariate;class iq method;
histogram score;
title 'means by iq and method';
proc univariate; class age method;
histogram score;
title 'means by age and method';
proc univariate;class age iq;
histogram score;
title 'means by age and iq';
run;
data train;set train;
lof=compress(iq||age||method);
proc glm;
class iq age method lof;
model score=iq age method lof;
title 'lack of fit glm. Refer to Table 28.8(b), p. 1197';
proc glm;class iq age method;
model score=iq age method;
title 'reduced model glm. Refer to Table 28.8(b). p. 1197';run;
/*Apple Sales Example. p. 1200*/
PROC IMPORT OUT= WORK.Apple 
            DATAFILE= "&location\chapter28.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH28TA11$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;class display;var sales;
title 'means by display';run;
proc means;class time;var sales;
title 'means by time';
proc means;class pattern;var sales;
title 'means by pattern';
proc glm;class pattern time display store;
model sales=pattern time display store(pattern);
output out=appleout r=rsales p=psales rstudent=rstud;
title 'latin square crossover design. reference page 1200';run;
proc capability data=appleout;var rsales;
qqplot rsales;
title 'probability plot. apple sales because we always do it';run;
symbol i=none v=dot;
proc gplot;
plot rsales*psales/vminor=0 hminor=0;
title "Residual Plot. apple sales because we always do it";run;