options nodate nonumber ps=55 ls=80;

title 'RPD -- Example 8.3 -- Algae Density: CO2 only, Both Reps';

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
1   .184
2   .664
3  1.553
4  1.910
5  2.585
6  3.009
7  3.403
8  3.892
9  4.367
10 4.551
11 4.656
12 4.754
13 4.842
14 4.969
;


proc iml;

use algae;
read all;

Y = density;

ONES = j(nrow(Y),1);

   /* Reduced model is a quadratic polynomial: */

X_R = ONES||day||day##2;

XTXINV_R = inv(X_R`*X_R);

BETAHAT_R = XTXINV_R*(X_R`*Y);

SSRES_R = Y`*Y - BETAHAT_R`*X_R`*Y;

print SSRES_R;


   /* Full model is a 13-degree polynomial: */

X_F = orpol(day,13);

   /* orpol produces an appropriate X matrix for polynomial regression of the
      specified degree;  the polynomials are actually orthogonalized, but this
      has no effect on SSRES */

XTXINV_F = inv(X_F`*X_F);

BETAHAT_F = XTXINV_F*(X_F`*Y);

SSRES_F = Y`*Y - BETAHAT_F`*X_F`*Y;

print SSRES_F;   /* the "pure error" sum of squares */


   /* To test for lack of fit ... */

SSLF = SSRES_R - SSRES_F;

F0 = (SSLF/(13-2))/(SSRES_F/(nrow(Y)-14));

PVALF0 = 1-probf(F0,13-2,nrow(Y)-14);

print SSLF F0 PVALF0;

quit;


/* PROC REG cannot easily be used for testing lack of fit;
   PROC GLM would also be somewhat awkward to use, in this case;
   PROC RSREG can be used easily, provided that the reduced
   model is quadratic: */

proc rsreg;
model density = day / lackfit;
run;


symbol interpol=rqclm95   /* regression, quadratic,
                             confidence limits (for) mean, 95% (limits) */
       co=blue            /* confidence limits printed in blue */
       value=dot;

   /* For a cubic regression fit, use:  interpol=rcclm95  */

axis offset=(5 pct, 5 pct);

proc gplot;
plot density*day / haxis=axis;
run;