options nodate nonumber ps=55 ls=80;

title 'RPD -- Example 9.2 -- Numerical Example';

data artificial;
input group y;
datalines;
1  8.90
1  8.76
2 11.78
2 12.07
3 14.50
3 12.48
4 16.79
4 16.57
;


proc iml;

use artificial;
read all;

   /* vector Y now contains the responses */


/* Cell Means Model: */

X1 = design(group);

   /* 'design' creates a matrix with an indicator variable column for each
      level of group */

XTXINV1 = inv(X1`*X1);

BETAHAT1 = XTXINV1*(X1`*Y);

print X1 XTXINV1 BETAHAT1;


/* Effects Model: Sum-to-Zero */

X2 = j(nrow(Y),1)||designf(group);

   /* 'designf' creates a matrix with columns for effects (of the levels of 
      group) under the sum-to-zero restriction, after substituting out the
      effect for the last level */

XTXINV2 = inv(X2`*X2);

BETAHAT2 = XTXINV2*(X2`*Y);

TAUHAT4 = - BETAHAT2[2] - BETAHAT2[3] - BETAHAT2[4];

print X2 XTXINV2 BETAHAT2 TAUHAT4;


/* Effects Model: Reference Cell */

X3 = j(nrow(Y),1)||X1[,1:3];
   
   /* 'X1[,1:3]' is the first three columns of X1 */

XTXINV3 = inv(X3`*X3);

BETAHAT3 = XTXINV3*(X3`*Y);

print X3 XTXINV3 BETAHAT3;

quit;


proc glm;
class group;
model y = group / solution;
   /* gives "reference cell model" estimates */
run;


proc glm;
class group;
model y = group / noint solution;
   /* gives "cell means model" estimates */
   /* Note: The listed ANOVA F-test is for H0: all means are ZERO. */
run;


/* PROC GLM does not automatically produce "effects model" estimates that 
   satisfy the sum-to-zero restriction */