|
DATA ANOVA;
INPUT Obs y group intercept g1 g2
g3 eff1 eff2;
CARDS;
1 57.607 1 1 1 0
0
-1 -1
2 72.872 1 1 1 0
0
-1 -1
3 59.590 1 1 1 0
0
-1 -1
4 62.096 1 1 1 0
0
-1 -1
5 35.548 1 1 1 0
0
-1 -1
6 49.867 1 1 1 0
0
-1 -1
7 53.850 1 1 1 0
0
-1 -1
8 42.697 1 1 1 0
0
-1 -1
9 57.699 1 1 1 0
0
-1 -1
10 55.889 1 1 1 0
0
-1 -1
11 96.261 2 1 0 1
0 1 0
12 112.327 2 1 0 1
0 1 0
13 89.650 2 1 0 1
0 1 0
14 95.550 2 1 0 1
0 1 0
15 106.094 2 1 0 1
0 1 0
16 95.218 2 1 0 1
0 1 0
17 73.611 2 1 0 1
0 1 0
18 92.811 2 1 0 1
0 1 0
19 84.821 2 1 0 1
0 1 0
20 101.658 2 1 0 1
0 1 0
21 -28.039 3 1 0 0 1
0 1
22 -30.638 3 1 0 0 1
0 1
23 -26.553 3 1 0 0 1
0 1
24 -21.488 3 1 0 0 1
0 1
25 -17.679 3 1 0 0 1
0 1
26 -26.164 3 1 0 0 1
0 1
27 -15.792 3 1 0 0 1
0 1
28 -22.087 3 1 0 0 1
0 1
29 -12.109 3 1 0 0 1
0 1
30 -57.962 3 1 0
0
1 0 1
;
RUN;
/*What do the various design matrices look like*/
/*Cell means coding*/
proc freq data=anova;
table g1*g2*g3 / list;
run;
/*Reference cell coding*/
proc freq data=anova;
table intercept*g2*g3 / list;
run;
/*Effect coding*/
proc freq data=anova;
table intercept*eff1*eff2 / list;
run;
/*Classic ANOVA coding*/
proc freq data=anova;
table intercept*g1*g2*g3 / list;
run;
/*Cell Means Coding -- each estimate is the mean for group x*/
proc glm data=anova;
model y=g1 g2 g3 / solution noint; /* notice there is no intercept in this model
*/
estimate 'g1 versus g2' g1 -1 g2 1 g3 0;
estimate 'g1 versus g3' g1 -1 g2 0 g3 1;
run;
quit;
/* HOW ARE STANDARD ERRORS COMPUTED */
proc iml;
use anova;
read all var {y} into
y;
read all var {g1 g2 g3}
into x;
close anova;
DF=NROW(Y)-NCOL(X);
I=I(30);
SSE=Y`*(I-X*INV(X`*X)*X`)*Y;*EXPANSION OF THIS EXPRESSION YIELDS: SUM(Y-E[Y])^2 -- THE SUM OF
SQUARES, WHERE E[Y]=X*(INV(X`*X)*X`*Y);
SIGMA=SSE/DF;
STD_ERROR=VECDIAG(ROOT(INV(X`*X)*SIGMA));
PRINT STD_ERROR; *COMPARE
THESE VALUES TO PROC glm OUTPUT, THEY ARE
IDENTICAL;
quit;
/*Reference Cell Coding -- Intercept is reference group mean
and each estimate is the difference between reference mean and mean
for group x*/
proc glm data=anova;
model y= g2
g3 / solution;
estimate 'm1' intercept 1 g2 0 g3 0;*Mean of the first group;
estimate 'm2' intercept 1 g2 1 g3 0;*Mean of the second group;
estimate 'm3' intercept 1 g2 0 g3 1;*Mean of the third group;
run;
quit;
/*Effect Coding -- Intercept is the mean of all group means
and each estimate is the difference between reference mean and mean
for group x*/
proc glm data=anova;
model y=
eff1 eff2 / solution;
estimate 'm1' intercept 1 eff1 -1 eff2 -1;*Mean of the first group;
estimate 'm2' intercept 1 eff1 1 eff2 0;*Mean of the second
group;
estimate 'm3' intercept 1 eff1 0 eff2 1;*Mean of the third
group;
/* So even in an effect coding scheme you can get cell mean results
*/
run;
quit;
/*Using the class statement results in the classic less than full
rank ANOVA coding with the last group as reference*/
proc glm data=anova;
class group;
model y=group / solution;
estimate 'g1 versus g2' group -1 1 0;
estimate 'g1 versus g3' group -1 0 1;
run;
quit;
|