Pharmacokinetics Estimation Example

pharmacokinetics.R

x <- c(0.27, 0.58, 1.02, 2.02, 3.62, 5.08, 7.07, 9.00, 12.15, 24.17)
y <- c(4.40, 6.90, 8.20, 7.80, 7.50, 6.20, 5.30, 4.90, 3.70, 1.05)

# function
f <- function(beta, x) {
  D <- 4.53
  V <- beta[1]
  k.a <- beta[2]
  k.e <- beta[3]
  D * k.a / (V * (k.a - k.e)) * (exp(-k.e * x) - exp(-k.a * x))
}

# gradient
g <- function(beta, x) {
  D <- 4.53
  V <- beta[1]
  k.a <- beta[2]
  k.e <- beta[3]
  # come up a few times
  e.k.a <- exp(-k.a * x)
  e.k.e <- exp(-k.e * x)
  e.diff <- e.k.e - e.k.a
  k.diff <- k.a - k.e
  # the gradient
  grad <- rep(0, 3)
  grad[1] <- -D * k.a / (V^2 * k.diff) * e.diff
  grad[2] <- D * k.a / (V * k.diff) * e.diff - D * k.a / (V * k.diff^2) * e.diff + D * k.a / (V * k.diff^2) * e.k.a * x
  grad[3] <- D * k.a / (V * k.diff^2) * e.diff - D * k.a / (V * k.diff) * e.k.e * x
  grad
}

# hessian
h <- function(beta, x) {
  D <- 4.53
  V <- beta[1]
  k.a <- beta[2]
  k.e <- beta[3]
  # come up a few times
  e.k.a <- exp(-k.a * x)
  e.k.e <- exp(-k.e * x)
  e.diff <- e.k.e - e.k.a
  k.diff <- k.a - k.e
  # the hessian
  hess <- matrix(0, nrow=3, ncol=3)
  hess[1,1] <- 2 * D * k.a / (V^3 * k.diff) * e.diff
  hess[1,2] <- -D / (V^2 * k.diff) * e.diff + D * k.a / (V^2 * k.diff^2) * e.diff - D * k.a / (V^2 * k.diff^2) * e.k.a * x
  hess[2,1] <- hess[1,2]
  hess[1,3] <- -D * k.a / (V^2 * k.diff^2) * e.diff + D * k.a / (V^2 * k.diff) * e.k.e * x
  hess[3,1] <- hess[1,3]
  hess[2,2] <- -D / (V^2 * k.diff^2) * e.diff + D / (V * k.diff) * e.k.a * x - D / (V * k.diff^2) * e.diff + 2 * D * k.a / (V * k.diff^3) * e.diff - D * k.a / (V * k.diff) * e.k.a * x + D / (V * k.diff) * e.k.a * x - D * k.a / (V * k.diff^2) * e.k.a * x - D * k.a / (V * k.diff^2) * e.k.a * x^2
  hess[3,3] <- 2 * D * k.a / (V * k.diff^3) * e.diff - D * k.a / (V * k.diff) * e.k.a * x - D * k.a / (V * k.diff^2) * e.k.a * x + D * k.a / (V * k.diff) * e.k.a * x^2
  hess[2,3] <- D / (V * k.diff^2) * e.diff - D / (V * k.diff) * e.k.e * x - 2 * D * k.a / (V * k.diff^3) * e.diff + D * k.a / (V * k.diff^2) * e.k.e * x + D * k.a / (V * k.diff^2) * e.k.a * x
  hess[3,2] <- hess[2,3]
  hess
}

# loss function
f.loss <- function(beta, x, y) {
  mean((y - f(beta, x))^2)
}

# gradient of loss function
g.loss <- function(beta, x, y) {
  grad <- rep(0, 3)
  n <- length(x)
  for (i in 1:n) {
    grad <- grad + 2 * (f(beta, x[i]) - y[i]) * g(beta, x[i])
  }
  grad / n
}

# hessian of loss function
h.loss <- function(beta, x, y) {
  hess <- matrix(0, nrow=3, ncol=3)
  n <- length(x)
  for (i in 1:n) {
    hess <- hess + 2 * tcrossprod(g(beta, x[i])) + 2 * (f(beta, x[i]) - y[i]) * h(beta, x[i])
  }
  hess / n
}

A.hat <- function(beta, x, y) {
  h.loss(beta, x, y)
}

B.hat <- function(beta, x, y) {
  n <- length(x)
  n * tcrossprod(g.loss(beta, x, y))
}

f.var <- function(beta, x, y, xi) {
  bcov <- beta.vcov(beta, x, y)
  grad <- g(beta, xi)
  ((bcov %*% grad)[,1] %*% grad)[1,1]
}

beta.vcov <- function(beta, x, y) {
  n <- length(x)
  ahat <- A.hat(beta, x, y)
  bhat <- B.hat(beta, x, y)
  ahat.inv <- solve(ahat)
  ahat.inv %*% bhat %*% ahat.inv / n
}

f.ci <- function(beta, x, y, xis) {
  n <- length(xis)
  vals <- f(beta, xis)
  fvar <- sapply(xis, function(xi) {f.var(beta, x, y, xi)})
  lower.mean <- vals + qnorm(0.025) * sqrt(fvar)
  upper.mean <- vals + qnorm(0.975) * sqrt(fvar)
  sigma.sq <- sum((y - f(beta, x))^2) / length(x)
  lower.error <- vals + qnorm(0.025) * sqrt(sigma.sq)
  upper.error <- vals + qnorm(0.975) * sqrt(sigma.sq)
  list(value=vals, lower.mean=lower.mean, upper.mean=upper.mean,
       lower.error=lower.error, upper.error=upper.error)
}


# attempt to optimize with Newton's method (something wrong with Hessian?)
beta <- c(exp(-0.8), exp(0.8), exp(-2.5))
for (it in 1:100) {
  grad <- g.loss(beta, x, y)
  hess <- h.loss(beta, x, y)
  beta <- beta - solve(hess, grad)
  #cat(beta, "\n")
  #cat(f.loss(beta, x, y), "\n")
}

# optimizing with just gradient and BFGS works well...
beta <- c(exp(-0.8), exp(0.8), exp(-2.5))
b <- optim(beta, f.loss, gr=g.loss, x=x, y=y, method="BFGS")
beta <- b$par

# 95% confidence interval on betas
b.cov <- beta.vcov(beta, x, y)
sds <- sqrt(diag(b.cov))
lowers <- round(beta + qnorm(0.025) * sds, 4)
uppers <- round(beta + qnorm(0.975) * sds, 4)
PrintCI <- function() {
  cat("95% Confidence Interval:\n")
  cat(paste("    V  : ", round(beta[1], 4), " (", lowers[1], ", ", uppers[1], ")\n", sep=""))
  cat(paste("    k_a: ", round(beta[2], 4), " (", lowers[2], ", ", uppers[2], ")\n", sep=""))
  cat(paste("    k_e: ", round(beta[3], 4), " (", lowers[3], ", ", uppers[3], ")\n", sep=""))
}
PrintCI()


# 95% confidence bands on function
xvals <- seq(0, 25, 0.01)
res <- f.ci(beta, x, y, xvals)
df <- data.frame(x=xvals,
                 y=res$value,
                 lower.error=res$lower.error,
                 upper.error=res$upper.error,
                 lower.mean=res$lower.mean,
                 upper.mean=res$upper.mean)

library(ggplot2)
quartz()
ggplot(data.frame(x=x,y=y), aes(x=x,y=y)) + geom_point() +
  geom_line(data=df, aes(x=x, y=y)) +
  geom_ribbon(data=df, aes(ymin=lower.mean, ymax=upper.mean), alpha=0.6) +
  geom_ribbon(data=df, aes(ymin=lower.error, ymax=upper.error), alpha=0.3)