options nodate nonumber ps=55 ls=80;

title 'RPD -- Example 8.4 -- Algae Density: CO2 only, Rep 1';

data algae;
input day density;
datalines;
1   .530
2  1.183
3  1.603
4  1.994
5  2.708
6  3.006
7  3.867
8  4.059
9  4.349
10 4.699
11 4.983
12 5.100
13 5.288
14 5.374
;


proc iml;

use algae;
read all;

Y = density;

   /* Cubic orthogonal polynomial fit: */

Q = orpol(day,3);

   /* Note:  Q is not quite the same as the matrix in the
      textbook because the columns of Q have length 1. */

QTQINV = inv(Q`*Q);   /* should be the identity matrix */

GAMMAHAT = QTQINV*(Q`*Y);

SSRES = Y`*Y - GAMMAHAT`*Q`*Y;

S2 = SSRES/(nrow(Q)-ncol(Q));

reset fuzz;   /* set to display very small values as zero */

print Q QTQINV GAMMAHAT S2;


   /* To translate back to the usual polynomial coefficient estimates ... */

X = j(nrow(Y),1)||day||day##2||day##3;

R = Q`*X;

BETAHAT = solve(R,GAMMAHAT);

print R BETAHAT;


   /* To test whether a quadratic model is adequate ... */

T0 = GAMMAHAT[4]/sqrt(QTQINV[4,4]*S2);

PVALT0 = 2*(1-probt(abs(T0),nrow(Q)-ncol(Q)));

print T0 PVALT0;

quit;