options nodate nonumber ps=55 ls=80;

title 'RPD -- Example 9.4 -- Listening-Reading Skills Test';

data listen;
input trt pretest posttest;
datalines;
1  89  87
1  82  86
1  88  94
1  94  96
2  89  84
2  90  94
2  91  97
2  92  93
3  89  96
3  99  97
3  84 100
3  87  98
;


proc iml;

use listen;
read all;

Y = posttest;

   /* Full model --- slopes and intercepts differ by trt: */

X1 = design(trt)||design(trt)#pretest;

   /* the operator # multiplies a matrix and a vector "elementwise",
      replicating the vector to match the size of the matrix */

print X1;

XTXINV1 = inv(X1`*X1);

BETAHAT1 = XTXINV1*(X1`*Y);

SSRES1 = Y`*Y - BETAHAT1`*X1`*Y;

S2_1 = SSRES1/(nrow(X1)-ncol(X1));

print SSRES1;


   /* Reduced model with same slopes, different intercepts: */

X2 = design(trt)||pretest;

print X2;

XTXINV2 = inv(X2`*X2);

BETAHAT2 = XTXINV2*(X2`*Y);

SSRES2 = Y`*Y - BETAHAT2`*X2`*Y;

S2_2 = SSRES2/(nrow(X2)-ncol(X2));

print SSRES2;


   /* Test for equal slopes: */

F0 = (SSRES2 - SSRES1)/(ncol(X1) - ncol(X2)) / S2_1;

PVALF0 = 1-probf(F0,ncol(X1)-ncol(X2),nrow(X1)-ncol(X1));

print F0 PVALF0;


   /* Further reduced model with same slopes and intercepts: */

X3 = j(nrow(Y),1)||pretest;

print X3;

XTXINV3 = inv(X3`*X3);

BETAHAT3 = XTXINV3*(X3`*Y);

SSRES3 = Y`*Y - BETAHAT3`*X3`*Y;

print SSRES3;


   /* Test for equal intercepts, assuming equal slopes: */

F0 = (SSRES3 - SSRES2)/(ncol(X2) - ncol(X3)) / S2_2;

   /* using s^2 from the equal-slopes model as the denominator, but could
      alternatively use s^2 from the full model (and change df below) */

PVALF0 = 1-probf(F0,ncol(X2)-ncol(X3),nrow(X2)-ncol(X2));

print F0 PVALF0;


   /* Joint test for equal intercepts and equal slopes: */

F0 = (SSRES3 - SSRES1)/(ncol(X1) - ncol(X3)) / S2_1;

PVALF0 = 1-probf(F0,ncol(X1)-ncol(X3),nrow(X1)-ncol(X1));

print F0 PVALF0;

quit;


proc glm;
class trt;
model posttest = trt pretest trt*pretest;
contrast 'equal slopes?' trt*pretest  1 -1  0,
                         trt*pretest  0  1 -1;
contrast 'equal slopes and intercepts?' trt          1 -1  0,
                                        trt          0  1 -1,
                                        trt*pretest  1 -1  0,
                                        trt*pretest  0  1 -1;
run;


proc glm;
class trt;
model posttest = trt pretest;
contrast 'equal intercepts (given equal slopes)?' trt  1 -1  0,
                                                  trt  0  1 -1;
run;