In smouksassi/coveffectsplot: Produce Forest Plots to Visualize Covariate Effects
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")
}
Here we illustrate how the ordinary differential equations (ODEs) model and the approach of varying one covariate at a time, can be expanded.
We link the same two-compartment PK model (from the PK Example vignette) to an indirect response pharmacodynamic (PD) model where the drug concentrations inhibit the rate constant of input (Kin).
The covariates model included several covariates effects on Clearance, Volume and Kin. The baseline PD value is controlled by the ratio of Kin/Kout.
In this vignette we do not go into a lot of details, as we assume that the user has read and run the code of the Introduction to coveffectsplot and PK Example vignettes and that the reader is familiar with PK/PD concepts.
At the end we show how we can add a table under a multiple parameters forest plot.
Specifying a PK/PD Model using mrgsolve
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)
Simulate Reference Subjects with BSV
We simulate at reference covariate values with between subject variability (BSV) and then we show a plot of the PK and PD profiles of five random subjects.
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())
Compute PD Parameters and Summarize BSV
Here we compute the PD baseline (where we start), nadir response (minimum response achieved) and the delta (difference) between the baseline and nadir. We then summarize and report the BSV around these parameters as ranges of 50 and 90% of patients. We then show a plot of the first 10 replicates as an example of the simulated PD profiles. Since the code is similar to the PK Example vignette it is not shown.
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
Construct ans Simulate at Combinations of Covariate of Interest
Similarly to the PK Example vignette we generate covariate combinations of interest and we simulate with uncertainty using an invented varcov matrix.
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
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"))
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
Compute PD Parameters and Distributions Plots
Similar to the above we compute the PD parameters, standardize by the median and provide a plot. Since the code is similar to the PK Example vignette it is not shown.
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"]
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
# 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()
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
# pdggridges+theme(legend.position = "none")
# ggsave("Figure_S_PD_3.png", device="png",type="cairo-png",
# width= 7, height = 5,dpi=72)
Summarize, add the BSV Ranges and Putting it all Together Using forest_plot
Here we show how a multiple parameters, multiple covariates and table can be done.
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), "]"))
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 = ))
smouksassi/coveffectsplot documentation built on Jan. 19, 2024, 8:24 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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") }
Here we illustrate how the ordinary differential equations (ODEs) model and the approach of varying one covariate at a time, can be expanded.
We link the same two-compartment PK model (from the PK Example vignette) to an indirect response pharmacodynamic (PD) model where the drug concentrations inhibit the rate constant of input (Kin).
The covariates model included several covariates effects on Clearance, Volume and Kin. The baseline PD value is controlled by the ratio of Kin/Kout.
In this vignette we do not go into a lot of details, as we assume that the user has read and run the code of the Introduction to coveffectsplot and PK Example vignettes and that the reader is familiar with PK/PD concepts.
At the end we show how we can add a table under a multiple parameters forest plot.
Specifying a PK/PD Model using mrgsolve
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)
Simulate Reference Subjects with BSV
We simulate at reference covariate values with between subject variability (BSV) and then we show a plot of the PK and PD profiles of five random subjects.
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())
Compute PD Parameters and Summarize BSV
Here we compute the PD baseline (where we start), nadir response (minimum response achieved) and the delta (difference) between the baseline and nadir. We then summarize and report the BSV around these parameters as ranges of 50 and 90% of patients. We then show a plot of the first 10 replicates as an example of the simulated PD profiles. Since the code is similar to the PK Example vignette it is not shown.
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
Construct ans Simulate at Combinations of Covariate of Interest
Similarly to the PK Example vignette we generate covariate combinations of interest and we simulate with uncertainty using an invented varcov matrix.
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
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"))
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
Compute PD Parameters and Distributions Plots
Similar to the above we compute the PD parameters, standardize by the median and provide a plot. Since the code is similar to the PK Example vignette it is not shown.
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"]
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
# 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()
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
# pdggridges+theme(legend.position = "none") # ggsave("Figure_S_PD_3.png", device="png",type="cairo-png", # width= 7, height = 5,dpi=72)
Summarize, add the BSV Ranges and Putting it all Together Using forest_plot
Here we show how a multiple parameters, multiple covariates and table can be done.
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), "]"))
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 = ))
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.