Binles VPC Examples

Population pharmacokinetic model of phenobarbital.
Hypothetical one‐compartment population pharmacokinetic model.

Overview

VPC1 - Population pharmacokinetic model of phenobarbital.

Optimized smoothing parameters
User specifying llambda and span parameters.
Optimized smoothing parameters with larger lambda optimization interval and 99% confidence interval for sim fits.

VPC2 - Hypothetical one‐compartment population pharmacokinetic model.

Stratified by GROUP (0 = No concomitant, 1 = Concomitant). Loess prediction corrected.
Stratified by GROUP (0 = No concomitant, 1 = Concomitant). User specified lambda by group.
Stratfied by GROUP + SEX (Sex added for example).

Population pharmacokinetic model of phenobarbital.

Setup

setwd("~/myPackages/Regression Approach to VPC")

#Model results - PsN
mod.res <- read.table(
    file="sdtab1",
    header=T, skip=1, stringsAsFactors=F
    )

mod.res <- subset(mod.res, subset=c(RES!=0))


#Simulated (stand alone simulation based on model above)
sim.pheno <- read.table(
    file="phenofull_sim.fit",
    header=F, stringsAsFactors=F
    )
#ID MDV TIME DV PRED REP
colnames(sim.pheno) <- c("ID", "MDV", "TIME", "DV", "PRED", "REP")

sim.pheno <- subset(sim.pheno, subset=c(MDV==0))

Optimized smoothing parameters

  vpc1 <- observed(mod.res, y = DV, x = TIME) %>%
  simulated(sim.pheno, y = DV) %>%
  predcorrect(pred = PRED) %>%
  binlessaugment(qpred = c(0.1, 0.5, 0.9), loess.ypc = TRUE) %>%
  binlessfit() %>%
  binlessvpcstats()


ggplot(vpc1$stats, aes(x=x)) +
   # facet_grid(~ GROUP) +
    geom_ribbon(aes(ymin=lo, ymax=hi, fill=qname, col=qname, group=qname), alpha=0.1, col=NA) +
    geom_line(aes(y=med, col=qname, group=qname)) +
    geom_line(aes(y=y, linetype=qname), size=1) +
    geom_point(aes(y = l.ypc, x = x), data = vpc1$obs)  +
    scale_colour_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_fill_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_linetype_manual(
        name="Observed Percentiles\n(black lines)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("dotted", "solid", "dashed"),
        labels=c("10%", "50%", "90%")) +
    guides(
        fill=guide_legend(order=2),
        colour=guide_legend(order=2),
        linetype=guide_legend(order=1)) +
    theme(
        legend.position="top",
        legend.key.width=grid::unit(2, "cm")) +
    labs(x="Time (h)", y="Concentration (ng/mL)")

Optimized llambda from vpc1

print(vpc1$llam.qpred)

##     q0.1     q0.5     q0.9 
## 4.308714 2.518344 3.565408

Optimized span from vpc1

print(vpc1$span)

## [1] 0.6072733

User specifying llambda and span parameters.

vpc1.1 <- observed(mod.res, y = DV, x = TIME) %>%
  simulated(sim.pheno, y = DV) %>%
  predcorrect(pred = PRED) %>%
  binlessaugment(qpred = c(0.1, 0.5, 0.9), loess.ypc = TRUE) %>%
  binlessfit(llam.quant = c(3.5, 1.5, 2.5), span = .75) %>%
  binlessvpcstats()


ggplot(vpc1.1$stats, aes(x=x)) +
   # facet_grid(~ GROUP) +
    geom_ribbon(aes(ymin=lo, ymax=hi, fill=qname, col=qname, group=qname), alpha=0.1, col=NA) +
    geom_line(aes(y=med, col=qname, group=qname)) +
    geom_line(aes(y=y, linetype=qname), size=1) +
    geom_point(aes(y = l.ypc, x = x), data = vpc1.1$obs)  +
    scale_colour_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_fill_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_linetype_manual(
        name="Observed Percentiles\n(black lines)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("dotted", "solid", "dashed"),
        labels=c("10%", "50%", "90%")) +
    guides(
        fill=guide_legend(order=2),
        colour=guide_legend(order=2),
        linetype=guide_legend(order=1)) +
    theme(
        legend.position="top",
        legend.key.width=grid::unit(2, "cm")) +
    labs(x="Time (h)", y="Concentration (ng/mL)")

Optimized smoothing parameters with larger lambda optimization interval and 99% confidence interval for sim fits.

  vpc1.2 <- observed(mod.res, y = DV, x = TIME) %>%
  simulated(sim.pheno, y = DV) %>%
  predcorrect(pred = PRED) %>%
  binlessaugment(qpred = c(0.1, 0.5, 0.9), interval = c(-2, 10), loess.ypc = TRUE) %>%
  binlessfit(conf.level = .99) %>%
  binlessvpcstats()


ggplot(vpc1.2$stats, aes(x=x)) +
    geom_ribbon(aes(ymin=lo, ymax=hi, fill=qname, col=qname, group=qname), alpha=0.1, col=NA) +
    geom_line(aes(y=med, col=qname, group=qname)) +
    geom_line(aes(y=y, linetype=qname), size=1) +
    geom_point(aes(y = l.ypc, x = x), data = vpc1.2$obs)  +
    scale_colour_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_fill_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_linetype_manual(
        name="Observed Percentiles\n(black lines)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("dotted", "solid", "dashed"),
        labels=c("10%", "50%", "90%")) +
    guides(
        fill=guide_legend(order=2),
        colour=guide_legend(order=2),
        linetype=guide_legend(order=1)) +
    theme(
        legend.position="top",
        legend.key.width=grid::unit(2, "cm")) +
    labs(x="Time (h)", y="Concentration (ng/mL)")

Hypothetical one‐compartment population pharmacokinetic model.

Setup

setwd("~/myPackages/Regression Approach to VPC")
est.PK <- read.table(
  file="hyp_sim_PK.fit",
  skip=1, header=T, stringsAsFactors=F
)
est.PK <- subset(est.PK, subset=c(MDV==0))

sims.PK <- read.table(
  file="hyp_sim_pk_sim.fit",
  header=F, stringsAsFactors=F
)
#ID TIME DV GROUP AMT MDV EVID PRED REP
colnames(sims.PK) <- c(
  "ID", "TIME", "DV", "GROUP", "AMT", "MDV", "EVID", "PRED", "REP"
)

sims.PK <- subset(sims.PK, subset=c(MDV==0))

Stratified by GROUP, loess prediction corrected.

vpc2 <- observed(est.PK, y = DV, x = TIME) %>%
  simulated(sims.PK, y = DV) %>%
  stratify(~ GROUP) %>%
  predcorrect(pred = PRED) %>%
  binlessaugment(qpred = c(0.1, 0.5, 0.9), loess.ypc = TRUE) %>%
  binlessfit() %>%
  binlessvpcstats()



ggplot(vpc2$stats, aes(x=x)) +
   facet_grid(vpc2$strat.formula) +
    geom_ribbon(aes(ymin=lo, ymax=hi, fill=qname, col=qname, group=qname), alpha=0.1, col=NA) +
    geom_line(aes(y=med, col=qname, group=qname)) +
    geom_line(aes(y=y, linetype=qname), size=1) +
    geom_point(aes(y = l.ypc, x = x), data = vpc2$obs)  +
    scale_colour_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_fill_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_linetype_manual(
        name="Observed Percentiles\n(black lines)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("dotted", "solid", "dashed"),
        labels=c("10%", "50%", "90%")) +
    guides(
        fill=guide_legend(order=2),
        colour=guide_legend(order=2),
        linetype=guide_legend(order=1)) +
    theme(
        legend.position="top",
        legend.key.width=grid::unit(2, "cm")) +
    labs(x="Time (h)", y="Concentration (ng/mL)")

Stratified by GROUP. User specified lambda. No loess prediction corrected fits.

newlambda <- data.table(
  group0 = c(2, 3, 4),
  group1 = c(2.5, 4, 5)
)

vpc2.1 <- observed(est.PK, y = DV, x = TIME) %>%
  simulated(sims.PK, y = DV) %>%
  stratify(~ GROUP) %>%
  binlessaugment(qpred = c(0.1, 0.5, 0.9)) %>%
  binlessfit(llam.quant = newlambda) %>%
  binlessvpcstats()



ggplot(vpc2.1$stats, aes(x=x)) +
   facet_grid(vpc2.1$strat.formula) +
    geom_ribbon(aes(ymin=lo, ymax=hi, fill=qname, col=qname, group=qname), alpha=0.1, col=NA) +
    geom_line(aes(y=med, col=qname, group=qname)) +
    geom_line(aes(y=y, linetype=qname), size=1) +
    geom_point(aes(y = y, x = x), data = vpc2.1$obs)  +
    scale_colour_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_fill_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("red", "blue", "red"),
        labels=c("10%", "50%", "90%")) +
    scale_linetype_manual(
        name="Observed Percentiles\n(black lines)",
        breaks=c("q0.1", "q0.5", "q0.9"),
        values=c("dotted", "solid", "dashed"),
        labels=c("10%", "50%", "90%")) +
    guides(
        fill=guide_legend(order=2),
        colour=guide_legend(order=2),
        linetype=guide_legend(order=1)) +
    theme(
        legend.position="top",
        legend.key.width=grid::unit(2, "cm")) +
    labs(x="Time (h)", y="Concentration (ng/mL)")

Stratified by GROUP and SEX - Optimized smoothing parameters. Default quantiles (.05, .5, .95). No loess prediction corrected.

sex <- sample(c("F","M"), length(unique(est.PK$ID)), replace=T)
newobs <- data.table(ID = unique(est.PK$ID), SEX = sex)
est.PK <- merge(est.PK, newobs, by="ID")


vpc2.2 <- observed(est.PK, y = DV, x = TIME) %>%
  simulated(sims.PK, y = DV) %>%
  stratify(~ GROUP + SEX) %>%
  binlessaugment() %>%
  binlessfit() %>%
  binlessvpcstats()



ggplot(vpc2.2$stats, aes(x=x)) +
   facet_grid(vpc2.2$strat.formula) +
    geom_ribbon(aes(ymin=lo, ymax=hi, fill=qname, col=qname, group=qname), alpha=0.1, col=NA) +
    geom_line(aes(y=med, col=qname, group=qname)) +
    geom_line(aes(y=y, linetype=qname), size=1) +
    geom_point(aes(y = y, x = x), data = vpc2.2$obs)  +
    scale_colour_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.05", "q0.5", "q0.95"),
        values=c("red", "blue", "red"),
        labels=c("5%", "50%", "95%")) +
    scale_fill_manual(
        name="Simulated Percentiles\nMedian (lines) 95% CI (areas)",
        breaks=c("q0.05", "q0.5", "q0.95"),
        values=c("red", "blue", "red"),
        labels=c("5%", "50%", "95%")) +
    scale_linetype_manual(
        name="Observed Percentiles\n(black lines)",
        breaks=c("q0.05", "q0.5", "q0.95"),
        values=c("dotted", "solid", "dashed"),
        labels=c("10%", "50%", "90%")) +
    guides(
        fill=guide_legend(order=2),
        colour=guide_legend(order=2),
        linetype=guide_legend(order=1)) +
    theme(
        legend.position="top",
        legend.key.width=grid::unit(2, "cm")) +
    labs(x="Time (h)", y="Concentration (ng/mL)")