/* Use the NLIN procedure in SAS to do
   the analysis for problem 3 on 
   assignment 9, Spring 2002  */

data set1;
  infile 'c:\courses\st511\snew\kinetics.dat';
  input y x;
  run;

proc print data=set1; run;

/* Get starting values  */

data set2; set set1;
  yinv = 1/y;
  xinv = 1/x;
  run;

proc reg data=set2 outest=set3;
  model yinv = xinv ;
  run;

proc print data=set3; run;


/* Add an extra case to the data file to 
   get an estimated mean for x=25
   and 95% coonfidence interval for the mean
   response and a 95% prediction interval */

data temp; 
  input y x;
  cards;
  . 25
  run;

data set1a; set set1 temp;  run;

/*  Use the NLIN procedure to fit the model
    Here we include formulas for first 
    partial derivatives and use the 
    Gauss-Newton procedure  */

proc nlin data=set1a method=gauss convergeparm=.000001
          maxiter=50 outest=setp ;
  parms b0 = 29.62 b1=13.45;
  model y = (b0*x)/(b1+x);
  der.bo = x/(b1+x);
  der.b1 = -(b0*x)/((b1+x)**2);
  output out=set4 residual=r predicted=yhat student=student
         stdi=stdi stdp=stdp l95=l95pred u95=u95pred 
         l95m=l95mean u95m=u95mean;
  run;

  proc print data=set4;  
    var y x r yhat l95pred u95pred l95mean u95mean;
    run;

/* Make residual plots  */

proc rank data=set4 normal=blom out=set4;
  var r; ranks q;
  run;

/* UNIX users can replace the following goptions
   line with

 filename graffile pipe 'lpr -Dpostscript';
 goptions gsfmode=replace gsfname=graffile
    targetdevice=ps300 rotate=landscape;
*/

goptions device=WIN target=WINPRTC rotate=landscape;


axis1 label=(f=swiss h=1.8  "Concentration")  
      value=(f=swiss h=1.8)  w=3.0  
      length= 6.5 in;

axis3 label=(f=swiss h=2.0 a=90 r=0 "Residuals")  
      value=(f=swiss h=1.8) w= 3.0  
      length = 5.0 in;

axis4 label=(f=swiss h=1.8  "q")
      order = -3 to 3 by 1 
      value=(f=swiss h=1.8)  w=3.0  
      length= 6.5 in;

 SYMBOL1 V=CIRCLE H=1.7 w=3 l=1 i=none ;
 SYMBOL2 V=DIAMOND H=1.7 w=3 l=3 i=join ;
 SYMBOL3 V=square H=1.7 w=3 l=9 i=join ;
 

 PROC GPLOT DATA=set4; 
   PLOT r*x / vaxis=axis3 haxis=axis1;
   TITLE1 ls=0.01in H=2.0 F=swiss "Residuals vs. Concentration";
	   footnote ls=0.01in;
   RUN;

 PROC GPLOT DATA=set4; 
   PLOT r*q / vaxis=axis3 haxis=axis4;
   TITLE1 ls=0.01in H=2.0 F=swiss "Normal Probability Plot";
	   footnote ls=0.01in;
   RUN;

 

;