%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\chapter21;
PROC IMPORT OUT= WORK.confidence 
            DATAFILE= "&location\chapter21.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH21TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value methfmt
 1='1=Utility   '
 2='2=Worry'
 3='3=Comparison';
data confidence;set confidence;
format method methfmt.;
proc print;
title "printout of data. Refer to Table 21.1, page 896";run;
proc sort data=confidence out=confidence;by method;
proc boxplot;
plot confid*method/boxconnect=mean;
title "boxplots of data";run;
symbol1 i=std2mtj v=none;
proc gplot;
plot confid*method;
title "observed standard errors around mean";run;
symbol1 i=std2mptj v=star;
proc gplot;
plot confid*method;
title "pooled (i.e., anova) standard errors around mean";run;
symbol1 i=std2tj v=none;
proc gplot data=confidence;
plot confid*method;
title "observed standard deviations around mean";run;
symbol1 i=join;
proc gplot;
plot confid*method=block;
title "Refer to Figure 21.2, page 896";run;

proc univariate;var confid;by method;
histogram confid/normal;
title "univariate histograms";run;
proc univariate;class method;var confid;
histogram confid/normal;
title "univariate histograms another way";run;

proc glm data=confidence;class method block;
model confid=method block;
output out=outconf p=yhat r=resid;
title "glm of method and block Refer to Table 21.3, page 899";run;
symbol1 i=none v=dot;
proc gplot;
plot resid*yhat;
title "Residual by Predicted plot. Refer to Figure 21.5a, page 903";run;
proc univariate;var resid;
qqplot resid/normal;
title "QQ Normal probability plot. Refer to Figure 21.5b. page 903";run;
proc glm data=confidence;
  class block method;
  model confid=block method;
  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='block' then call symput('ssa', ss);
  if hypothesistype=1 and source='method' then call symput('ssb', ss);
  if hypothesistype=1 and source='block' then call symput('dfa', df);
  if hypothesistype=1 and source='method' then call symput('dfb', df);
run;
proc means data=confidence;var confid;
output out=meant mean=meant;
title 'calculate grand mean';run;
data _null;set meant;
call symput('meant',meant);run;
proc means data=confidence;class block;var confid;
output out=meanr mean=meanr;
title "calculate means by block";run;
proc means data=confidence;class method;var confid;
output out=means mean=means;
title "calculate means by method";run;
proc sort data=means out=means;by method;
proc sort data=confidence out=confidence;by method;
data allmeans;merge confidence means;by method;
proc sort data=meanr out=meanr;by block;
proc sort data=allmeans out=allmeans;by block;
data allmeant;merge allmeans meanr;by block;
meant=&meant;
if block>. and method>.;
drop _type_ _freq_;

ss=(meanr - meant)*(means - meant)*confid;
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;
  msa = &ssa/(&dfb+1);
  msb = &ssb/(&dfa+1);
  ssab = (&total*&total) / ( msa*msb );
  ssrem = &overall - &ssa - &ssb - ssab;
  f = ssab/( ssrem/((&dfa+1)*(&dfb+1) - (&dfa+1) - (&dfb+1)) );
  p_value = 1- cdf('F',f,  1, (&dfa+1)*(&dfb+1) - (&dfa+1) - (&dfb+1) );
run;
proc print data=final;
title "Tukey additivity test. Refer to page 904";run;
run;
proc glm data=confidence;
  class block method;
  model confid=block method;
lsmeans method/adjust=tukey pdiff cl alpha=.05;
title "Tukey follow up tests. Refer to Page 905";run;
data table21_4;
input time distract sex;
cards;
12 1 1
 8 1 1
 7 1 1
 5 1 1
 3 1 2
 9 1 2
 5 1 2
 9 1 2
14 2 1
16 2 1
15 2 1
13 2 1
11 2 2
 9 2 2
10 2 2
14 2 2 
;
proc print;
title "More than one replicate per block. Refer to Table 21.4, page 907";run;
proc glm;class distract sex;
model time=distract|sex;
title "More than one replicate in each block. Refer to Figure 21.6, page 908";run;
  /* Calculating efficiency (page 911-912)*/
proc print data=model;run;
proc print data=overall;run;
data _null_;set overall;
file print;
if source='Error' then call symput('mse', ms);
if source='Error' then call symput('dfe', df);
if source='Corrected Total' then call symput('totdf',df);
data efficient;file print;
eff=(&ssa+(&dfa+1)*&dfb*&mse)/
  (&totdf*&mse);
  dfb=&dfa+&dfe;
  dfa=&dfe;
eff2=eff*((dfa+1)*(dfb+3))/((dfa+3)*(dfb+1));
put 'raw efficiency is ' eff;
put 'adjusted efficiency is ' eff2;
title "Calculating Raw and adjusted efficiency- Refer to Page 912";
run;

/* power analyses on page 910-911*/
data means;
input factora depvar;
cards;
0 0
1 1.5
2 3
;
proc glmpower data=means;
class factora ;
model depvar=factora;
power
stddev=2
ntotal=30
power=.;
plot x=n min=5 max=120;
title 'simple power analysis with plot';
title2 'note that ntotal represents the full sample size';run;

/*One way anovas can also be done using Proc Power*/
proc power;
ONEWAYANOVA
std=2
groupmeans= 0|1.5|3
npergroup= 7 8 9 10 11 12 15 18 21 24 27 30
power=.;
plot x=n min=6 max=60;
title "three-group anova results";
title2 "note that npergroup is the number of required blocks";run;