SAS codes for some of the examples and Exercises taken from the book :
Design and Analysis of Experiments, Douglas C. Montgomery
fig2.6
proc iml;
n=20;
x=1:n;
ind=1:5;
y=j(5,n,0);
z=j(n,5,0);
do i=1 to 5;
y[i,]=pdf('chisquared',x,i);
end;
z=y`;
cn='c1':'c5';
create dat1 from z[colname=cn];
append from z;
close dat1;
create dat2 var{ind x};
append;
close dat2;
run;
data dat3;
merge dat1 dat2;
run;
proc print data=dat3;
run;
symbol1 i=join v=dot;
proc gplot data=dat3;
plot (c1-c5)*x/overlay;
run;
Table-2.1
/* Tension bond strength data
j modified unmodified
mortar(1) mortar(2)
------------------------------------- */
data tension;
input modstatus tensval @@;
datalines;
1 16.85 2 17.50 1 16.40 2 17.63
1 17.21 2 18.25 1 16.35 2 18.00
1 16.52 2 17.86 1 17.04 2 17.75
1 16.96 2 18.22 1 17.15 2 17.90
1 16.59 2 17.96 1 16.57 2 18.15
;
proc sort data=tension;
by modstatus;
run;
/* The following procedure will do lots of stuff for me.. if you want only
boxplot then you better use the PROC BOXPLOT statement.. I have shown in the
next code segment.. */
proc univariate data=tension plots;
by modstatus; /* whenever you are doing some analysis using "BY" statement
you need to sort the data first "BY" that variable*/
var tensval;
run;
/* for page-24*/
proc boxplot data=tension;
plot tensval*modstatus;
run;
/*for small output at the bottom of page-35 */
proc means data=tension mean var std n;
by modstatus;
var tensval;
run;
/* for doing a 2-sample t-test similar on page-38
The following code produces some outputs.. which gives the 95% CI for mean and
standard deviations of both the groups.. and it also performs the test for
equality of two variances.. (from the assumptions.. if the two samples have
almost equal variance then only we can do the pooled variance type of test..) */
proc ttest data=tension;
class modstatus;
var tensval;
run;
/* following code performs the test of equality of two variances LEVENE test.
you can also specify BARTLETT */
proc glm data=tension;
class modstatus;
model tensval=modstatus;
means modstatus / hovtest=levene ;
run;
/* for normal probability plots on page-39 */
proc univariate data=tension;
by modstatus;
probplot tensval/normal(mu=est sigma=est);
run;
/* table-2.5 data for the hardness testing experiment.*/
data hardness;
input sn tip1 tip2;
datalines;
1 7 6
2 3 3
3 3 5
4 4 3
5 8 8
6 3 2
7 2 4
8 9 9
9 5 4
10 4 5
;
proc ttest data=hardness;
paired tip1*tip2;
run;
/* if I need to modify the data to do some tests I can do the following to get
the suitable data...
data new1(keep=tip grp);
set hardness;
grp=1;
tip=tip1;
output;
grp=2;
tip=tip2;
output;
run;
the output will be a dataset with a single variable tip with a class variable
called "grp" which is 1 for TIP-1 and 2 for TIP-2. Suppose I have the dataset
in the following form :
sn tip tip_type
1 283 1
2 374 2
4 549 2 ......
and you want to do the pairwise comparison or something with tip1 & tip2
then you can use the following code to get the data on tip1 & tip2 seperately.
data tens1;
set tension;
if modstatus=1 ;tval1=tensval;
keep tval1;
data tens2;
set tension;
if modstatus=2 ; tval2=tensval;
keep tval2;
data tens;
merge tens1 tens2;
run;
*/
/* for doing the paired t-test shown on the page-49. The following program gives
the output which is similar to the output given in the table 2.6 */
proc ttest data=hardness;
paired tip1*tip2;
run;
/* Exercise-2.5 carbonated beverage data*/
data bever;
input days @@;
datalines;
108 124 124 106 115 138 163 159 134 139
;
/* the following code does everything except for the fact that
I needed to test Ho:=120 vs H1: > 120... but I am testing. Ho vs ~Ho */
proc ttest data=bever alpha=0.05 H0=120;
var days;
run;
/* the following will give you the required QQplot for testing normality.*/
proc univariate data=bever;
qqplot days/normal(mu=est sigma=est);
run;
/*here is the bad news.. there is no inbuilt function which does a one-sided
t-test for you.. you have to write a small code for doing this test.. eg.*/
PROC MEANS NOPRINT DATA=sound;
VAR spl;
OUTPUT OUT=temp MEAN=xbar STD=sd N=n
RUN;
DATA temp2;
SET temp;
KEEP xbar mu sd n t pgreater pless ptwoside;
INPUT mu;
t = (xbar-mu)/(sd/sqrt(n));
df = n - 1;
pgreater = 1 - probt(t,df);
pless = probt(t,df);
ptwoside = 2*MIN(1-ABS(probt(t,df)),ABS(probt(t,df)));
cards;
75 /* this is the mu(0) for which you are testing..*/
;
PROC PRINT;
RUN;
/* Exercise2.7 electronic instrument */
data elect;
input times @@;
datalines;
159 280 101 212 224 379 179 264
222 362 168 250 149 260 485 170
;
proc univariate data=elect;
probplot times/normal(mu=est sigma=est);
run;
/* test for normality of the data is not very perfect from this
normal probability plot */
Exercise-2.9
/* I did not enter the actual data given in the example but I am trying to
do the same problem with a different data used in exercise-2.7 */
data elect;
input grp hours @@;
datalines;
1 159 1 280 1 101 1 212 1 224 1 379 1 179 1 264
2 222 2 362 2 168 2 250 2 149 2 260 2 485 2 170
;
proc sort data=elect;
by grp;
run;
/*
proc univariate data=elect plots;
by grp;
var hours;
run;
*/
proc ttest data=elect;
class grp;
var hours;
run;
/* I am afraid that through SAS you cannot do the Z-test.. (or test for normal
distribution using some standard PROC) you will have to write your own code for
that test.. */
Exercise-2.11
/* for testing the equality of two variances.. you can use PROC TTEST */
proc ttest data=flares;
class type;
var burntime;
run;
/* the output will contain a test result called "FOLDED F" test.. and this is
nothing really big or mysterious..
Folded F = max( S1^2,S2^2)
-------------------
min( S1^2,S2^2)
this test, tests for the equality of the two variances.. vs not equal.*/
Exercise-2.12
/* C2F6 flow rate data...*/
data flowrate;
input rate etch @@;
etch2=etch*etch; // create new variable related to observed data
datalines;
1 2.7 1 4.6 1 2.6 1 3.0 1 3.2 1 3.8
2 4.6 2 3.4 2 2.9 2 3.5 2 4.1 2 5.1
;
proc sort data=flowrate;
by rate;
run;
/* following code is going to work for both the equality of avg etch uniformity
test and the equality of variability test.. at 95% level..part(a)-part(c)*/
proc ttest data=flowrate;
class rate;
var etch;
run;
/*following is the code for doing the boxplot on the two types.part(d)*/
proc boxplot data=flowrate;
plot etch*rate;
run;
Exercise-2.14
data flow;
input etch @@;
datalines;
* ** ** ******
* ********** *
;
/* this gives the answer for test of Ho:SD=1 ,this gives the 95% CI for
sigma^2... and the part of normality assumption is : if the data is normal then
only we can say that (n-1)S^2/sig^2 ~ Chisqr(n-1) */
proc ttest data=flow;
var etch;
run;
/* part(d) gives the normal probability plot */
proc univariate data=flow;
probplot etch/normal(mu=est sigma=est);
run;
Exercise-2.15
/* diameter of a ball bearing */
data ball;
input cali1 cali2;
datalines;
0.265 0.264
0.265 0.265
0.266 0.264
0.267 0.266
;
data new1(keep=cali diam);
set ball;
cali=1;
diam=cali1;
output;
cali=2;
diam=cali2;
output;
run;
proc print data=new1;
run;
proc sort data=new1;
by cali;
run;
proc ttest data=new1;
class cali;
var diam;
run;
Table-3.1
/* table 3.1 Tensile strength experiment */
data strength;
input wt obse1-obse5;
datalines;
15 07 12 14 19 07
20 07 17 18 25 10
25 15 12 18 22 11
30 11 18 19 19 15
35 09 18 19 23 11
;
data new1(keep=wt obse ID);
set strength;
wt=15;
obse=obse1;
output;
wt=20;
obse=obse2;
output;
wt=25;
obse=obse3;
output;
wt=30;
obse=obse4;
output;
wt=35;
obse=obse5;
output;
run;
/* now I am trying to do the Boxplot shown in the fig-3.1
Note : for using the PROC BOXPLOT you need to have the data sorted by the
clurstering variable.. */
proc sort data=new1;
by wt;
run;
/* this is the code for procing Box-plot
proc boxplot data=new1;
plot obse*wt;
run;
*/
/* the following code will give the scatter plot shown in fig-3.2
proc plot data=new1;
plot obse*wt;
run;
*/
/* here is the one-way analysis of the data shown above.. this is mentioned
as example 3.1 in the book at page-71. PROC ANOVA could also be used.. but I
think in PROC GLM I have more options and more freedom to do lots of stuff*/
proc glm data=new1;
class wt;
model obse=wt/ss3 solution clparm ;
output out=resout student=std_residual r=raw_residual p=predicted ;
proc print;
run;
/*if you want the parameter estimates.. you have to specify the
option "CLPARM" in the model statement as I have shown... */
/* MODEL ADEQUECY CHECKING
following are some plots which tests for the validity of normality
of the data which I wanted to fit assuming the normal distribution */
proc univariate data=resout;
var raw_residual;
probplot raw_residual/normal(mu=est sigma=est);
run;
proc plot data=resout;
plot raw_residual*time="*"; /* I do not have a time variable right now.*/
plot raw_residual*predicted="*";
run;
/* now we are interested in some contract estimation.... page-94
/* the following code produces the output shown in table-3.11 */
proc glm data=new1;
class wt;
model obse=wt;
contrast 'C1' wt 0 0 0 1 -1;
contrast 'C2' wt 1 0 1 -1 -1;
contrast 'C3' wt 1 0 -1 0 0;
contrast 'C4' wt 1 -4 1 1 1;
run;
/* Here is the code for producing output on the page-99 . */
proc glm data=new1;
class wt;
model obse=wt;
lsmeans wt /pdiff adjust=tukey tdiff; /* page-98 pdiff-> pairwise difference
tdiff-> T-statistics for all
these pairwaise difference */
run;
/* Example-3.10 page-103
in the above explained method you can write whatever method you like to get
adjusted for the pairwise contrast..
For comparing treatment means with some CONTROL say the 5th. just do the
following changes in the above code..
proc glm data=new1;
class wt;
model obse=wt;
lsmeans wt /pdiff=control('35') adjust=bon tdiff; /* page-98 pdiff-> pairwise difference
tdiff-> T-statistics for all
these pairwaise difference */
run;
/*
knowing all these methods.. producing output similar to page-105 is not a big
deal.. in fact most of the output has already been produced.. except for Press
statistics.. */
Table-3.7
/* peak discharge data */
data peak;
input meth obs @@;
datalines;
1 0.34 1 0.12 1 1.23 1 0.70 1 1.75 1 0.12
2 0.91 2 2.94 2 2.14 2 2.36 2 2.86 2 4.55
3 6.31 3 8.37 3 9.75 3 6.09 3 9.82 3 7.24
4 17.15 4 11.82 4 10.95 4 17.20 4 14.35 4 16.82
;
/* following is the day to produce mean median and all those summary statistics
of the dataset in the manner you wanted... */
proc means data=peak;
class meth;
var obs;
output out=meanout mean=avgobs var=varin p50=mediann;
run;
symbol1 i=join v=dot;
proc gplot data=meanout;
plot varin*meth;
run;
/*
proc boxplot data=peak;
plot obs*meth;
run;
*/
/* the following code will give you the result shown in the table-3.10 at the
bottom of page 85.*/
proc glm data=peak;
class meth;
model obs=meth / ss1 ss3 solution clparm;
lsmeans meth / pdiff adjust=bon; /* this gives the pairwise difference
testing for the method described by
"adjust" command.. */
means meth/tukey scheffe bon; /* this statement tries to classify the
variables in to small groups.. */
output out=resout p=predict r=resid;
run;
/* for testing the normal distribution of the data.. look at the residual plot
and this shows that the data is not normal... or better to say it needs some
transformation.. */
options pagesize=30;
proc plot data=resout;
plot resid*predict="*";
run;
/* so I am going to try the square root transformation in the data */
data newpeak(keep=meth obs obsnew);
set peak;
obsnew=sqrt(obs);
output;
run;
/* after transformation the result of GLM is as follows : */
proc glm data=newpeak;
class meth;
model obsnew=meth / ss1 ss3 solution clparm;
output out=resoutnew p=predictnew r=residnew;
run;
proc plot data=resoutnew;
plot residnew*predictnew="*";
run;
/* these results shows that doing transformation was much helpful.. if you look
at the the Error Sum Squares in the two models then you will able to see the
difference .. and also the residual plots clarifies this distinction..*/
Exercise-3.1
/* Tensile strength data Exercise-3.1*/
data tech;
input tec stren @@;
datalines;
1 3129 2 3200 3 2800 4 2600
1 3000 2 3300 3 2900 4 2700
1 2865 2 2975 3 2985 4 2600
1 2890 2 3150 3 3050 4 2765
;
/* part(a) the result of this analysis shows significant effect of mixing
techiniques on the tensile strength.. with p-value=0.0005*/
proc glm data=tech;
class tec;
model stren=tec/p;
run;
/* part(b) the first part of the code gives the average measurement of the
tensile strength and then second part of the code plots those means...*/
proc means data=tech;
class tec;
var stren;
output out=meanout mean=avgstr var=varin;
run;
symbol1 i=join v=dot;
proc gplot data=meanout;
plot avgstr*tec;
run;
/* part(c) Using Fisher's LSD method for pairwise comparison of means. I guess
that Least significant difference method is same as Bonferroni.*/
proc glm data=tech;
class tec;
model stren=tec;
lsmeans tec / pdiff adjust=bon tdiff;
means tec / lsd;
output out=resout r=residual p=predicted;
run;
/* note that in case of equal sample size, option "t" is same as "LSD" */
/* part(d)-normal probability plot of the raw_residual */
proc univariate data=resout;
var residual;
probplot residual/normal(mu=est sigma=est);
run;
/*looking at the probability plot I can say that the data seemed to be more or
less normal.. and residual vs predicted plot is also not revealing any evidence
against normality.*/
proc gplot data=resout;
plot residual*predicted="*";
run;
Exercise-3.3
/* for Exercise-3.3 we have to calculate the CI for means.. it can be done in
two ways.. one way to do that is to use LSMEANS statement and the other way is
to use the MEANS statement.. both with proper options are shown below. If you
want (1-alpha) CI , then you should specify the(alpha=0.02,say) in the option*/
proc glm data=tech;
class tec;
model stren=tec/p;
means tec/ t clm;
lsmeans tec / cl;
run;
Exercise-3.6
/* manufacturer of television sets */
data tv;
input coat cond @@;
datalines;
1 143 2 152 3 134 4 129
1 141 2 149 3 136 4 127
1 150 2 137 3 132 4 132
1 146 2 143 3 127 4 129
;
/*since the data is not sorted ,before doing further analysis we should
sort the data in the appropriate manner. */
proc sort data=tv;
by coat;
run;
/* part(a): The model statement gives P-value 0.0003 and therefore the effects
are significant. part(b):this gives the parameter estimates. */
proc glm data=tv;
class coat;
model cond=coat/ss1 ss3 solution clparm;
run;
/* for getting the 99% CI for testing the difference between coating type1-type4 is
shown by the following statement.("estimate") */
proc glm data=tv;
class coat;
model cond=coat/alpha=0.01;
estimate 'coat 1-4' coat 1 0 0 -1;
run;
/* part(d) : LSD estimate for all the pairs */
proc glm data=tv;
class coat;
model cond=coat;
means coat/lsd;
run;
/*part(e) ; basically I have to produce profile plot */
proc means data=tv;
class coat;
var cond;
output out=meanout mean=avg var=varin;
run;
symbol1 i=join v=dot;
proc gplot data=meanout;
plot avg*coat;
run;
/* for doing Exercise-3.7 we have to get hold of the residuals first
and then we can plot those to have a look for Model adequecy test.*/
proc glm data=tv;
class coat;
model cond=coat;
output out=resout r=residual p=predicted;
proc print;
run;
proc gplot data=resout;
plot residual*predicted="*";
run;
Exercise-3.10
/* different types of circuits */
data circuit1;
input ctype @@;
datalines;
1 2 3
;
data circuit2;
input rtime1-rtime3;
datalines;
9 20 6
12 21 5
10 23 8
8 17 16
15 30 7
;
data circuit;
merge circuit1 circuit2;
run;
data cir(keep=rtime ctype);
set circuit;
ctype=1;
rtime=rtime1;
output;
ctype=2;
rtime=rtime2;
output;
ctype=3;
rtime=rtime3;
output;
run;
/* before start working with the dataset I should sort it.*/
proc sort data=cir;
by ctype;
run;
/* part(a) , usual model checkup.. for part(b) all you have to do
is the addition of line "LSMEANS" below the model statement. */
proc glm data=cir;
class ctype;
model rtime=ctype;
lsmeans ctype/pdiff adjust=tukey cl;
run;
/*part(c) : here you have to make the profile plot... steps are
(1) first get the means(classwise) and then
(2) plot those means wrt ctype (the class variable) */
proc means data=cir;
class ctype;
var rtime;
output out=meanout mean=avg var=varin;
run;
symbol1 i=join v=dot c=red;
symbol2 i=join v=dot c=blue;
proc gplot data=meanout;
plot avg*ctype varin*ctype/overlay ;
run;
quit;
/* part(e) can be answered from the plot above done.. just look at the minimum
avg and that ctype should be taken..
For the residual analysis.. all you hae to do is the similar stuff..(following)*/
proc glm data=cir;
class ctype;
model rtime=ctype;
output out=resout r=residual p=predicted;
run;
proc gplot data=resout;
plot residual*predicted="*";
run;
Exercise-3.16
/*Exercise-3.16, Effectiveness of insulation material.*/
data failtime;
input mat time @@;
datalines;
1 110 1 157 1 194 1 178
2 1 2 2 2 4 2 18
3 880 3 1256 3 5276 3 4355
4 495 4 7040 4 5307 4 10050
5 7 5 5 5 29 5 2
;
/* before starting the analysis let me just sort the data "by" class variable*/
proc sort data=failtime;
by mat;
run;
/*part(a) standard question !! where all the five materials have the same
effect on failure time or not */
proc glm data=failtime;
class mat;
model time=mat;
output out=resout r=residual p=predicted;
run;
/* Assuming that the ANOVA assumptions are satisfied.. the p-value says different materials
have different effect.
But for the sake of completeness we should check for the
(1) equality of variance condition
(2) normality of residuals(probability plot)[obtained from above PROC statement]
(3) residual vs predicted plot... */
proc glm data=failtime;
class mat;
model time=mat;
means mat /hovtest=bartlett;
run;
proc univariate data=resout;
var residual;
probplot residual/normal(mu=est sigma=est);
run;
proc plot data=resout;
plot residual*predicted="*";
run;
/* as a result of the above analysis part(b) we saw that the res*pred plot and
the bartlett's both suggested that doing ANOVA with the raw data doesn't
make much sense...
So looking at this some data transformation in needed. */
data new1;
set failtime;
logtime=log(time);
output;
run;
/* now doing the similar analysis with the transformed data gives a better result.
and it suggests that we can stick to this.. here Bartlett's test of equality of
variance also supports the use of ANOVA */
proc glm data=new1;
class mat;
model logtime=mat;
means mat /hovtest=bartlett;
output out=resout1 r=resid1 p=predic1;
run;
proc univariate data=resout1;
var resid1;
probplot resid1/normal(mu=est sigma=est);
run;
proc plot data=resout1;
plot resid1*predic1="*";
run;
/*if you look at the Probplot of the residual then, you will be able to
see a great difference in the two analysis.. one more thing that in the
untransferred data the residual plot shows the range on (-6000+) whereas
after transformation it is (-2+)..*/
Table-4.1
/* Hardness testing Experiment */
data hardness;
input tiptype couptype hardns @@;
datalines;
1 1 9.3 1 2 9.4 1 3 9.6 1 4 10.0
3 1 9.2 3 2 9.4 3 3 9.5 3 4 9.7
2 1 9.4 2 2 9.3 2 3 9.8 2 4 9.9
4 1 9.7 4 2 9.6 4 3 10.0 4 4 10.2
;
/*before using this data, it has to be sorted */
proc sort data=hardness;
by tiptype couptype;
run;
/*Example-3.1 page 131 (this is the skeleton)*/
proc glm data=hardness noprint;
class tiptype couptype;
model hardns= tiptype couptype;
run;
/* for the result shown in page-134 fig-4.2 following is the code.*/
proc glm data=hardness noprint;
class tiptype couptype;
model hardns = tiptype couptype;
lsmeans tiptype / pdiff adjust=t tdiff cl; /* here option "t" is same as "LSD" */
output out=resout p=predicted r=residual h=leverage;
run;
/*now to produce a profile plot similar to shown on page-135 */
proc means data=hardness noprint;
class tiptype;
var hardns;
output out=meanout mean=avg var=varin;
run;
/*
symbol1 i=join c=red v=dot;
symbol2 i=join c=green v=dot;
proc gplot data=meanout;
plot avg*tiptype varin*tiptype /overlay;
plot residual*predicted;
run;
*/
/* this plot shows that the variance of the four tip type are almost the same
and therefore ANOVA assumption for equality of variance is Ok. Idea of anova
is that the variance within blocks should be small too. */
proc means data=hardness noprint;
class couptype;
var hardns;
output out=minout mean=cavg var=cvarin;
run;
/*
symbol1 i=join v="m";
symbol2 i=join v="v";
proc gplot data=minout;
plot cavg*couptype cvarin*couptype /overlay;
run;
*/
/* but this is not enough, you should check for the normality
of the residuals. */
proc univariate data=resout noprint;;
var residual;
probplot residual/normal(mu=est sigma=est);
run;
/* the above residual plot doesn't show any evidence against normality,
also the following shows no evidence against normality.*/
proc plot data=resout;
plot residual*(couptype tiptype predicted)="*"/vref=0;
run;
Table-4.9
/* Rocket propellant problem */
data rocket;
input batch oper trt $ obsn @@;
datalines;
1 1 a 24 3 1 c 18 5 1 e 22
1 2 b 20 3 2 d 38 5 2 a 30
1 3 c 19 3 3 e 26 5 3 b 20
1 4 d 24 3 4 a 27 5 4 c 29
1 5 e 24 3 5 b 21 5 5 d 31
2 1 b 17 4 1 d 26
2 2 c 24 4 2 e 31
2 3 d 30 4 3 a 26
2 4 e 27 4 4 b 23
2 5 a 36 4 5 c 22
;
proc sort data=rocket;
by batch oper trt;
run;
proc glm data=rocket noprint;
class batch oper trt;
model obsn = trt oper batch;
output out=resout r=residual p=predicted;
run;
proc means data=rocket noprint;
class trt;
var obsn;
output out=meanout mean=trtavg var=trtvar;
run;
symbol1 i=join v=dot;
symbol2 i=join v=star;
proc gplot data=meanout;
plot trtavg*trt trtvar*trt/overlay;
run;
proc plot data=resout;
plot residual*(predicted trt batch oper)="*"/vref=0;
run;
Table-4.22
/*table-4.22 BIBD example for catalyst experiment.*/
data catal;
input trt block recntime @@;
datalines;
1 1 73 3 1 73
1 2 74 3 2 75
1 4 71 3 3 68
2 2 75 4 1 75
2 3 67 4 3 72
2 4 72 4 4 75
;
/*before doing anything we should sort the data */
proc sort data=catal;
by trt block;
run;
proc glm data=catal noprint;
class trt block;
model recntime = trt block;
run;
/* the result of the above procedure is not exactly the stuff we are looking for.
look at the type-3 sum of squares.. this gives me the SS corresponding to every
factor adjusted for rest. but this is not what I want !! Our anova table should
look like :
sources of variation
---------------------
SS trt(adjusted for block)
SS block
SS Error
SS total corrected
----------------------- so look at the model SS in the overall ANOVA and then take the
SS trt(adjusted) from the type-3 SS and then subtract it to get the SS BLOCK.. eg in this
setup SS trt(adjust) = 22.75, SS block =55.0
OR, do the analysis as "model recntime=block trt" then take the type-1 SS for blocks
and tyoe-3 SS for trt..this will give the required result.. */
proc means data=catal noprint;
class trt;
var recntime;
output out= meanout mean=avgtime var=vartime;
run;
/* the following code gives an idea of the average profile of the response variable
with respect to the treatment..
footnote1;
footnote2 h=1 j=l 'Red-> Average time, Blue-> Variance of time';
symbol1 i=join c=red v=dot;
symbol2 i=join c=blue v=dot;
proc gplot data=meanout;
plot avgtime*trt vartime*trt/overlay;
run;
*/
/*for producing the output similar on page-160, ie table-4.26 you can use the following:*/
proc glm data=catal;
class trt block;
model recntime = trt block;
lsmeans trt / pdiff tdiff adjust=tukey cl;
means trt/ tukey;
run;
Exercise-4.1
/* tensile strength Exercixe-4.1*/
data tensile;
input chem bolt stren @@;
datalines;
1 1 73 3 1 75
1 2 68 3 2 68
1 3 74 3 3 78
1 4 71 3 4 73
1 5 67 3 5 68
2 1 73 4 1 73
2 2 67 4 2 71
2 3 75 4 3 75
2 4 74 4 4 75
2 5 70 4 5 69
;
proc sort data=tensile;
by chem bolt;
run;
/*I am trying to look at the data to get a rough idea about the interaction
between trt and blocking variable. For that you need to do this kind of plot. */
symbol i=join v=dot;
proc gplot data=tensile;
plot stren*chem=bolt;
run;
/* now comes the actual analysis of the data given */
proc glm data=tensile;
class bolt chem;
model stren=bolt chem;
lsmeans chem/pdiff tdiff adjust=tukey cl;
means chem / t clm;
output out=resout r=residual p=predict;
run;
/*
proc plot data=resout;
plot residual*(predict chem)="*"/vref=0;
run;
proc univariate data=resout;
var residual;
probplot residual/normal(mu=est sigma=est);
run;
*/
Exercise4.5
/* Nozzle design vs Shape factor Exercise-4.5*/
data nozzle;
input nozdes vel shape @@;
datalines;
1 11.73 0.78 1 14.37 0.80 1 16.59 0.81 1 20.43 0.75 1 23.46 0.77 1 28.74 0.78
2 11.73 0.85 2 14.37 0.85 2 16.59 0.92 2 20.43 0.86 2 23.46 0.81 2 28.74 0.83
3 11.73 0.93 3 14.37 0.92 3 16.59 0.95 3 20.43 0.89 3 23.46 0.89 3 28.74 0.83
4 11.73 1.14 4 14.37 0.97 4 16.59 0.98 4 20.43 0.88 4 23.46 0.86 4 28.74 0.83
5 11.73 0.97 5 14.37 0.86 5 16.59 0.78 5 20.43 0.76 5 23.46 0.76 5 28.74 0.75
;
proc sort data=nozzle;
by nozdes vel;
run;
/*part(a) following is the code which performs the general analysis and
generates a scatter plot of the data */
proc glm data=nozzle noprint;
class nozdes vel;
model shape=nozdes vel;
lsmeans nozdes/pdiff tdiff adjust=tukey cl; /*this gives the pairwise contrast estimate
and its CI with t-test for equality of the pair.*/
means nozdes / tukey alpha=0.01; /*this will give some idea about the grouping of response
variable, more clear explanation can be obtained from
the package "Guided data analysis" of SAS.*/
output out=resout r=resid p=predict; /*this will give me residual which
will be used in the analysis of part(b) */
run;
/* gives the scatter plot of the data with clearly shown the different nozzles. */
symbol i=join v=dot h=1;
proc gplot data=nozzle;
plot shape*vel=nozdes;
run;
/* part(b): residual plots.. */
proc univariate data=resout noprint;
probplot resid/normal(mu=est sigma=est);
run;
proc plot data=resout;
plot resid*(predict nozdes)="*"/vref=0;
run;
/*note :residual plots doesnot reveal any significant violation of normality.
part(c): plot of avg shape for each nozzle design is as follows */
proc means data=nozzle noprint;
by nozdes;
var shape;
output out=meanout mean=avgshape var=varshape;
run;
symbol i=join v=dot h=1;
proc gplot data=meanout;
plot avgshape*nozdes;
run;
Exercise4.27
/* Exercise-4.27 mileage performance for Gasoline Additives */
data gas;
input additive car mileage @@;
datalines;
1 2 17 1 3 14 1 4 13 1 5 12
2 1 14 2 2 14 2 4 13 2 5 10
3 1 12 3 3 13 3 4 12 3 5 9
4 1 13 4 2 11 4 3 11 4 4 12
5 1 11 5 2 12 5 3 10 5 5 8
;
proc sort data=gas;
by additive car;
run;
/* following is the basic analysis of the data.. over all model fit and
the sum of squares.. note that take the block SS from the type-1 SS and
trt SS from type-3 SS.. bcos you need trt(adjust) and block(unadjust)*/
proc glm data=gas noprint;
class car additive ;
model mileage=car additive ;
lsmeans additive/pdiff tdiff adjust=tukey cl;
estimate 't1-t2' additive 1 -1 0 0 0;
estimate 't1+t2-2t3' additive 1 1 -2 0 0;
estimate 't1+t2+t3-3t4' additive 1 1 1 -3 0;
estimate 't1+t2+t3+t4-4t5' additive 1 1 1 1 -4;
contrast 't1-t2' additive 1 -1 0 0 0;
contrast 't1+t2-2t3' additive 1 1 -2 0 0;
contrast 't1+t2+t3-3t4' additive 1 1 1 -3 0;
contrast 't1+t2+t3+t4-4t5' additive 1 1 1 1 -4;
output out=resout r=resid p=predict;
run;
/* above code also provides the estimate and SS of a set of orthogonal contrasts.
which was demanded in Exercise-4.28
But I think we should not stop over here, we should have a look at the profile
plots and residual plots.. to check the normality conditions. */
proc means data=gas;
class additive;
var mileage;
output out=meanout mean=avgmile var=varmile;
run;
symbol1 i=join h=1 v=dot; symbol2 i=join h=1 v=dot;
proc gplot data=meanout;
plot avgmile*additive varmile*additive/overlay;
run;
proc univariate data=resout noprint;
probplot resid/normal(mu=est sigma=est);
run;
proc plot data=resout;
plot resid*(predict additive)="*"/vref=0;
run;
Table-5.1
/* Example-5.3.1 Life data of the Battery design example (3^2 factorial expt) */
data life;
input mtype temp time @@;
datalines;
1 15 130 1 15 155 1 15 74 1 15 180 1 70 34 1 70 40 1 70 80 1 70 75
2 15 150 2 15 188 2 15 159 2 15 126 2 70 136 2 70 122 2 70 106 2 70 115
3 15 138 3 15 110 3 15 168 3 15 160 3 70 174 3 70 120 3 70 150 3 70 139
1 125 20 1 125 70 1 125 82 1 125 58 2 125 25 2 125 70 2 125 58 2 125 45
3 125 96 3 125 104 3 125 82 3 125 60
;
proc sort data=life;
by mtype temp;
run;
/* in the following code, look for the type-3 analysis of interaction !! if
this shows that the interaction is zero then you should remove interaction and
then fit the model again. but here interaction is not zero so you have to keep the
interaction in the model.*/
proc glm data=life noprint;
class mtype temp;
model time = mtype|temp/solution clparm p;
output out=resout r=residual p=predicted;
run;
/* now I am trying to do the profile plots to check for the interaction
proc means data=life noprint;
class mtype temp;
var time;
output out=meanout mean=timeavg var=timevar;
run;
symbol i=join v=dot h=1;
proc gplot data=meanout;
plot timeavg*temp=mtype;
run;
above code is for producing the plot similar to page-182 figure-5.9 This also
suggests that the interaction is present */
/* Now we can go for Model adequacy checking.. get the residual from the
above "proc glm" code , some residuals plots for adequacy of model*/
proc univariate data=resout;
var residual;
probplot residual/normal(mu=est sigma=est);
run;
proc plot data=resout;
plot residual*(predicted mtype temp)="*"/vref=0;
run;
Table-5.13
/* Example-5.3, The soft-drink bottling problem (eg of 3*2*2 factorial expt)*/
data dev;
input carb press speed height @@;
datalines;
10 25 200 -3 10 25 200 -1 10 25 250 -1 10 25 250 0
10 30 200 -1 10 30 200 0 10 30 250 1 10 30 250 1
12 25 200 0 12 25 200 1 12 25 250 2 12 25 250 1
12 30 200 2 12 30 200 3 12 30 250 6 12 30 250 5
14 25 200 5 14 25 200 4 14 25 250 7 14 25 250 6
14 30 200 7 14 30 200 9 14 30 250 10 14 30 250 11
;
proc sort data=dev;
by carb press speed;
run;
proc glm data=dev ;
class carb press speed;
model height=carb|press|speed;
output out=resout1 r=resid1 p=predict1;
run;
/*looking at the ANOVA we can say that carb*speed press*speed and the
3 factor inter are close to zero. If you wanna get the parameter estimates
you can use "solution clparm" in the model statement as an option.
So I will fit the following model : */
proc glm data=dev;
class carb press speed;
model height= carb|press speed;
output out=resout2 r=resid2 p=predict2;
run;
/*now if you wanna compare the two models, look at the adjusted R^2 value then
you can see that this is a better model, you can also look at MSE. I would like
to look at the two residual plots and see the pattern.*/
proc univariate data=resout1 noprint;
probplot resid1/normal(mu=est sigma=est);
run;
proc plot data=resout1;
plot resid1*predict1="*"/vref=0;
run;
proc univariate data=resout2 noprint;
probplot resid2/normal(mu=est sigma=est);
run;
proc plot data=resout2;
plot resid2*predict2="*"/vref=0;
run;
/* if you look at the two probability plots of residuals then forsure you will say
that the second model is better than the first one.. */
/*
proc means data=dev noprint;
class carb press speed;
var height;
output out=meanout mean=avght;
run;
proc sort data=meanout;
by speed;
run;
symbol i=join h=1 v=dot;
proc gplot data=meanout;
by speed;
plot avght*carb=press;
run;
*/
/* these three plots done above shows that there is negligible
interaction between speed and pressure.*/
Example-5.4
/* Example-5.4 Fitting response curves.. we will work with the same dataset
as in table-5.1(page-176) here : mtype-> qualitative & temp -> quantitative*/
data life;
input mtype temp time @@;
datalines;
1 15 130 1 15 155 1 15 74 1 15 180 1 70 34 1 70 40 1 70 80 1 70 75
2 15 150 2 15 188 2 15 159 2 15 126 2 70 136 2 70 122 2 70 106 2 70 115
3 15 138 3 15 110 3 15 168 3 15 160 3 70 174 3 70 120 3 70 150 3 70 139
1 125 20 1 125 70 1 125 82 1 125 58 2 125 25 2 125 70 2 125 58 2 125 45
3 125 96 3 125 104 3 125 82 3 125 60
;
proc sort data=life;
by mtype temp;
run;
proc glm data=life;
class mtype;
model time=temp|mtype temp*temp|mtype/solution clparm;
run;
Example-5.5
/* Example-5.5 (data shown on table-5.16) data for tool life expt.*/
data tool;
input angle speed life @@;
datalines;
15 125 -2 15 125 -1 15 150 -3 15 150 0 15 175 2 15 175 3
20 125 0 20 125 2 20 150 1 20 150 3 20 175 4 20 175 6
25 125 -1 25 125 0 25 150 5 25 150 6 25 175 0 25 175 -1
;
proc sort data=tool;
by angle speed;
run;
/* both the variables angle and speed are quantitative.. */
proc glm data=tool noprint;
model life= angle|angle|speed|speed/ss1;
output out=resout p=predict;
run;
/* I dont know how to get the estimate of the main effects given in the book
I mean SS due to A , B, A^2, B^2 ... but my code gives the rest of the table.*/
/*now I will try to produce contour plot shown in figure-5.19 on page-206
and figure-5.20 surface plot on same page..
step-1 : first make a grid of speed and angle
step-2 : then plot those points on the grid; */
proc g3grid data=resout out=gridresout;
grid speed*angle=predict /naxis1=41 naxis2=41 spline;
run;
proc gcontour data=gridresout;
plot speed*angle=predict;
run;
proc g3d data=gridresout;
plot speed*angle=predict/grid zmin=-4;
run;
Example-5.6
/* example-5.6, given on page 207, data in table-5.19 (Intensity level at target detection)*/
data intens1;
input level @@;
datalines;
90 86 96 84 100 92 92 81 102 87 106 90 105 97 96 80
114 93 112 91 108 95 98 83
;
data intes2;
do i=1 to 3;
do j=1 to 4;
do k=1 to 2;
clutter=i;
operator=j;
filter=k;
output;
end;
end;
end;
keep filter operator clutter;
data intensity;
merge intes2 intens1;
run;
/* this is an example of blocking in factorial experiment setup, code
for the output shown in table-5.20 */
proc glm data=intensity;
class filter operator clutter;
model level= filter|clutter operator;
run;
/*looking at the anova result interaction is almost zero... lets look at
the profile plot for interaction.*/
proc means data=intensity noprint;
class filter clutter;
var level;
output out=meanout mean=avglevel;
run;
symbol i=join v=dot;
proc gplot data=meanout;
plot avglevel*clutter=filter;
run;
/* If you look at the profile plot then also you can see that though the two
lines are not very parallel there is very less evidence against interaction.*/
Exercise-5.1
/*Exercise-5.1 .... yield of some chemical process.*/
data chem2;
input yield @@;
datalines;
90.4 90.7 90.2 90.2 90.6 90.4 90.1 90.5 89.9
90.3 90.6 90.1 90.5 90.8 90.4 90.7 90.9 90.1
;
data chem1;
do i=150,160,170;
do k=1 to 2;
do j=200,215,230;
temp=i;
rep=k;
press=j;
output;
end;
end;
end;
keep temp rep press;
data chem;
merge chem1 chem2;
run;
/* part(a) some general analysis of the data */
proc glm data=chem noprint;
class temp rep press;
model yield= temp|press rep/solution clparm;
output out=resout1 r=residual1 p=predicted1;
run;
/* it shows that interaction is almost zero... lets verify this from the profile plot.
following is the code for doing interaction plot */
proc means data=chem;
class press temp;
var yield;
output out=meanout mean=avgyield;
run;
symbol i=join v=dot h=1;
proc gplot data=meanout;
plot avgyield*temp=press;
run;
/*since the interaction is almost zero we will fit a model without interaction*/
proc glm data=chem noprint;
class temp rep press;
model yield= temp press rep/solution clparm;
output out=resout2 r=residual2 p=predicted2;
run;
/* after fitting this model, we can see that the MSE has been reduced a bit.
it may be the case that R-adjusted has also increased. */
/* part(b), following is the code for residual plots.. */
proc univariate data=resout1 noprint;
var residual1;
probplot residual1/normal(mu=est sigma=est);
run;
proc plot data=resout1;
plot residual1*(predicted1 temp press)="*"/vref=0;
run;
proc univariate data=resout2 noprint;
var residual2;
probplot residual2/normal(mu=est sigma=est);
run;
proc plot data=resout2;
plot residual2*(predicted2 temp press)="*"/vref=0;
run;
/* for answering part(c), we will have to look at the profile plots,
if I am interested in the maximum yield then we would go for
temp =170 and press=215 */
Exercise-5.8
/* Exercise-5.8 light output data of an oscillation tube.*/
data tube2;
input light @@;
datalines;
580 1090 1392 568 1087 1380 570 1085 1386
550 1070 1328 530 1035 1312 579 1000 1299
546 1045 867 575 1053 904 599 1066 889
;
data tube1;
do i =1 to 3;
do k=1 to 3;
do j=100,125,150;
ctype=i;rep=k;temp=j;
output;
end;
end;
end;
keep ctype rep temp;
data tube;
merge tube1 tube2;
run;
/* part(a) says to test if there is any interaction or not ?
so for that we will do the test for interaction=0 and the
profile plot to see the visual evidence for interaction.*/
proc means data=tube;
class ctype temp;
var light;
output out=meanout mean=avglight;
run;
symbol i=join v=dot h=1;
proc gplot data=meanout;
plot avglight*ctype=temp;
run;
proc glm data=tube;
class ctype rep temp;
model light=ctype|temp rep/solution clparm;
output out=resout1 r=resid1 p=predict1;
run;
/* NOTE : both the profile plot and the test shows that there is very
strong evidence against interaction=0. so the model we have fit is the right model.
so if you go for a reduced model without interaction then I guess MSE would increase. */
proc glm data=tube;
class ctype rep temp;
model light=ctype temp rep/solution clparm;
output out=resout2 r=resid2 p=predict2;
run;
/* in fact the MSE has increased a lot... so ofcourse the reduced model is worse.
Now lets see some residual plots.. to check for the adequacy of the model.*/
proc univariate data=resout1 noprint;
probplot resid1/normal(mu=est sigma=est);
run;
proc univariate data=resout2 noprint;
probplot resid2/normal(mu=est sigma=est);
run;
/* the above conclusion about the bigger model is justified by the normal
probability plot and residual vs predicted plots for both the models..*/
proc plot data=resout1;
plot resid1*predict1="*"/vref=0;
run;
proc plot data=resout2;
plot resid2*predict2="*"/vref=0;
run;
Exercise-5.16
/* Exercise-5.16 page-215, strength of paper data */
data paper2;
input stren @@;
datalines;
196.6 197.7 199.8 198.4 199.6 200.6 196.0 196.0 199.4 198.6 200.4 200.9
198.5 196.0 198.4 197.5 198.7 199.6 197.2 196.9 197.6 198.1 198.0 199.0
197.5 195.6 197.4 197.6 197.0 198.5 196.6 196.2 198.1 198.4 197.8 199.8
;
data paper1;
do i=2,4,8;
do k=1,2;
do j=3,4;
do l=400,500,650;
press=l; cooktime=j;rep=k;conc=i;
output;
end;
end;
end;
end;
keep press cooktime rep conc;
data paper;
merge paper1 paper2;
run;
/*part(a) start with the preliminary ANOVA of the model.*/
proc glm data=paper;
class cooktime rep press conc;
model stren= cooktime|press|conc rep;
output out=resout1 r=residual1 p=predicted1;
run;
proc glm data=paper;
class cooktime rep press conc;
model stren= cooktime press|conc rep;
output out=resout2 r=residual2 p=predicted2;
run;
/* look at the profile plots to have an idea about the interaction among factors..*/
proc means data=paper noprint;
class cooktime press conc;
var stren;
output out=meanout mean=avgstren;
run;
proc sort data=meanout;
by cooktime;
run;
symbol i=join v=dot h=1;
proc gplot data=meanout;
by cooktime;
plot avgstren*press=conc;
run;
/* if you look at the profile plot then for the last part(c) of the question
the answer would be cooktime=4 conc=2 & press=650, this will produce the
maximum paper strength.*/
/*Now if we look at the residual plots for the two models*/
data resout;
merge resout1 resout2;
run;
proc plot data=resout;
plot residual1*predicted1="*"/vref=0 vaxis=-2 to 2;
plot residual2*predicted2="*"/vref=0 vaxis=-2 to 2;
run;
proc univariate data=resout noprint;
probplot residual1/normal(mu=est sigma=est);
probplot residual2/normal(mu=est sigma=est);
run;
/* both the plots do not show any evidence against normality,
So I would suggest to consider the model with all main effect
and all 2 factor interactions.. */
Table-6.4
/* this data is a part of the data shown in table-5.13 , here we
have carbonation at two levels : 10 & 12*/
data dev1;
input rep carb press speed height @@;
datalines;
1 10 25 200 -3 1 10 25 200 -1 1 10 25 250 -1 1 10 25 250 0
2 10 30 200 -1 2 10 30 200 0 2 10 30 250 1 2 10 30 250 1
1 12 25 200 0 1 12 25 200 1 1 12 25 250 2 1 12 25 250 1
2 12 30 200 2 2 12 30 200 3 2 12 30 250 6 2 12 30 250 5
;
/* carb -> 10, 12
press -> 25, 30
speed -> 200, 250
and then the observation variable "height". eg of 2*2*2 factorial experiment.
Following is the code for producing outputs on page -234 in the table-6.5 and table-6.6*/
proc glm data=dev1 noprint;
class carb press speed rep;
model height=carb|press|speed/ss1;
/* following are the effect estimate of coded factors. x'y/(2^[k-1])
and the reg coeff are x'y/(2^k) */
estimate 'A(carb)' carb -1 1/divisor=2;
estimate 'B(press)' press -1 1/divisor=2;
estimate 'C(speed)' speed -1 1/divisor=2;
estimate 'AB' carb*press 1 -1 -1 1/divisor=4;
estimate 'BC' press*speed 1 -1 -1 1/divisor=4;
estimate 'AC' carb*speed 1 -1 -1 1/divisor=4;
estimate 'ABC' carb*press*speed -1 1 1 -1 1 -1 -1 1/divisor=8;
output out=resout1 p=predicted1 r=residual1;
run;
/*looking at the p-values we have decided to fit this model.. */
proc glm data=dev1 ;
class carb press speed rep; /* note if you want to get the actual estimates
to write the predicted equation then remove the
class statement from the code this will give you
the actual estimates: page-237-238 */
model height=carb|press speed/ss1 solution;
/* following are the effect estimate of coded factors. x'y/(2^[k-1])
and the reg coeff are x'y/(2^k) */
estimate 'A(carb)' carb -1 1/divisor=1;
estimate 'B(press)' press -1 1/divisor=1;
estimate 'C(speed)' speed -1 1/divisor=1;
estimate 'AB' carb*press 1 -1 -1 1/divisor=2;
output out=resout2 p=predicted2 r=residual2;
run;
/* I have checked with the predicted value given in the table-6.7 on page-238
this is same as the predicted value saved in "resout2".
Now, I am trying to do the contour plots and the surface plots shown in figure-6.7
page-236. */
proc g3grid data=resout2 out=gridout;
grid press*carb=predicted2/naxis1=21 naxis2=21 spline;
run;
proc gcontour data=gridout;
plot press*carb=predicted2/grid;
run;
/* now I would like to see the use of PROC FACTEX, i guess what I had to do
above for getting the estimates of main effect and itneractions...can be
done in an easier way..*/
proc factex;
factors carb press speed;
model estimate=(carb|press|speed);
output out=newdes designrep=2
carb nvals=(10 12)
press nvals=(25 30)
speed nvals=(200 250);
run;
proc print data=newdes;
run;
/* I am sorry, at this point of time I am not able to find the way out for avoiding
the estimate statements.. hopefully I will find some wayout..*/
Example-6.2
/* Example-6.2 (table 6.10) Pilot plant filtration rate experiment. */
proc factex;
factors D C B A;
output out=outdes;
run;
data filter;
input frate @@;
datalines;
45 71 48 65 68 60 80 65 43 100 45 104 75 86 70 96
;
data plant;
merge outdes filter;
run;
/* trying to produce the output in table-6.12 page 249*/
proc glm data=plant noprint outstat=statout; /* "outstat" contains the SS of all model parameters
which can be used to produce Cp-Plot*/
class A B C D;
model frate=A|B|C|D/ss1;
output out=resout1 r=residual1 p=predicted1;
/*following are the effect estimates and for reg coeff divide all of them by 2*/
estimate 'A' A -1 1/divisor=1;
estimate 'B' B -1 1/divisor=1;
estimate 'C' C -1 1/divisor=1;
estimate 'D' D -1 1/divisor=1;
estimate 'AB' A*B 1 -1 -1 1/divisor=2;
estimate 'AC' A*C 1 -1 -1 1/divisor=2;
estimate 'AD' A*D 1 -1 -1 1/divisor=2;
estimate 'BC' B*C 1 -1 -1 1/divisor=2;
estimate 'BD' B*D 1 -1 -1 1/divisor=2;
estimate 'CD' C*D 1 -1 -1 1/divisor=2;
run;
/* note that Error df is zero and therefore F-values and t-values are not calculated
so we have to use some other technique to know about the important factors..
one method is Cp-plot , basically plotting the SS for each of the factors.. */
proc glm data=plant noprint outstat=statout; /* "outstat" contains the SS of all model parameters
which can be used to produce Cp-Plot*/
model frate=A|B|C|D/ss1;
run;
proc plot data=statout ;
plot SS*_SOURCE_=_TYPE_;
run;
/* looking at the plot, I will keep the factors with large SS, this idea allows
me to fit a model with-> A,AC,AD,C,D. This is equivalent to getting the normal
probability plot of effect estimates */
/* now getting profile plots for main effects and interactions*/
proc means data=plant noprint;
class A;
var frate;
output out=meanouta mean=avga;
run;
proc means data=plant noprint;
class C;
var frate;
output out=meanoutc mean=avgc;
run;
proc means data=plant noprint;
class D;
var frate;
output out=meanoutd mean=avgd;
run;
data meanout;
merge meanouta meanoutc meanoutd;
run;
symbol i=join h=1 v=dot;
axis1 label=('text for x-axis');
axis2 label=('text for y-axis');
proc gplot data=meanout;
plot avga*A avgc*C avgd*D/overlay haxis=axis1 vaxis=axis2;
run;
/* following is the code for producing the output similar on page-250 (table-6.13)*/
proc glm data=plant noprint;
class A B C D;
model frate= A|C|D;
run;
/* so assuming that the model selected from the Cp-plot study is right the predicted
values are as follows */
proc glm data=plant;
class A B C D;
model frate= A C D A*C A*D;
output out=resout2 r=residual2 p=predicted2;
/* this will produce the output same as shown on page-251*/
proc print;
run;
/* to check if the residuals are normal or not...fig-6.13 page-252 */
proc univariate data=resout2;
probplot residual2/normal(mu=est sigma=est);
run;
/* similarly you can plot the response surface and the two contour plots
shown on page-253 , Fixing D=1 gives me the following model.. so I have
get the predicted values with respect to this model and then get the corresponding
contour plot for that model.... I have shown for one case, similarly we can get the
other one... */
proc glm data=plant;
class A B C D;
model frate= A C A*C;
output out=resout3 r=residual3 p=predicted3;
run;
proc g3grid data=resout3 out=gridca;
grid C*A=predicted3/naxis1=41 naxis2=41 spline;
run;
proc gcontour data=gridca;
plot C*A=predicted3/grid autolabel;
run;
/*Now I am going to try, how to do the Half-Normal Plot ???
This is nothing but the normal probability plot for the absolute
value of effects estimate.... If I am not able to do this, then also
you dont need to worry, the same thing can be done using the plot I
have suggested above.. */
Example-6.5
/* Example-6.5 page-265, data in table-6.16 The Oxide thickness expt.*/
proc factex;
factors D C B A ;
output out=desout1 designrep=4;
run;
data desout2;
input thickness @@;
datalines;
378 415 380 450 375 391 384 426 381 416 371 445 377 391 375 430
376 416 379 446 371 390 385 433 381 420 372 448 377 391 376 430
379 416 382 449 373 388 386 430 375 412 371 443 379 386 376 428
379 417 383 447 369 391 385 431 383 412 370 448 379 400 377 428
;
data oxide;
merge desout1 desout2;
run;
proc glm data=oxide outstat=ssout;
class A B C D;
model thickness=A|B|C|D/ss1;
output out=resout1 p=predicted1 r=residual1;
estimate 'A' A -1 1/divisor=1;
estimate 'B' B -1 1/divisor=1;
estimate 'D' C -1 1/divisor=1;
estimate 'D' D -1 1/divisor=1;
estimate 'AB' A*B 1 -1 -1 1/divisor=2;
estimate 'AC' A*C 1 -1 -1 1/divisor=2;
estimate 'AD' A*D 1 -1 -1 1/divisor=2;
estimate 'BC' B*C 1 -1 -1 1/divisor=2;
estimate 'BD' B*D 1 -1 -1 1/divisor=2;
estimate 'CD' C*D 1 -1 -1 1/divisor=2;
run;
/*
proc plot data=ssout;
plot SS*_SOURCE_=_TYPE_;
run;
looking at this plot we decide to keep A,AB,AC,B & C, now refitting
the model with these factors..*/
proc glm data=oxide outstat=ssout;
class A B C D;
model thickness=A|B A|C/ss1 solution clparm;
run;
/* to get the regression coefficients for the model remove the
class statement from the above code, and then you can get it.*/
Exercise-6.1
/* Exercise-6.1 page 276...
Study of Cutting speed(A), Tool geometry(B) & Cutting Angle(C) on the
life of a machine tool */
proc factex;
factors C B A;
output out=mach1 designrep=3;
run;
data mach2;
input life @@;
datalines;
22 32 35 55 44 40 60 39
31 43 34 47 45 37 50 41
25 29 50 46 38 36 54 47
;
data mach;
merge mach1 mach2;
run;
/* part(a): estimating all the effects...*/
proc glm data=mach outstat=ssout;
class A B C;
model life=A|C|B/ss1;
estimate 'A' A -1 1/divisor=1;
estimate 'B' B -1 1/divisor=1;
estimate 'C' C -1 1/divisor=1;
estimate 'AB' A*B 1 -1 -1 1/divisor=2;
estimate 'AC' A*C 1 -1 -1 1/divisor=2;
estimate 'BC' B*C 1 -1 -1 1/divisor=2;
estimate 'ABC' A*B*C -1 1 1 -1 1 -1 -1 1/divisor=4;
output out=resout1 r=residual1 p=predicted1;
run;
/* part(b):
proc plot data=ssout;
plot SS*_SOURCE_=_TYPE_;
run;
looking at the plot we decide to reduce the model to C,B,AC. Exactly
the same result will be obtained from the probability plot or half normal
probability plot of the effect estimates.. */
proc glm data=mach ;
class A B C;
model life=C B A*C/ss1;
run;
/* part(c): following will give you the estimates of the parameters in the regression model. */
proc glm data=mach noprint;
model life= C B A*C/ss1;
output out=resout2 r=residual2 p=predicted2;
run;
/* part(d): analysising the residuals of this reduced model..*/
proc univariate data=resout1 noprint;
probplot residual1/normal(mu=est sigma=est);
run;
proc univariate data=resout2 noprint;
probplot residual2/normal(mu=est sigma=est);
run;
/* normal probability plots do not reveal any trivial evidence against
normality.. now we should check for other */
proc plot data=resout1;
plot residual1*predicted1="*"/vref=0 vaxis=-20 to 20;
run;
proc plot data=resout2 ;
plot residual2*predicted2="*"/vref=0 vaxis=-20 to 20;
run;
/* I guess these residual vs predicted plot shows that the second model
is better and there is no evidence against normality in either case. */
/* part(d): now I have to do the interaction plots and main effect plots*/
proc means data=mach noprint;
class A;
var life;
output out=outa mean=avga;
run;
proc means data=mach noprint;
class B;
var life;
output out=outb mean=avgb;
run;
proc means data=mach noprint;
class C;
var life;
output out=outc mean=avgc;
run;
data outmean;
merge outa outb outc;
run;
/* always plot the main effects together.. so that you can see which
one in more prominent..*/
symbol i=join h=1 v=dot;
legend1 label=none;
proc gplot data=outmean;
plot avga*A avgb*B avgc*C/overlay legend=legend1;
run;
proc means data=mach noprint;
class A B;
var life;
output out=outab mean=avgab;
run;
symbol i=join h=1 v=dot;
proc gplot data=outab;
plot avgab*A=B;
run;
proc means data=mach noprint;
class B C;
var life;
output out=outbc mean=avgbc;
run;
symbol i=join h=1 v=dot;
proc gplot data=outbc;
plot avgbc*B=C;
run;
proc means data=mach noprint;
class A C;
var life;
output out=outac mean=avgac;
run;
symbol i=join h=1 v=dot;
proc gplot data=outac;
plot avgac*A=C;
run;
/* now i am doing the exercise-6.2, since this is related
to the model in previuos question I will do it in this code itself.*/
proc g3grid data=resout2 out=gridout;
grid A*C=predicted2/naxis1=31 naxis2=31 spline;
run;
proc gcontour data=gridout;
plot A*C=predicted2/grid autolabel;
run;
proc g3grid data=resout2 out=gridout;
grid B*C=predicted2/naxis1=31 naxis2=31 spline;
run;
proc gcontour data=gridout;
plot B*C=predicted2/grid autolabel;
run;
proc g3grid data=resout2 out=gridout;
grid B*A=predicted2/naxis1=31 naxis2=31 spline;
run;
proc gcontour data=gridout;
plot B*A=predicted2/grid autolabel;
run;
Exercise-6.6
/* Exercise-6.6, this question is the modification of question 6.1 */
%adxgen
%adxff
%adxcc
%adxinit
%adxffd(design,3,8,1)
%adxadcen(design,t1 t2 t3,4,block);
data mach1;
input life @@;
datalines;
22 32 35 55 44 40 60 39 36 40 43 45
;
data machfin;
merge design mach1;
proc print;
run;
Exercise-6.7
/* Exercise-6.7, yield of a chemical experiment...*/
proc factex;
factors rep D C B A;
output out=part1
rep nvals=(1,2);
run;
data part2;
input yield @@;
datalines;
90 74 81 83 77 81 88 73 98 72 87 85 99 79 87 80
93 78 85 80 78 80 82 70 95 76 83 86 90 75 84 80
;
data chempr;
merge part1 part2;
run;
/* part(a): Estimation of all the effects... note that there are 32 obsn
and 16 factors. so it is possible to estimate the error this time..
ods output Estimates =estime;
ods listing close;
*/
proc glm data=chempr outstat=ssout noprint;
class rep D C B A;
model yield= A|B|C|D/ss1;
output out=resout1 r=residual1 p=predicted1;
estimate 'A' A -1 1/divisor=1;
estimate 'B' B -1 1/divisor=1;
estimate 'C' C -1 1/divisor=1;
estimate 'D' D -1 1/divisor=1;
estimate 'AB' A*B 1 -1 -1 1/divisor=2;
estimate 'AC' A*C 1 -1 -1 1/divisor=2;
estimate 'AD' A*D 1 -1 -1 1/divisor=2;
estimate 'BC' B*C 1 -1 -1 1/divisor=2;
estimate 'BD' B*D 1 -1 -1 1/divisor=2;
estimate 'CD' C*D 1 -1 -1 1/divisor=2;
run;
/*
ods listing;
proc print data=estime;
run;
proc univariate data=estime noprint;
probplot Estimate/normal(mu=est sigma=est);
run;
/* Lets find this out... get hold of the estimate values somehow..good work..
look at the above lines using "ods" this is the way to save most of the outputs..
in fact you can save anyout that is there in the "RESULT" as a tree structure..
Reference : http://www.sfu.ca/sasdoc/sashtml/stat/chap15/sect3.htm
look at "Using ODS with the SAS Explorer" */
/* above written "estimate" statements will give the effect estimates..
now to get an idea about which factors are relvant for the model ?
Below is the method to see the similar conclusion using effect SS.
proc plot data=ssout;
plot SS*_SOURCE_=_TYPE_;
run;
/* Looking at the plot, I am considering the following factors in my model
A, AB, ABC, ABD & D */
/* Note : Examples before this are ok.. becoz fortunately.. the reduced models are such that
it contained the linear as well as the higher factors.. but this time we have only higher
factors but not the linear factor.... and this is a problem.. */
proc glm data=chempr;
class rep A B C D;
model yield= A D A*B A*B*C A*B*D/ ss1 ;
output out=resout2 r=residual2 p=predicted2;
run;
/* I dont know how to fit this model... here I want the anova of the form
source df
A 1
D 1
AB 1
ABC 1
ABD 1 ..... but what I am getting is not at all wanted.. we need
to work on this to remove the error...*/
/*
data resout;
merge resout1 resout2;
run;
proc univariate data=resout noprint;
probplot residual1/normal(mu=est sigma=est);
probplot residual2/normal(mu=est sigma=est);
run;
proc plot data=resout;
plot residual1*predicted1="*"/vref=0 vaxis=-7.5 to 7.5;
plot residual2*predicted2="*"/vref=0 vaxis=-7.5 to 7.5;
run;
*/
Exercise-6.23
/* Exercise-6.23 Study of yield with four factors.. */
proc factex;
factors D C B A;
output out=part1;
/* D nvals=(225 250)
C nvals=(60 80)
B nvals=(14 18)
A nvals=(2.5 3); */
run;
data part2;
input yield @@;
datalines;
12 18 13 16 17 15 20 15 10 25 13 24 19 21 17 23
;
data proce;
merge part1 part2;
run;
/* doing part(a): first I will get the estimates at then I will make
the normal probability plot of these estimates..*/
ods output Estimates=estim;
ods listing close;
proc glm data=proce outstat=ssout;
class A B C D;
model yield= A|B|C|D/ss1;
output out=resout1 r=resid1 p=predict1;
estimate 'A' A -1 1/divisor=1;
estimate 'B' B -1 1/divisor=1;
estimate 'C' C -1 1/divisor=1;
estimate 'D' D -1 1/divisor=1;
estimate 'AB' A*B 1 -1 -1 1/divisor=2;
estimate 'AC' A*C 1 -1 -1 1/divisor=2;
estimate 'AD' A*D 1 -1 -1 1/divisor=2;
estimate 'BC' B*C 1 -1 -1 1/divisor=2;
estimate 'BD' B*D 1 -1 -1 1/divisor=2;
estimate 'CD' C*D 1 -1 -1 1/divisor=2;
estimate 'ABC' A*B*C -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ABD' A*B*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ACD' A*C*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'BCD' B*C*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ABCD' A*B*C*D 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1/divisor=8;
run;
ods listing;
/*
proc print data=estim;
run;
data estim1;
set estim;
absestim=abs(Estimate);
run;
title 'Normal Probability plot';
symbol v=square h=1;
proc univariate data=estim noprint;
probplot Estimate/normal(mu=est sigma=est) ;
inset normal;
run;
/* Well, I have done the normal probability plot of effects but the problem
is that I dont have the labels ... so I will stick to the SS plots
proc print data=ssout;
run;
proc plot data=ssout;
plot SS*_SOURCE_=_TYPE_;
run;
/* looking at the SS of the effects we can say what are the important factors
in the model.. also it is clear from the normal probability plot of the effect estimates.
and the factors are : A,C,D,AC AD. This is verified from the profile plots */
proc glm data=proce;
class A B C D;
model yield= A|C A|D/ss1;
output out=resout2 r=resid2 p=predict2;
run;
/* part(C): the required regression model is obtained from the
regression coefficient estimates given by the following code...*/
proc glm data=proce;
model yield= A|C A|D;
run;
/* part(d): Analysis of the residual is as follows... */
proc plot data=resout2;
plot resid2*(predict2 A C D)="*"/vref=0;
run;
proc univariate data=resout2 noprint;
probplot resid2/normal(mu=est sigma=est);
inset normal;
run;
/* looking at the residual plots.. I can say that there is no strong
evidence against normality of the residuals..*/
/* If I ask that which combination is the best for the yield. then we have to
do the profile plots..*/
proc means data=proce noprint;
class A;
var yield;
output out=outa mean=avga;
run;
proc means data=proce noprint;
class B;
var yield;
output out=outb mean=avgb;
run;
proc means data=proce noprint;
class C;
var yield;
output out=outc mean=avgc;
run;
proc means data=proce noprint;
class D;
var yield;
output out=outd mean=avgd;
run;
data meanout;
merge outa outb outc outd;
run;
title 'Profile plots';
symbol i=join h=1 v=dot;
legend1 label=none;
proc gplot data=meanout;
plot avga*A avgb*B avgc*C avgd*D/overlay legend=legend1;
run;
/* Exercise-6.24,to get the response surface of the model in part(c)... */
proc g3grid data=resout2 out=gridout;
grid A*C=predict2/naxis1=31 naxis2=31 spline;
run;
proc gcontour data=gridout;
plot A*C=predict2/grid;
run;
proc g3grid data=resout2 out=gridout;
grid A*D=predict2/naxis1=31 naxis2=31 spline;
run;
proc gcontour data=gridout;
plot A*D=predict2/grid;
run;
proc g3grid data=resout2 out=gridout;
grid D*C=predict2/naxis1=31 naxis2=31 spline;
run;
proc gcontour data=gridout;
plot D*C=predict2/grid;
run;
Exercise-6.31
/* Exercise-6.31.... Resistivity on a silicon wafer */
proc factex;
factors D C B A;
output out=part1;
run;
data part2;
input resis @@;
datalines;
1.92 11.28 1.09 5.75 2.13 9.53 1.003 5.35 1.60 11.73 1.16 4.68 2.16 9.11 1.07 5.30
;
data silic;
merge part1 part2;
run;
data silic2;
set silic;
logresis=log(resis);
output;
run;
/* part(a): Estimation of main effects and interactions, normal probability plot
and profile plots are as follows */
proc means data=silic;
class A;
var resis;
output out=outa mean=avga;
run;
proc means data=silic;
class B;
var resis;
output out=outb mean=avgb;
run;
proc means data=silic;
class C;
var resis;
output out=outc mean=avgc;
run;
proc means data=silic;
class D;
var resis;
output out=outd mean=avgd;
run;
data meanout;
merge outa outb outc outd;
run;
symbol i=join h=1 v=dot;
legend1 label=none;
proc gplot data=meanout;
plot avga*A avgb*B avgc*C avgd*D/overlay legend=legend1;
run;
/* if you look at this plot(above).. then you can see that the
main effect of C & D are almost negligible */
proc means data=silic;
class A B;
var resis;
output out=outab mean=avgab;
run;
symbol i=join v=dot h=1;
proc gplot data=outab;
plot avgab*A=B;
run;
proc means data=silic;
class A C;
var resis;
output out=outac mean=avgac;
run;
symbol i=join v=dot h=1;
proc gplot data=outac;
plot avgac*A=C;
run;
proc means data=silic;
class A D;
var resis;
output out=outad mean=avgad;
run;
symbol i=join v=dot h=1;
proc gplot data=outad;
plot avgad*A=D;
run;
proc means data=silic;
class B C;
var resis;
output out=outbc mean=avgbc;
run;
symbol i=join v=dot h=1;
proc gplot data=outbc;
plot avgbc*B=C;
run;
proc means data=silic;
class B D;
var resis;
output out=outbd mean=avgbd;
run;
symbol i=join v=dot h=1;
proc gplot data=outbd;
plot avgbd*B=D;
run;
proc means data=silic;
class C D;
var resis;
output out=outcd mean=avgcd;
run;
symbol i=join v=dot h=1;
proc gplot data=outcd;
plot avgcd*C=D;
run;
quit;
/* up to this point i was doing all the interaction plot and the
main effect plots.. to get a visual look about the interactions..*/
/* Above interaction plot suggests that AD BD & CD are almost zero.. */
/* here SS and its corresponding P-value will give me the idea that
which one is playing important role in the model */
ods output Estimates=estimet;
proc glm data=silic outstat=ssout ;
class A B C D;
model resis=A|B|C|D/ss1;
output out=resout1 r=resid1 p=predict1;
estimate 'A' A -1 1/divisor=1;
estimate 'B' B -1 1/divisor=1;
estimate 'C' C -1 1/divisor=1;
estimate 'D' D -1 1/divisor=1;
estimate 'AB' A*B 1 -1 -1 1/divisor=2;
estimate 'AC' A*C 1 -1 -1 1/divisor=2;
estimate 'AD' A*D 1 -1 -1 1/divisor=2;
estimate 'BC' B*C 1 -1 -1 1/divisor=2;
estimate 'BD' B*D 1 -1 -1 1/divisor=2;
estimate 'CD' C*D 1 -1 -1 1/divisor=2;
estimate 'ABC' A*B*C -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ABD' A*B*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ACD' A*C*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'BCD' B*C*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ABCD' A*B*C*D 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1/divisor=8;
run;
ods listing;
proc univariate data=estimet noprint;
probplot Estimate/normal(mu=est sigma=est);
inset normal;
run;
proc plot data=ssout;
plot SS*_SOURCE_=_TYPE_;
run;
/* SS plot shows me that A B and AB are the only significant factors in
this model.. and rest are very close to zero.. and the normal plot of
estimates also shows the same thing.. so our reduced model will be
resis= A B AB */
/* part(b): model adequacy checking for the reduced model..*/
proc glm data=silic;
class A B C D;
model resis=A|B/ss1;
output out=resout2 r=resid2 p=predict2;
run;
proc univariate data=resout2 noprint;
probplot resid2/normal(mu=est sigma=est);
inset normal;
run;
proc plot data=resout2;
plot resid2*(predict2 A B)="*"/vref=0;
run;
/* these plots suggest that the model is not a good fit ,the
normality assumptions are violated..*/
/*part(c): transformation of datat---> log(y)... */
Exercise-6.31(d)
proc factex;
factors D C B A;
output out=part1;
run;
data part2;
input resis @@;
datalines;
1.92 11.28 1.09 5.75 2.13 9.53 1.003 5.35 1.60 11.73 1.16 4.68 2.16 9.11 1.07 5.30
;
data silic0;
merge part1 part2;
run;
data silic;
set silic0;
logresis=log(resis);
output;
run;
/* part(a): Estimation of main effects and interactions, normal probability plot
and profile plots are as follows */
proc means data=silic;
class A;
var logresis;
output out=outa mean=avga;
run;
proc means data=silic;
class B;
var logresis;
output out=outb mean=avgb;
run;
proc means data=silic;
class C;
var logresis;
output out=outc mean=avgc;
run;
proc means data=silic;
class D;
var logresis;
output out=outd mean=avgd;
run;
data meanout;
merge outa outb outc outd;
run;
symbol i=join h=1 v=dot;
legend1 label=none;
proc gplot data=meanout;
plot avga*A avgb*B avgc*C avgd*D/overlay legend=legend1;
run;
/* if you look at this plot(above).. then you can see that the
main effect of C & D are almost negligible */
proc means data=silic;
class A B;
var logresis;
output out=outab mean=avgab;
run;
symbol i=join v=dot h=1;
proc gplot data=outab;
plot avgab*A=B;
run;
proc means data=silic;
class A C;
var logresis;
output out=outac mean=avgac;
run;
symbol i=join v=dot h=1;
proc gplot data=outac;
plot avgac*A=C;
run;
proc means data=silic;
class A D;
var logresis;
output out=outad mean=avgad;
run;
symbol i=join v=dot h=1;
proc gplot data=outad;
plot avgad*A=D;
run;
proc means data=silic;
class B C;
var logresis;
output out=outbc mean=avgbc;
run;
symbol i=join v=dot h=1;
proc gplot data=outbc;
plot avgbc*B=C;
run;
proc means data=silic;
class B D;
var logresis;
output out=outbd mean=avgbd;
run;
symbol i=join v=dot h=1;
proc gplot data=outbd;
plot avgbd*B=D;
run;
proc means data=silic;
class C D;
var logresis;
output out=outcd mean=avgcd;
run;
symbol i=join v=dot h=1;
proc gplot data=outcd;
plot avgcd*C=D;
run;
quit;
/* up to this point i was doing all the interaction plot and the
main effect plots.. to get a visual look about the interactions..*/
/* Above interaction plot suggests that AD BD & CD are almost zero.. */
ods output Estimate=estim;
/* here SS and its corresponding P-value will give me the idea that
which one is playing important role in the model */
proc glm data=silic outstat=ssout ;
class A B C D;
model logresis=A|B|C|D/ss1;
output out=resout1 r=resid1 p=predict1;
estimate 'A' A -1 1/divisor=1;
estimate 'B' B -1 1/divisor=1;
estimate 'C' C -1 1/divisor=1;
estimate 'D' D -1 1/divisor=1;
estimate 'AB' A*B 1 -1 -1 1/divisor=2;
estimate 'AC' A*C 1 -1 -1 1/divisor=2;
estimate 'AD' A*D 1 -1 -1 1/divisor=2;
estimate 'BC' B*C 1 -1 -1 1/divisor=2;
estimate 'BD' B*D 1 -1 -1 1/divisor=2;
estimate 'CD' C*D 1 -1 -1 1/divisor=2;
estimate 'ABC' A*B*C -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ABD' A*B*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ACD' A*C*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'BCD' B*C*D -1 1 1 -1 1 -1 -1 1/divisor=4;
estimate 'ABCD' A*B*C*D 1 -1 -1 1 -1 1 1 -1 -1 1 1 -1 1 -1 -1 1/divisor=8;
run;
ods listing;
proc print data=estim;
run;
proc print data=ssout;
run;
proc univariate data=estim noprint;
probplot Estimate/normal(mu=est sigma=est);
inset normal;
run;
proc plot data=ssout;
plot SS*_SOURCE_=_TYPE_;
run;
/* SS plot shows me that A B and AB are the only significant factors in
this model.. and rest are very close to zero.. and the normal plot of
estimates also shows the same thing.. so our reduced model will be
logresis= A B AB */
/* part(b): model adequacy checking for the reduced model..*/
proc glm data=silic;
class A B C D;
model logresis=A B/ss1;
output out=resout2 r=resid2 p=predict2;
run;
proc univariate data=resout2 noprint;
probplot resid2/normal(mu=est sigma=est);
inset normal;
run;
proc plot data=resout2;
plot resid2*(predict2 A B)="*"/vref=0;
run;
/* these plots suggest that the model is not a good fit ,the
normality assumptions are violated..*/
/*part(c): transformation of datat---> log(y)... */
/* After taking the log() transformation of the data, I saw that the
interaction between A B totally gone... and some other inteaction
came up... like the model with (A B CD AC) would be perfect from
the profile plot...
from normal prob plot of estimates... or the estimate effect estimates..
it is clear that we should take some effects like : A C D AC AD ...
so in this case... after taking the log transformation the ANOVA result
and the result given by effect estimates vary... */
Example-7.1
/* Example-7.1 , data is shown in table-7.1, chemical process Expt.*/
options ps=60;
proc factex;
factors speed feed angle;
blocks nblocks=2;
model
output out=blocd speed nvals=(300 500)
feed nvals=(20 30 )
angle nvals=(6 8 )
blockname=machine cvals=('A' 'B');
run;
proc print data=blocd;
run;
proc print data=part1;
run;
data part2;
input yield @@;
datalines;
28 36 18 31 25 32 19 30 27 32 23 29
;
data chem;
merge part1 part2;
proc print;
run;
Example-7.2
/* Example-7.2,page-293(similar situation as in the example-6.2)
factors are :
temperature -> A
pressure -> B
concentration -> C
stirring rate -> D
response -> filteration rate
and the experiment is done in 2 blocks each of size 8 and the
factor confounded with blocks is : ABCD
proc factex;
factors D C B A;
output out=part1;
proc print;
run;
*/
proc factex;
factors D C B A;
size design=16;
model res=max;
blocks nblocks=2;
output out=part1;
run;
proc sort data=part1;
by BLOCK;
run;
data part2;
input rate @@;
datalines;
71 48 68 65 43 104 86 70 25 45 40 60 80 25 55 76
;
data filter;
merge part1 part2;
run;
/* now to produce outputs like shown in table-7.6, this will
give me the SS of the required model and the break-up of
all possible interactions.. For getting the effect estimates
you need to write all the "estimate" statements properly.*/
proc glm data=filter outstat=ssout;
class A B C D BLOCK;
model rate=A|B|C|D@3 BLOCK/ss1;
run;
quit;
/* following commands give me the file of regression estimates for
main effects and interaction effects...*/
ods output ParameterEstimates=estim;
ods listing close;
proc glm data=filter;
model rate=A|B|C|D@3 /ss1;
run;
ods listing;
proc print data=estim;
run;
/* this will show me which effects are significant in the model
we could also get the estimate effects from the main model and
then do the normal probability plot */
proc plot data=ssout;
plot SS*_SOURCE_=_TYPE_;
run;
/* this plot says that the factors important are :
A AC AD C D . Now to produce an output similar to table-7.7*/
proc glm data=filter noprint;
class A B C D BLOCK;
model rate= BLOCK A|C|D@2/ss1;
run;
Example-7.3
/* Example-7.3 ,page-300, Partial confounding scheme.. */
proc factex;
factors C B A;
size design=8;
model estimate=(A|B|C@2); /*forcing to confound ABC*/
blocks nblocks=2;
output out=rep1 blockname=block ;
run;
data rep11;
set rep1;
do i=1 to 8;
rep =1;
end;
output;
keep rep block A B C;
proc factex;
factors C B A;
size design=8;
model estimate=(A B C A*C B*C A*B*C); /* forcing to confound AB */
blocks nblocks=2;
output out=rep2 blockname=block ;
run;
data rep21;
set rep2;
do i=1 to 8;
rep =2;
end;
output;
keep rep block A B C;
/*this is the way to combine the data sets..*/
data part1;
set rep11 rep21;
run;
data part2;
input height @@;
datalines;
-3 2 2 1 0 -1 -1 6 1 0 1 1 -1 3 0 5
;
data fillh;
merge part1 part2;
proc print;
run;
/* following is the code which gives the required SS taken properly.. which is
same as the output shown in table-7.3...*/
proc glm data=fillh;
class rep block A B C;
model height=rep block(rep) A|B|C /ss1 ;
run;
/* note : "block(rep)" implies block within replication .. and
SS(ABC) & SS(AB) have been calulated from one replication only..*/
Example-8.1
/* Example-8.1, the same example with filteration rate data.. as shown in example-6.2 */
proc factex;
factors C B A D;
size fraction=2;
model res=4;
examine confounding aliasing;
output out=part1;
run;
data part2;
input rate @@;
datalines;
45 100 45 65 75 60 80 96
;
data filter;
merge part1 part2;
run;
proc glm data=filter;
class A B C D;
model rate=A|B|C|D /ss1;
estimate 'A+BCD' A -1 1 B*C*D -1 1 1 -1 1 -1 -1 1;
run;
/* the above estimate statement is not correct.. I dont know yet, about how
to get the estimates of main effects and interaction in this case.
but here is an alternate way I am using a different section of SAS for this.. */
%adxgen
%adxff
%adxinit
%adxffa(filter,rate,C B A D,4);
/* the above code gives the required output.. though this is 1/2 of the outputs
from the result given in table-8.4. this gives the normal probability plot
of estimates which shows that except B & AB all are important factors in the model.
which is also justified by the ANOVA output...*/
Example-8.2
/*Example-8.2 , design for process improvement... */
proc factex;
factors D C B A E;
size fraction=2;
model res=5;
examine confounding aliasing;
output out=part1;
run;
data part2;
input yield @@;
datalines;
8 9 34 52 16 22 45 60 6 10 30 50 15 21 44 63
;
data proce;
merge part1 part2;
run;
/* this will give me the SS kind of stuff.. */
proc glm data=proce;
class A B C D E;
model yield=A|B|C|D|E/ss1;
run;
/* the following code will give me the estimates and the normal
probability plot of the estimates..*/
%adxgen
%adxff
%adxinit
%adxffa(proce,yield,A B C D E,5);
/* for getting the regression coefficients and SS .. use the following code..
so for this particular case you dont need the above "proc glm" statement.*/
proc glm data=proce;
model yield=A|B|C|D|E;
run;
/* looking at the normal probability plot of estimates and the SS the new
reduced model is as follows..table-8.7*/
proc glm data=proce;
model yield=A|B C;
output out=resout r=resid p=predict;
run;
/* this is the plot done in figure-8.7 page-313*/
symbol v=dot h=1;
proc univariate data=resout;
probplot resid/normal(mu=est sigma=est);
inset normal;
run;
Example-8.3
/* Example-8.3 the same data as in example-8.1, with a different design sequence..*/
/* this is the design created by factex...*/
proc factex;
factors C B A D;
size fraction=2;
model res=4;
examine A(3) C;
output out=part1;
run;
/* this is the way to get the other half of the fraction..*/
data part1;
set part1;
D=-D;
output;
keep A B C D;
data part2;
input rate @@;
datalines;
43 71 48 104 68 86 70 65
;
data filter;
merge part1 part2;
proc print;
run;
%adxgen
%adxff
%adxinit
%adxffa(filter,rate,A B C D,4);
/* note this is half of the estimates given in the table on page-317.
one more thing, that the output of these adx codes are not the main effects
or interactions.. they are the l(A) & l'(A)'s resply.. */
Example-8.4
/* Example-8.4 , table-8.10 Injection molding experiment...*/
/* the default defining fraction taken is as follows :
proc factex;
factors D C B A E F;
size fraction=4;
model res=max;
examine A(3) C;
output out=deslay;
run;
/* but if we want the actual design given to us then we have to follow the
code I have written below...*/
proc factex;
factors D C B A;
output out=part1;
run;
data part1;
set part1;
E=A*B*C;
F=B*C*D;
output;
keep A B C D E F;
data part2;
input shrink @@;
datalines;
6 10 32 60 4 15 26 60 8 12 34 60 16 5 37 52
;
data inject;
merge part1 part2;
run;
/* to produce the output shown in table-8.11 on page-321*/
proc glm data=inject;
model shrink=A|B|C|D|E|F;
run;
%adxgen
%adxff
%adxinit
%adxffa(inject,shrink,A B C D E F,4);
/* looking at this normal probability plot...and SS from ANOVA,
the reduced model would be y= A B AB*/
proc glm data=inject;
model shrink=A|B;
output out=resout r=resid p=predict;
run;
/* the above will give you the model written on page 322*/
symbol v=dot h=1;
proc univariate data=resout noprint;
probplot resid/normal(mu=est sigma=est);
inset normal;
run;
proc plot data=resout;
plot resid*predict="*"/vref=0;
run;
Example-8.6
/* Example-8.6 CNC machine expt....page-332 data on table-8.16*/
proc factex;
factors E D C B A;
output out=part1;
run;
data part1;
set part1;
F=A*B*C;
G=A*B*D;
H=B*C*D*E;
output;
keep A B C D E F G H;
data part2;
input dev @@;
datalines;
2.76 6.18 2.43 4.01 2.48 5.91 2.39 3.35 4.40 4.10 3.22 3.78 5.32 3.87 3.03
2.95 2.64 5.50 2.24 4.28 2.57 5.37 2.11 4.18 3.96 3.27 3.41 4.30 4.44 3.65 4.41 3.40
;
data part2;
set part2;
logdev=log(dev);
output;
keep logdev;
data cnc;
merge part1 part2;
proc print;
run;
/* to get the sum of squares as shown in table-8.17 we have to do the GLM stuff..*/
proc glm data=cnc ;
model logdev=A|B|C|D|E|F|G|H/ss1;
run;
/* for getting the estimates... regression estimates... and the normal probability plots.
to get the actuals effect estimates... I mean the estimated aliased effect.. just do the
2*Estimate vector of the following output... */
%adxgen
%adxff
%adxinit
%adxffa(cnc,logdev,A B C D E F G H,4);
/* looking at the above plot and the ANOVA output we came up to the result
that the reduced model is :::: Oops... one thing is missing.. BLOCKS..
so you have to get the block data also in the model and then do it..
output shown in table-8.18*/
proc glm data=cnc;
model logdev=A B D A*D/ss1;
output out=resout r=resid p=predict;
run;
symbol v=dot h=1; /* plot in figure-8.20.....*/
proc univariate data=resout noprint;
probplot resid/normal(mu=est sigma=est);
inset normal;
run;
proc plot data=resout;
plot resid*predict="*"/vref=0;
run;
/* to get the plot in figure-8.21... the profile plot...we have
to get the means.. and then plot it..*/
proc means data=cnc noprint;
class A D;
var logdev;
output out=meanout mean=avgad;
run;
symbol i=join h=1 v=dot;
proc gplot data=meanout;
plot avgad*A=D;
run;
Exercise-8.28
/* Exercise-8.28, 10 factors data on sand casting...*/
proc factex;
factors d c b a;
output out=part1;
run;
data part1;
set part1;
e=c*d;
f=b*d;
g=b*c;
h=a*c;
j=a*b;
k=a*b*c;
output;
keep a b c d e f g h j k;
data part2;
input p arcp fandt;
datalines;
.958 1.364 1.363
1.00 1.571 1.555
.977 1.419 1.417
.775 1.077 1.076
.958 1.364 1.363
.958 1.364 1.363
.813 1.124 1.123
.906 1.259 1.259
.679 .969 .968
.781 1.081 1.083
1.00 1.571 1.556
.896 1.241 1.242
.958 1.364 1.363
.818 1.130 1.130
.841 1.161 1.160
.955 1.357 1.356
;
data sand;
merge part1 part2;
run;
/* part(a) & part(b): this gives the estimate effects and the normal
probability plot of these estimated effects...this also gives the aliasing
structure upto two factor interaction...*/
%adxgen
%adxff
%adxinit
%adxffa(sand,p,a b c d e f g h j k,4);
proc glm data=sand noprint outstat=ssout;
class a b c d e f g h j k;
model p=a|b|c|d|e|f|g|h|j|k@2/ss1;
run;
proc plot data=ssout;
plot SS*_SOURCE_=_TYPE_;
run;
proc glm data=sand;
class k f;
model p=k f /ss1;
output out=resout r=resid p=predict;
run;
proc plot data=resout;
plot resid*predict="*"/vref=0;
run;
proc univariate data=resout noprint;
probplot resid/normal(mu=est sigma=est);
inset normal;
run;
Example-9.1
/*Example-9.1 Syrup loss data (table-9.1)
notations : press-> C, nozz->A speed->B*/
proc factex;
factors C A B/nlev=3;
output out=part1
C nvals=(10 15 20)
A nvals=(1 2 3)
B nvals=(100 120 140)
designrep=2;
run;
proc sort data=part1;
by C;
run;
data part2;
input sloss @@;
datalines;
-35 -45 -40 17 -65 20 -39 -55 15 -25 -60 15 24 -58 4 -35 -67 -30
110 -10 80 55 -55 110 90 -28 110 75 30 54 120 -44 44 113 -26 135
4 -40 31 -23 -64 -20 -30 -61 54 5 -30 36 -5 -62 -31 -55 -52 4
;
data syrup;
merge part1 part2;
run;
/* to produce the result similar to the table-9.2, the usual anova */
proc glm data=syrup;
class A B C;
model sloss=A|B|C/ss1;
run;
/* now to produce the images that are shown on page-369..The interaction
plots..AB AC & BC . I will plot it for AB and the code for AC & BC will be
exactly similar..
proc means data=syrup noprint;
class A B;
var sloss;
output out=meanab mean=avgab;
run;
symbol i=join v=dot h=1;
proc gplot data=meanab;
plot avgab*A=B;
run;
*/
/* now to get the regression equations as given on the bottom of the page-369,
in table-9.3, the coeffcients in the model will give the regression eqn.. */
proc sort data=syrup;
by A;
run;
proc glm data=syrup noprint;
by A;
model sloss=B C B*B C*C B*C/ss1;
output out=resout r=resid p=predict;
proc print;
run;
/* now to show the contour plots shown in figure-9.5 on page-371 */
proc g3grid data=resout out=gridout;
by A;
grid C*B=predict/naxis1=31 naxis2=31 spline;
run;
proc gcontour data=gridout;
by A;
plot C*B=predict/grid autolabel;
run;
/* note: in some of the previous questions we had to do the
contour plot for one level fixed and I could not do that at that
point of time.. so here it goes... you can do that in this manner.*/
Example-9.2
/*Example-9.2, fabricated data to illustrate
the confounding scheme in 3^2 factorial expt..
proc factex;
factors A B /nlev=3;
size design=9;
model res=max;
blocks nblocks=3;
examine C A;
output out=deslay blockname=block
A nvals=(0 1 2)
B nvals=(0 1 2);
run;
*/
data manip;
input block A B val @@;
datalines;
1 0 0 4 1 1 1 -4 1 2 2 0
2 1 0 -2 2 2 1 1 2 0 2 8
3 0 1 5 3 1 2 -5 3 2 0 0
;
proc glm data=manip;
class A B block;
model val=block A|B/ss1;
run;
/* this above code will give you the output shown on page 375, table-9.4*/
Exercise-9.6
/* Exercise-9.6,Chemical process expt. */
proc factex;
factors temp press/nlev=3;
output out=part1 designrep=2
temp nvals=(80 90 100)
press nvals=(100 120 140);
run;
proc sort data=part1;
by temp press;
run;
data part2;
input yield @@;
datalines;
47.58 48.77 64.97 69.22 80.92 72.60
51.86 82.43 88.47 84.23 93.95 88.54
71.18 92.77 96.57 88.72 76.58 83.04
;
data chem;
merge part1 part2;
run;
/* part(a): usual ANOVA to get an over all idea about the factors..
I will extract the residuals to plot for part(b) to get an idea about
the model adequacies.. */
proc glm data=chem noprint;
class press temp;
model yield=temp|press/ss1;
output out=resout r=resid p=predict;
run;
/* the p-values in the ANOVA suggests that interaction between AB can be
ignored.. but we will have to check that exactly which order of the interaction
should be ignored..
part(b): residual analysis...
proc univariate data=resout noprint;
probplot resid/normal(mu=est sigma=est);
run;
proc plot data=resout;
plot resid*predict="*"/vref=0;
run;
/* the normal probability plot of residuals and residual vs predicted plots are
not bad and doesn;t give any evidence against normality...
part(c) : verification of the model... */
proc glm data=chem;
class press temp;
model yield=press|temp;
estimate 'temp' temp -1 0 1/divisor=2;
estimate 'press' press -1 0 1/divisor=2;
estimate 'temp^2' temp 1 -2 1/divisor=2;
run;
/*part(d): the following code will give you the required model ... */
proc glm data=chem;
model yield= press temp press*press temp*temp press*temp/ss1;
run;
/*part(e): Contour plot...*/
proc g3grid data=resout out=gridout;
grid press*temp=predict/naxis1=41 naxis2=41 spline;
run;
proc gcontour data=gridout;
plot press*temp=predict/grid autolabel;
run;
/*part(f): looking at the profile plot.. or interaction plot..*/
proc means data=chem noprint;
class press temp;
var yield;
output out=meanout mean=avg;
run;
symbol i=join v=dot h=1;
proc gplot data=meanout;
plot avg*temp=press;
run;
/* looked at both the plots and concluded the same result...
recommended combination would be... temp-90, press-140 OR temp-100 press-120*/
Example-11.1
/*Example-11.1,Operating conditions..
First-order Response Surface Model...*/
proc factex;
factors x1 x2;
output out=part1;
run;
%adxgen
%adxff
%adxcc
%adxinit
%adxadcen(part1,x1 x2,5);
data part2;
input y @@;
datalines;
39.3 40.0 40.9 41.5 40.3 40.5 40.7 40.2 40.6
;
data time;
merge part1 part2;
run;
/* If I use the "PROC RSREG" procedure it will fit the second order
response surface model by default, but what I want here is the "first order
response surfac model"...So, how to do that ??? */
proc rsreg data=time;
model y= x1 x2;
run;
/* to get the first-order model fit... all you have to do is the
following.. get the output from the second order fit of the model
eg. from the output of the above code .. I am getting the model
fit for the first order model fit...(table-11.2, page-433)
Regression DF Sum of Squares
(model) Linear 2 2.825000
(Residual)
(pure-Quadratic) 1 0.002722
(Interaction-Crossproduct) 1 0.002500
(Pure error - Total Error) 4 0.172000
(total) ---------- get from the proc glm statement..-----
following is the predicted regression equation ---------->
Parameter
Estimate
Standard from Coded
Parameter DF Estimate Error t Value Pr > |t| Data
Intercept 1 40.460000 0.092736 436.29 <.0001 40.460000
x1 1 0.775000 0.103682 7.47 0.0017 0.775000
x2 1 0.325000 0.103682 3.13 0.0350 0.325000
You can also use PROC GLM to get such results.. eg. */
proc glm data=http://www.geocities.com/prima001/comps/time;
model y= x1|x2 x1*x1 x2*x2/ss1;
run;
/* in this case you have the whole thing clear and nice...
Source DF Squares Mean Square F Value Pr > F
(model SS will
be obtained from
the sum of these two)
x1 1 2.40250000 2.40250000
x2 1 0.42250000 0.42250000
(interaction) x1*x2 1 0.00250000 0.00250000
(quadratic will be sum
of these two terms)
x1*x1 1 0.00272222 0.00272222
x2*x2 0 0.00000000 .
(Pure error)Error 4 0.17200000 0.04300000
(total) Corrected Total 8 3.00222222
following will give me the regression equation coefficients.....
Parameter Estimate Error t Value Pr > |t|
Intercept 40.46000000 0.09273618 436.29 <.0001
x1 0.77500000 0.10368221 7.47 0.0017
x2 0.32500000 0.10368221 3.13 0.0350
*/
quit;
Example-11.2
/*Example-11.2,Central composite design for
example-11.2, Second order Response Surface model.*/
data ccd;
input zi1 zi2 yield @@;
datalines;
80 170 76.5
80 180 77
90 170 78
90 180 79.5
85 175 79.9
85 175 80.3
85 175 80.0
85 175 79.7
85 175 79.8
92.07 175 78.4
77.93 175 75.6
85 182.07 78.5
85 167.93 77
;
/* trying to produce a similar output as shown in table-11.7 on page-443*/
proc rsreg data=ccd;
model yield = zi1 zi2/lackfit;
ridge max;
run;
proc glm data=ccd;
model yield= zi1 zi2 zi1*zi1 zi2*zi2 zi1*zi2/ss1;
run;
Example-11.3
/*Example-11.3,A three component mixture..
The (3,2) Simplex Lattice design for the yarn Elongation problem*/
data yarn;
input x1 x2 x3 eval @@;
datalines;
1 0 0 11.7
.5 .5 0 15.3
0 1 0 9.4
0 .5 .5 10.5
0 0 1 16.4
.5 0 .5 16.9
;
/* trying to produce the output shown on page-477,table-11.13 */
proc glm data=yarn noprint;
model eval=x1|x2|x3@2/ss1 noint;
output out=resout r=resid p=predict;
run;
/* Now I am trying to make the contour plot shown at the bottom of the page-477
figure-11.34...I dont know how to do such plots in SAS.. following is just the
ordinary contour plots...*/
proc g3grid data=resout out=gout;
grid x1*x2=predict/naxis1=31 naxis2=31 ;
run;
proc gcontour data=gout;
plot x1*x2=predict/grid;
run;
quit;
Example-12.2
/*Example-12.2, A Measurement systems Capability Study..*/
data part1;
do i=1 to 20;
do j=1 to 3;
do k=1,2;
part = i;
oper = j;
output;
end;
end;
end;
keep part oper;
data part2;
input meas @@;
datalines;
21 20 20 20 19 21 24 23 24 24 23 24 20 21 19 21 20 22 27 27 28 26 27 28
19 18 19 18 18 21 23 21 24 21 23 22 22 21 22 24 22 20 19 17 18 20 19 18
24 23 25 23 24 24 25 23 26 25 24 25 21 20 20 20 21 20 18 19 17 19 18 19
23 25 25 25 25 25 24 24 23 25 24 25 29 30 30 28 31 30 26 26 25 26 25 27
20 20 19 20 20 20 19 21 19 19 21 23 25 26 25 24 25 25 19 19 18 17 19 17
;
data syms;
merge part1 part2;
run;
/* the following code will give you the output required for
the table-12.4 on page-521 */
proc mixed data=syms method=type1 noinfo;
class part oper;
model meas=;
random part|oper;
run;
/* since the p-value for the interaction term is very large, ie there
is very less evidence of interaction.. so we will fit a reduced model
for getting the output similar to the output shown on page-522,table-12.5
proc mixed data=syms method=type1 noinfo;
class part oper;
model meas=;
random part oper;
run;
/* Example-12.3, it says that if I have the operator as FIXED factor
and "part" as a random factor... then how to analyse the data...this
example refers to a REstricted model... and the restriction is that
sum of all the interaction over the sum of fixed effects is zero.*/
proc mixed data=syms method=type1 cl noinfo;
class part oper;
model meas=oper/cl;
random part part*oper ;
run;
/* the above code gives the output which is given on table-12.7, page 528..
This refers to example-12.4 The unrestricted model...
Note : "cl" option in the proc mixed statement gives the confidence interval for
variance component and there are many methods for doing this.. default method is REML*/
Table-13.14
/*table-13.14 ,page-574,Tensile strength data..Split plot design*/
data part1;
do i =1 to 4;
do j=1 to 3;
do k=1 to 3;
temp=i; rep=j; method=k;
output;
end;
end;
end;
keep temp rep method;
data part2;
input stren @@;
datalines;
30 34 29 28 31 31 31 35 32 35 41 26 32 36 30 37 40 34
37 38 33 40 42 32 41 39 39 36 42 36 41 40 40 40 44 45
;
data paper;
merge part1 part2;
run;
/* the following code produced the output shown in table-13.16.*/
proc glm data=paper;
class temp rep method;
model stren=rep|method|temp/ss1;
test h=method e=rep*method;
test h=temp e=rep*temp;
test h=temp*method e=rep*method*temp;
run;
/* we need to writes lots of test statement in an appropriate way to get
the appropriate p-values for the appropriate tests.. but in the PROC MIXED
all you have to do is to specify the model properly, I mean random and fixed
effects properly... then it will give you the whole thing on its own...*/
proc mixed data=paper noinfo cl method=type1;
class temp rep method;
model stren = rep method|temp;
random rep*method|rep*temp;
run;
Example-13.1
/* Example-13.1 ,Coded purity data... page-560,table-13.3 Split plot design.
Small problem !! I can;t use proc factex to construct the mixed level designs..*/
data part1;
do k=1 to 3;
do i =1 to 3;
do j =1 to 4;
rep=k;
supp=i;
batch=j;
output;
end;
end;
end;
keep rep supp batch;
data part2;
input code @@;
datalines;
1 -2 -2 1 1 0 -1 0 2 -2 1 3
-1 -3 0 4 -2 4 0 3 4 0 -1 2
0 -4 1 0 -3 2 -2 2 0 2 2 1
;
data purity;
merge part1 part2;
run;
/* now I am trying to produce the output shown on page-561, table13.4
some basic ANOVA for split plot design..*/
proc glm data=purity;
class supp batch;
model code=supp batch(supp)/ss1;
test h=supp e=batch(supp);
output out=resout r=resid p=predict;
/*proc print; you can check that this is the output given on page 563*/
run;
/* this will produce the output shown on page 562, table-13.6*/
proc mixed data=purity method=type1 noinfo cl;
class supp batch;
model code=supp;
random batch(supp);
run;
/* following is the code to produce output shown in figure-13.3 , Of course
I would like to check for the normality of the residuals..*/
proc plot data=resout;
plot resid*(predict supp)="*"/vref=0;
run;
proc univariate data=resout noprint;
probplot resid/normal(mu=est sigma=est);
inset normal;
run;
Example-13.2
/* Design with both nested and factorial factors..
Example-13.2, Assembly time data */
data part1;
do k=1 to 3;
do r=1,2;
do i=1,2;
do j=1 to 4;
fix = k; rep=r; layout=i; oper=j;
output;
end;
end;
end;
end;
keep rep layout oper fix;
data part2;
input time @@;
datalines;
22 23 28 25 26 27 28 24 24 24 29 23 28 25 25 23
30 29 30 27 29 30 24 28 27 28 32 25 28 27 23 30
25 24 27 26 27 26 24 28 21 22 25 23 25 24 27 27
;
data assem;
merge part1 part2;
run;
/* table on page-13.11 ANOVA output..*/
proc glm data=assem;
class fix layout oper;
model time= fix|layout fix|oper(layout)/ss1;
test h=fix e=fix*oper(layout);
test h=layout e=oper(layout);
test h=fix*layout e=fix*oper(layout);
output out=resout1 r=resid1 p=predict1;
run;
/* code for the output similar to the output shown in table-13.13 page-572
Note: till this point of time I have no idea about how to do the analysis
for "RESTRICTED MODEL" */
proc mixed data=assem noinfo method=type1 cl;
class fix oper layout;
model time= fix|layout;
random oper(layout) fix*oper(layout);
run;
Exercise-13.20
/* Exercise-13.20 Pigment dispersion in paint.
whole plot -> mix
split-plot -> application method. */
data part1;
do i=1 to 3;
do j=1 to 3;
do k=1 to 4;
day=i; method=j; mix=k;
output;
end;
end;
end;
keep day method mix;
data part2;
input disp @@;
datalines;
64.5 66.3 74.1 66.5 68.3 69.5 73.8 70.0 70.3 73.1 78.0 72.3
65.2 65.0 73.8 64.8 69.2 70.3 74.5 68.3 71.2 72.8 79.1 71.5
66.2 66.5 72.3 67.7 69.0 69.0 75.4 68.6 70.8 74.2 80.1 72.4
;
data paint;
merge part1 part2;
run;
proc mixed data=paint method=type3 cl noinfo;
class mix method day;
model disp=day|method;
random day*mix|mix*method;
run;
proc glm data=paint;
class mix method day;
model disp=mix|method|day/ss3;
test h=day e=mix*day;
test h=method e=method*mix;
test h=method*day e=method*mix*day;
output out=resout2 r=resid2 p=predict2;
run;
proc plot data=resout2;
plot resid2*predict2="*"/vref=0;
run;