%let location=C:\Documents and Settings\woodph\My Documents\glm2\kutnersolutions\chapter25;
PROC IMPORT OUT= WORK.ratings 
            DATAFILE= "&location\chapter25.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH25TA01$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "data means. Refer to Table 25.1, page 1036";run;
proc means;class officer candidat;var rating;
title "data listing. Refer to Table 25.1, page 1036";run;
symbol1 v=dot;
proc gplot;
plot officer*rating/vaxis=1 to 5 by 1 haxis=40 to 90 by 10 hminor=0 vminor=0 noframe;
title 'Reference figure 25.2 page 1036';run;
proc glm;class officer;
model rating=officer;
random officer/test;
means officer/hovtest=bf clm;
title 'random effects anova. Refer to Table 25.2, page 1037';run;
proc glm;class officer;
model rating=officer;
title 'fixed effects anova output for comparison- note that the ss are the same';run;
/*First Method: random mean and confidence intervals- using Proc Mixed*/
proc mixed method=reml asycov cl covtest alpha=.1;
class officer;
model rating=/cl alpha=.1;
random officer;
title 'upper and lower confidence intervals using Proc Mixed';
title2 "refer to page 1039";run;

/* Calculating the confidence interval for the grand mean
 It doesn't seem that one can do this in SAS via any program statements. 
 You more or less have to go out and get the numbers in the book as
 described in the book. I'll annotate what I'm doing as I go along*/

*Tell SAS to write the anova results to a data set called anova
 and to write the fit statistics for the model to a data set called fit;

ods listing close;
ods output OverallANOVA=anova FitStatistics=fit;
proc glm;class officer;
model rating=officer;
random officer;
run;quit;
ods listing;

data _null_;set anova;
*We're going to need the mean square treatment (in this case the mean square associated
 with officer rating) and the total degrees of freedom for the model- I'm writing these
 to two macro variables, mstr and df;
if source = 'Model' then call symput( 'mstr', ms);
  if source = 'Corrected Total' then call symput('df', df);
data _null_;
  set fit;
*Now, we'll output the overall dependent mean (because that's the statistic
 we're trying to get the confidence interval on) we'll call that macro variable "mean";
  call symput('mean',  DepMean);
run;

data temp;
  t = tinv( .95,4);*we need the critical t associated with the interval we want.;
  df = &df+1;*total df is r x n, see page 1038, formula 25.12;
  s2 = &mstr/df;*formula 25.13;
  Lower = &mean - t*sqrt(s2);*formula 25.15;
  Upper = &mean + t*sqrt(s2);
run;
proc print data=temp;
title "calculated confidence interval on mu, refer to pp. 1038-1039";
run;
/*How to calculate the 90% confidence interval of the ICC*/
*As before, we need the anova summary table;
ods listing close;
ods output overallanova=anova;
proc glm data=ratings;
   class officer;
   model rating = officer;
run;
quit;
ods listing;
data _null_;
  set anova;
  if source='Model' then do;
          call symput('mstr', ms);*capture the ms treatment (officers);
		  call symput('dfmodel', df);*capture the df;
		  end;
  if source='Error' then do;
          call symput('mse', ms);*capture the mean square error;
          call symput('dferr', df);*capture error df;
		  end;
run;
options macrogen;
data temp;
  lower_f = finv( .95, &dfmodel, &dferr);
  upper_f = finv( .05, &dfmodel, &dferr);
  mstr =&mstr;
  mse = &mse;
  L = (1/4)*( (&mstr/&mse)*(1/lower_f) - 1 );*formula 25.18a;
  U = (1/4)*( (&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 "Confidence interval on ICC. Refer to page 1040-1041";
run;
/* MLS procedure*/
data interval;set ratings;
dfe=&dferr;dfm=&dfmodel;
mstr=&mstr;mse=&mse;
c1=1/&dfmodel;
c2=-1/&dfmodel;
lhat=(&mstr-&mse)/&dfmodel;
f1=finv(.95,&dfmodel,9999);
f2=finv(.95,&dferr,9999);
f3=finv(.95,9999,&dfmodel);
f4=finv(.95,9999,&dferr);
f5=finv(.95,&dfmodel,&dferr);
f6=finv(.95,&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;
cl90l=lhat-hl;
cl90u=lhat+hu;
sqrt90l=sqrt(cl90l);
sqrt90u=sqrt(cl90u);
run;
proc print data=interval (obs=1);
title "MLS procedure- refer to Table 25.3 & numbers on pp. 1046-1047";run;
/*Sheffield Foods Example*/
PROC IMPORT OUT= WORK.fat 
            DATAFILE= "&location\chapter25.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH25TA12$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;class method lab;var fat;
title "Sheffield Foods Example. Refer to Table 25.12, page 1071";run;
proc sort data=fat out=fat;by method;
proc gplot;by method;
plot lab*fat/vaxis=0 to 4 by 1 haxis=2.5 to 5.5 vminor=0 hminor=0 noframe;
title 'reference figure 25.3 lab plots by method';run;

proc means data=fat;class method lab;
var fat;
output out=meand mean=meanfat;
data meand;set meand;
if method>. and lab>.;
symbol1 i=join v=dot;
proc gplot data=meand;
plot meanfat*lab=method;
title 'reference figure 25.4 estimated treatment means';run;

data effect;set fat;
if method=1 then effect1=1;
else if method=2 then effect1=-1;


proc mixed method=ml asycov;
class lab;
model fat=effect1/cl;
random lab effect1*lab;
title 'proc mixed results for estimate and se';
run;

proc glm;class lab effect1;
model fat=effect1|lab/solution;
random lab effect1*lab/test q;
title 'unbalanced analysis, but done in GLM';run;