## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( dev = "png", 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(egg) library(data.table) library(patchwork) library(GGally) 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) } 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") } derive.exposure <- function(time, CP) { n <- length(time) x <- c( Cmax = max(CP), Clast = CP[n], AUC = sum(diff(time) * (CP[-1] + CP[-n])) / 2 ) data.table(paramname=names(x), paramvalue=x) } cor2cov <- function (cor, sd) { if (missing(sd)) { sd <- diag(cor) } diag(cor) <- 1 n <- nrow(cor) diag(sd, n) %*% cor %*% diag(sd, n) } stat_sum_df <- function(fun, geom="point", ...) { stat_summary(fun.data = fun, geom=geom, ...) } summary_conc <- function(CP) { x <- c( Cmed = median(CP), Clow = quantile(CP, probs = 0.05), Cup = quantile(CP, probs = 0.95) ) data.table(paramname=names(x), paramvalue=x) } nbsvsubjects <- 1000 nsim <- 10 # uncertainty replicates for vignette you might want a higher number ## ----pkmodel, collapse=TRUE--------------------------------------------------- codepkmodelcov <- ' $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) $PARAM @annotated // reference values for covariate WT : 85 : Weight (kg) SEX : 0 : Sex (0=Female, 1=Male) ALB : 45 : Albumin (g/L) $PKMODEL cmt="GUT CENT PER", depot=TRUE, trans=11 $MAIN double CLi = CL * pow((ALB/45.0), CLALB)* (SEX == 1.0 ? (1.0+CLSEX) : 1.0)* pow((WT/85.0), CLWT)*exp(nCL); double V2i = V * (SEX == 1.0 ? (1.0+VSEX) : 1.0)* pow((WT/85.0), VWT)*exp(nVC); double KAi = KA; double V3i = Vp *pow((WT/85.0), 1); double Qi = Qp *pow((WT/85.0), 0.75); $OMEGA @annotated @block nCL : 0.09 : ETA on CL nVC : 0.01 0.09 : ETA on Vc $TABLE double CP = CENT/V2i; $CAPTURE CP KAi CLi V2i V3i Qi WT SEX ALB ' modcovsim <- mcode("codepkmodelcov", codepkmodelcov) partab <- setDT(modcovsim@annot$data)[block=="PARAM", .(name, descr, unit)] partab <- merge(partab, melt(setDT(modcovsim@param@data), meas=patterns("*"), var="name")) knitr::kable(partab) ## ----pksimulation, fig.width=7, message=FALSE--------------------------------- idata <- data.table(ID=1:nbsvsubjects, WT=85, SEX=0, ALB=45) ev1 <- ev(time = 0, amt = 100, cmt = 1) data.dose <- ev(ev1) data.dose <- setDT(as.data.frame(data.dose)) data.all <- data.table(idata, data.dose) outputsim <- modcovsim %>% data_set(data.all) %>% mrgsim(end = 24, delta = 0.25) %>% as.data.frame %>% as.data.table outputsim$SEX <- factor(outputsim$SEX, labels="Female") # Only plot a random sample of N=500 set.seed(678549) plotdata <- outputsim[ID %in% sample(unique(ID), 500)] # New facet label names for dose variable albumin.labs <- c("albumin: 45 ng/mL") names(albumin.labs) <- c("45") wt.labs <- c("weight: 85 kg") names(wt.labs) <- c("85") p1 <- ggplot(plotdata, aes(time, CP, group = ID)) + geom_line(alpha = 0.2, size = 0.1) + facet_grid(~ WT + ALB + SEX, labeller = labeller(ALB = albumin.labs, WT = wt.labs)) + labs(y = "Plasma Concentrations", x = "Time (h)") p2 <- ggplot(plotdata, aes(time, CP, group = ID)) + geom_line(alpha = 0.2, size = 0.1) + facet_grid(~ WT + ALB + SEX, labeller = labeller(ALB = albumin.labs, WT = wt.labs)) + scale_y_log10() + labs(y = "Plasma~Concentrations\n(logarithmic scale)", x = "Time (h)") egg::ggarrange(p1, p2, ncol = 2) ## ----computenca, fig.width=7, message=FALSE----------------------------------- derive.exposure <- function(time, CP) { n <- length(time) x <- c( Cmax = max(CP), Clast = CP[n], AUC = sum(diff(time) * (CP[-1] + CP[-n])) / 2 ) data.table(paramname=names(x), paramvalue=x) } refbsv <- outputsim[, derive.exposure(time, CP), by=.(ID, WT, SEX, ALB)] p3 <- ggplot(refbsv, 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("white", "#FF000050", "#FF0000A0", "#FF0000A0", "#FF000050", "white"), 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", axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank()) + labs(x="PK Parameters", y="") + scale_x_log10() + coord_cartesian(expand=FALSE) # Obtain the standardized parameter value by dividing by the median. refbsv[, stdparamvalue := paramvalue/median(paramvalue), by=paramname] p4 <- ggplot(refbsv, aes( x = stdparamvalue, 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("white", "#FF000050", "#FF0000A0", "#FF0000A0", "#FF000050", "white"), 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", axis.text.y = element_blank(), axis.ticks.y = element_blank(), axis.title.y = element_blank()) + labs(x="Standardized PK Parameters", y="") + scale_x_log10() + coord_cartesian(expand=FALSE, xlim = c(0.3,3)) p3 + p4 ## ----computebsvpk , fig.width=7 , message=FALSE------------------------------- 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, SEX = 0) covdatasim$SEX<- ifelse(covdatasim$SEX==0,1,0) covdatasim$SEX <- as.factor(covdatasim$SEX ) covdatasim$SEX <- factor(covdatasim$SEX,labels = c("Female","Male")) covdatasimpairs <- covdatasim covdatasimpairs$Weight <- covdatasimpairs$WT covdatasimpairs$Sex <- covdatasimpairs$SEX covdatasimpairs$Albumin <- covdatasimpairs$ALB ggpairsplot <- GGally::ggpairs(covdatasimpairs, columns = c("Weight","Sex","Albumin"),mapping = aes(colour=SEX), diag= list( continuous = GGally::wrap("densityDiag", alpha = 0.3,colour=NA), discrete = GGally::wrap("barDiag", alpha =0.3, position = "dodge2") ), lower = list( continuous = GGally::wrap("points", alpha = 0.2, size = 2), combo = GGally::wrap("facethist", alpha = 0.2, position = "dodge2") ), upper = list( continuous = GGally::wrap("cor", size = 4.75, align_percent = 0.5), combo = GGally::wrap("box_no_facet", alpha =0.3), discrete = GGally::wrap("facetbar", alpha = 0.3, position = "dodge2") ) ) covdatasim$SEX <- as.numeric(covdatasim$SEX)-1 ggpairsplot + theme_bw(base_size = 12) + theme(axis.text = element_text(size=9)) ## ----fig.width=7 ,message=FALSE, fig.height=5--------------------------------- idata <- data.table::copy(covdatasim) idata$covname <- NULL ev1 <- ev(time=0, amt=100, cmt=1) data.dose <- as.data.frame(ev1) data.all <- data.table(idata, data.dose) outcovcomb<- modcovsim %>% data_set(data.all) %>% zero_re() %>% mrgsim(end=24, delta=0.25) %>% as.data.frame %>% as.data.table outcovcomb$SEX <- as.factor(outcovcomb$SEX ) outcovcomb$SEX <- factor(outcovcomb$SEX, labels=c("Female", "Male")) stat_sum_df <- function(fun, geom="ribbon", ...) { stat_summary(fun.data = fun, geom = geom, ...) } stat_sum_df_line <- function(fun, geom="line", ...) { stat_summary(fun.data = fun, geom = geom, ...) } fwt <- function(x, xcat, which, what, from, to, ...) { what <- sub("WT", "\nWeight", what) sprintf("%s %s [%s to %s[", which, what, signif_pad(from, 3, FALSE), signif_pad(to, 3, FALSE)) } f <- function(x, xcat, which, what, from, to, ...) { what <- sub("ALB", "\nAlbumin", what) sprintf("%s %s [%s to %s[", which, what, signif_pad(from, 3, FALSE), signif_pad(to, 3, FALSE)) } plotlines<- ggplot(outcovcomb, aes(time,CP,col=SEX ) )+ geom_line(aes(group=ID),alpha=0.1,size=0.1)+ facet_grid(table1::eqcut(ALB,2,f) ~ table1::eqcut(WT,2,fwt),labeller = label_value)+ labs(colour="Sex",caption ="Simulation without Uncertainty\nFull Covariate Distribution\nwithout BSV/Uncertainty", x = "Time (h)", y="Plasma Concentrations")+ coord_cartesian(ylim=c(0,3.5)) plotranges<- ggplot(outcovcomb, aes(time,CP,col=SEX,fill=SEX ) )+ stat_sum_df(fun="median_hilow",alpha=0.2, mapping = aes(group=interaction(table1::eqcut(WT,2,fwt), SEX, table1::eqcut(ALB,2,f)) ), colour = "transparent")+ stat_sum_df_line(fun="median_hilow",size =2, mapping = aes(linetype = SEX, group=interaction(table1::eqcut(WT,2,fwt), SEX,table1::eqcut(ALB,2,f))))+ facet_grid(table1::eqcut(ALB,2,f) ~ table1::eqcut(WT,2,fwt), labeller = label_value)+ labs(linetype="Sex",colour="Sex",fill="Sex", caption ="Simulation with Full Covariate Distribution with BSV 95% (Covariate Effects + BSV) Percentiles", x = "Time (h)", y="Plasma Concentrations")+ coord_cartesian(ylim=c(0,3.5)) plotranges ## ----fig.width=7-------------------------------------------------------------- theta <- unclass(as.list(param(modcovsim))) theta[c("WT", "SEX", "ALB")] <- 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) ## ----fig.width=7-------------------------------------------------------------- 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) ## ----fig.width=7,fig.height=4, fig.height=5----------------------------------- idata <- copy(covdatasim) ev1 <- ev(time=0, amt=100, cmt=1) data.dose <- as.data.frame(ev1) iter_sims <- NULL for(i in 1:nsim) { # you might want to resample your covariate database here # e.g. subsample from a large pool of patient # include uncertainty on your covariate distribution data.all <- data.table(idata, data.dose, sim_parameters[i]) out <- modcovsim %>% data_set(data.all) %>% #zero_re() %>% #omat(rxode2::cvPost(2000, matrix(c(0.09,0.01,0.01,0.09), 2, 2), #type = "invWishart")) %>% # unc on bsv uncomment and increase nsim for CPT:PSP paper mrgsim(start=0, end=24, delta=0.25) %>% as.data.frame %>% as.data.table out[, rep := i] iter_sims <- rbind(iter_sims, out) } f <- function(x, xcat, which, what, from, to, ...) { what <- sub("ALB", "\nAlbumin", what) sprintf("%s %s [%s to %s[", which, what, signif_pad(from, 3, FALSE), signif_pad(to, 3, FALSE)) } fwt <- function(x, xcat, which, what, from, to, ...) { what <- sub("WT", "\nWeight", what) sprintf("%s %s [%s to %s[", which, what, signif_pad(from, 3, FALSE), signif_pad(to, 3, FALSE)) } iter_sims_summary_all <- iter_sims %>% mutate(WT=table1::eqcut(WT,2,fwt),ALB=table1::eqcut(ALB,2,f)) %>% group_by(time,WT,ALB,SEX)%>% summarize( P50= median(CP) , P05 = quantile(CP,0.05), P95= quantile(CP,0.95)) iter_sims_summary_all$SEX <- as.factor(iter_sims_summary_all$SEX ) iter_sims_summary_all$SEX <- factor(iter_sims_summary_all$SEX,labels = c("Female","Male")) legendlabel<- "Median\n5%-95%" plotrangesunc<- ggplot(iter_sims_summary_all, aes(time,P50,col=SEX,fill=SEX,group=SEX,linetype=SEX) )+ geom_ribbon(aes(ymin=P05,ymax=P95),alpha=0.3,linetype=0)+ geom_line(size=1)+ facet_grid(ALB ~ WT, labeller = label_value)+ labs(linetype=legendlabel,colour=legendlabel,fill=legendlabel, caption ="Simulation with joint correlated covariate distributions with uncertainty and between subject variability", x = "Time (h)", y="Plasma Concentrations")+ coord_cartesian(ylim=c(0,3.5)) plotrangesunc+ theme(axis.title.y = element_text(size=12))+ coord_cartesian(ylim=c(0,4)) ## ----fig.width=7, include=TRUE, message=FALSE--------------------------------- out.df.parameters <- iter_sims[, derive.exposure(time, CP), by=.(rep, ID, WT, SEX, ALB)] refvalues <- out.df.parameters[,.(medparam = median(paramvalue)), by=.(paramname,rep)] ## ----fig.width=7,fig.height=5 ,message=FALSE---------------------------------- setkey(out.df.parameters, paramname, rep) out.df.parameters <- merge(out.df.parameters,refvalues) out.df.parameters[, paramvaluestd := paramvalue/medparam] out.df.parameters[, SEXCAT := ifelse( SEX==0,"Female","Male")] out.df.parameters[, REF := "All Subjects"] out.df.parameters[, WTCAT4 := table1::eqcut( out.df.parameters$WT,4,varlabel = "Weight")] out.df.parameters[, ALBCAT3 := table1::eqcut( out.df.parameters$ALB,3,varlabel = "Albumin")] nca.summaries.long <- melt(out.df.parameters, measure=c("REF","WTCAT4","ALBCAT3","SEXCAT"), value.name = "covvalue",variable.name ="covname" ) nca.summaries.long$covvalue <- as.factor( nca.summaries.long$covvalue) nca.summaries.long$covvalue <- reorder(nca.summaries.long$covvalue,nca.summaries.long$paramvalue) nca.summaries.long$covvalue <- factor(nca.summaries.long$covvalue, levels =c( "1st tertile of Albumin: [31.0,44.0)" , "2nd tertile of Albumin: [44.0,46.0)" , "3rd tertile of Albumin: [46.0,54.0]" , "Male" , "Female" , "All Subjects" , "1st quartile of Weight: [40.6,71.3)" , "2nd quartile of Weight: [71.3,85.0)" , "3rd quartile of Weight: [85.0,98.2)" ,"4th quartile of Weight: [98.2,222]" )) nca.summaries.long$covvalue2 <- factor(nca.summaries.long$covvalue, labels =c( "T1: [31.0,44.0)" , "T2: [44.0,46.0)" , "T3: [46.0,54.0]" , "Male" , "Female" , "All Subjects" , "Q1: [40.6,71.3)" , "Q2: [71.3,85.0)" , "Q3: [85.0,98.2)" , "Q4: [98.2,222]" )) nca.summaries.long$covname<- as.factor(nca.summaries.long$covname) nca.summaries.long$covname<- factor(nca.summaries.long$covname, levels =c("WTCAT4","SEXCAT","ALBCAT3","REF"), labels = c("Weight","Sex","Albumin","REF")) func <- function(bob) c(min(bob), median(bob), max(bob)) boxplotMV<- ggplot(nca.summaries.long , aes(x=covvalue2 , y=paramvalue ))+ facet_grid ( paramname ~covname, scales="free", labeller=label_parsed, switch="both",space="free_x") + geom_boxplot(outlier.shape = NA) + theme_bw(base_size = 12)+ theme(axis.title=element_blank(), strip.placement = "outside", axis.text.x = element_text(angle=20,vjust = 1, hjust = 1, face = "bold"), strip.text.y.left = element_text(angle= 0,vjust = 1, hjust = 1,face = "bold"))+ scale_y_continuous(breaks = scales::pretty_breaks(n=4) ) boxplotMV ## ----fig.width=7, fig.height=4 ,message=FALSE--------------------------------- ggridgesplot<- ggplot(nca.summaries.long, aes(x=paramvaluestd,y=covvalue, fill=factor(..quantile..), height=..ndensity..))+ facet_grid(covname~paramname,scales="free_y")+ 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.01,scale=0.9, quantiles = c(0.05,0.5, 0.95))+ geom_vline( aes(xintercept = 1),size = 1)+ 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(data=data.frame (xintercept=1), aes(xintercept =xintercept ),size = 1)+ theme_bw()+ theme(legend.position = "none")+ labs(x="Effects Of Covariates on PK Parameter",y="")+ scale_x_continuous(breaks=c(0.5,0.8,1/0.8,1/0.5,1/0.25),trans ="log" )+ coord_cartesian(xlim=c(0.25,3)) ggridgesplot ## ----fig.width=7, fig.height=6------------------------------------------------ coveffectsdatacovrep <- nca.summaries.long %>% dplyr::group_by(paramname,covname,covvalue) %>% dplyr::summarize( mid= median(paramvaluestd), lower= quantile(paramvaluestd,0.05), upper = quantile(paramvaluestd,0.95)) %>% dplyr::filter(!is.na(mid)) coveffectsdatacovrepbsv <- coveffectsdatacovrep[coveffectsdatacovrep$covname=="REF",] coveffectsdatacovrepbsv$covname <- "BSV" coveffectsdatacovrepbsv$covvalue <- "90% of patients" coveffectsdatacovrepbsv$label <- "90% of patients" coveffectsdatacovrepbsv$lower <- bsvranges$P05 coveffectsdatacovrepbsv$upper <- bsvranges$P95 coveffectsdatacovrepbsv2 <- coveffectsdatacovrep[coveffectsdatacovrep$covname=="REF",] coveffectsdatacovrepbsv2$covname <- "BSV" coveffectsdatacovrepbsv2$covvalue <- "50% of patients" coveffectsdatacovrepbsv2$label <- "50% of patients" coveffectsdatacovrepbsv2$lower <- bsvranges$P25 coveffectsdatacovrepbsv2$upper <- bsvranges$P75 coveffectsdatacovrepbsv<- rbind(coveffectsdatacovrep,coveffectsdatacovrepbsv2, coveffectsdatacovrepbsv) 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 <- gsub(": ", ":\n", coveffectsdatacovrepbsv$label) coveffectsdatacovrepbsv$covname <-factor(as.factor(coveffectsdatacovrepbsv$covname ), levels = c("Weight","Sex","Albumin","REF","BSV")) coveffectsdatacovrepbsv$label <- factor(coveffectsdatacovrepbsv$label, levels =c( "1st tertile of Albumin:\n[31.0,44.0)" , "2nd tertile of Albumin:\n[44.0,46.0)" , "3rd tertile of Albumin:\n[46.0,54.0]" , "Male", "Female" , "All Subjects","90% of patients","50% of patients" , "1st quartile of Weight:\n[40.6,71.3)" , "2nd quartile of Weight:\n[71.3,85.0)" , "3rd quartile of Weight:\n[85.0,98.2)" ,"4th quartile of Weight:\n[98.2,222]" )) coveffectsdatacovrepbsv$label <- factor(coveffectsdatacovrepbsv$label, labels =c("T1:\n[31.0,44.0)" , "T2:\n[44.0,46.0)" , "T3:\n[46.0,54.0]" , "Male", "Female" , "All Subjects","90% of patients","50% of patients" , "Q1:\n[40.6,71.3)" , "Q2:\n[71.3,85.0)" , "Q3:\n[85.0,98.2)" , "Q4:\n[98.2,222]" )) interval_legend_text <- "Median (points)\n90% intervals (horizontal lines) of joint effects: covariate distributions, uncertainty and between subject variability" interval_bsv_text <- "BSV (points)\nPrediction Intervals (horizontal lines)" ref_legend_text <- "Reference (vertical line)\nClinically relevant limits\n(gray area)" area_legend_text <- "Reference (vertical line)\nClinically relevant limits\n(gray area)" #emf("Figure_PKdist_forest.emf",width= 15, height = 7.5) png("./coveffectsplot_full.png",width =9.5 ,height = 8,units = "in",res=72) forest_plot(coveffectsdatacovrepbsv[coveffectsdatacovrepbsv$covname!="REF"& coveffectsdatacovrepbsv$covname!="BSV", ], ref_area = c(0.8, 1/0.8),x_range = c(0.4,3), strip_placement = "outside",base_size = 18, y_label_text_size = 9,x_label_text_size = 10, xlabel = "Fold Change Relative to Reference", ref_legend_text =ref_legend_text, area_legend_text =ref_legend_text , interval_legend_text = interval_legend_text, plot_title = "", interval_bsv_text = interval_bsv_text, facet_formula = "covname~paramname", facet_switch = "y", table_facet_switch = "both", reserve_table_xaxis_label_space = FALSE, facet_scales = "free_y", facet_space = "free", paramname_shape = FALSE, table_position = "below", show_table_yaxis_tick_label = TRUE, table_text_size= 4, plot_table_ratio = 1, show_table_facet_strip = "both", logxscale = TRUE, major_x_ticks = c(0.5,0.8,1/0.8,1/0.5), return_list = FALSE) dev.off()