nlme, augPred, IV infusion

All,

I'm using nlme to estimate a nonlinear mixed-efffects model for pharmacokinetic data resulting from an IV infusion.   The fixed-effects form is a hard-coded equation used for IV infusions (open compartment, linear elimination).  The estimates seem to be okay, and a manual plot of the fixed effects curve is also okay, but the usual plot of the augPred does not seem to work.  Specifically, the fixed effects curve is not constant across subjects and the curve doesn't correspond to the estimated fixed effects; the EB estimated curves are also off.   Any suggestions much appreciated. Reproducible example below.

Cheers,
David



--
David Afshartous, Ph.D.
Research Associate Professor
PSTAT?: ASA Accredited Professional Statistician
Department of Biostatistics
Vanderbilt University Medical Center


library(nlme)
conc <- c(32.04, 50.82, 40.2, 53.12, 38.37, 61.36, 143.8, 27.42, 51.13, 37.14, 27.99, 136.4, 34.36, 135, 156.9, 52.92, 62.57, 56.02, 98.66, 76.36, 64.07, 47.92, 23.99, 93.42, 43.22, 41.45, 107.6, 79.74, 81.21, 49.09, 47.94, 108, 29.04, 32.56, 64.4, 28, 67.39, 60.24, 38.68, 33.4, 104.3, 48.14, 24.29, 10.7, 12.8, 156.8, 119.4, 51.23, 15.61, 16.75, 26.26, 58.16, 65.2, 57.65, 105.1, 104.2, 111.4, 125.1, 96, 94.1, 272.2, 153.8, 87.51, 74.75, 115.6, 253, 75.22, 229.2, 295.5, 129.9, 210.3, 124.8, 153.7, 142.9, 117.8, 99.37, 106.5, 162.6, 65.8, 105.6, 138.3, 151.7, 220.9, 93.46, 117.9, 152.7, 98.63, 130.3, 115.7, 66.49, 109, 90.64, 68.2, 55.67, 149.4, 100.5, 73.34, 48.35, 75.8, 172.6, 157.6, 129.6, 56.63, 77.57, 89.74, 204.2, 145.8, 107.7, 63.74, 45.59, 48.06, 63.76, 104.8, 44.84, 49.91, 177.5, 91.75, 150.6, 69.72, 48.21, 53.71, 16.51, 83.85, 43.39, 98.15, 65.3, 70.5, 45.15, 52.6, 81.9, 25.08, 32.37, 45.75, 61.8, 47, 34.78, 32.02, 85.34, 55.24, 17.73, 15.34, 14.53, 130.1, 96.41, 78.6, 16.93, 40.98, 68.38, 68.93, 68, 54.83)
Tinf.ind <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1)
Dose <- c(3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 3000, 2000, 3000, 2000, 2000, 3000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 3000, 2000, 3000, 2000, 2000, 3000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 3000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 3000, 2000, 3000, 2000, 2000, 3000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 3000)
Time <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 8, 8, 12, 8, 8, 8, 12, 8, 8, 8, 8, 12, 8, 8, 8, 8, 6, 8, 8, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 8, 6, 6, 6, 6)
ID <- c(19, 18, 25, 29, 39, 24, 53, 43, 17, 8, 28, 52, 10, 51, 54, 32, 49, 34, 42, 35, 36, 45, 21, 46, 11, 20, 33, 40, 50, 14, 27, 41, 15, 31, 26, 4, 23, 13, 5, 2, 38, 16, 6, 1, 7, 47, 44, 30, 3, 9, 12, 48, 37, 22, 19, 18, 25, 29, 39, 24, 53, 43, 17, 8, 28, 52, 10, 51, 54, 32, 49, 34, 42, 35, 36, 45, 21, 46, 11, 20, 33, 40, 50, 14, 27, 41, 15, 31, 26, 4, 23, 13, 5, 2, 38, 16, 6, 1, 7, 47, 44, 30, 3, 9, 12, 48, 37, 22, 18, 25, 29, 24, 53, 17, 28, 52, 51, 54, 32, 34, 45, 21, 46, 11, 33, 40, 50, 14, 27, 41, 15, 31, 26, 23, 13, 5, 2, 38, 16, 6, 1, 7, 47, 44, 30, 3, 9, 12, 48, 37, 22)
dat <- data.frame(ID, conc, Time, Tinf.ind, Dose)
dat <- groupedData(conc ~ Time | ID, data = dat,
       labels = list(x = "Time (hours)", y = "Concentration"))

## uses hard-coded functional form for IV infusion where Tinf.ind is indicator variable that
## determines which part of the function corresponds to the data point, pre versus post infusion:
mod1  <- nlme(conc ~ (1 - Tinf.ind) * (0 + (Dose/0.5) * (1/ (k*V)) * (1 - exp(- k * Time) ) )
           + (Tinf.ind) * (0 + (Dose/0.5) * (1/ (k*V)) * (1 - exp(- k * 0.5) ) * exp(-k* (Time - 0.5))),
      fixed = c(k + V  ~ 1),
     random = pdDiag(V~ 1),
     start = c(.3, 30),
     data = dat)
summary(mod1)

k <- fixef(mod1)[1]; V <- fixef(mod1)[2];
## manual plot of resulting fixed effects curve:
IV.function <- function(t)
  {
    (1 - (t >=.5 )) * (D/Tinf)*(1/(k*V)) * (1 - exp(-k * t)) + (t >=.5)*(D/Tinf)*(1/(k*V))*(1 - exp(-k*Tinf))*exp(-k*(t-Tinf))
  }
curve(IV.function, 0, 8)

## this does not seem to work:
plot(augPred(mod1, level=0:1),
     main = "IV infusion, nonlinear mixed-effects model")