%let location=C:\Documents and Settings\woodph\My Documents\glmwood\kutnersolutions\Chapter16;
PROC IMPORT OUT= WORK.kenton 
            DATAFILE= "&location\Chapter16.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH16TA04$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc print;
title "Kenton Food Company example. Refer to Table 16.1, page 686";run;
proc means sum mean n maxdec=1;class package;var cases;
title "Summary Statistics: Kenton Food Company example. Refer to Table 16.1, page 686";run;
proc means sum mean n maxdec=1;var cases;
title "Summary Statistics: Kenton Food Company example. Refer to Table 16.1, page 686";run;
axis1 offset=(0.5 cm) minor=none;
axis2 order=5 to 35 by 5;
symbol1 c=black v=circle;
symbol2 c=red v=circle;
symbol3 c=yellow v=circle;
symbol4 c=blue v=circle;
symbol5 c=green v=circle;
proc gplot;
plot cases*package=store/haxis=axis1 vaxis=axis2;
title "Scatter plot of number of cases sold by package design. Kenton foods- Refer ot Figure 16.3, page 686";run;
/*Maximum likelihood estimators for anova*/
/*You can get an idea of the multivariate normal likelihood by looking at the marginal likelihoods:*/
proc means data=kenton;var cases;class package;
output out=mkenton mean=mcases;
data mkenton;set mkenton;if package>.;
proc sort data=mkenton out=mkenton;by package;
proc sort data=kenton out=kenton;by package;
data kenton;merge kenton mkenton;by package;
data Like;set kenton;
Q=(cases-mcases)**2.;*Eq. 16.15;
proc means sum;var q;
title "Refer to Equation 16.15 p. 687";run;
proc means sum;class package;var q;
title "Refer to Eq. 16.15a p. 687";run;
data mldata;set kenton;
do lval=10 to 33;
 q=(cases-lval)**2.;
 output;
end;
proc means data=mldata;class package lval;var q;
output out=marg_ml sum=sigma_q;
data marg_ml;set marg_ml;if package>. and lval>.;
proc print data=marg_ml;run;
proc means data=marg_ml;run;

symbol1 v=dot i=spline l=1 color=green width=2;
symbol2 v=circle i=spline l=2 color=blue width=2;
symbol3 v=triangle i=spline l=3 color=red width=2;
symbol4 v=square i=spline l=4 color=black width=2;
axis1 minor=none;
proc gplot data=marg_ml;
plot sigma_q*lval=package/noframe vaxis=axis1 haxis=axis1 href=(14.6 13.4 19.5 27.2);
title "marginal likelihood values for Kenton foods example. Refer to page 688";run;
proc glm data=kenton;class package;
model cases=package;
output out=kenton r=resids;
proc sort data=kenton out=kenton;by package;
proc transpose data=kenton out=kentont prefix=store;by package;var resids;run;
data kentont;set kentont;drop _name_;
label package="Package Design";
total=sum(of store1-store5);
proc print;format _numeric_ 4.1;
title "refer to Table 16.2, page 689";run;

data kenton;set kenton;
id=package*10+store;
ccases=cases;
proc standard data=kenton out=kenton mean=0;var ccases;
symbol1 i=none v=none;
symbol2 i=needle;
proc gplot;
plot cases*id;
plot2 ccases*id;
title "Refer to 16.4a Partitioning of Total Variance to cell means and mean deviations";
run;
/*I'm stumped as to how to do the other two simply ;-)
  One option is to fill in your own lines given this plot*/
symbol1 i=none v=dot;
symbol2 i=none ;
proc gplot;
plot mcases*id;
plot2 ccases*id/vref=0;
title "Refer to 16.4bc Partitioning of Total Variance to cell means and mean deviations";
run;
proc glm data=kenton;class package;
model cases=package;
lsmeans package/cl alpha=.05;
title "Oneway Anova summary table. Refer to Figure 16.5, page 695; 699";run;
proc univariate data=kenton;var cases;class package;
histogram cases/normal;inset mean="Mean" std="St.Dev.";
title "Rough diagram corresponding to Figure 16.6";
title2 "Distribution based on cell-specific standard deviation, whereas 16.6 assumes MSE";run;
data critf;
file print;
critf=finv(1-.05, 3, 15);
put "Critical F is" critf=;
probf=1-probf(18.59,3,15);
put "p-level for observed F" probf=;
title "Critical F for test in GLM above. Refer to page 699";run;
proc means data=kenton;var ccases;
title "Weighted means for Kenton foods example. (data for incentive pay example aren't available";
title2 "refer to page 703";run;
proc print data=kenton;var package store e1 e2 e3;
title "Regression data for model 16.79. Refer to Table 16.4, page 707";run;
proc reg data=kenton;
model cases=e1 e2 e3/stb;
title "Fitted regression. Refer to Table 16.4 b and c";run;
data kenton;set kenton;
if package=1 then we1=1;else we1=0;
if package=2 then we2=1;else we2=0;
if package=3 then we3=1;else we3=0;
if package=4 then do;we1=-1;we2=-1;we3=-.8;end;
proc corr;var e1 e2 e3 we1 we2 we3;
title "weighted means. Notice how unequal weighting causes means of weighted contrasts to be zero";run;
proc reg;
model cases=we1 we2 we3/stb;
title "weighted contrasts. Note differences with table 16.4 as well as identities";
title2 "refer to page 710";run;
data kenton;set kenton;
if package=1 then d1=1;else d1=0;
if package=2 then d2=1;else d2=0;
if package=3 then d3=1;else d3=0;
if package=4 then d4=1;else d4=0;
drop _name_ _type_ _freq_;
proc reg;
model cases=d1 d2 d3 d4/noint;
equal: test d1=d2=d3=d4;
title "Cell means model. Refer to page 712";run;
/*It is, to my knowledge, not possible to do the simulations of the overall F described on page 715,
  but one can, using LSMEANS, generate analogous simulations for the pairwise comparisons of means
  which illustrates teh general idea*/
data example;
trt=1;y=3;output;
trt=1;y=7;output;
trt=2;y=8;output;
trt=2;y=10;output;
proc glm;class trt;
model y=trt;
lsmeans trt/adjust=simulate(nsamp=6 report) pdiff;
title "First simulation example. Refer to page 713";run;
data quality;
input method defect;
cards;
1 1.1
1 .5
1 -2.1
2 4.2
2 3.7
2 .8
3 3.2
3 2.8 
3 6.3
;
proc glm;class method;
model defect=method;
lsmeans method/pdiff adjust=simulate(report);
title "Second simulation example. Refer to page 715- note that this simulation is for the pairwise comparisons, not the overall F";run;
/*Power p. 716. It's easiest to use GLM power to do these study planning exercises, but following the math
  in the book shows you how it was done*/
data kentonm;
cases=12.5;package=1;cellwgt=5;output;
cases=13;package=2;cellwgt=5;output;
cases=18;package=3;cellwgt=4;output;
cases=21;package=4;cellwgt=5;output;
proc glmpower;
class package;
weight cellwgt;
model cases=package;

 contrast "1&2 vs 3&4" package 1 1 -1 -1;
 contrast "1 vs. 4" package 1 1 1 -3;
 contrast "All Pairwise" package 1 -1 0 0,
                         package      1 0 -1 0,
						 package 	   1 0 0 -1,
						 package 	   0 1 -1 0,
						 package 	   0 1 0 -1,
						 package 	   0 0 1 -1;
	
power
         stddev = 3.5
         alpha  = 0.025 .05
         ntotal = .
         power  = 0.8;
 plot x=power min=.5 max=.95;
 title "Some power plot Scenarios for Kenton Foods";
title2 "refer to 717-718";run;
/*I also demonstrate various power for contrasts in the Chapter17.sas data set if you're more interested in those
  comparisons*/
data Tire;
      input tire $ wear1 wear2 wear3 wear4 CellWgt;
      datalines;
tire1 24 27 25 26 1
tire2 25 25 26 23 1
tire3 25 25 25 28 1
tire4 25 25 26.5 23.5 1
;
   run;
/*pp. 760-761 give the ways to estimate the width of the Scheffe CI for planned comparisons.
  This is straightforward. I find it useful, though, to show power calculations for planned contrasts
  in GLM Power, as this is likely to be slightly more useful and you can adjust the experimentwise alpha
  in this if you wish. Four different scenarios are taken from Chapter 16, p. 721*/
proc glmpower;
   class tire;
      model wear1 wear2 wear3 wear4=tire;
      weight CellWgt;
      contrast "High vs. Low price" tire 1 -1 -1 1;
      contrast "Local vs. National" tire 1 1 -3 1;
	  contrast "All Pairwise" tire 1 -1 0 0,
	                          tire      1 0 -1 0,
							tire 	   1 0 0 -1,
							tire 	   0 1 -1 0,
							tire 	   0 1 0 -1,
							tire 	   0 0 1 -1;
	contrast "tire 1 vs. 4" tire 1 0 0 -1;
power
         stddev = 2 3
         alpha  = 0.025
         ntotal = .
         power  = 0.8;
title "Example power analyses of proposed snow tire study. Refer to discusssion on pages 719-722";run;