/*  Fit using SAS PROC NLIMIXED  */

options ps=55 ls=80 nodate;

/*  input the data */

data arg; infile 'argconc.dat';
  input obsno indiv dose time conc;
run;

/*  create variables for defining the PK model */

data arg; set arg;
  tinf=240;
  t1=1;
  if time>tinf then t1=0;
  t2=tinf*(1-t1)+t1*time;
run;

/* put method=firo in the proc nlmixed statement to get first order
   approximation.  Default is adaptive Gaussian quadrature. Achieving
   convergence is not straightforward - we tried different optimization
   methods and starting values, using more abcissae than the default (=61),
   and requiring a more stringent convergence criterion than the default.
   See the documentation for details.

.  Here, we show using the Trust Region optimization method rather than
   the default Dual Quasi-Newton method.  */

/* first, fix theta at 0.22 as estimated based on individual estimates */

proc nlmixed data=arg method=gauss tech=trureg qmax=101 qtol=1e-5;
  parms beta1=-5.5 beta2=-2.0 s2b1=0.14 cb12=0.006 s2b2=0.006 sig=23.0;
  cl=exp(beta1+b1);
  v=exp(beta2+b2);
  pred=(dose/cl)*(1-exp(-cl*t2/v))*exp(-cl*(1-t1)*(time-tinf)/v);
  var=(sig**2)*(pred**(2*0.22));
  model conc ~ normal(pred,var);
  random b1 b2 ~ normal([0,0],[s2b1,cb12,s2b2])
     subject=indiv;
run;

/* allow theta to be estimated */

proc nlmixed data=arg method=gauss tech=trureg qmax=101 qtol=1e-5;
  parms beta1=-5.5 beta2=-2.0 s2b1=0.14 cb12=0.006 s2b2=0.006 sig=23.0
     theta=0.22;
  cl=exp(beta1+b1);
  v=exp(beta2+b2);
  pred=(dose/cl)*(1-exp(-cl*t2/v))*exp(-cl*(1-t1)*(time-tinf)/v);
  var=(sig**2)*(pred**(2*theta));
  model conc ~ normal(pred,var);
  random b1 b2 ~ normal([0,0],[s2b1,cb12,s2b2])
     subject=indiv;
run;