%let location=C:\Documents and Settings\woodph\My Documents\glm2\kutnersolutions\chapter22;
PROC IMPORT OUT= WORK.cracker 
            DATAFILE= "&location\chapter22.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH22TA02$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;run;
proc sort data=cracker out=cracker;by cond;
proc univariate;var sales priors;class cond;
histogram sales priors/normal;
title "histograms by condition of covariate and dependent variable";
title "Not shown in text, but useful information";run;
symbol1 repeat=3 v=dot i=rl;
axis1 order=20 to 50 by 10;
axis2 order=15 to 35 by 5;
proc format; 
value condfmt 1='1=Treatment 1'
2='2=Treatment 2'
3='3=Treatment 3';
data cracker;set cracker;
format cond condfmt.;
proc gplot;
plot sales*priors=cond/haxis=axis2 vaxis=axis1 noframe vminor=0 hminor=0;
title "Visual of equality of slopes for ancova";
title2 "Refer to Figure 22.5, page 927";run;
/*The analyses performing analysis of covariance are presented in two forms. The first, which is
presented in the text, shows how ancova may be performed by suitably coded variables and then conducted
using regression. This approach is useful on many software platforms which don't have a GLM module.
After that, the simpler code for doing analysis of covariance using GLM is shown, which is less work, because
the user doesn't have to code all those manifest variables*/
proc print;var cond subno sales priors x i1 i2 i1x i2x;
Title "Regression variables for Single-Factor Covariance Analysis- Refer to Table 22.2, page 927";run;
/*I decided to show the test for equality of slopes, since that is a usual first step. It'd discussed
after the analysis of covariance work, but you can flip back to the earlier page*/
proc reg;
model sales=priors i1 i2 i1x i2x/stb;
eqslopes: test i1x=i2x=0;
title "ancova via regression- test for equality of slopes";
title2 "Refer to page 933 & formula 22.24";run;
proc reg;
model sales=i1 i2 priors/stb covb;
title 'ancova via reg';
class: test i1=i2=0;
eqeff: test i1=i2;
title "Analysis of Covariance via coded regression variables";
title2 "Refer to Table 22.3, page 928";run;

proc glm;class cond;
model sales=priors|cond;
title "test of equality of slopes via GLM";
title2 "Refer to page 933 and formula 22.24";run;
proc glm;class cond;
model sales=priors cond;
lsmeans cond/adjust=scheffe cl pdiff;
output out=resids r=rsales p=psales;
means cond/tukey;
title 'Analysis of covariance via GLM';
title2 "Refer to Table 22.3, page 928";
title3 "Least suqares means and 95% confidence regions. Refer to pages 930-931";
title4 "Raw Means for the same data";run;
proc univariate;class cond;var rsales;
histogram rsales/normal;
title "Refer to Figure 22.6a, page 928";
proc univariate;var rsales;
qqplot rsales;
title "Normal Probability Plot of residuals. Refer to Figure 22.6b, page 928";run;
proc sort data=resids out=resids;by cond;
proc capability;var rsales;by cond;
histogram rsales;
title "Same histogram in Proc Capability";run;
proc capability;var rsales;
qqplot rsales;
title "Same qqplot in Proc Capability";run;
PROC IMPORT OUT= WORK.flower 
            DATAFILE= "&location\chapter22.xls" 
            DBMS=EXCEL REPLACE;
     SHEET="CH22TA06$"; 
     GETNAMES=YES;
     MIXED=NO;
     SCANTEXT=YES;
     USEDATE=YES;
     SCANTIME=YES;
RUN;
proc means;
title "Summary data for Salable Flowers";run;
proc means;class variety moisture;var y x;
title "Summary information by Class. Refer to Table 22.6, page 935";run;
proc sort data=flower out=flower;by variety moisture;
proc print;var y x;by variety moisture;
title "Data listing by class, Refer to Table 22.6, page 935";run;
proc format; 
value varmoisf
11="A1B1"
12="A1B2"
21="A2B1"
22="A2B2"
other="Miscoded";

data flower;set flower;
label varmois="Variety 1st digit Moisture 2nd digit";
varmois=variety*10+moisture;
format varmois varmoisf.;
symbol1 repeat=3 v=dot i=rl;
axis1 order=40 to 100 by 10;
axis2 order=0 to 20 by 5;
proc gplot;
plot y*x=varmois/vaxis=axis1 haxis=axis2 noframe vminor=0 hminor=0;
title "Scatterplot of Salable Flours. Refer to Figure 22.7, page 935";run;
proc univariate;class variety moisture;
histogram y x/normal;
title "univariate plots within cell";run;
proc glm;class variety moisture;
model y=x|variety|moisture;
title "test of slopes";run;
proc glm;class variety moisture;
model y=x variety|moisture/solution e;
title "Salable Flowers ancova. Refer to Table 22.7, refer to page 936";run;
/*Now, if you want the analysis presented in the text, you need to center the covariate to zero
and dummy code variety and moisture as follows:*/
proc standard data=flower out=flower mean=0;var x;
data flower;set flower;
if variety=2 then variety=-1;
if moisture=2 then moisture=-1;
proc glm;
model y=x variety moisture variety*moisture/solution e;
title "Ancova based on centered covariate and dummy coded information. Refer to Table 22.7a, page 936";run;
proc glm;class variety moisture;
model y=x variety moisture;
lsmeans variety moisture/pdiff adjust=bonf cl alpha=.025;
title "Least squares means comparisons- Bonferroni adjustment. Refer to page 937";
title2 "note the alpha of .025 necessary to get an experimentwise error rate of .05";run;
proc glm;class variety moisture;
model y=x variety moisture;
output out=flowero r=ry p=py;
proc univariate;class variety moisture;var ry;
histogram ry/normal;
title "Diagnostics referred to on p. 934";
proc univariate;var ry;
qqplot ry;
title "Normal Probability Plot of residuals. Referred to on page 934";run;

proc glm data=flower;*class variety moisture;
model y=x variety moisture;
lsmeans variety moisture/at x=0 out=slicedm;
estimate 'lsmean1' intercept 1 x 0 variety 1 moisture 1;
estimate 'lsmean2' intercept 1 x 0 variety 1 moisture -1;
estimate 'lsmean3' intercept 1 x 0 variety -1 moisture 1;
estimate 'lsmean4' intercept 1 x 0 variety -1 moisture -1;
title "using the at command in lsmeans for general covariate levels."
title2 "Estimate to calculate individual lsmeans. Refer to Figure 22.8, page 936";run;
proc print data=slicedm;
title "main effect lsmeans at covariate=0";run;
proc print data=flowerc;run;
proc glm data=flowerc;class variety moisture subno;
model y=x variety moisture subno;
title "Analogous randomized block design for flower data- see disscussion on page 938";
title "other analyses are possible";run;
data flowerb;set flower;
blockedx=x;
proc rank groups=5 data=flowerb out=flowerb;var blockedx;
proc glm;class blockedx moisture variety;
model y=blockedx moisture variety;
title "blocked main effects model. Refer to page 939";
proc glm;class blockedx moisture variety;
model y=blockedx moisture variety blockedx*moisture blockedx*variety;
title "blocking with interactions. Refer to page 939. Other analyses are possible";run;
data diff;set cracker;
salesd=sales-priors;
proc glm;class cond;
model salesd=cond;
title "difference score as an alternaive to covariance analysis. See pp. 939-940 (cracker data)";run;