/*3-D Polynomial Response Surface, Test for Lack of Fit**/

/*------------------------------------------------------------------*
|  Use moisture and temperature to predict germination percentage   |
| Data uses 3 temperatures and 3 moistures in all combinations.     |
| Middle moisture and temperature used multiple times so we can     |
| do a lack of fit test.   There are a total of 12 observations.    |
| We consider the following quadratic surface                       |
|    mu(m,t)=beta0+beta1*m+beta2*m^2+beta3*t+beta4*t^2+beta5*m*t    | 
*-----------------------------------------------------------------*/
Data response;
  input  germpct moist temp trt @@;
     *xtemp = temp+.0000*ranuni(1928766);*  <-- no jittering;
     xtemp = temp+.0001*ranuni(1928766);*  <-- imperceptible jittering
                  of temp avoids one point just above another
                  (SAS would not plot such points in G3D) ;
   m2 = moist*moist; t2=temp*temp; mt = moist*temp;
   svar = 1.5  ;   *<---- A size variable for graphing;
   cvar = 'green';  *<---- A color variable ;
cards;
  14   90 110 1    42  60 110 2      30  30 110 3

  42   90  90 4    85  60  90 5      51  30  90 6
                   92  60  90 5
                   84  60  90 5
                   78  60  90 5

  46   90  70 7    38  60  70 8      17  30  70 9
    ;

Data grid;  ** create a grid of temperature and moisture values for
               the floor of our plot;
 do moist = 30 to 90 by 2;  do temp = 70 to 110 by 2;
   m2 = moist*moist; t2=temp*temp; mt = moist*temp;
   svar = 0.5; cvar = 'cyan';
   output; end; end;

data all; set response grid;  proc print data=all(obs=20);  /* note obs= option */
  title "First 20 observations in extended data set";

proc reg; 
  title "PROC REG to fit quadratic surface";
  model germpct = moist temp t2 m2 mt;
  output out=out1 predicted = G_hat;  /*G_hat=predicted germpct */

Data out1; set out1; if germpct = . then germpct = G_hat;
  If xtemp ne . then temp = xtemp;  *<-- jittered version of temp;
  drop m2 t2 mt;
  ** Note germpct now has actual data followed by predictions over
     the moisture, temp grid we created;

proc print data=out1 (Obs = 30);
  title "Observed values and regression predictions";

proc g3d; 
 title "Observed data and fitted response surface";
scatter moist*temp=germpct / shape = "balloon" size = svar
   color = cvar noneedle;   
/* noneedle omits vertical lines connecting floor to datapoints */
/* omitting this option leads to a pincushion-like  plot */
quit;

proc glm data=response; class trt; model germpct = trt;
 title "Getting pure error SS - can you obtain lack of fit F test?";
 title2 "Make sure the df make sense to you";

proc rsreg data=response; 
model germpct = moist temp;
*model germpct = moist temp/lackfit;  
/* pure-error should be based only on the replicated design point,  right? */
run;
/* in the output, find the estimated combo of MOIST and TEMP at which GERMPCT is maximized */