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

data set1;
  infile 'c:\courses\st511\snew\pnonlin.dat';
  input y x1 x2;
  run;

proc print data=set1; run;

/* Get starting values  */

data set2; set set1;
  logy = log(y);
  logx1 = log(x1);
  logx2 = log(x2);
  run;

proc reg data=set2 outest=set3;
  model logy = logx1 logx2 ;
  run;

proc print data=set3; run;


/* Add an extra case to the data file to 
   get an estimated mean for x1=50 and x2=60
   and 95% coonfidence interval for the mean
   response and a 95% prediction interval */

data temp; 
  input y x1 x2;
  cards;
  . 50 60 
  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 = 9.591 b1=0.51485 b2=0.29845;
  model y = b0*(x1**b1)*(x2**b2);
  der.bo = (x1**b1)*(x2**b2);
  der.b1 = b0*(x1**b1)*(x2**b2)*log(x1);
  der.b2 = b0*(x1**b1)*(x2**b2)*log(x2);
  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 x1 x2 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  "Temperature(C)")  
      value=(f=swiss h=1.8)  w=3.0  
      length= 6.5 in;

axis2 label=(f=swiss h=1.8  "Pressure(psi)")  
      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*x1 / vaxis=axis3 haxis=axis1;
   TITLE1 ls=0.01in H=2.0 F=swiss "Residuals vs. Temperature";
	   footnote ls=0.01in;
   RUN;

PROC GPLOT DATA=set4; 
   PLOT r*x2 / vaxis=axis3 haxis=axis2;
   TITLE1 ls=0.01in H=2.0 F=swiss "Residuals vs. Pressure";
	   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;

 

/* Use the estimated parameters and the covariance
   matrix produced by the outest=  option in the
   nlin procedure to obtain an approximate t-test
   for beta1=beta2  */

proc print data=setp; run;

data setpp; set setp;
  if(_TYPE_="FINAL" or _TYPE_="COVB");
  keep b0 b1 b2;
  run;

proc print data=setpp; run; 

proc iml;
  start b1vsb2;

use setpp;
read all into x;

use set1;
read all into w;
nn=nrow(w);

np=ncol(x);
nr=nrow(x);
b = x[1, ];
vcov = x[2:nr, ];

c={0 1 -1};

t=c*t(b)/sqrt(c*vcov*t(c));
at=abs(t);
df=nn-np;
pvalue=(1-probt(at,df))*2;
print t df pvalue;

finish b1vsb2;

run b1vsb2; 




/* Another way to test the hypothesis that beta1=beta2
   is to fit a model with the contraints that
   beta1=beta2 and perform an approximate f-test

    f=(SSE(restricted)-SSE(unrestricted))/(df)/(MSE(unrestricted))
*/

data set1; set set1;
  x12=x1*x2;
  run;

proc nlin data=set1 method=gauss convergeparm=.000001
          maxiter=50 outest=setpr;
  parms b0 = 10. b1=0.4 ;
  model y = b0*(x12**b1);
  der.bo = (x12**b1);
  der.b1 = b0*(x12**b1)*log(x12);
  run;


proc print data=setpr;  run;


data setpr2; set setpr;
 if(_TYPE_="FINAL");
    sser=_SSE_;
	keep sser;
 run;

data setp2; set setp;
  if(_TYPE_="FINAL");
   sse=_SSE_;
   keep sse;
   run;

data both; merge setp2 setpr2; run;

proc iml;
  start ff;
  use both;
  read all into sse;
  f = (sse[1,2]-sse[1,1])/(sse[1,1]/15);
  pvalue = 1-probf(f,1,15);
  print sse f pvalue;
  finish;
  run;

run ff;