%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\chapter26;
PROC IMPORT OUT= WORK.learn 
            DATAFILE= "&location\chapter26.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH26TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value schlfmt
1='1=Atlanta      '
2='2=Chicago      '
3='3=San Francisco';run;
data learn;set learn;
format school schlfmt.;run;
proc means;class school teacher;var learning;
title 'data summary';run;
symbol1 v=dot;
symbol2 v=circle;
proc gplot;
plot school*learning=teacher/vminor=0 hminor=0 haxis=0 to 30 by 10 vaxis=0 to 4 by 1;
title 'gplot of the data (reference 26.3, page 1096)';run;
proc glm data=learn;class teacher school;
model learning=school teacher(school)/solution;
lsmeans teacher(school)/slice=school cl alpha=.1;
title 'ss(teacher) via slice, reference 26.4, page 1097';run;

proc glm data=learn;class teacher school;
model learning=school teacher(school)/solution;
contrast 'Instructors in Atlanta' teacher(school) 1 -1 0 0 0 0;
contrast 'Instructors in Chicago' teacher(school) 0 0 1 -1 0 0;
contrast 'Instructors in San Fran' teacher(school) 0 0 0 0 1 -1;
output out=resids p=plearn r=rlearn;
title 'nested effects with glm (reference 26.4, page 1097)';run;
proc capability data=resids noprint;
  qqplot rlearn;
title 'qq plot for overall model reference 26.4, page 1100';run;
proc gplot data=resids;
plot school*rlearn;
title 'residual dot plot- reference 26.4, page 1100';run;
proc glm data=learn;class teacher school;
model learning=school teacher(school)/solution;
lsmeans school/adjust=tukey pdiff cl alpha=.05;
title 'estimation of school means reference p. 1101';run;
proc glm data=learn;class teacher school;
model learning=school teacher(school)/solution;
lsmeans school/adjust=tukey pdiff cl alpha=.1;
title 'estimation of school means reference 90% CI p. 1101';run;
proc glm data=learn;class teacher school;
model learning=school teacher(school);
lsmeans teacher(school)/alpha=.03333 slice=school cl pdiff;
title 'estimation of intructor differences, reference page 1102';run;
/*Considering unbalanced case in table 26.6*/
PROC IMPORT OUT= WORK.learn 
            DATAFILE= "&location\chapter26.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH26TA06$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;class school teacher;var learning;
title 'unbalanced example data summary See data on page 105';run;
proc gplot;
plot school*learning=teacher;
title 'gplot of the data compare to full data in fig. 26.3, page 1096';run;
proc glm data=learn;class teacher school;
model learning=school teacher(school)/solution;
lsmeans teacher(school)/slice=school cl alpha=.1;
output out=resids p=plearn r=rlearn;
title 'ss(teacher) via slice, reference 26.7, page 1106';run;

proc capability data=resids noprint;
  qqplot rlearn;
title 'qq plot for overall model unbalanced data compare to reference 26.4, page 1100';run;
proc gplot data=resids;
plot school*rlearn;
title 'residual dot plot unbalanced data compare to reference 26.4, page 1100';run;
proc glm data=learn;class teacher school;
model learning=school teacher(school)/solution;
lsmeans school/adjust=tukey pdiff cl alpha=.05;
title 'estimation of school means unbalanced data compare to reference p. 1101';run;
proc glm data=learn;class teacher school;
model learning=school teacher(school)/solution;
lsmeans school/adjust=tukey pdiff cl alpha=.1;
title 'estimation of school means unbalanced data compare to reference 90% CI p. 1101';run;
proc glm data=learn;class teacher school;
model learning=school teacher(school);
lsmeans teacher(school)/alpha=.03333 slice=school cl pdiff;
title 'estimation of intructor differences unbalanced data compare to reference page 1102';run;
/*Bread Crustiness, page 1109*/
PROC IMPORT OUT= WORK.bread 
            DATAFILE= "&location\chapter26.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH26TA09$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc format;
value tempfmt
1='1=Low   '
2='2=Medium'
3='3=High';
data bread;set bread;
format temp tempfmt.;
proc means mean;class temp batch obs;
title 'data design Refer to Table 26.9, page 1109';run;
proc gplot;
plot temp*crust=batch/vaxis=0 to 4 by 1 vminor=0 haxis=0 to 20 by 5 hminor=0;
title 'gplot of data reference 26.5, page 1109';run;
proc glm;class temp batch;
model crust=temp batch(temp);
random batch(temp);
  temperature_effect: test h=temp e=batch(temp);
  lsmeans temp / e=batch(temp) pdiff cl;
title 'glm analysis reference 26.10 and text page 1110';run;
proc mixed data=bread;
  class temp batch;
  model crust = temp;
  random batch(temp)/cl;
  lsmeans temp / cl pdiff;
title "glm analysis via Proc Mixed. Reference 26.10 and text page 1110";
run;
/*Group Decision Making Example*/
PROC IMPORT OUT= WORK.gdecide 
            DATAFILE= "&location\chapter26.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH26TA12$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means mean;class nation grpsize observe subno;var ninterac;
proc format;
value obsfmt
 1="1=US Teams Observer 1"
 2="2=US Teams Observer 2"
 3="3=US Teams Observer 3"
 4="4=US Teams Observer 4";
 value gfmt
 1="1=Four Members"
 2="2=Eight Members";
 value nationfmt
 1="1=US Teams     "
 2="2=Foreign Teams";
data gdecide;set gdecide;
observen=nation*(nation-1)+observe;
format observen obsfmt. grpsize gfmt. nation nationfmt.;
label grpsize="Group Size"
ninterac="Number of Interactions";
title "Data for Crossed-Nested Three-Factor Study- Group Decision Making Example";run;
proc gplot;
plot grpsize*ninterac=observen/vaxis=0 to 2 by 1 haxis=0 to 30 by 10 vminor=0 hminor=0;
title "Dot Plots for Cross-Nested Design. Refer to Figure 26.6, page 1117";run;
proc glm;class nation grpsize observe;
model ninterac=nation grpsize observe(nation) nation*grpsize grpsize*observe(nation);
test h=nation e=observe(nation);
test h=grpsize e=grpsize*observe(nation);
test h=nation*grpsize e=grpsize*observe(nation);
random observe(nation)/test;
lsmeans nation/cl adjust=t pdiff tdiff e=observe(nation);
title "Cross-Nested Design Experiment- Group Decision Making. Refer to Figure 26.7, page 1117";run;
/*Health Awareness example. Problem 26.9*/
PROC IMPORT OUT= WORK.health 
            DATAFILE= "&location\chapter26.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH26PR09$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;class state city;
var aware;
title 'design means';run;
proc glm;class state city;
model aware=state city(state);
output out=resids p=paware r=raware;
title 'residuals of the model 29.a';
run;
proc capability data=resids noprint;var raware;
  qqplot raware;
title 'qq plot for overall model reference 26.9a';run;
data resids;set resids;
nresid=raware;
proc rank normal=blom data=resids out=resids;var nresid;
proc corr;var nresid raware;
title 'correlation between normal and actual residuals 26.9a';run;
proc gplot;
plot state*raware;
title 'residuals by state';run;
proc gplot;
plot state*raware=city;
title 'residual plots by city within state';run;
proc glm;class state city;
model aware=state city(state);
title 'test of state differences 26.10ab';run; 
proc glm data=health;class state city;
model aware=state city(state);
lsmeans state/ alpha=.01 cl;
lsmeans state/ alpha=.1 cl pdiff;
lsmeans state/ alpha=.1 cl pdiff adjust=tukey;
lsmeans city(state)/slice=state alpha=.05 cl pdiff;
title '26.11, 95% confidence interval for mu 11';run;

proc glm data=health;class state city;
model aware=state city(state);
random city(state);
lsmeans state / cl pdiff alpha=.1 adjust=tukey;
title '26.12';run;
proc glm;class state city;
model aware=state city(state);
random state city(state);
test h=state e=city(state);
title '26.13';run;

proc mixed data=health method=reml;
  class state city;
  model aware = /cl alpha=.01;
  random state city(state)/cl;  
run;

ods listing close;
ods output OverallANOVA=anova FitStatistics=fit modelanova=model;
proc glm data=health;class state city;
model aware=state city(state);
random state city(state);
run;quit;
ods listing;
proc print data=anova;run;
proc print data=fit;run;
proc print data=model;run;
data _null_;set model;
if hypothesistype=3 and source='state' then call symput('mstr',ms);
if hypothesistype=3 and source='state' then call symput('dfmodel',df);
if hypothesistype=3 and source='city(state)' then call symput('mse',ms);
if hypothesistype=3 and source='city(state)' then call symput('dferr',df);

run;
options macrogen;
data temp;
  lower_f = finv( .95, &dfmodel, &dferr);
  upper_f = finv( .05, &dfmodel, &dferr);
  mstr =&mstr;
  mse = &mse;
  L = (1/15)*( (&mstr/&mse)*(1/lower_f) - 1 );*formula 25.18a;
  U = (1/15)*( (&mstr/&mse)*(1/upper_f) - 1 );*formula 25.18b;
  L_star = L/(1+L);*formula 25.19;
  U_star = U/(1+U);
run;
proc print data=temp;
title 'icc and confidence interval for the driver effect';
run;
data _null_;set model;
if hypothesistype=3 and source='state' then call symput('mstr',ms);
if hypothesistype=3 and source='state' then call symput('dfmodel',df);
if hypothesistype=3 and source='city(state)' then call symput('mse',ms);
if hypothesistype=3 and source='city(state)' then call symput('dferr',df);

run;
options macrogen;
data temp;
  lower_f = finv( .95, &dfmodel, &dferr);
  upper_f = finv( .05, &dfmodel, &dferr);
  mstr =&mstr;
  mse = &mse;
  L = (1/15)*( (&mstr/&mse)*(1/lower_f) - 1 );*formula 25.18a;
  U = (1/15)*( (&mstr/&mse)*(1/upper_f) - 1 );*formula 25.18b;
  L_star = L/(1+L);*formula 25.19;
  U_star = U/(1+U);
run;
proc print data=temp;
title 'icc and confidence interval';
run;
data _null_;set model;
if hypothesistype=3 and source='state' then call symput('mstr',ms);
if hypothesistype=3 and source='state' then call symput('dfmodel',df);
if hypothesistype=3 and source='city(state)' then call symput('mse',ms);
if hypothesistype=3 and source='city(state)' then call symput('dferr',df);

run;

/* MLS procedure*/
data interval;
dfe=&dferr;dfm=&dfmodel;
mstr=&mstr;mse=&mse;
c1=1/15;
c2=-1/15;
lhat=(c1*&mstr+c2*&mse);
f1=finv(.995,&dfmodel,9999);
f2=finv(.995,&dferr,9999);
f3=finv(.995,9999,&dfmodel);
f4=finv(.995,9999,&dferr);
f5=finv(.995,&dfmodel,&dferr);
f6=finv(.995,&dferr,&dfmodel);
g1=1-(1/f1);
g2=1-(1/f2);
g3=((f5-1)**2-(g1*f5)**2-(f4-1)**2)/f5;
g4=f6*(((f6-1)/f6)**2-(((f3-1)/f6)**2)-g2**2);
hl=(((g1*c1*&mstr)**2)+(((f4-1)*c2*&mse)**2)-g3*c1*c2*&mse*&mstr)**.5;
hu=(((f3-1)*c1*&mstr)**2+((g2*c2*&mse)**2)-(g4*c1*c2*&mse*&mstr))**.5;
cl95l=lhat-hl;
cl95u=lhat+hu;
sqrt95l=sqrt(cl95l);
sqrt95u=sqrt(cl95u);
output;
run;
proc print data=interval;
title 'MLS interval for the driver effect';run;

proc mixed data=health method=reml asycov cl covtest alpha=.1;
class city state;
model aware=/ddfm=satterth cl alpha=.05 intercept;
random state city(state)/cl alpha=.05;
title 'upper and lower confidence intervals using Proc Mixed';
run;
data _null_;set model;
if hypothesistype=3 and source='state' then call symput('mstr',ms);
if hypothesistype=3 and source='state' then call symput('dfmodel',df);
if hypothesistype=3 and source='city(state)' then call symput('mse',ms);
if hypothesistype=3 and source='city(state)' then call symput('dferr',df);

run;

data temp;
mstr=&mstr;mse=&mse;
dftr=&dfmodel;dfe=&dferr;
  s_mu2 = ( &mstr - &mse)/15;
  dfround = round( (15*s_mu2)**2/( (&mstr**2)/&dfmodel + (&mse**2)/&dferr ), 1);
  dfexact =(15*s_mu2)**2/( (&mstr**2)/&dfmodel + (&mse**2)/&dferr );
  lower_chi = cinv( .975, dfround);
  upper_chi = cinv( .025, dfround);
  lower2 = (dfexact*s_mu2)/lower_chi;
  upper2 = (dfexact*s_mu2)/upper_chi;
  lower = sqrt(lower2);
  upper = sqrt(upper2);
  output;
run;
proc print data=temp;
title 'Satterthwaite procedure';
run;