options nodate nonumber ps=55 ls=80;

title 'RPD -- Example 8.2 -- 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
;


symbol value=dot
       interpol=none;

proc gplot;
plot density*day;
run;


proc iml;

use algae;
read all;

Y = density;

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

   /* Try a cubic polynomial fit: */

X = ONES||day||day##2||day##3;   /* ## is elementwise power */

XTXINV = inv(X`*X);

BETAHAT = XTXINV*(X`*Y);

SSRES = Y`*Y - BETAHAT`*X`*Y;

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

SE_BETAHAT = sqrt(vecdiag(XTXINV)#S2);   /* estimated standard errors */

print X BETAHAT SE_BETAHAT;


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

T0 = BETAHAT[4]/SE_BETAHAT[4];

PVALT0 = 2*(1-probt(abs(T0),nrow(X)-ncol(X)));   /* probt(t,d) = c.d.f. at t 
                                                    of t(d)-distribution */

print T0 PVALT0;


   /* To test whether a straight-line model is adequate ... */

X_R = ONES||day;

XTXINV_R = inv(X_R`*X_R);

BETAHAT_R = XTXINV_R*(X_R`*Y);

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

F0 = ((SSRES_R-SSRES)/2)/S2;

PVALF0 = 1-probf(F0,2,nrow(X)-ncol(X));

print F0 PVALF0;

quit;


proc glm;
model density = day day*day day*day*day / solution;
contrast 'quadratic adequate?' day*day*day 1;
contrast 'simple linear adequate?' day*day 1, day*day*day 1;
run;


/* To use PROC REG, the polynomial terms must be in the data set: */

data algae2;
set algae;      /* read records from existing data set 'algae' */
day2 = day**2;
day3 = day**3;
run;


proc reg;
model density = day day2 day3;
test day3 = 0;
test day2 = 0, day3 = 0;
run;