# PK Example

## Specifying a PK Model using mrgsolve

Here we illustrate the approach with a two-compartment PK model defined with an ODE and covariates on the PK parameters Clearance and Volume.

codepkmodelcov <- '
$PARAM CL = 4, V=10 , KA=0.5, Vp =50, Qp= 10, CLALB = -0.8, CLSEX = 0.2, CLWT = 1, VSEX = 0.07, VWT = 1, WT=85, SEX=0, AGE=45, ALB =45$CMT GUT CENT PER
$MAIN 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 KAi = KA; double Vpi = Vp *pow((WT/70.0), 1); double Qpi = Qp *pow((WT/70.0), 0.75); double Keli = CLi/Vi ; double Kpti = Qpi/Vi ; double Ktpi = Qpi/Vpi ;$OMEGA
0.3
0.01 0.3

$ODE dxdt_GUT = -KAi*GUT; dxdt_CENT = KAi*GUT-Keli*CENT-Kpti*CENT+Ktpi*PER; dxdt_PER = Kpti*CENT-Ktpi*PER;$TABLE
double CP   = CENT/ Vi;
double CPER =  PER/Vpi;
$CAPTURE CP CPER KAi CLi Vi Vpi Qpi ' modcovsim <- mcode("codepkmodelcov", codepkmodelcov) #> Building codepkmodelcov ... done. ### Simulate Reference Subjects with BSV We simulate with between subject variability (BSV) and provide a plot of concentrations on linear and log linear scales. idata <- tibble( ID = 1:1000, WT = 85, ALB = 45, SEX = 0 ) ev1 <- ev(time = 0, amt = 100, cmt = 1) data.dose <- ev(ev1) data.dose <- as.data.frame(data.dose) data.all <- merge(idata, data.dose) outputsim <- modcovsim %>% data_set(data.all) %>% carry.out(WT, AGE, SEX, ALB, CLi) %>% mrgsim(end = 24, delta = 0.25) outputsim <- as.data.frame(outputsim) outputsim <- outputsim %>% arrange(ID, time) outputsim$SEX <- as.factor(outputsim$SEX) set.seed(678549) p1 <- ggplot(data = outputsim[outputsim$ID %in% sample(unique(outputsim$ID), 500), ], aes(time, CP, group = ID)) + geom_line(alpha = 0.2, size = 0.1) + facet_grid(ALB ~ WT+SEX, labeller = label_both) + labs(y = "Plasma Concentrations", color = "Sex", x = "Time (h)") set.seed(678549) p2 <- ggplot(data = outputsim[outputsim$ID %in% sample(unique(outputsim$ID), 500), ], aes(time, CP, group = ID)) + geom_line(alpha = 0.2, size = 0.1) + facet_grid(ALB ~ WT+SEX, labeller = label_both) + scale_y_log10() + labs(y = expression(Log[10]~~Plasma~~Concentrations), color = "Sex", x = "Time (h)") p1+p2 #> Warning: Transformation introduced infinite values in continuous y-axis ### Compute NCA PK Parameters, Plot and Summarize BSV We then compute simple NCA parameters, provide a plot of the parameters as well as of the standardized ones. We also summarize and report the BSV as ranges of 50 and 90% of patients for each NCA parameter. NCATYPICAL <- outputsim %>% group_by(ID, ALB, WT, SEX) %>% summarise ( Cmax = max(CP, na.rm = TRUE), Clast = CP[n()], AUC = sum(diff(time) * na.omit(lead(CP) + CP)) / 2 ) %>% gather(paramname, paramvalue, Cmax, Clast, AUC) p3 <- ggplot(NCATYPICAL, aes(x=paramvalue,y=paramname,fill=factor(..quantile..),height=..ndensity..))+ facet_wrap(~paramname,scales="free",ncol=1)+ 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.25,0.5,0.75, 0.95)) + scale_fill_manual( name = "Probability", values = c("#FF0000A0", "#0000FFA0", "white","white", "#0000FFA0","#FF0000A0"), labels = c("(0, 0.05]", "(0.05, 0.25]", "(0.25, 0.5]","(0.5, 0.75]", "(0.75, 0.95]","(0.95, 1]") )+ theme_bw()+ labs(x="NCA parameters",y="")+ scale_x_log10()+ coord_cartesian(expand = FALSE) NCATYPICALREF <- NCATYPICAL%>% group_by (paramname) %>% mutate(medparam = median(paramvalue), paramvalue = paramvalue / medparam) p4 <- ggplot(NCATYPICALREF, aes(x=paramvalue,y=paramname,fill=factor(..quantile..),height=..ndensity..))+ facet_wrap(~paramname,scales="free_y",ncol=1)+ 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.25,0.5,0.75, 0.95)) + scale_fill_manual( name = "Probability", values = c("#FF0000A0","#0000FFA0", "white","white", "#0000FFA0","#FF0000A0"), labels = c("(0, 0.05]", "(0.05, 0.25]", "(0.25, 0.5]","(0.5, 0.75]", "(0.75, 0.95]","(0.95, 1]") )+ theme_bw()+ theme(legend.position = "none")+ labs(x="Standardized NCA parameters",y="")+ scale_x_log10()+ coord_cartesian(expand = FALSE) p3+p4  BSVRANGES<- NCATYPICALREF %>% summarize( P05 = quantile(paramvalue, 0.05), P25 = quantile(paramvalue, 0.25), P50 = quantile(paramvalue, 0.5), P75 = quantile(paramvalue, 0.75), P95 = quantile(paramvalue, 0.95) ) BSVRANGES paramname P05 P25 P50 P75 P95 AUC 0.5248993 0.7945283 1 1.215652 1.486990 Clast 0.2432740 0.6147814 1 1.477597 2.190552 Cmax 0.7851306 0.9229614 1 1.068844 1.150875 ## Construct Combinations of Covariate of Interest Based on our observed covariate data, we compute percentiles of interest that we will use to simulate data at. Common practice is to compute the 5,25,75,95 percentiles (the median being the reference). In some cases, we might want to explore the min, max or other extreme case scenarios. Care should be taken as this approach might generate unrealistic combination of covariates that can never appear in a real patient. expand.modelframe (written by Benjamin Rich) is a function defined in the setup section of the vignette and can be found in the source code. reference.values <- data.frame(WT = 85, ALB = 45, SEX = 0 ) covcomb <- expand.modelframe( WT = c(56,72,85,98,128), #P05,P25,P50,P75,P95 ALB = c(40,45,50),#P05,P50,P95 SEX = c(0,1),#Refernce is for SEX =0 rv = reference.values) covcomb$ID <- 1:nrow(covcomb)
covcomb <- covcomb[!duplicated(
paste(covcomb$WT,covcomb$AGE,covcomb$ALB,covcomb$SEX)),]
#remove duplicate cov combinations unless you want the ref to be repeated for each covariate type
covcomb
WT ALB SEX covname ID
1 56 45 0 WT 1
2 72 45 0 WT 2
3 85 45 0 WT 3
4 98 45 0 WT 4
5 128 45 0 WT 5
6 85 40 0 ALB 6
8 85 50 0 ALB 8
10 85 45 1 SEX 10

### Simulation at Unique Combinations of Covariates

For this first step, we simulate without uncertainty and without BSVzero_re() at unique combination of covariates generated in the previous section. We then plot the results to visualize the effects.

idata <- as.data.frame(covcomb)
ev1 <- ev(time=0,amt=100, cmt=1)
data.dose <- ev(ev1)
data.dose<-as.data.frame(data.dose)
data.all<-merge(idata,data.dose)
outcovcomb<- modcovsim %>%
data_set(data.all) %>%
carry.out(WT,AGE,SEX,ALB,CLi) %>%
zero_re() %>%
mrgsim(end=24, delta=1)
outcovcomb<-as.data.frame(outcovcomb)
outcovcomb <- outcovcomb %>%
arrange(ID,time,WT)
outcovcomb$SEX <- as.factor(outcovcomb$SEX )

ggplot(outcovcomb, aes(time,CP,col=factor(WT),linetype=factor(SEX)) )+
geom_line(aes(group=interaction(ID)),alpha=1,size=1.5)+
facet_grid(ALB~ WT,labeller = label_both)+
labs(linetype="Sex",colour="Weight",caption ="Simulation without Uncertainty\nwithout BSV",
x = "Time (h)", y="Plasma Concentrations")+
coord_cartesian(ylim=c(0,6))

### Adding Uncertainty from a Varcov Matrix

We will invent a varcov matrix by assuming 15% relative standard errors and correlations of 0.2 across the board. We then simulate a 100 set of parameters using a multivariate normal (kept at 100 for the vignette, use more replicates for a real project). Also, unless the model was written in a way to allow unconstrained parameter values, care should be taken to make sure the simulated parameters are valid and make sense. When available, use the set of parameters from a bootstrap run.

thmeans <- c(4,10,0.5,50,10,-0.8,0.2,1,0.07,1)
thvariances<- (thmeans*0.15)^2
thecorrelations <- matrix(ncol=length(thmeans),nrow=length(thmeans))
diag(thecorrelations)<- 1
thecorrelations[lower.tri(thecorrelations, diag = FALSE)]<- 0.2
thecorrelations[upper.tri(thecorrelations, diag = FALSE)]<- 0.2
thevarcovmatrix<- diag(sqrt(thvariances))%*%thecorrelations%*%diag(sqrt(thvariances))
as.data.frame(thevarcovmatrix)
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10
0.36000 0.18000 0.0090000 0.90000 0.18000 0.014400 3.6e-03 0.018000 0.0012600 0.018000
0.18000 2.25000 0.0225000 2.25000 0.45000 0.036000 9.0e-03 0.045000 0.0031500 0.045000
0.00900 0.02250 0.0056250 0.11250 0.02250 0.001800 4.5e-04 0.002250 0.0001575 0.002250
0.90000 2.25000 0.1125000 56.25000 2.25000 0.180000 4.5e-02 0.225000 0.0157500 0.225000
0.18000 0.45000 0.0225000 2.25000 2.25000 0.036000 9.0e-03 0.045000 0.0031500 0.045000
0.01440 0.03600 0.0018000 0.18000 0.03600 0.014400 7.2e-04 0.003600 0.0002520 0.003600
0.00360 0.00900 0.0004500 0.04500 0.00900 0.000720 9.0e-04 0.000900 0.0000630 0.000900
0.01800 0.04500 0.0022500 0.22500 0.04500 0.003600 9.0e-04 0.022500 0.0003150 0.004500
0.00126 0.00315 0.0001575 0.01575 0.00315 0.000252 6.3e-05 0.000315 0.0001102 0.000315
0.01800 0.04500 0.0022500 0.22500 0.04500 0.003600 9.0e-04 0.004500 0.0003150 0.022500

### Generating Sets of Parameters

Here we generate the sim_parameters dataset usingmvrnormand then incorporate the uncertainty by simulating using a different set of parameters (row) for each replicate.

nsim <- 100 # for vignette to make it run faster otherwise increase to 1000
sim_parameters <- MASS::mvrnorm(n = nsim, mu=as.numeric(thmeans),
Sigma=thevarcovmatrix, empirical = TRUE)
# library(mvtnorm) rmvnorm is another option that can be explored

colnames(sim_parameters) <- colnames(thevarcovmatrix) <- c("CL","V","KA",
"Vp","Qp",
"CLALB","CLSEX","CLWT",
"VSEX","VWT")

sim_parameters<- as.data.frame(sim_parameters)
head(sim_parameters)
CL V KA Vp Qp CLALB CLSEX CLWT VSEX VWT
3.660692 9.464078 0.4700449 30.53759 8.209243 -0.7920392 0.2030781 0.9514649 0.0665269 0.9091075
3.703555 10.647515 0.4572856 38.79312 10.542653 -0.8030430 0.2153613 0.6727316 0.0779609 1.0290469
3.214785 10.564932 0.4140212 46.86862 11.130283 -0.6993439 0.2139758 1.1733537 0.0581141 0.8185978
2.980833 10.525490 0.5887814 52.56325 12.077388 -0.5773783 0.2552327 0.7902157 0.0691301 0.9814391
3.176026 10.811632 0.5823588 48.35893 10.266314 -0.7019336 0.2068864 1.0969736 0.0827653 1.0836803
3.491516 6.886621 0.3800583 53.31957 10.082666 -0.7863597 0.1704918 0.8668090 0.0665500 0.8479502

### Simulation with Uncertainty

We illustrate how you can iterate over a set of parameters value using a for loop. We then overlay the previous simulation without uncertainty on the one with uncertainty to visualize the effect of adding the uncertainty.

iter_sims <- NULL
for(i in 1:nsim) {
covcomb$covname<- NULL covcomb$ID <- 1:nrow(covcomb)
idata <- as.data.frame(covcomb)
ev1 <- ev(time=0,amt=100, cmt=1)
data.dose <- ev(ev1)
data.dose<-as.data.frame(data.dose)
data.all<-merge(idata,data.dose)
data.all$amt<- 100 data.all$CL <- as.numeric(sim_parameters[i,1])
data.all$V <- as.numeric(sim_parameters[i,2]) data.all$KA <- as.numeric(sim_parameters[i,3])
data.all$Vp <- as.numeric(sim_parameters[i,4]) data.all$Qp <- as.numeric(sim_parameters[i,5])
data.all$CLALB <- as.numeric(sim_parameters[i,6]) data.all$CLSEX <- as.numeric(sim_parameters[i,7])
data.all$CLWT <- as.numeric(sim_parameters[i,8]) data.all$VSEX <- as.numeric(sim_parameters[i,9])
data.all$VWT <- as.numeric(sim_parameters[i,10]) out <- modcovsim %>% data_set(data.all) %>% carry.out(ALB,WT,SEX,amt) %>% zero_re() %>% mrgsim(start=0,end=24,delta=0.5) dfsimunc <- as.data.frame(out%>% mutate(rep = i) ) iter_sims <- rbind(iter_sims,dfsimunc) } ggplot(iter_sims, aes(time,CP,col=factor(WT),linetype=factor(SEX) ) )+ geom_line(aes(group=interaction(ID,rep)),alpha=0.3,size=0.3)+ geom_line(data=outcovcomb,aes(group=interaction(ID)),color="black")+ facet_grid(ALB~ WT,labeller = label_both)+ labs(linetype="No Uncertainty\nSex", colour="Uncertainty Replicates\nWeight", caption ="Simulation with Uncertainty\nwithout BSV" , x="Time (h)", y = "Plasma Concentrations")+ coord_cartesian(ylim=c(0,6)) ### Compute NCA PK Parameters and Boxplots Similar to the above we compute the NCA parameters, standardize by the median and provide a plot. Since the code is similar it is not shown refer to the source of the vignette to look at it. We illustrate some data manipulation to construct more informative labels that will help in the plotting. out.df.univariatecov.nca <- out.df.univariatecov.nca %>% ungroup() %>% dplyr::mutate( covname = case_when( ID== 1 ~ "Weight", ID== 2 ~ "Weight", ID== 3 ~ "REF", ID== 4 ~ "Weight", ID== 5 ~ "Weight", ID== 6 ~ "Albumin", ID== 7 ~ "Albumin", ID== 8 ~ "Sex" ), covvalue =case_when( ID== 1 ~ paste(WT,"kg"), ID== 2 ~ paste(WT,"kg"), ID== 3 ~ "85 kg-Female-45 g/L", ID== 4 ~ paste(WT,"kg"), ID== 5 ~ paste(WT,"kg"), ID== 6 ~ paste(ALB,"g/L"), ID== 7 ~ paste(ALB,"g/L"), ID== 8 ~ "Male" ) ) out.df.univariatecov.nca.long <- gather(out.df.univariatecov.nca, paramname, paramvalue, Cmax:AUC) out.df.univariatecov.nca.long$covname <-factor(as.factor(out.df.univariatecov.nca$covname ), levels = c("Weight","Sex","Albumin", "REF") ) out.df.univariatecov.nca.long$covvalue <-factor(as.factor(out.df.univariatecov.nca$covvalue ), levels = c("56 kg", "72 kg", "98 kg", "128 kg", "Male", "40 g/L", "50 g/L", "85 kg-Female-45 g/L" ) ) # the order of y ticks is controlled by the factor levels. ggplot(out.df.univariatecov.nca.long, aes(paramname,paramvalue,col=factor(paramname),shape=factor(SEX) ))+ geom_point(alpha=0.1,size=1)+ geom_hline(yintercept = 1)+ geom_boxplot()+ facet_grid(ALB~ WT,labeller = label_both)+ theme(axis.text.x=element_text(angle=30))+ labs(y="PK Parameter Value",x="",shape="Sex",col="") ### Alternative View of the Data: Distributions and Intervals Here we provide an alternative visual summary of the PK parameters. It shows the distribution, quantiles of interest and a simpler median, 90% pointinterval summary right below it. It isolates each covariate effects in one panel keeping the reference on its own. It is exactly the same data as the boxplots including AUC only. Which visual presentation do you prefer? Which one enables you to clearly see and compare the covariate effects? ggplot(out.df.univariatecov.nca.long, aes(x=paramvalue,y=covvalue,fill=factor(..quantile..),height=..ndensity..))+ facet_grid(covname~paramname,scales="free",space="free")+ annotate( "rect", xmin = 0.8, xmax = 1.25, 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("#FF0000A0", "white","white", "#0000FFA0"), labels = c("(0, 0.05]", "(0.05, 0.5]","(0.5, 0.95]", "(0.95, 1]") )+ stat_summaryh(aes(x=paramvalue,y=covvalue),fun.data="median_hilow_h", color="steelblue", geom="pointrangeh",size=0.5,position=position_nudge(y=-0.1), fun.args = list(conf.int=0.9),inherit.aes = FALSE)+ 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() ### Adding the BSV and Using forest_plot To contrast the covariate effects with random unexplained variability we add to the data the BSV intervals computed in an earlier section. We then do some data manipulation and formatting to produce a plot from the package function forest_plot. It can be overwhelming to understand what each argument does. For that reason we provide a shinyapp that will help you in setting all the options you might need. To simplify we will only keep AUC before revisiting more than one parameter plots at the end. out.df.univariatecov.nca.long$covname<- as.character(out.df.univariatecov.nca.long$covname) out.df.univariatecov.nca.long$covvalue<- as.character(out.df.univariatecov.nca.long$covvalue) coveffectsdatacovrep <- out.df.univariatecov.nca.long %>% dplyr::group_by(paramname,covname,covvalue) %>% dplyr::summarize( mid= median(paramvalue), lower= quantile(paramvalue,0.05), upper = quantile(paramvalue,0.95)) 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 <-factor(as.factor(coveffectsdatacovrepbsv$label ), levels = c("56 kg", "72 kg", "98 kg", "128 kg", "Male", "40 g/L", "50 g/L", "85 kg-Female-45 g/L", "90% of patients", "50% of patients" ) ) coveffectsdatacovrepbsv$covname <-factor(as.factor(coveffectsdatacovrepbsv$covname ), levels = c("Weight","Sex", "Albumin", "REF", "BSV") ) interval_legend_text = "Median (points)\n90% CI (horizontal lines)" interval_bsv_text = "BSV (points)\nPrediction Intervals (horizontal lines)" png("./coveffectsplot3.png",width =9 ,height = 6,units = "in",res=72) coveffectsplot::forest_plot(coveffectsdatacovrepbsv[coveffectsdatacovrepbsv$paramname=="AUC",],
ref_area = c(0.8, 1/0.8),
x_range = c(0.5,2),
strip_placement = "inside",
base_size = 18,
y_label_text_size = 12,
xlabel = "Fold Change Relative to Reference",
ref_legend_text =
"Reference (vertical line)\nClinically relevant limits\n(gray area)",
area_legend_text =
"Reference (vertical line)\nClinically relevant limits\n(gray area)",
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 = "free",
paramname_shape = FALSE,
table_position = "right",
table_text_size=4,
plot_table_ratio = 3,
show_table_facet_strip = "none",
logxscale = TRUE,
major_x_ticks = c(0.5,0.8,1/0.8,1/0.5),
return_list = FALSE)
#> Scale for 'x' is already present. Adding another scale for 'x', which will
#> replace the existing scale.
dev.off()
#> png
#>   2

## Customization of the Plots

In this section, we first show the forest_plot built-in theme, then how you to get the ggplots as a list for further editing using ggplot code. Alternatively, you can save the data as a csv and launch the app using this command run_interactiveforestplot() for point and click editing.

### Using theme_benrich along Additional Options

This is achieved by setting theme_benrich = TRUE, specifying that you want no legend legend_position = "none". With this theme active you can controltable_title text and table_title_size.


png("./coveffectsplot4.png",width =9 ,height = 6,units = "in",res=72)
coveffectsplot::forest_plot(coveffectsdatacovrepbsv
[coveffectsdatacovrepbsv$paramname=="AUC",], ref_area = c(0.8, 1/0.8), x_range = c(0.5,2), xlabel = "Fold Change Relative to Reference", x_label_text_size= 10, facet_formula = "covname~paramname", theme_benrich = TRUE, table_title_size = 15, table_title = "Median [90% CI]", interval_legend_text = interval_legend_text, interval_bsv_text = interval_bsv_text, legend_position = "none", strip_placement = "outside", base_size = 12, facet_switch = "y", facet_scales = "free_y", facet_space = "free", paramname_shape = FALSE, table_position = "right", table_text_size=4, plot_table_ratio = 3, show_table_facet_strip = "none", logxscale = TRUE, major_x_ticks = c(0.25,0.5,0.8,1/0.8,1/0.5,1/0.25), return_list = FALSE) dev.off() #> png #> 2 ### Returning a List of ggplots This is achieved setting return_list = TRUE and saving it into an object. The list will contain two objects the first being the main plot and the second the table. We show how you can modify the look of the plots using ggplot code. png("./coveffectsplot0.png",width =9 ,height = 6,units = "in",res=72) plotlists <- coveffectsplot::forest_plot(coveffectsdatacovrepbsv [coveffectsdatacovrepbsv$paramname=="AUC",],
ref_area = c(0.8, 1/0.8),
xlabel = "Fold Change Relative to Reference",
ref_legend_text = "Reference (vertical line)\nClinically relevant limits\n(gray area)",
area_legend_text = "Reference (vertical line)\nClinically relevant limits\n(gray area)",
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 = "free",
paramname_shape = FALSE,
table_position = "right",
table_text_size=4,
plot_table_ratio = 4,
show_table_facet_strip = "none",
logxscale = TRUE,
major_x_ticks = c(0.25,0.5,0.8,1/0.8,1/0.5,1/0.25),
return_list = TRUE)
plotlists
#> [[1]]
#>
#> [[2]]
dev.off()
#> png
#>   2
main_plot <- plotlists[[1]] + theme(
panel.spacing=unit(0, "pt"),
panel.grid=element_blank(),
panel.grid.minor=element_blank(),
legend.position="bottom",
strip.placement.y="outside",
strip.background.y=element_blank(),
strip.text.y=element_text(
hjust=1,
vjust=1,
face="bold",
size=rel(1)),
legend.text = element_text(size=rel(0.5)),
plot.margin = margin(t=0,r=0,b=0,l=5,unit="pt"))

table_plot <- plotlists[[2]] + theme(
panel.border=element_blank(),
panel.spacing=unit(0, "pt"),
strip.background.y=element_blank(),
legend.text = element_text(size=rel(0.5)),
plot.margin = margin(t=0,r=5,b=0,l=0,unit="pt"))

png("./coveffectsplot5.png",width =8.5 ,height = 6,units = "in",res=72)
egg::ggarrange(
main_plot,
table_plot,
nrow = 1,
widths = c(3, 1)
)
dev.off()
#> png
#>   2

## Plots with Multiple PK Parameters

You can also have plots with more than one PK parameter. You may want to facet by parameter, or to use different shape by parameter.

### Facet by Parameter

This is achieved by setting paramname_shape = FALSE and facet_formula = "covname~paramname" we also suppress the table by using table_position = "none".

png("./coveffectsplot6.png",width =9.5 ,height = 6,units = "in",res=72)
forest_plot(coveffectsdatacovrepbsv,
ref_area = c(0.8, 1/0.8),
x_range = c(0.5,2),
xlabel = "Fold Change Relative to Reference",
facet_formula = "covname~paramname",
interval_legend_text = interval_legend_text,
interval_bsv_text = interval_bsv_text,
facet_switch = "y",
facet_scales = "free_y",
facet_space = "free",
paramname_shape = FALSE,
table_position = "none",
table_text_size=4,
base_size = 9,
plot_table_ratio = 4,
show_table_facet_strip = "none",
logxscale = TRUE,
major_x_ticks = c(0.5,0.8,1/0.8,1/0.5),
x_label_text_size = 10,
return_list = FALSE)
dev.off()
#> png
#>   2

### Shape by Parameter

This is achieved by setting paramname_shape = TRUE we also illustrate how you can use legend_order to choose the legend ordering and few other options.

png("./coveffectsplot7.png",width =9.5 ,height = 6,units = "in",res=72)
forest_plot(coveffectsdatacovrepbsv
[coveffectsdatacovrepbsv\$paramname!="AUC",],
ref_area = c(0.8, 1/0.8),
x_range = c(0.35,1/0.35),
xlabel = "Fold Change Relative to Reference",
ref_legend_text = "Reference\nClinically relevant limits\n(0.8-1.25)",
area_legend_text = "Reference\nClinically relevant limits\n(0.8-1.25)",
interval_legend_text = "Median\n90% CI",
interval_bsv_text = "BSV\nPrediction Intervals",
facet_formula = "covname~.",
paramname_shape = TRUE,
legend_order =c("shape","pointinterval","ref", "area"),
legend_shape_reverse = TRUE,
bsv_col = scales::muted("red"),
interval_col = scales::muted("blue"),
facet_switch = "y",
facet_scales = "free_y",
facet_space = "free",
table_position = "none",
table_text_size=4,
base_size = 9,
plot_table_ratio = 4,
show_table_facet_strip = "none",
logxscale = TRUE,
major_x_ticks = c(0.5,0.8,1/0.8,1/0.5),
legend_space_x_mult = 0.01,
legend_position = "right",
return_list = FALSE)
dev.off()
#> png
#>   2

While we covered varying one at a time covariate value (marginal effects), we can use observed or simulated distribution of correlated covariates and simulate joint covariate effects.