/*  SAS code for analzing the ice crsytal
    growth data stored in the file

      c:\courses\st511\sas\crystal.dat   */


data set1;
   infile 'c:\courses\st511\sas\crystals.dat';
   input time length;
   run;

/* Fit a straight line model  */

proc reg data=set1;
  model length = time / ss1 ss2;
  output out=set2 residual=resid predicted=yhat;
  run;

/* Plot the observations and the estimated line  */

 goptions cback=white colors=(black)
    targetdevice=ps300 rotate=landscape;

   axis1 label=(f=swiss h=2.5)
      value=(f=swiss h=2.0)  w=3.0  length= 6.5 in;

   axis2 label=(f=swiss h=2.2 a=270 r=90)
      value=(f=swiss h=2.0) w= 3.0  length = 5.5 in;

       SYMBOL1 V=circle H=2.0 w=3 l=3 i=none ;
       SYMBOL2 V=none H=2.0 w=3 l=1 i=spline ;
    

 proc gplot data=set2;
   plot (length yhat)*time / overlay vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Ice Crystal Growth";
       label length='Axial Length ';
       label time = 'Time(seconds) ';
 run;


/* Construct a normal probability plot
   and other residual plots  */
   
proc rank data=set2 normal=blom out=set2;
  var resid;
  ranks q;
  run;

  axis1 label=(f=swiss h=2.5)
      value=(f=swiss h=2.0)  w=3.0  length= 6.5 in;

   axis2 label=(f=swiss h=2.2 a=270 r=90)
      value=(f=swiss h=2.0) w= 3.0  length = 5.5 in;

       SYMBOL1 V=circle H=2.0 w=3 l=1 i=none ;
 
 proc gplot data=set2;
   plot resid*time / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Residuals vs. Time";
       label resid='Residuals ';
       label time = 'Time(seconds) ';
 run;
 
 proc gplot data=set2;
   plot resid*q / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Normal Probability Plot";
       label resid='Residuals ';
       label q = 'Standard Normal Quantiles ';
 run;


/* Fit a quadratic model  */

data set1; set set1;
  time2=time*time;
  run;

proc reg data=set1;
  model length = time time2/ ss1 ss2;
  output out=set3 residual=resid predicted=yhat;
  run;

/* Plot the observations and the estimated line  */

 goptions cback=white colors=(black)
    targetdevice=ps300 rotate=landscape;

   axis1 label=(f=swiss h=2.5)
      value=(f=swiss h=2.0)  w=3.0  length= 6.5 in;

   axis2 label=(f=swiss h=2.2 a=270 r=90)
      value=(f=swiss h=2.0) w= 3.0  length = 5.5 in;

       SYMBOL1 V=circle H=2.0 w=3 l=3 i=none ;
       SYMBOL2 V=none H=2.0 w=3 l=1 i=spline ;
    

 proc gplot data=set3;
   plot (length yhat)*time / overlay vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Ice Crystal Growth";
       label length='Axial Length ';
       label time = 'Time(seconds) ';
 run;


/* Construct a normal probability plot
   and other residual plots  */
   
proc rank data=set3 normal=blom out=set3;
  var resid;
  ranks q;
  run;

  axis1 label=(f=swiss h=2.5)
      value=(f=swiss h=2.0)  w=3.0  length= 6.5 in;

   axis2 label=(f=swiss h=2.2 a=270 r=90)
      value=(f=swiss h=2.0) w= 3.0  length = 5.5 in;

       SYMBOL1 V=circle H=2.0 w=3 l=1 i=none ;
 
 proc gplot data=set3;
   plot resid*time / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Residuals vs. Time";
       label resid='Residuals ';
       label time = 'Time(seconds) ';
 run;
 
 proc gplot data=set3;
   plot resid*q / vaxis=axis2 haxis=axis1;
   title  H=3.0 F=swiss "Normal Probability Plot";
       label resid='Residuals ';
       label q = 'Standard Normal Quantiles ';
 run;


data set1; set set1;
   time2 = time*time;
   run;

proc reg data=set1;
  model length = time time2/ ss1 ss2;
  output out=set2 residual=resid predicted=yhat;
  run;


/* Peform a goodness of fit test 
   for the quadratic model  */


data set1; set set1;
   group=time;
   run;

proc glm data=set1;
  class group;
  model length = time time2 group/ ss1 ss2;
  run;