/*  SAS code for fitting a non-linear model to the
    weight loss data from Venables & Ripley  */

data set1;
  infile 'wtloss2.dat';
  input case time y;
  label  y = ' Weight (kg)'
         time = 'Time (days)';
  run;


proc print data=set1;
  run;


/*  Use SASGRAPH to plot the data  */



/* goptions device=WIN target=WINPRTC rotate=landscape; */

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

  axis1 label = (h=3 r=0 a=90 f=swiss 'Weight (kg)')
        value = (h=2.0 f=swiss)
        length = 5 in  width=8.0 style=1;
        
  axis2 label = (h=2.5 f=swiss 'Time (days)')
        value = (h=2 f=swiss)
        length = 5.5 in  width=8.0 style=1
        order = 0 to 250 by 50;

  symbol1 v=circle i=none h=2 width=5;
  
proc gplot data=set1;
  plot y*time / vaxis=axis1 haxis=axis2;
  title h=2.5 f=swiss c=black 'Weight Loss from an Obese Patient';
  run;
  
  
/*  Fit a model using the Gauss-Newton algorithm
    where formulas for derivatives are given.
    A grid search is used to find starting values.  */

proc nlin data=set1 method=gauss;
  parms b0 = 50 to 100 by 5
        b1 = 75 to 125 by 5
        b2 = 100 to 200 by 10;
  temp = 2**(-(time/b2));
  model y = b0 + b1*temp;
  bounds b0>1e-20,  b1>1e-20, b2>1e-20;
  der.b0 = 1;
  der.b1 = temp;
  der.b2 = temp*(b1*time*log(2))/(b2**2);
  output out=set2 predicted=yhat residual=r student=sr
          l95 u95 l95m u95m stdp=stdp stdr=stdr;
  run;


/* Plot the studentized residuals  */

proc univariate data=set2 normal;
  var sr;  run;

proc rank data=set2 normal=blom out=set2;
  var sr; ranks q;
  run;



  axis1 label = (h=3 r=0 a=90 f=swiss ' ')
        value = (h=2.0 f=swiss)
        length = 5 in  width=8.0 style=1;

  axis2 label = (h=2.5 f=swiss 'Standard Normal Quantiles')
        value = (h=2 f=swiss)
        length = 5.5 in  width=8.0 style=1;

  axis3 label = (h=2.5 f=swiss 'Time (days)')
        value = (h=2 f=swiss)
        length = 5.5 in  width=8.0 style=1;

  symbol1 v=circle i=none h=2 width=5;
  symbol2 v=none i=spline width=8 line=1;

proc gplot data=set2;
  plot sr*q / vaxis=axis1 haxis=axis3;
  title h=3.0 f=swiss c=black 'Studentized Residuals';
  run;

proc gplot data=set2;
  plot sr*q / vaxis=axis1 haxis=axis2;
  title h=3.5 f=swiss c=black 'Studentized Residuals';
  run;


  axis4 label = (h=3 r=0 a=90 f=swiss 'Weight (kg)')
        value = (h=2.0 f=swiss)
        length = 5 in  width=8.0 style=1;

proc gplot data=set2;
  plot (y yhat)*time / overlay vaxis=axis4 haxis=axis2;
  title h=2.5 f=swiss 'Weight Loss from an Obese Patient';
  run;