***************************
* bleach.sas (lab 7) ***
***************************;
options pageno=1 nodate;
ods graphics on;
ods listing gpath="c:\temp";
ods pdf file="mypath\bleach.pdf";
data bleach;
input amtblch staintyp $ time;
cards;
3 ink 3600
3 ink 3920
3 ink 3340
3 ink 3173
3 jam 495
3 jam 236
3 jam 515
3 jam 573
3 jam 555
3 tomato 733
3 tomato 525
3 tomato 793
3 tomato 1026
3 tomato 510
5 ink 2029
5 ink 2271
5 ink 2156
5 ink 2493
5 ink 2805
5 jam 428
5 jam 432
5 jam 335
5 jam 288
5 tomato 880
5 tomato 759
5 tomato 1138
5 tomato 780
5 tomato 1625
7 ink 3660
7 ink 4105
7 ink 4545
7 ink 3569
7 ink 3342
7 jam 410
7 jam 225
7 jam 437
7 jam 350
7 jam 140
7 tomato 347
7 tomato 584
7 tomato 781
;
run;
/* Fit two-way anova model, compute treatment (joint) means and
marginal means, and output residuals for diagnostics*/
title 'Model 1: 2-way ANOVA model fit to data on original scale';
proc glm data=bleach plots=diagnostics;
class amtblch staintyp;
model time=amtblch staintyp amtblch*staintyp/ss1 ss2 ss3;
lsmeans amtblch staintyp amtblch*staintyp/out=outmns;
* output dataset in the above statement useful for producing profile plot without
using ods graphics (e.g., usinf proc plot or proc gplot) but not used here;
output out=resout p=preds student=instdres;
run;
title 'Model 1: 2-way ANOVA model fit to data on original scale';
title2 'Refit to get alternate profile plot';
proc glm data=bleach;
class staintyp amtblch;
model time=amtblch staintyp amtblch*staintyp/ss3;
ods exclude nobs classlevels; *cut down on output;
run;
title2;
/*
The code commented out below can be used to get profile plots without
using ods graphics, but we'll take advantgae of ods graphics here and
not run this code.
Produce profile plots of both types. In these low-resolution plots
you should connect the points that have the same plotting symbol to
see the profiles.
title 'Profile plots to depict interaction from model 1';
proc plot data=outmns;
where amtblch ne . and staintyp ne '';
plot lsmean*staintyp=amtblch;
plot lsmean*amtblch=staintyp;
run;
*/
/* Produce residuals vs fitteds plot to detect nonconstant variance.
Uncomment if not using ods graphics.
title '(Internally studentized) residuals vs fitteds from Model 1';
proc plot data=resout;
plot instdres*preds;
run;
*/
/* Compute treatment means and sd's so that we can regress log(sd) on
log(mean) to determine the appropriate transformation. The option NOPRINT
suppresses the output. The means and SDs are written to a new data set called mnsout */
proc means data=bleach noprint;
by amtblch staintyp;
var time;
output out=mnsout mean=timemean std=timesd;
run;
/* Compute log of treatmenat means and sd's */
data regdata;
set mnsout;
logtimmn=log(timemean);
logtimsd=log(timesd);
run;
ods graphics off;
/* Regress log(sd) on log(mean) */
title 'Regression of log(SD) on log(mean) to determine transformation';
proc reg data=regdata;
model logtimsd=logtimmn;
run;
ods graphics on;
/* The above regression gave alphahat=.53, so lambdahat=1-.53=.47
which is quite close to .5, so we choose a square root transformation.
Now compute square root of orginal response*/
data trans;
set bleach;
sqrttime=sqrt(time);
run;
/* and re-analyzed sqrt(time) rather than time, with the same two-way layout
model. Again, compute joint means and marginal means, and compute residuals
for model diagnostics. In addition, I've also included four contrasts.
The first and second of these examine the linearity of the mean response
(which is now the square root of time until the stain is removed) and
amount of bleach for stains of type 1 only.
The third and fourth contrasts examine the linearity of the mean response
and amount of bleach for stains of type 2 only. Finally, I've commented
out the fifth and sixth contrasts and left the coefficients of those
contrasts unspecified. Its your job to determine what those contrast
coefficients should be to examine the linearity of the mean response
versus amount of bleach for stains of type 3. */
title 'Model 2: 2-way ANOVA model after taking sqrt transformation';
proc glm data=trans plots=diagnostics;
class amtblch staintyp;
model sqrttime=amtblch staintyp amtblch*staintyp/ss1 ss2 ss3;
lsmeans amtblch staintyp amtblch*staintyp/out=outmns2;
output out=resout2 p=preds student=instdres;
contrast 'linear bleach, stain 1' amtblch -1 0 1
amtblch*staintyp -1 0 0 0 0 0 1 0 0;
contrast 'nonlinear bleach, stain 1' amtblch 1 -2 1
amtblch*staintyp 1 0 0 -2 0 0 1 0 0;
contrast 'linear bleach, stain 2' amtblch -1 0 1
amtblch*staintyp 0 -1 0 0 0 0 0 1 0;
contrast 'nonlinear bleach, stain 2' amtblch 1 -2 1
amtblch*staintyp 0 1 0 0 -2 0 0 1 0;
/*
contrast 'linear bleach, stain 3' amtblch ? ? ?
amtblch*staintyp ? ? ? ? ? ? ? ? ?;
contrast 'nonlinear bleach, stain 3' amtblch ? ? ?
amtblch*staintyp ? ? ? ? ? ? ? ? ?;
*/
run;
title 'Model 2: 2-way ANOVA model fit to data after sqrt transformation';
title2 'Refit to get alternate profile plot';
proc glm data=trans;
class staintyp amtblch;
model sqrttime=amtblch staintyp amtblch*staintyp/ss3;
ods exclude nobs classlevels; *cut down on output;
run;
title2;
/* If not using ODS graphics uncomment the code below to get the profile
plots using proc plot
* Produce profile plots ;
title 'Profile plots depicting interaction from Model 2';
proc plot data=outmns2;
where amtblch ne . and staintyp ne '';
plot lsmean*staintyp=amtblch;
plot lsmean*amtblch=staintyp;
run;
* Produce residual vs predicted plot (looks much better after transformation);
title '(Internally studentized) residuals vs fitteds from Model 2';
proc plot data=resout2;
plot instdres*preds;
run;
*/
ods graphics off;
* In the call to PROC GLM below I refit model 2 using a cell
means parameterization: y_ijk=mu_ij+e_ijk
Note that i refers to the levels of factor A and j refers to the levels
of factor B and the order of the mu_ij's in the CONTRAST statement is
mu_11, mu_12, mu_13, mu_21, mu_22,...
Which factor is treated as factor A is controlled by the order that the variables
appear on the CLASS statement. Below I have made staintype factor A by putting it first
on the CLASS statement.
This makes specifying the contrasts acrss levels of amtblch for a given staintyp
a bit easier. E.g., treating staintyp as factor A, the quadratic contrast coefficients
across levels of amtblch for staintypes 1, 2 and 3 are
1 -2 1 0 0 0 0 0 0,
0 0 0 1 -2 1 0 0 0, and
0 0 0 0 0 0 1 -2 1, respectively, whereas if we had treated amtblch as factor
A they would have been
1 0 0 -2 0 0 1 0 0,
0 1 0 0 -2 0 0 1 0,
0 0 1 0 0 -2 0 0 1, respectively, which are a bit less convenient to specify.;
title 'Model 2: 2-way ANOVA model fit to data after sqrt transformation';
title2 'Refit using cell-means parameterization to make contrasts easier';
proc glm data=trans;
class staintyp amtblch;
model sqrttime=amtblch*staintyp/ss3 noint;
contrast 'linear bleach, stain 1' amtblch*staintyp -1 0 1 0 0 0 0 0 0;
contrast 'nonlinear bleach, stain 1' amtblch*staintyp 1 -2 1 0 0 0 0 0 0;
contrast 'linear bleach, stain 2' amtblch*staintyp 0 0 0 -1 0 1 0 0 0;
contrast 'nonlinear bleach, stain 2' amtblch*staintyp 0 0 0 1 -2 1 0 0 0;
contrast 'linear bleach, stain 3' amtblch*staintyp 0 0 0 0 0 0 -1 0 1;
contrast 'nonlinear bleach, stain 3' amtblch*staintyp 0 0 0 0 0 0 1 -2 1;
ods exclude nobs classlevels; *cut down on output;
run;
title2;
ods pdf close;