nlme, augPred, IV infusion

PS: slight correction to below: the fixed effects curve will not be constant since there will be a dependence on covariates, but the overall shape from plotting the augPred is still way off

--
David Afshartous, Ph.D.
Research Associate Professor
PSTAT?: ASA Accredited Professional Statistician
Department of Biostatistics
Vanderbilt University Medical Center
________________________________________
From: r-sig-mixed-models-bounces at r-project.org [r-sig-mixed-models-bounces at r-project.org] on behalf of Afshartous, David [d.afshartous at vanderbilt.edu]
Sent: Saturday, February 09, 2013 4:14 PM
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] 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")



        [[alternative HTML version deleted]]


_______________________________________________
R-sig-mixed-models at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models