/* Code to fit integrated form of Michaelis-Menton equ */

/* Philip Dixon, Dec 2000 - May 2001 */

/* the code below is not (at least yet) set up as a macro */
/*  if there is interest in this, let me know and I will convert it */

/* the current version assumes: */

/*    the X variable is named TIME */
/*    the Y variable is named CONC */
/*    groups are marked by the variable named POT */

/* the data set is assumed to be sorted BY POT TIME */
/*   that is, observations are arranged in order of increasing X */
/*   this speeds up the integration considerably, because previous */
/*    computations can be reused */

/* in various places, comments are used to indicate optional code */
/*   that can be implemented by uncommenting the indicated statements */


/* the code produces two data sets: */
/*   coeff: estimated coefficients for each pot */
/*   resids: residuals and predicted values for each pot */
    
/* the data step at the end of the file processes the coeff */
/*   data set into a more usuable form with */
/*   estimates and their se's for each coefficient and pot */
    
proc nlin method = dud outest=coeff;

   by pot;

   retain c t;

   deltat = 1;           /* time increment for integrator */

/* uncomment the next two lines to fix the starting concentration */
/*  c0 = 370;            /* Fixed starting conc */
/*  parms i = 0.25 0.5 0.75 1.0 km= 50 100 200 cmin=5 30 45;  */
   
/* current version estimates starting concentration */   
   parms c0= 300 350 400  i = 0.25 0.5 0.75 1.0 km= 50 100 200 cmin=5 30 45;

   bounds i > 0,   cmin < km < c0, cmin > 0;

   if (_obs_ = 1) and (c = .) then do;
      c = c0;
          t = 0;
          end;

   do until (t >= time);
      dcdt = i*(c-cmin) / (km + c-cmin);
          c = c - dcdt*deltat;
          t = t + deltat;

          end;

   if (t > time) then do;        
         /* if overshoot needed time, then back off */
       c = c + dcdt*(t-time);
       t = time;
       end;

   model conc = c;

/* the current version assumes constant error variance */
/* uncomment the following line if Var(error) is proportional */
/*   to the fitted concentration  */
/*   _weight_ = 1/c;         /* variance proportional to conc */

   output out=resids p=Chat r=resid;

run;


data coeff2;
  set adams.coeff;
  retain c0hat ihat kmhat cminhat c0se ise kmse cminse;
  by pot;

  if _type_ = 'FINAL' then do;
     c0hat = c0;
     ihat = i;
     kmhat = km;
     cminhat = cmin;
     end;
  if _type_ = 'COVB' then do;
     if _NAME_ = 'c0' then c0se = sqrt(c0);
     if _NAME_ = 'i' then ise = sqrt(i);
     if _NAME_ = 'km' then kmse = sqrt(km);
     if _NAME_ = 'cmin' then cminse = sqrt(cmin);
     end;

  if last.pot then output;

  keep pot c0hat -- cminse _SSE_;

proc print data=coeff2;
  title 'Estimates and ses from integrated fits';
  var pot _sse_ c0hat c0se ihat ise kmhat kmse cminhat cminse;
 run;