/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 17.8                                              */
/*   TITLE: SAS Tasks 17.1 & 17.2 - Nonlinear fitting                 */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "SAS Task 17.1 & 17.2 - Nonlinear fitting";

/*
/  Data from Table 17.1, amount of nitrogen in US lakes.
/     x1 = avg influent N concentration
/     x2 = water retention time
/     y  = mean annual N concentration
/---------------------------------------------------------------------*/
data Lakes;
   input x1 x2 y;
   obs = _N_;
datalines;
 5.548  0.137  2.590
 4.896  2.499  3.770
 1.964  0.419  1.270
 3.586  1.699  1.445
 3.824  0.605  3.290
 3.111  0.677  0.930
 3.607  0.159  1.600
 3.557  1.699  1.250
 2.989  0.340  3.450
18.053  2.899  1.096
 3.773  0.082  1.745
 1.253  0.425  1.060
 2.094  0.444  0.890
 2.726  0.225  2.755
 1.758  0.241  1.515
 5.011  0.099  4.770
 2.455  0.644  2.220
 0.913  0.266  0.590
 0.890  0.351  0.530
 2.468  0.027  1.910
 4.168  0.030  4.010
 4.810  3.400  1.745
34.319  1.499  1.965
 1.531  0.351  2.555
 1.481  0.082  0.770
 2.239  0.518  0.720
 4.204  0.471  1.730
 3.463  0.036  2.860
 1.727  0.721  0.760
;

/*
/  Fitting all data.  Save the parameter estimates for examining
/  later.
/---------------------------------------------------------------------*/
title2 "SAS Task 17.1 - Fitting all data";
proc nlin data=Lakes outest=ParmAll(where=(_TYPE_="FINAL")) hougaard;
   parameters t1 = 1 t2 = 1;
   model y = x1 / (1 + t1 * x2**t2);
run;
title2;

/*
/  Fitting all data except for the observations with extreme outliers
/  in x1.  Save the parameter estimates for examining later.
/---------------------------------------------------------------------*/
title2 "SAS Task 17.1 - Dropping outliers in x1";
proc nlin data=Lakes outest=ParmSub(where=(_TYPE_="FINAL"));
   where ^((obs = 10) | (obs = 23));
   parameters t1 = 1 t2 = 1;
   model y = x1 / (1 + t1*x2**t2);
run;
title2;

/*
/  Combine and compare parameter estimates.
/---------------------------------------------------------------------*/
title2 "SAS Task 17.1 - Comparing parameter estimates, SSE";
data ParmAll; set ParmAll; Source = "All data    ";
data ParmSub; set ParmSub; Source = "Drop 10 & 23";
data Parm; set ParmAll ParmSub;
proc print data=Parm noobs;
   var Source t1 t2 _SSE_;
run;
title2;



title2 "SAS Task 17.2 - Estimating maximum concentration and time to it";
data Theophylline;
   input Time Concentration;
   Time = round(60*Time)/60;
   cards;
 0.166 10.1
 0.333 14.8
 0.5   19.9
 0.666 22.1
 1     20.8
 1.5   20.3
 2     19.7
 2.5   18.9
 3     17.3
 4     16.1
 5     15.0
 6     14.2
 8     13.2
10     12.3
12     10.8
24      6.5
30      4.6
48      1.7
;
proc nlmixed data=Theophylline;
   parameters t1=1 t2=1 t3=1;
   model Concentration
            ~ normal(t3*(exp(-t1*Time) - exp(-t2*Time)),s2);
   Estimate "MaxTime" (log(t2) - log(t1))/(t2 - t1);
   Estimate "MaxConc" 
      t3*(  exp(-t1*((log(t2) - log(t1))/(t2 - t1)))
          - exp(-t2*((log(t2) - log(t1))/(t2 - t1))));
run;
title2;

title1;