inst/doc/PKPD_Example.R

In coveffectsplot: Produce Forest Plots to Visualize Covariate Effects

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  message =FALSE,
  warning =FALSE,
  fig.width = 7,
  comment = "#>"
)
if (capabilities(("cairo"))) {
  knitr::opts_chunk$set(dev.args = list(png = list(type = "cairo")))
}
options(rmarkdown.html_vignette.check_title = FALSE)
library(coveffectsplot)
library(mrgsolve)
library(ggplot2)
library(ggstance)
library(ggridges)
library(tidyr)
library(dplyr)
library(table1)
library(patchwork)
library(ggh4x)
library(data.table)
theme_set(theme_bw())
#utility function to simulate varying one covariate at a time keeping the rest at the reference
expand.modelframe <- function(..., rv, covcol="covname") {
  args <- list(...)
  df <- lapply(args, function(x) x[[1]])
  df[names(rv)] <- rv
  res <- lapply(seq_along(rv), function(i) {
    df[[covcol]] <- names(rv)[i]
    df[[names(rv)[i]]] <- args[[names(rv)[i]]]
    as.data.frame(df)
  })
  do.call(rbind, res)
}
cor2cov <- function (cor, sd) 
{
    if (missing(sd)) {
        sd <- diag(cor)
    }
    diag(cor) <- 1
    n <- nrow(cor)
    diag(sd, n) %*% cor %*% diag(sd, n)
}
nbsvsubjects <- 100
nsim <- 100 # uncertainty replicates for vignette you might want a higher number
round_pad <- function(x, digits = 2, round5up = TRUE) {
  eps <- if (round5up) x * (10^(-(digits + 3))) else 0
  formatC(round(x + eps, digits), digits = digits, format = "f", flag = "0")
}

## ----pkpdmodel, collapse=TRUE-------------------------------------------------
codepkpdmodelcov <- '
$PARAM @annotated
KA     : 0.5   : Absorption rate constant Ka (1/h)
CL     : 4     : Clearance CL (L/h)
V      : 10    : Central volume Vc (L)
Vp     : 50    : Peripheral volume Vp (L)
Qp     : 10    : Intercompartmental clearance Q (L/h)
CLALB  : -0.8  : Ablumin on CL (ref. 45 g/L)
CLSEX  : 0.2   : Sex on CL (ref. Female)
CLWT   : 1     : Weight on CL (ref. 85 kg)
VSEX   : 0.07  : Sex on Vc (ref. Female)
VWT    : 1     : Weight on Vc (ref. 85 kg)
KIN    : 3     : Zero-order Rate constant of biomarker production (amount/h)
KOUT   : 0.06  : First-order Rate constant of biomarker loss (1/h)
IC50   : 3     : Drug concentration producing 50% of maximum inhibition
IMAX   : 0.999 : Maximum Inhibition Response
gamma  : 0.55  : Sigmoidicity factor of the sigmoid Emax equation
KINWT  : 0.4   : Weight on KIN (ref. 85 kg)
KINAGE : -0.08 : Age on KIN (ref. 40 years)
KINHLTY: 1.5   : Weight on CL (ref. 85 kg)

$PARAM @annotated // reference values for covariate
WT     :  85    : Weight (kg)
SEX    :  0     : Sex (0=Female, 1=Male)
ALB    :  45    : Albumin (g/L)
AGE    :  40    : Age (years)
HEALTHY:  0     : Health Status (0=Diseased, 1=Healthy)

$CMT GUT CENT PER RESP
$GLOBAL
#define CP   (CENT/Vi)
#define CPER (PER/Vpi)
#define INH  (IMAX*pow(CP,gamma)/(pow(IC50,gamma)+pow(CP,gamma)))
#define PDRESP RESP

$MAIN
double KAi = KA;
double Vpi = Vp *pow((WT/70.0),    1);
double Qpi = Qp *pow((WT/70.0), 0.75);
double CLi = CL *
    pow((ALB/45.0), CLALB)*
    (SEX == 1.0 ? (1.0+CLSEX) : 1.0)*
    pow((WT/85.0), CLWT)*exp(ETA(1)); 
double Vi = V *
    (SEX == 1.0 ? (1.0+VSEX) : 1.0)*
    pow((WT/85.0), VWT)*exp(ETA(2));  
double KINi = KIN *
  pow((AGE/40), KINAGE)*
  (HEALTHY == 1.0 ? KINHLTY : 1.0)*
  pow((WT/85.0), KINWT)*exp(ETA(3));
double RESP_0 = KINi/KOUT;

$OMEGA
0.09 
0.01 0.09
$OMEGA
0.25

$ODE
dxdt_GUT    = -KAi *GUT;
dxdt_CENT   =  KAi *GUT  - (CLi+Qpi)*CP  + Qpi*CPER;
dxdt_PER    =                   Qpi*CP   - Qpi*CPER;
dxdt_RESP   =  KINi*(1-INH) - KOUT*RESP;

$CAPTURE CP PDRESP KAi CLi Vi Vpi Qpi WT SEX ALB AGE HEALTHY
'
modpkpdsim <- mcode("codepkpdmodelcov", codepkpdmodelcov)
partab <- setDT(modpkpdsim@annot$data)[block=="PARAM", .(name, descr, unit)]
partab <- merge(partab, melt(setDT(modpkpdsim@param@data), meas=patterns("*"), var="name"))
knitr::kable(partab)


## ----pkpdsimulation, fig.width=7,fig.height=4, message=FALSE------------------
idata <- data.table(ID=1:nbsvsubjects, WT=85, SEX=0, ALB=45, AGE=40, HEALTHY = 0)
ev1 <- ev(time = 0, amt = 100, cmt = 1, ii = 24, addl = 20)
data.dose <- ev(ev1)
data.dose <- setDT(as.data.frame(data.dose))
data.all <- data.table(idata, data.dose)

set.seed(678549)
outputpkpdsim <- modpkpdsim %>%
  data_set(data.all) %>%
  mrgsim(end = 28*24, delta = 0.25) %>%
  as.data.frame %>%
  as.data.table

outputpkpdsim$HEALTHY <- as.factor(outputpkpdsim$HEALTHY)

yvar_names <- c(
  'CP'="Plasma Concentrations",
  'RESP'="PD Values"
)
set.seed(678549)
outputpkpdsimlong <- outputpkpdsim[outputpkpdsim$ID %in%
sample(unique(outputpkpdsim$ID), 5), ] %>% 
  gather(key,value,CP,RESP)

ggplot(data =outputpkpdsimlong ,
       aes(time, value, group = ID)) +
  geom_line(alpha = 0.8, size = 0.3) +
  facet_grid(key ~ID,scales="free_y",switch="y",
             labeller = labeller(key=yvar_names)) +
  labs(y = "", color = "Sex", x = "Time (h)")+
  theme(strip.placement = "outside",
        axis.title.y=element_blank())

## ----computenca , fig.width=7,  message=FALSE---------------------------------
derive.exposure <- function(time, PDRESP) {
  x <- c(
    nadir = min(PDRESP, na.rm = TRUE),
    baselinepd = PDRESP[1L],
    deltapd = PDRESP[1L]-min(PDRESP, na.rm = TRUE)
  )
  data.table(paramname=names(x), paramvalue=x)
}
refbsv <- outputpkpdsim[, derive.exposure(time, PDRESP),
                        by=.(ID, WT, SEX, ALB, AGE, HEALTHY)]

refbsv[, stdparamvalue := paramvalue/median(paramvalue), by=paramname]

bsvranges <- refbsv[,list(
    P05 = quantile(stdparamvalue, 0.05),
    P25 = quantile(stdparamvalue, 0.25),
    P50 = quantile(stdparamvalue, 0.5),
    P75 = quantile(stdparamvalue, 0.75),
    P95 = quantile(stdparamvalue, 0.95)), by = paramname]
bsvranges

## ----covcomb , fig.width=7----------------------------------------------------
reference.values <- data.frame(WT = 85, ALB = 45, AGE = 40, SEX = 0, HEALTHY = 0)   
covcomb <- expand.modelframe(
  WT  = c(56,128), 
  AGE = c(20,60),
  ALB = c(40,50),
  SEX = c(1),#Refernce is for SEX =0
  HEALTHY = c(1),#Refernce is for HEALTHY =0
  rv = reference.values)

# Add the reference
covcomb <- rbind(covcomb, data.table(reference.values, covname="REF"))
covcomb$ID <- 1:nrow(covcomb)

covcomb

## ---- fig.width=7 ,message=FALSE, include=FALSE-------------------------------
idata <- data.table::copy(covcomb)
idata$covname <-  NULL
ev1 <- ev(time=0, amt=100, cmt=1, ii = 24, addl = 20)
data.dose <- as.data.frame(ev1)
data.all <- data.table(idata, data.dose)

outcovcomb<- modpkpdsim %>%
  data_set(data.all) %>%
  zero_re() %>% 
  mrgsim(start=0,end=24*28,delta=0.25)%>%
  as.data.frame %>%
  as.data.table

outcovcomb$SEX <- as.factor(outcovcomb$SEX )
outcovcomb$SEX <- factor(outcovcomb$SEX, labels=c("Female", "Male"))

outcovcomb$HEALTHY <- as.factor(outcovcomb$HEALTHY )


theta <- unclass(as.list(param(modpkpdsim)))
theta[c("WT", "SEX", "ALB","AGE","HEALTHY")] <- NULL
theta <- unlist(theta)
as.data.frame(t(theta))

varcov <- cor2cov(
  matrix(0.2, nrow=length(theta), ncol=length(theta)),
  sd=theta*0.25)
rownames(varcov) <- colnames(varcov) <- names(theta)
as.data.frame(varcov)

set.seed(678549)
# mvtnorm::rmvnorm is another option that can be explored
sim_parameters <- MASS::mvrnorm(nsim, theta, varcov, empirical=T) %>% as.data.table
head(sim_parameters)

idata <- data.table::copy(covcomb)
idata$covname <-  NULL
ev1 <- ev(time=0, amt=100, cmt=1, ii = 24, addl = 20)
data.dose <- as.data.frame(ev1)

iter_sims <- NULL
for(i in 1:nsim) {
  data.all  <- data.table(idata, data.dose, sim_parameters[i])
  out <- modpkpdsim %>%
    data_set(data.all) %>%
    zero_re() %>% 
    mrgsim(start=0, end=28*24, delta=0.25) %>%
    as.data.frame %>%
    as.data.table
  out[, rep := i]
  iter_sims <- rbind(iter_sims, out)
}
iter_sims$SEX <- as.factor(iter_sims$SEX )
iter_sims$SEX <- factor(iter_sims$SEX, labels=c("Female", "Male"))


## ---- fig.width=7, fig.height=6, message=FALSE, warning=FALSE-----------------
albumin.labs <- c("albumin: 40 ng/mL","albumin: 45 ng/mL","albumin: 50 ng/mL")
names(albumin.labs) <- c("40","45","50")
wt.labs <- c("weight: 85 kg","weight: 56 kg","weight: 128 kg")
names(wt.labs) <- c("85","56","128")

age.labs <- c("age: 20 years","age: 40 years","age: 60 years")
names(age.labs) <- c("20","40","60")

pdprofiles <- ggplot(iter_sims[iter_sims$rep<=10,], aes(time/24,PDRESP,col=factor(WT),linetype=factor(HEALTHY) ) )+
  geom_line(aes(group=interaction(ID,rep)),alpha=0.3,size=0.3)+
  geom_line(data=outcovcomb,aes(group=interaction(ID)),color="black")+
  facet_nested(ALB+SEX~ AGE+WT,  labeller = 
                 labeller( WT = wt.labs,
                           ALB = albumin.labs,
                           AGE = age.labs))+
  labs(linetype="Black Lines\nNo Uncertainty\nHealthy Status",
       colour="Colored Lines\nUncertainty\nReplicates\n(1 to 10)\nWeight (kg)",
       caption ="Simulation\nwith Uncertainty without BSV" ,
       x="Days", y = "PD Values")+
  guides(colour = guide_legend(override.aes = list(alpha = 1)))

pdprofiles

## ---- fig.width=7,fig.height=6, include=FALSE, message=FALSE------------------

out.df.univariatecov.nca <- iter_sims[, derive.exposure(time, PDRESP),
                                      by=.(rep, ID, WT, SEX, ALB, AGE, HEALTHY)]
out.df.univariatecov.nca

refvalues <- out.df.univariatecov.nca[
  ALB==45 & WT==85 & SEX=="Female"& AGE==40 & HEALTHY==0,
  .(medparam = median(paramvalue)), by=paramname]
refvalues

covcomb <- as.data.table(covcomb)

covcomb[covname=="WT",  covvalue := paste(WT,"kg")]
covcomb[covname=="ALB", covvalue := paste(ALB,"g/L")]
covcomb[covname=="AGE", covvalue := paste(AGE,"years")]
covcomb[covname=="SEX", covvalue := "Male"]
covcomb[covname=="HEALTHY", covvalue := "Diseased"]
covcomb[covname=="REF", covvalue := "85 kg-Female-45 g/L-40 years-healthy"]
covcomb
covcomb[covname=="REF", covvalue := "85 kg-Female\n45 g/L-40 years\nhealthy"]
covcomb <- as.data.table(covcomb)
out.df.univariatecov.nca <- merge(
  out.df.univariatecov.nca,
  covcomb[, .(ID, covname, covvalue)])

setkey(out.df.univariatecov.nca, paramname)
setkey(refvalues, paramname)
out.df.univariatecov.nca <- merge(out.df.univariatecov.nca,refvalues)

out.df.univariatecov.nca[, paramvaluestd := paramvalue/medparam]

out.df.univariatecov.nca$covvalue <-factor(as.factor(out.df.univariatecov.nca$covvalue ),
                                               levels =  c("56 kg",
                                                           "85 kg",
                                                           "128 kg",
                                                           "Male",
                                                           "40 g/L",
                                                           "50 g/L",
                                                           "20 years",
                                                           "60 years",
                                                           "Diseased",
"85 kg-Female\n45 g/L-40 years\nhealthy")
)
out.df.univariatecov.nca$covname2 <- as.factor(out.df.univariatecov.nca$covname2)
out.df.univariatecov.nca$covname2 <- factor(out.df.univariatecov.nca$covname2,
                                            levels= c( "Weight", "Sex",
                                                    "Albumin","Age", "Healthy",
                                                    "Reference")
                                            )
boxplotdat <- out.df.univariatecov.nca
boxplotdat[covname=="WT",  covname2 := "Weight"]
boxplotdat[covname=="ALB", covname2 := "Albumin"]
boxplotdat[covname=="SEX", covname2 := "Sex"]
boxplotdat[covname=="AGE", covname2 := "Age"]
boxplotdat[covname=="HEALTHY", covname2 := "Healthy"]
boxplotdat[covname=="REF", covname2 := "Reference"]


## ---- fig.width=7, fig.height=5, message=FALSE--------------------------------
boxplotpd <- ggplot(boxplotdat,
       aes(x=covvalue ,y=paramvalue))+
  facet_grid(paramname ~covname2,scales="free",switch="both",
             labeller = label_parsed)+
    geom_boxplot()+
  theme(axis.title = element_blank(),strip.placement = "outside")+
  labs(y="PD Parameter Values",x="Covariate Value")
boxplotpd

## ---- fig.width=7 ,message=FALSE, include=FALSE-------------------------------
#  pdprofiles<- pdprofiles+theme(axis.title.y = element_text(size=15))+
#    guides(colour=guide_legend(override.aes = list(alpha=1,size=0.5)),
#           linetype=guide_legend(override.aes = list(size=0.5)))
# pdprofiles
# ggsave("pd3.png", device="png",type="cairo-png",width= 7, height = 5,dpi=72)
#  boxplotpd
#  ggsave("pd4.png", device="png",type="cairo-png",width= 7, height = 4,dpi=2*72)
#   png("Figure_S_PD_2.png", type="cairo-png",width= 2*7*72, height =5*72)
#  egg::ggarrange(pdprofiles,boxplotpd,nrow=1)
#  dev.off()

## ---- fig.width=7,fig.height=5,message=FALSE----------------------------------
pdggridges<- ggplot(out.df.univariatecov.nca,
       aes(x=paramvaluestd,y=covvalue,fill=factor(..quantile..),height=..ndensity..))+
  facet_grid(covname2~paramname,scales="free_y",space="free")+
  annotate( "rect",
            xmin = 0.5,
            xmax = 2,
            ymin = -Inf,
            ymax = Inf,
            fill = "gray",alpha=0.4
  )+
  stat_density_ridges(
    geom = "density_ridges_gradient", calc_ecdf = TRUE,
    quantile_lines = TRUE, rel_min_height = 0.001,scale=0.9,
    quantiles = c(0.05,0.5, 0.95)) +
  scale_fill_manual(
    name = "Probability", values = c("white", "#0000FFA0","#0000FFA0", "white"),
    labels = c("(0, 0.05]", "(0.05, 0.5]","(0.5, 0.95]", "(0.95, 1]")
  )+
  geom_vline( aes(xintercept = 1),size = 1)+
  theme_bw()+
  labs(x="Effects Relative to Parameter Reference value",y="")+
  scale_x_continuous(breaks=c(0.25,0.5,0.8,1/0.8,1/0.5,1/0.25))+
  scale_x_log10()
pdggridges

## ---- fig.width=7 ,message=FALSE, include=FALSE-------------------------------
# pdggridges+theme(legend.position = "none")
#  ggsave("Figure_S_PD_3.png", device="png",type="cairo-png",
#         width= 7, height = 5,dpi=72)


## ---- fig.width=7, fig.height=7 ,message=FALSE--------------------------------
coveffectsdatacovrep <- out.df.univariatecov.nca %>% 
  dplyr::group_by(paramname,covname,covvalue) %>% 
  dplyr::summarize(
    mid= median(paramvaluestd),
    lower= quantile(paramvaluestd,0.05),
    upper = quantile(paramvaluestd,0.95))

coveffectsdatacovreplabel<-   coveffectsdatacovrep %>%
  mutate(
    label= covvalue,
    LABEL = paste0(format(round(mid,2), nsmall = 2),
                   " [", format(round(lower,2), nsmall = 2), "-",
                   format(round(upper,2), nsmall = 2), "]"))


## ---- fig.width=7, fig.height=7 ,message=FALSE--------------------------------
setkey(bsvranges, paramname)
coveffectsdatacovrepbsv <- coveffectsdatacovrep[coveffectsdatacovrep$covname=="REF",]
coveffectsdatacovrepbsv$covname <- "BSV"
coveffectsdatacovrepbsv$covvalue <- "50% of patients"
coveffectsdatacovrepbsv$label <-    "50% of patients"
coveffectsdatacovrepbsv$lower <- bsvranges$P25
coveffectsdatacovrepbsv$upper <- bsvranges$P75

coveffectsdatacovrepbsv2 <- coveffectsdatacovrep[coveffectsdatacovrep$covname=="REF",]
coveffectsdatacovrepbsv2$covname <- "BSV"
coveffectsdatacovrepbsv2$covvalue <- "90% of patients"
coveffectsdatacovrepbsv2$label <-    "90% of patients"
coveffectsdatacovrepbsv2$lower <- bsvranges$P05
coveffectsdatacovrepbsv2$upper <- bsvranges$P95
coveffectsdatacovrepbsv<- rbind(coveffectsdatacovrep,coveffectsdatacovrepbsv,coveffectsdatacovrepbsv2)
coveffectsdatacovrepbsv <- coveffectsdatacovrepbsv %>% 
  mutate(
    label= covvalue,
    LABEL = paste0(format(round(mid,2), nsmall = 2),
                   " [", format(round(lower,2), nsmall = 2), "-",
                   format(round(upper,2), nsmall = 2), "]"))
coveffectsdatacovrepbsv<- as.data.frame(coveffectsdatacovrepbsv)
coveffectsdatacovrepbsv$label <- as.factor(coveffectsdatacovrepbsv$covvalue )
coveffectsdatacovrepbsv$label <- reorder(coveffectsdatacovrepbsv$label,
                                         coveffectsdatacovrepbsv$lower)


coveffectsdatacovrepbsv$covname <-factor(as.factor(coveffectsdatacovrepbsv$covname ),levels =c("WT","SEX","ALB","AGE","HEALTHY", "REF", "BSV"),
labels=  c("Weight","Sex","Albumin","Age","Healthy", "Reference", "BSV")
)

interval_legend_text <- "Median (points)\n90% CI (horizontal lines)"
interval_bsv_text <- "BSV (points)\nPrediction Intervals (horizontal lines)"
ref_legend_text <- "Reference\n(vertical line)\nClinically relevant limits\n(gray area)"
area_legend_text <- "Reference\n(vertical line)\nClinically relevant limits\n(gray area)"
png("./Figure_S_PD_4.png",width =12 ,height = 9,units = "in",res=72)
coveffectsplot::forest_plot(coveffectsdatacovrepbsv,
                            ref_area = c(0.5, 1/0.5),
                            x_range = c(0.25,4),
                            strip_placement = "outside",
                            base_size = 16,
                            y_label_text_size = 12,
                            y_label_text_width = 50,
                            xlabel = "Fold Change Relative to Reference",
                            ref_legend_text = ref_legend_text,
                            area_legend_text =area_legend_text,
                            interval_legend_text = interval_legend_text,
                            interval_bsv_text = interval_bsv_text,
                            facet_formula = "covname~paramname",
                            facet_switch = "y",
                            facet_scales = "free_y",
                            facet_space = "fixed",
                            paramname_shape = FALSE,
                            table_position = "below",
                            table_text_size=4,
                            plot_table_ratio = 1,
                            table_facet_switch = "both",
                            show_table_facet_strip = "both",
                            show_table_yaxis_tick_label = TRUE,
                            logxscale = TRUE,
                            major_x_ticks          = c(0.5, 1,  1/0.5),
                            major_x_labels         = c("1/2", "1", "2"),
                            table_margin = c(0,5.5,0,0),
                            plot_margin =c(0,5.5,0,0),
                            reserve_table_xaxis_label_space = FALSE,
                            return_list = FALSE)
dev.off()
# consider returning a list and editing the y axis label line breaks height
# theme(axis.text.y = element_text(lineheight = ))