Nothing
R/braidReport.R
In braidReports: Visualize Combined Action Response Surfaces and Report BRAID Analyses
Defines functions
basicBraidAnalysis
runBraidAnalysis.default
runBraidAnalysis.formula
runBraidAnalysis
recastFillScale
recastPositionScale
makeBraidReport
Documented in
basicBraidAnalysis
makeBraidReport
runBraidAnalysis
runBraidAnalysis.default
runBraidAnalysis.formula
#' Render a BRAID Report
#'
#' Produces a one page report depicting the results of a full BRAID analysis
#' for a single combination.
#'
#' @param analysis An object of class `braidAnalysis` produced by the
#' [runBraidAnalysis()] or [basicBraidAnalysis()] functions
#' @param compounds A length-2 character vector containing the names of the
#' two compounds tested in the combination
#' @param levels Two levels at which the IAE should be evaluated
#' @param limits Two values representing the maximal achievable concentrations
#' for the compounds tested, used to esitmate the IAE
#' @param control A named list of additional control parameters adjusting the
#' appearance of the resulting report
#'
#' @details
#' This function attempts, however foolhardily, to encompass many of the
#' details, plots, and values that the user might wish to report for a complete
#' BRAID analysis of a given drug combination. All reports are built for a
#' single 8.5-by-11 inch page, either in landscape or potrait orientation, but
#' reports can be customized to contain more or less information. Here is a
#' full list of what *can* appear in the BRAID report:
#'
#' * A raw and smoothed plot of the actual measured response data; the raw plot
#' is built using [geom_braid()] or, for irregularly laid out data,
#' [geom_braid_glass()], while the smoothed data is built using
#' [geom_braid_smooth()]. (Included in all layouts)
#' * A plot of residual errors and a smoothed surface plot of the predicted
#' additive surface based on the dose response behavior of the individual
#' compounds alone. In cases of pronounced interaction, this will differ
#' significantly from the best-fit BRAID plots. (Included in the dense layout)
#' * A plot of residual errors and a smoothed surface plot of the best-fitting
#' BRAID surface. (Included in the all layouts)
#' * A table of the best-fitting BRAID response surface parameters (Included in
#' all layouts)
#' * A table of estimated IAE values at the specified effect levels (Included in
#' all layouts)
#' * Two tables of the dose required of one drug to produce a desired effect
#' level (the first value in `levels`) in the presence of several doses of the
#' other drug; used to gauge the degree of potentiation. (Included in the
#' standard and dense layouts)
#' * Two plots depicting the predicted dose response of one drug in the presence
#' of various levels of the other, also used to gauge potentiation. (Included
#' in the standard and dense layouts)
#'
#' So the resulting report page can contain from six (simple layout) to twelve
#' (dense layout) elements depicting different aspects of the BRAID analysis.
#'
#' The precise appearance of the report page is controlled by various elements
#' of the `control` parameter. Though the default value of the parameter is
#' an empty list, several fields will be filled in if they are unspecified. The
#' full set of possible control options is:
#'
#' * `abbs`: A pair of character strings specifying the abbreviations of the
#' tested compounds. By default, the abbreviations consist of the firs three
#' characters of each compound's name, but for some compound names this is not
#' an appropriate abbreviation Abbreviations are used in many axis labels and
#' tables
#' * `units`: If included, a single string or pair of strings specifying the
#' dose units for the two drugs tested, included in axis labels and tables. If
#' left unspecified, units will not be included
#' * `xscale`: Either a character string specifying a named transformation
#' object (e.g "log2", "log10", "sqrt") or a `ggplot2` continuous x-position
#' scale generated by [ggplot2::scale_x_continuous()] or related functions. This
#' scale will be applied to the x-dimension of all surface plots and the
#' x-dimension of the first potentiation plot. If a `name` is specified for the
#' scale, this will be the x-axis label; otherwise other defaults will be used.
#' Default value for this control parameter is "log10".
#' * `yscale`: Either a character string specifying a named transformation
#' object (e.g "log2", "log10", "sqrt") or a `ggplot2` continuous y-position
#' scale generated by [ggplot2::scale_y_continuous()] or related functions. This
#' scale will be applied to the y-dimension of all surface plots and the
#' x-dimension of the second potentiation plot. If a `name` is specified for the
#' scale, this will be the y-axis label; otherwise other defaults will be used.
#' Default value for this control parameter is "log10".
#' * `fillscale`: If included, continuous fill scale object generated by one of
#' several `ggplot2` continuous fill scale functions. This fill scale will
#' control the fill appearance for all *effect* surface plots; fill colors
#' in residual error plots will use a different color palette. In addition,
#' any names, labels, breaks, transformations, etc. included in this scale will
#' also be applied to the y-axis of both potentiation plots, as those also
#' represent the modeled effect. If unspecified, will be set to a brewer
#' continuous color scale with palette "RdYlBu".
#' * `colorscale`: If included, a discrete color scale object generated by one
#' of several `ggplot2` discrete color scale functions. This color scale
#' controls the colors chosen for the curves in the two potentiation plots. If
#' left unspecified, will default to the output of
#' [ggplot2::scale_color_discrete()]
#' * `xname`: A string specifying the desired x-axis labels in surface plots.
#' Will be overridden if control parameter `xscale` is a scale object with a
#' non-empty `name` attribute. If unspecified, defaults the abbreviation of the
#' first compound followed by the units if included.
#' * `yname`: A string specifying the desired y-axis labels in surface plots.
#' Will be overridden if control parameter `yscale` is a scale object with a
#' non-empty `name` attribute. If unspecified, defaults the abbreviation of the
#' second compound followed by the units if included.
#' * `effectname`: The name of the modeled effect variable. Could be "Growth"
#' or "Survival" or "Activity". The default value is simply "Effect"
#' * `title`: A string containing the overall title of the report page. If left
#' unspecified, will simply be the first compound "vs." the second
#' * `contourcolor`: Contours of the smoothed surfaces at the levels specified
#' by the parameter `levels` are included in all smoothed plots. By default,
#' they are black, but specifying this control parameter will change that color
#' * `irregular`: If `TRUE`, the data are not assumed to lie on a regular grid
#' in the plotted, and will be visualized with [geom_braid_glass()] rather than
#' [geom_braid()]
#' * `swidth`: A numeric value specifying how widely the smooth surfaces should
#' be smoothed, passed as the `width` aesthetic to [geom_braid_smooth()]
#' * `sheight`: A numeric value specifying how far in the y-direction the smooth
#' surfaces should be smoothed, passed as the `height` aesthetic to
#' [geom_braid_smooth()]
#' * `npoints`: The density of points used in the smoothed plots. See
#' [geom_braid_smooth()] for details
#' * `leveltext`: A pair of strings indicating how the two IAE levels should be
#' displayed in tables. In some cases, the precise number at which the IAE is
#' calculated does not reflect the level that the user wishes to express. So
#' a user might want to refer to a relative survival value of 0.1 as IAE90 (for
#' 90% killing); or a log2-fold growth inhibition as IAE50 (for 50% inhibition);
#' passing "50" or "90" as the `leveltext` in such cases will produce the
#' desired appearance. If left unspecified, labels will simply be the string
#' representations of the paramter `levels`
#' * `metrics`: A named numeric or named character vector specifying additional
#' metrics for the combination; they will be added to the table containing the
#' calculated IAE values. Numeric values will be rounded to three significant
#' figures; string values will be included as is. Names will be parsed by
#' default, so the string "CI\[90\]" will be displayed with the "90" as a
#' subscript
#' * `surftheme`: A `ggplot2` `theme` object specifying any additional theme
#' adjustments to add to all response surface plots
#' * `curvetheme`: A `ggplot2` `theme` object specifying any additional theme
#' adjustments to add to both potentiation curve plots
#' * `layout`: The specific layout to be used; determines which report elements
#' are included. Can be "simple", "standard" (the default), or "dense"
#' * `orientation`: The expected orientantion of the rendered page; can be
#' "portrait" (the default) or "landscape"
#'
#'
#' @return A graphical object containing all plots and tables, arranged
#' according to the desired format. The resulting object is optimized for a
#' single page, either portrait or landscape as specified in `control`
#' @export
#'
#' @examples
#' surface <- synergisticExample
#' analysis <- runBraidAnalysis(measure~concA+concB, surface,
#' defaults=c(0,1), getCIs=FALSE)
#'
#' report <- makeBraidReport(analysis,c("A Drug","B Drug"),
#' levels=c(0.5, 0.9),limits=c(5,5))
#' print(report)
#'
#' control <- list(abbs=c("A","B"),units=c("\u00B5M"),leveltext=c("50","90"),
#' xscale=scale_x_log10(breaks=c(0.1,0.5,2,10),
#' labels=as.character),
#' fillscale=scale_fill_viridis_c(option="A"),
#' colorscale=scale_color_brewer(palette="Set1"),
#' title="Example Analysis")
#' nextReport <- makeBraidReport(analysis,c("A Drug","B Drug"),
#' levels=c(0.5, 0.9),limits=c(5,5),
#' control=control)
#' print(nextReport)
makeBraidReport <- function(analysis,compounds,levels,limits,control=list()) {
newerggplot <- utils::packageVersion("ggplot2")>="3.5.0"
if (!inherits(analysis,"braidAnalysis")) {
stop("Parameter 'analysis' must be an object of class 'braidAnalysis'.")
}
# Control options:
#
# xscale, yscale, fillscale, colorscale
# xname, yname, effectname, title
# contourcolor
# abbs, units
# irregular, swidth, sheight, npoints
# leveltext, metrics
# layout, orientation
defaultControls <- list(
xscale = "log10",
yscale = "log10",
contourcolor = "black",
npoints = 50,
layout = "standard",
orientation = "portrait"
)
defaultControls[names(control)] <- control
control <- defaultControls
if (!is.null(control$base_size)) {
base_size <- control$base_size
} else if (control$layout=="standard") {
base_size <- 8
} else if (control$layout=="dense") {
base_size <- 7
} else if (control$layout=="simple") {
base_size <- 9
} else { stop(paste0("Unrecognized layout options '",control$layout,"'.")) }
braidPar <- stats::coef(analysis$braidFit)
braidDir <- ifelse (braidPar[["Ef"]]>=braidPar[["E0"]],1,-1)
# build additive parameter vector
additivePar <- braidPar
if (is.null(analysis$hillFit1) || is.null(analysis$hillFit2)) {
additiveDir <- braidDir
} else {
hfit1 <- analysis$hillFit1
hfit2 <- analysis$hillFit2
hpar1 <- stats::coef(analysis$hillFit1)
hpar2 <- stats::coef(analysis$hillFit2)
hDir1 <- ifelse(hpar1[["Ef"]]>=hpar1[["E0"]],1,-1)
hDir2 <- ifelse(hpar2[["Ef"]]>=hpar2[["E0"]],1,-1)
if (stats::var(hfit1$fitted.values)==0 | stats::var(hfit1$act)==0) { cor1 <- 0 }
else { cor1 <- stats::cor(hfit1$act,hfit1$fitted.values) }
if (stats::var(hfit2$fitted.values)==0 | stats::var(hfit2$act)==0) { cor2 <- 0 }
else { cor2 <- stats::cor(hfit2$act,hfit2$fitted.values) }
if (cor1>=cor2) {
additiveDir <- hDir1
additivePar[["E0"]] <- hpar1[["E0"]]
} else {
additiveDir <- hDir2
additivePar[["E0"]] <- hpar2[["E0"]]
}
}
additivePar[["kappa"]] <- 0
if (!is.null(analysis$hillFit1)) {
hpar1 <- stats::coef(analysis$hillFit1)
if (sign(additiveDir*(hpar1[["Ef"]]-additivePar[["E0"]]))>=0) {
additivePar[c("IDMA","na","EfA")] <- hpar1[c("IDM","n","Ef")]
}
}
if (!is.null(analysis$hillFit2)) {
hpar2 <- stats::coef(analysis$hillFit2)
if (sign(additiveDir*(hpar2[["Ef"]]-additivePar[["E0"]]))>=0) {
additivePar[c("IDMB","nb","EfB")] <- hpar2[c("IDM","n","Ef")]
}
}
additivePar[["Ef"]] <- additiveDir*max(additiveDir*additivePar[c("EfA","EfB")])
mainDF <- data.frame(conc1=analysis$concs[,1],conc2=analysis$concs[,2],act=analysis$act)
mainDF$bfit <- evalBraidModel(mainDF$conc1,mainDF$conc2,braidPar)
mainDF$afit <- evalBraidModel(mainDF$conc1,mainDF$conc2,additivePar)
mainDF$berror <- mainDF$act-mainDF$bfit
mainDF$aerror <- mainDF$act-mainDF$afit
r2add <- stats::cov.wt(cbind(mainDF$act,mainDF$afit),wt=analysis$weights,cor=TRUE)$cor[1,2]
r2brd <- stats::cov.wt(cbind(mainDF$act,mainDF$bfit),wt=analysis$weights,cor=TRUE)$cor[1,2]
r2add <- signif(r2add^2,3)
r2brd <- signif(r2brd^2,3)
surflabels <- list()
surfscales <- list()
if (!is.null(control$units) && length(control$units)==1) {
control$units <- c(control$units,control$units)
}
if (is.null(control$xname)) {
if (is.null(control$abbs)) { surflabels$x <- substring(compounds[[1]],1,3) }
else ( surflabels$x <- control$abbs[[1]])
if (!is.null(control$units)) {
surflabels$x <- paste0(surflabels$x," (",control$units[[1]],")")
}
} else { surflabels$x <- control$xname }
if (is.null(control$yname)) {
if (is.null(control$abbs)) { surflabels$y <- substring(compounds[[2]],1,3) }
else ( surflabels$y <- control$abbs[[2]])
if (!is.null(control$units)) {
surflabels$y <- paste0(surflabels$y," (",control$units[[2]],")")
}
} else { surflabels$y <- control$yname }
if (is.null(control$effectname)) { surflabels$fill <- "Effect" }
else { surflabels$fill <- control$effectname }
if (is.character(control$xscale)) {
if (newerggplot) {
surfscales$x <- scale_x_continuous(transform=control$xscale)
} else {
surfscales$x <- scale_x_continuous(trans=control$xscale)
}
} else if (inherits(control$xscale,"ScaleContinuousPosition")) {
surfscales$x <- control$xscale
} else {
stop("Control parameter \"xscale\" must be a transform string or a ggplot2 position scale,")
}
if (is.character(control$yscale)) {
if (newerggplot) {
surfscales$y <- scale_y_continuous(transform=control$yscale)
} else {
surfscales$y <- scale_y_continuous(trans=control$yscale)
}
} else if (inherits(control$yscale,"ScaleContinuousPosition")) {
surfscales$y <- control$yscale
} else {
stop("Control parameter \"yscale\" must be a transform string or a ggplot2 position scale,")
}
if (is.null(control$fillscale)) {
surfscales$fill <- scale_fill_distiller(palette="RdYlBu",direction=-braidDir)
} else if (inherits(control$fillscale,"ScaleContinuous")) {
surfscales$fill <- control$fillscale
} else {
stop("Control parameter \"fillscale\" must be ggplot2 fill scale,")
}
valrange <- range(c(mainDF$act,mainDF$afit,mainDF$bfit))
errrange <- range(c(mainDF$aerror,mainDF$berror))
blankdf <- data.frame(conc1=stats::median(mainDF$conc1,na.rm=TRUE),
conc2=stats::median(mainDF$conc2,na.rm=TRUE),
act=valrange,bfit=valrange,afit=valrange,
berror=errrange,aerror=errrange)
if (is.null(control$irregular)) { irregular <- FALSE }
else { irregular <- control$irregular }
if (is.null(control$swidth)) { swidth <- NULL }
else { swidth <- control$swidth }
if (is.null(control$sheight)) { sheight <- NULL }
else { sheight <- control$sheight }
if (is.null(sheight)) {
if (is.null(swidth)) {
smooth_aes <- aes()
} else {
smooth_aes <- aes(width=swidth)
}
} else {
if (is.null(swidth)) {
smooth_aes <- aes(height=sheight)
} else {
smooth_aes <- aes(width=swidth,height=sheight)
}
}
crel <- is.finite(log(mainDF$conc1)) & is.finite(log(mainDF$conc2))
if (length(limits)==1) { limits <- c(limits,limits) }
limits <- limits[1:2]
limitdf <- data.frame(conc1=c(),conc2=c())
if (any(crel)) {
crng1 <- range(mainDF$conc1[crel])
crng2 <- range(mainDF$conc2[crel])
if (!any(limits<c(crng1[[1]],crng2[[1]])) &&
any(limits<=c(crng1[[2]],crng2[[2]]))) {
checkvec <- limits<=c(crng1[[2]],crng2[[2]])
if (all(checkvec)) {
limitdf <- data.frame(conc1=c(0,limits[[1]],limits[[1]]),
conc2=c(limits[[2]],limits[[2]],0))
} else if (checkvec[[1]]) {
limitdf <- data.frame(conc1=c(limits[[1]],limits[[1]]),
conc2=c(0,Inf))
} else {
limitdf <- data.frame(conc1=c(0,Inf),
conc2=c(limits[[2]],limits[[2]]))
}
}
}
plots <- list()
plots[[1]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="act"))+geom_blank(data=blankdf)
plots[[2]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="act",z="act"))+geom_blank(data=blankdf)
plots[[3]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="aerror"))+ geom_blank(data=blankdf)
plots[[4]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="afit",z="afit"))+geom_blank(data=blankdf)
plots[[5]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="berror"))+geom_blank(data=blankdf)
plots[[6]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="bfit",z="bfit"))+geom_blank(data=blankdf)
if (irregular) {
plots[[1]] <- plots[[1]]+geom_braid_glass()
plots[[3]] <- plots[[3]]+geom_braid_glass()
plots[[5]] <- plots[[5]]+geom_braid_glass()
} else {
plots[[1]] <- plots[[1]]+geom_braid()
plots[[3]] <- plots[[3]]+geom_braid()
plots[[5]] <- plots[[5]]+geom_braid()
}
plots[[2]] <- plots[[2]]+
geom_braid_smooth(smooth_aes,npoints=control$npoints)+
geom_braid_contour(smooth_aes,npoints=control$npoints,breaks=levels,colour=control$contourcolor)+
annotate("contour",x=limitdf$conc1,y=limitdf$conc2,colour=control$contourcolor,linetype=2)
plots[[4]] <- plots[[4]]+
geom_braid_smooth(smooth_aes,npoints=control$npoints)+
geom_braid_contour(smooth_aes,npoints=control$npoints,breaks=levels,colour=control$contourcolor)+
annotate("contour",x=limitdf$conc1,y=limitdf$conc2,colour=control$contourcolor,linetype=2)
plots[[6]] <- plots[[6]]+
geom_braid_smooth(smooth_aes,npoints=control$npoints)+
geom_braid_contour(smooth_aes,npoints=control$npoints,breaks=levels,colour=control$contourcolor)+
annotate("contour",x=limitdf$conc1,y=limitdf$conc2,colour=control$contourcolor,linetype=2)
plots[[1]] <- plots[[1]]+
labs(x=surflabels$x,y=surflabels$y,fill=surflabels$fill,title="Measured Data")+
surfscales$x+surfscales$y+surfscales$fill
plots[[2]] <- plots[[2]]+
labs(x=surflabels$x,y=surflabels$y,fill=surflabels$fill,title="Smoothed Data")+
surfscales$x+surfscales$y+surfscales$fill
plots[[3]] <- plots[[3]]+
labs(x=surflabels$x,y=surflabels$y,title="Additive Error")+
surfscales$x+surfscales$y+scale_fill_gradient2("Error")
plots[[4]] <- plots[[4]]+
labs(x=surflabels$x,y=surflabels$y,fill=surflabels$fill,
title=paste0("Additive Fit (R\u00B2=",r2add,")"))+
surfscales$x+surfscales$y+surfscales$fill
plots[[5]] <- plots[[5]]+
labs(x=surflabels$x,y=surflabels$y,title="BRAID Error")+
surfscales$x+surfscales$y+scale_fill_gradient2("Error")
plots[[6]] <- plots[[6]]+
labs(x=surflabels$x,y=surflabels$y,fill=surflabels$fill,
title=paste0("BRAID Fit (R\u00B2=",r2brd,")"))+
surfscales$x+surfscales$y+surfscales$fill
for (i in 1:6) {
plots[[i]] <- plots[[i]]+
theme(text=element_text(size=base_size),
legend.key.size=grid::unit(0.8,"lines"))
if (!is.null(control$surftheme) && inherits(control$surftheme,"theme")) {
plots[[i]] <- plots[[i]]+control$surftheme
}
}
paramdf <- data.frame(name=c("ID['M,A']","ID['M,B']","n[a]","n[b]","kappa",
"E[0]","E['f,A']","E['f,B']","E[f]"),
value=signif(stats::coef(analysis$braidFit),3))
bmodel <- analysis$braidFit$model
if (!is.null(analysis$braidFit$ciMat)) {
paramdf$lo <- NA_real_
paramdf$hi <- NA_real_
paramdf$lo[bmodel] <- signif(pmin(analysis$braidFit$ciMat[,1],paramdf$value[bmodel]),3)
paramdf$hi[bmodel] <- signif(pmax(analysis$braidFit$ciMat[,2],paramdf$value[bmodel]),3)
paramdf$text <- ifelse(is.na(paramdf$lo),
as.character(paramdf$value),
paste0(as.character(paramdf$value)," (",
as.character(paramdf$lo)," \u2013 ",
as.character(paramdf$hi),")"))
} else {
paramdf$text <- as.character(paramdf$value)
}
if (!any((7:8)%in%bmodel)) {
paramdf <- paramdf[c(1:6,9),]
} else if (!(9%in%bmodel)) {
paramdf <- paramdf[1:8,]
}
paramvec <- paramdf$text
names(paramvec) <- paramdf$name
partab <- parameterTable(paramvec,parse=TRUE,title="Best-Fit BRAID",base_size=base_size)
if (is.null(control$leveltext)) {
leveltext <- as.character(levels)
} else {
leveltext <- control$leveltext
}
metricdf <- data.frame(name=paste0("IAE['",leveltext,"']"),
value=signif(estimateIAE(braidPar,levels,limits),3))
if (!is.null(analysis$braidFit$ciMat)) {
iaemat <- calcBraidConfInt(
analysis$braidFit,
function(par) estimateIAE(par,levels,limits)
)
metricdf$lo <- signif(iaemat[,1],3)
metricdf$hi <- signif(iaemat[,3],3)
} else {
metricdf$lo <- NA_real_
metricdf$hi <- NA_real_
}
metricdf$text <- ifelse(is.na(metricdf$lo),
as.character(metricdf$value),
paste0(as.character(metricdf$value)," (",
as.character(metricdf$lo)," \u2013 ",
as.character(metricdf$hi),")"))
if (!is.null(control$metrics)) {
if (is.null(names(control$metrics))) {
warning("Control value 'metrics' must be a named vector. Values ignored.")
} else if (is.numeric(control$metrics)) {
metricdf <- rbind(metricdf,data.frame(name=names(control$metrics),
value=as.numeric(control$metrics),
lo=NA_real_,hi=NA_real_,
text=as.character(control$metrics)))
} else if (is.character(control$metrics)) {
metricdf <- rbind(metricdf,data.frame(name=names(control$metrics),
value=NA_real_,
lo=NA_real_,hi=NA_real_,
text=control$metrics))
} else {
warning("Control value 'metrics' must either be numeric or a character vector. Values ignored.")
}
}
metricvec <- metricdf$text
names(metricvec) <- metricdf$name
metrictab <- parameterTable(metricvec,parse=TRUE,title="Summary Metrics",base_size=base_size)
spacetab1 <- gtable::gtable(widths=unit(10,"mm"),heights=unit(10,"mm"),name="spacer")
spacetab2 <- gtable::gtable(widths=unit(10,"mm"),heights=unit(5,"mm"),name="spacer")
if (control$layout != "simple") {
pls <- suppressWarnings(layer_scales(plots[[1]]))
plsx <- pls$x
xbreaks <- plsx$break_info()$major_source
xbreaks <- plsx$trans$inverse(xbreaks)
xvector <- plsx$range$range
xvector <- seq(xvector[[1]],xvector[[2]],length=control$npoints)
xvector <- plsx$trans$inverse(xvector)
xvector[xvector<0] <- 0
xvector <- unique(xvector)
plsy <- pls$y
ybreaks <- plsy$break_info()$major_source
ybreaks <- plsy$trans$inverse(ybreaks)
yvector <- plsy$range$range
yvector <- seq(yvector[[1]],yvector[[2]],length=control$npoints)
yvector <- plsy$trans$inverse(yvector)
yvector[yvector<0] <- 0
yvector <- unique(yvector)
xbreaks <- sort(unique(c(0,xbreaks[
xbreaks>=min(analysis$concs[,1],na.rm=TRUE) &
xbreaks<=max(analysis$concs[,1],na.rm=TRUE)
])))
if (length(xbreaks)==6) { xbreaks <- xbreaks[c(1,3,4,5)] }
else if (length(xbreaks)>6) { xbreaks <- xbreaks[seq(1,length(xbreaks),by=2)] }
ybreaks <- sort(unique(c(0,ybreaks[
ybreaks>=min(analysis$concs[,2],na.rm=TRUE) &
ybreaks<=max(analysis$concs[,2],na.rm=TRUE)
])))
if (length(ybreaks)==6) { ybreaks <- ybreaks[c(1,3,4,5)] }
else if (length(ybreaks)>6) { ybreaks <- ybreaks[seq(1,length(ybreaks),by=2)] }
iaedf1 <- data.frame(name=ybreaks,
value=signif(invertBraidModel(DB=ybreaks,
effect=levels[[1]],
bpar=braidPar),3))
if (!is.null(analysis$braidFit$ciMat)) {
iaeci1 <- calcBraidConfInt(
analysis$braidFit,
function(p) invertBraidModel(DB=ybreaks,effect=levels[[1]],bpar=p)
)
iaedf1$lo <- signif(iaeci1[,1],3)
iaedf1$hi <- signif(iaeci1[,3],3)
} else {
iaedf1$lo <- NA_real_
iaedf1$hi <- NA_real_
}
iaedf1$text <- ifelse(is.na(iaedf1$lo),
as.character(iaedf1$value),
paste0(as.character(iaedf1$value)," (",
as.character(iaedf1$lo)," \u2013 ",
as.character(iaedf1$hi),")"))
iaedf1 <- iaedf1[,c("name","text")]
iaetitle1 <- paste0("IC['",leveltext[[1]],"']~of~'",compounds[[1]],"'")
iaetab1 <- dataFrameTable(iaedf1,title=iaetitle1,names=c(surflabels$y,"Estimate"),parse="title",base_size=base_size)
iaedf2 <- data.frame(name=xbreaks,
value=signif(invertBraidModel(DA=xbreaks,
effect=levels[[1]],
bpar=braidPar),3))
if (!is.null(analysis$braidFit$ciMat)) {
iaeci2 <- calcBraidConfInt(
analysis$braidFit,
function(p) invertBraidModel(DA=xbreaks,effect=levels[[1]],bpar=p)
)
iaedf2$lo <- signif(iaeci2[,1],3)
iaedf2$hi <- signif(iaeci2[,3],3)
} else {
iaedf2$lo <- NA_real_
iaedf2$hi <- NA_real_
}
iaedf2$text <- ifelse(is.na(iaedf2$lo),
as.character(iaedf2$value),
paste0(as.character(iaedf2$value)," (",
as.character(iaedf2$lo)," \u2013 ",
as.character(iaedf2$hi),")"))
iaedf2 <- iaedf2[,c("name","text")]
iaetitle2 <- paste0("IC['",leveltext[[1]],"']~of~'",compounds[[2]],"'")
iaetab2 <- dataFrameTable(iaedf2,title=iaetitle2,names=c(surflabels$x,"Estimate"),parse="title",base_size=base_size)
plotscales <- list()
plotscales$x1 <- surfscales$x
plotscales$x2 <- recastPositionScale(surfscales$y,"x")
plotscales$y <- recastPositionScale(surfscales$fill,"y")
if (is.null(control$colorscale)) {
plotscales$color <- scale_color_discrete()
} else {
plotscales$color <- control$colorscale
}
plotscales$fill <- recastFillScale(plotscales$color)
curvedf1 <- data.frame(
conc1 = rep(xvector,times=length(ybreaks)),
conc2 = rep(ybreaks,each=length(xvector))
)
curvedf1$act <- evalBraidModel(curvedf1$conc1,curvedf1$conc2,braidPar)
if (!is.null(analysis$braidFit$ciMat)) {
curveci1 <- calcBraidConfInt(
analysis$braidFit,
function(p) evalBraidModel(curvedf1$conc1,curvedf1$conc2,bpar=p)
)
curvedf1$act_lo <- curveci1[,1]
curvedf1$act_hi <- curveci1[,3]
}
curvedf1$conc2 <- factor(curvedf1$conc2)
curveplot1 <- ggplot(curvedf1,aes_string(x="conc1",y="act"))+
geom_hline(yintercept=levels,linetype=2,colour="black")
if (!is.null(analysis$braidFit$ciMat)) {
curveplot1 <- curveplot1 +
geom_ribbon(aes_string(fill="conc2",ymin="act_lo",ymax="act_hi"),alpha=0.3)
}
curveplot1 <- curveplot1 +
geom_line(aes_string(colour="conc2"))+
labs(x=surflabels$x,y=surflabels$fill,color=surflabels$y,
title=paste0("Potentiation of ",compounds[[1]]))+
plotscales$x1+plotscales$y+plotscales$color
if (!is.null(analysis$braidFit$ciMat)) {
curveplot1 <- curveplot1 +
labs(fill=surflabels$y)+
plotscales$fill
}
curveplot1 <- curveplot1 +
theme(text=element_text(size=base_size),
legend.key.size=grid::unit(0.8,"lines"))
curvedf2 <- data.frame(
conc1 = rep(xbreaks,each=length(yvector)),
conc2 = rep(yvector,times=length(xbreaks))
)
curvedf2$act <- evalBraidModel(curvedf2$conc1,curvedf2$conc2,braidPar)
if (!is.null(analysis$braidFit$ciMat)) {
curveci2 <- calcBraidConfInt(
analysis$braidFit,
function(p) evalBraidModel(curvedf2$conc1,curvedf2$conc2,bpar=p)
)
curvedf2$act_lo <- curveci2[,1]
curvedf2$act_hi <- curveci2[,3]
}
curvedf2$conc1 <- factor(curvedf2$conc1)
curveplot2 <- ggplot(curvedf2,aes_string(x="conc2",y="act"))+
geom_hline(yintercept=levels,linetype=2,colour="black")
if (!is.null(analysis$braidFit$ciMat)) {
curveplot2 <- curveplot2 +
geom_ribbon(aes_string(fill="conc1",ymin="act_lo",ymax="act_hi"),alpha=0.3)
}
curveplot2 <- curveplot2 +
geom_line(aes_string(colour="conc1"))+
labs(x=surflabels$y,y=surflabels$fill,color=surflabels$x,
title=paste0("Potentiation of ",compounds[[2]]))+
plotscales$x2+plotscales$y+plotscales$color
if (!is.null(analysis$braidFit$ciMat)) {
curveplot2 <- curveplot2 +
labs(fill=surflabels$x)+
plotscales$fill
}
curveplot2 <- curveplot2 +
theme(text=element_text(size=base_size),
legend.key.size=grid::unit(0.8,"lines"))
if (!is.null(control$curvetheme) && inherits(control$curvetheme,"theme")) {
curveplot1 <- curveplot1 + control$curvetheme
curveplot2 <- curveplot2 + control$curvetheme
}
}
if (is.null(control$title)) {
control$title = paste0(compounds[[1]]," vs. ",compounds[[2]])
}
titleGrob <- grid::textGrob(control$title,gp=grid::gpar(fontsize=14,fontface="bold"))
if (control$layout=="standard") {
if (control$orientation=="portrait") {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots[c(1,2,5,6)],ncol=2))
g2 <- rbind(partab,spacetab1,metrictab,size="max")
# g3 <- rbind(iaetab1,spacetab2,iaetab2,size="max")
# g4 <- cowplot::plot_grid(curveplot1,curveplot2,nrow=1)
g3 <- cowplot::plot_grid(iaetab1,iaetab2,curveplot1,curveplot2,ncol=2,rel_heights=c(2,3))
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(g1,g2,rel_widths=c(2,1)),
g3,
# cowplot::plot_grid(g3,g4,rel_widths=c(2,5)),
ncol=1,rel_heights=c(1,11,11)
)
} else {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots[c(1,2,5,6)],ncol=2))
g2 <- cowplot::plot_grid(partab,metrictab)
g3 <- cowplot::plot_grid(iaetab1,iaetab2)
g4 <- cowplot::plot_grid(curveplot1,curveplot2,ncol=1)
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(
cowplot::plot_grid(g1,g2,ncol=1,rel_heights=c(2,1)),
cowplot::plot_grid(g3,g4,ncol=1,rel_heights=c(2,5)),
nrow=1,rel_widths=c(5,3)
),
ncol=1,rel_heights=c(1,16)
)
}
} else if (control$layout=="dense") {
if (control$orientation=="portrait") {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots,ncol=2))
g2 <- rbind(partab,spacetab1,metrictab,size="max")
g3 <- rbind(iaetab1,spacetab2,iaetab2,size="max")
g4 <- cowplot::plot_grid(curveplot1,curveplot2,nrow=1)
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(g1,g2,rel_widths=c(2,1)),
cowplot::plot_grid(g3,g4,rel_widths=c(2,5)),
ncol=1,rel_heights=c(1,16,6)
)
} else {
for (i in 1:6) {
plots[[i]] <- plots[[i]]+theme(legend.position="bottom")
}
curveplot1 <- curveplot1+theme(legend.position="bottom")
curveplot2 <- curveplot2+theme(legend.position="bottom")
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots,ncol=2))
g2 <- rbind(partab,spacetab2,metrictab,spacetab2,
iaetab1,spacetab2,iaetab2,size="max")
g4 <- cowplot::plot_grid(curveplot1,curveplot2,ncol=1)
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(g1,g2,g4,nrow=1,rel_widths=c(6,3,4)),
ncol=1,rel_heights=c(1,16)
)
}
} else if (control$layout=="simple") {
if (control$orientation=="portrait") {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=c(plots[c(1,2,5,6)],list(partab),list(metrictab)),ncol=2))
cowplot::plot_grid(
titleGrob,
g1,
ncol=1,rel_heights=c(1,20)
)
} else {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots[c(1,2,5,6)],ncol=2))
g2 <- rbind(partab,spacetab1,metrictab,size="max")
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(g1,g2,rel_widths=c(5,2)),
ncol=1,rel_heights=c(1,16)
)
}
}
}
recastPositionScale <- function(oldScale,dimension="x") {
if (utils::packageVersion("ggplot2")>="3.5.0") {
if (dimension=="x") {
newScale <- scale_x_continuous(
name=oldScale$name,
breaks=oldScale$breaks,
minor_breaks=oldScale$minor_breaks,
n.breaks=oldScale$n.breaks,
labels=oldScale$labels,
limits=oldScale$limits,
expand=oldScale$expand,
oob=oldScale$oob,
na.value=oldScale$na.value,
transform=oldScale$trans
)
} else {
newScale <- scale_y_continuous(
name=oldScale$name,
breaks=oldScale$breaks,
minor_breaks=oldScale$minor_breaks,
n.breaks=oldScale$n.breaks,
labels=oldScale$labels,
limits=oldScale$limits,
expand=oldScale$expand,
oob=oldScale$oob,
na.value=oldScale$na.value,
transform=oldScale$trans
)
}
} else {
if (dimension=="x") {
newScale <- scale_x_continuous(
name=oldScale$name,
breaks=oldScale$breaks,
minor_breaks=oldScale$minor_breaks,
n.breaks=oldScale$n.breaks,
labels=oldScale$labels,
limits=oldScale$limits,
expand=oldScale$expand,
oob=oldScale$oob,
na.value=oldScale$na.value,
trans=oldScale$trans
)
} else {
newScale <- scale_y_continuous(
name=oldScale$name,
breaks=oldScale$breaks,
minor_breaks=oldScale$minor_breaks,
n.breaks=oldScale$n.breaks,
labels=oldScale$labels,
limits=oldScale$limits,
expand=oldScale$expand,
oob=oldScale$oob,
na.value=oldScale$na.value,
trans=oldScale$trans
)
}
}
}
recastFillScale <- function(colorscale) {
discrete_scale(
"fill",
palette=colorscale$palette,
name=colorscale$name,
breaks=colorscale$breaks,
labels=colorscale$labels,
limits=colorscale$limits,
expand=colorscale$expand,
na.translate=colorscale$na.translate,
na.value=colorscale$na.value,
drop=colorscale$drop,
guide=colorscale$guide,
position=colorscale$position
)
}
#' BRAID Surface Analysis
#'
#' Performs a convenient pre-built set of BRAID and dose-response analysis
#' tasks
#'
#' @inheritParams braidrm::findBestBraid
#' @param ... Additional parameters to be passed to [braidrm::findBestBraid()]
#'
#' @return An object of class `braidAnalysis`, containing the following values:
#' * `concs`: a width-two array containing the two tested doses for each
#' measurement
#' * `act`: a numeric vector with as many values as `concs` has rows,
#' containing the measured values for each measurement
#' * `weights`: a numeric vector of weights, the same length as `act`,
#' specifying the weight given to each measurement in fitting. All weights are
#' 1 by default
#' * `braidFit`: a fit object of class `braidrm` containing the best-fit BRAID
#' surface according to the given constraints
#' * `hillFit1`: If the given data contains measurements of the first drug in
#' isolation, those measurements are fit using [basicdrm::findBestHillModel];
#' the results of this analysis are stored as an object of class `hillrm` as
#' `hillFit1`. If no such measurements are found, this will be `NULL`
#' * `hillFit2`: the corresponding fit for measurements of the second drug
#' alone, if they are included; `NULL` otherwise
#'
#' @export
#' @examples
#' surface <- synergisticExample
#'
#' analysis <- runBraidAnalysis(measure~concA+concB, surface, defaults=c(0,1))
#'
#' names(analysis)
runBraidAnalysis <- function(formula,data,defaults,weights=NULL,start=NULL,
direction=0,lower=NULL,upper=NULL,useBIC=TRUE,...) UseMethod("runBraidAnalysis")
#' @export
#' @rdname runBraidAnalysis
runBraidAnalysis.formula <- function(formula,data,defaults,weights=NULL,start=NULL,
direction=0,lower=NULL,upper=NULL,useBIC=TRUE,...) {
mf <- stats::model.frame(formula=formula, data=data)
concs <- stats::model.matrix(attr(mf, "terms"), data=mf)
tms <- attr(concs,"assign")
for (i in seq(length(tms),1,by=-1)) {
if (tms[i]==0) { concs <- concs[,-i] }
}
act <- stats::model.response(mf)
weights <- eval(substitute(weights),data)
bfit <- runBraidAnalysis.default(concs,act,defaults,weights,start,
direction,lower,upper,useBIC,...)
bfit$call <- match.call()
return(bfit)
}
#' @export
#' @rdname runBraidAnalysis
runBraidAnalysis.default <- function(formula,data,defaults,weights=NULL,start=NULL,
direction=0,lower=NULL,upper=NULL,useBIC=TRUE,...) {
concs <- formula
act <- data
if (is.null(weights)) { weights <- rep(1,length(act)) }
else if (length(weights)==1) { weights <- rep(weights,length(act)) }
if (any(is.finite(log(concs[,1])) & concs[,2]==0)) {
conc_r1 <- concs[concs[,2]==0,1]
act_r1 <- act[concs[,2]==0]
weights_r1 <- weights[concs[,2]==0]
if (is.null(start)) { start_r1 <- NULL }
else { start_r1 <- c(start[c(1,3)],defaults) }
if (is.null(lower)) { lower_r1 <- NULL }
else { upper_r1 <- lower[c(1,3,6,7)] }
if (is.null(lower)) { upper_r1 <- NULL }
else { upper_r1 <- upper[c(1,3,6,7)] }
hfit1 <- basicdrm::findBestHillModel(conc_r1,act_r1,defaults,weights_r1,start_r1,
direction,lower_r1,upper_r1,useBIC)
} else {
hfit1 <- NULL
}
if (any(is.finite(log(concs[,2])) & concs[,1]==0)) {
conc_r2 <- concs[concs[,1]==0,2]
act_r2 <- act[concs[,1]==0]
weights_r2 <- weights[concs[,1]==0]
if (is.null(start)) { start_r2 <- NULL }
else { start_r2 <- c(start[c(2,4)],defaults) }
if (is.null(lower)) { lower_r2 <- NULL }
else { upper_r2 <- lower[c(2,4,6,8)] }
if (is.null(lower)) { upper_r2 <- NULL }
else { upper_r2 <- upper[c(2,4,6,8)] }
hfit2 <- basicdrm::findBestHillModel(conc_r2,act_r2,defaults,weights_r2,start_r2,
direction,lower_r2,upper_r2,useBIC)
} else {
hfit2 <- NULL
}
bfit <- findBestBraid(concs,act,defaults,weights=weights,start=start,
direction=direction,lower=lower,upper=upper,useBIC=useBIC,...)
structure(
list(concs=concs,act=act,weights=weights,
braidFit=bfit,hillFit1=hfit1,hillFit2=hfit2),
class="braidAnalysis"
)
}
#' Basic BRAID Analysis Conversion
#'
#' @param bfit A BRAID fit object of class `braidrm`
#'
#' @details
#' While we strongly recommend using the [runBraidAnalysis()] function for a
#' more complete treatment of a combination; there may be circumstances in
#' which is necessary or preferable to use an existing BRAID fit object (of
#' class `braidrm`). Thsi function takes such a fit and wraps it in a minimal
#' `braidAnalysis` object which can then be passed to [makeBraidReport()]
#'
#' @return A BRAID analysis object of class `braidanalysis` containin only the
#' results from the given BRAID fit
#' @export
#' @examples
#' surface <- antagonisticExample
#' fit <- braidrm(measure ~ concA + concB, surface, model="kappa2")
#'
#' analysis <- basicBraidAnalysis(fit)
#'
#' names(analysis)
basicBraidAnalysis <- function(bfit) {
structure(
list(concs=bfit$concs,act=bfit$act,weights=bfit$weights,
braidFit=bfit,hillFit1=NULL,hillFit2=NULL),
class="braidAnalysis"
)
}
Try the braidReports package in your browser
Any scripts or data that you put into this service are public.
braidReports documentation built on Sept. 30, 2024, 1:06 a.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.
Defines functions basicBraidAnalysis runBraidAnalysis.default runBraidAnalysis.formula runBraidAnalysis recastFillScale recastPositionScale makeBraidReport
Documented in basicBraidAnalysis makeBraidReport runBraidAnalysis runBraidAnalysis.default runBraidAnalysis.formula
#' Render a BRAID Report
#'
#' Produces a one page report depicting the results of a full BRAID analysis
#' for a single combination.
#'
#' @param analysis An object of class `braidAnalysis` produced by the
#' [runBraidAnalysis()] or [basicBraidAnalysis()] functions
#' @param compounds A length-2 character vector containing the names of the
#' two compounds tested in the combination
#' @param levels Two levels at which the IAE should be evaluated
#' @param limits Two values representing the maximal achievable concentrations
#' for the compounds tested, used to esitmate the IAE
#' @param control A named list of additional control parameters adjusting the
#' appearance of the resulting report
#'
#' @details
#' This function attempts, however foolhardily, to encompass many of the
#' details, plots, and values that the user might wish to report for a complete
#' BRAID analysis of a given drug combination. All reports are built for a
#' single 8.5-by-11 inch page, either in landscape or potrait orientation, but
#' reports can be customized to contain more or less information. Here is a
#' full list of what *can* appear in the BRAID report:
#'
#' * A raw and smoothed plot of the actual measured response data; the raw plot
#' is built using [geom_braid()] or, for irregularly laid out data,
#' [geom_braid_glass()], while the smoothed data is built using
#' [geom_braid_smooth()]. (Included in all layouts)
#' * A plot of residual errors and a smoothed surface plot of the predicted
#' additive surface based on the dose response behavior of the individual
#' compounds alone. In cases of pronounced interaction, this will differ
#' significantly from the best-fit BRAID plots. (Included in the dense layout)
#' * A plot of residual errors and a smoothed surface plot of the best-fitting
#' BRAID surface. (Included in the all layouts)
#' * A table of the best-fitting BRAID response surface parameters (Included in
#' all layouts)
#' * A table of estimated IAE values at the specified effect levels (Included in
#' all layouts)
#' * Two tables of the dose required of one drug to produce a desired effect
#' level (the first value in `levels`) in the presence of several doses of the
#' other drug; used to gauge the degree of potentiation. (Included in the
#' standard and dense layouts)
#' * Two plots depicting the predicted dose response of one drug in the presence
#' of various levels of the other, also used to gauge potentiation. (Included
#' in the standard and dense layouts)
#'
#' So the resulting report page can contain from six (simple layout) to twelve
#' (dense layout) elements depicting different aspects of the BRAID analysis.
#'
#' The precise appearance of the report page is controlled by various elements
#' of the `control` parameter. Though the default value of the parameter is
#' an empty list, several fields will be filled in if they are unspecified. The
#' full set of possible control options is:
#'
#' * `abbs`: A pair of character strings specifying the abbreviations of the
#' tested compounds. By default, the abbreviations consist of the firs three
#' characters of each compound's name, but for some compound names this is not
#' an appropriate abbreviation Abbreviations are used in many axis labels and
#' tables
#' * `units`: If included, a single string or pair of strings specifying the
#' dose units for the two drugs tested, included in axis labels and tables. If
#' left unspecified, units will not be included
#' * `xscale`: Either a character string specifying a named transformation
#' object (e.g "log2", "log10", "sqrt") or a `ggplot2` continuous x-position
#' scale generated by [ggplot2::scale_x_continuous()] or related functions. This
#' scale will be applied to the x-dimension of all surface plots and the
#' x-dimension of the first potentiation plot. If a `name` is specified for the
#' scale, this will be the x-axis label; otherwise other defaults will be used.
#' Default value for this control parameter is "log10".
#' * `yscale`: Either a character string specifying a named transformation
#' object (e.g "log2", "log10", "sqrt") or a `ggplot2` continuous y-position
#' scale generated by [ggplot2::scale_y_continuous()] or related functions. This
#' scale will be applied to the y-dimension of all surface plots and the
#' x-dimension of the second potentiation plot. If a `name` is specified for the
#' scale, this will be the y-axis label; otherwise other defaults will be used.
#' Default value for this control parameter is "log10".
#' * `fillscale`: If included, continuous fill scale object generated by one of
#' several `ggplot2` continuous fill scale functions. This fill scale will
#' control the fill appearance for all *effect* surface plots; fill colors
#' in residual error plots will use a different color palette. In addition,
#' any names, labels, breaks, transformations, etc. included in this scale will
#' also be applied to the y-axis of both potentiation plots, as those also
#' represent the modeled effect. If unspecified, will be set to a brewer
#' continuous color scale with palette "RdYlBu".
#' * `colorscale`: If included, a discrete color scale object generated by one
#' of several `ggplot2` discrete color scale functions. This color scale
#' controls the colors chosen for the curves in the two potentiation plots. If
#' left unspecified, will default to the output of
#' [ggplot2::scale_color_discrete()]
#' * `xname`: A string specifying the desired x-axis labels in surface plots.
#' Will be overridden if control parameter `xscale` is a scale object with a
#' non-empty `name` attribute. If unspecified, defaults the abbreviation of the
#' first compound followed by the units if included.
#' * `yname`: A string specifying the desired y-axis labels in surface plots.
#' Will be overridden if control parameter `yscale` is a scale object with a
#' non-empty `name` attribute. If unspecified, defaults the abbreviation of the
#' second compound followed by the units if included.
#' * `effectname`: The name of the modeled effect variable. Could be "Growth"
#' or "Survival" or "Activity". The default value is simply "Effect"
#' * `title`: A string containing the overall title of the report page. If left
#' unspecified, will simply be the first compound "vs." the second
#' * `contourcolor`: Contours of the smoothed surfaces at the levels specified
#' by the parameter `levels` are included in all smoothed plots. By default,
#' they are black, but specifying this control parameter will change that color
#' * `irregular`: If `TRUE`, the data are not assumed to lie on a regular grid
#' in the plotted, and will be visualized with [geom_braid_glass()] rather than
#' [geom_braid()]
#' * `swidth`: A numeric value specifying how widely the smooth surfaces should
#' be smoothed, passed as the `width` aesthetic to [geom_braid_smooth()]
#' * `sheight`: A numeric value specifying how far in the y-direction the smooth
#' surfaces should be smoothed, passed as the `height` aesthetic to
#' [geom_braid_smooth()]
#' * `npoints`: The density of points used in the smoothed plots. See
#' [geom_braid_smooth()] for details
#' * `leveltext`: A pair of strings indicating how the two IAE levels should be
#' displayed in tables. In some cases, the precise number at which the IAE is
#' calculated does not reflect the level that the user wishes to express. So
#' a user might want to refer to a relative survival value of 0.1 as IAE90 (for
#' 90% killing); or a log2-fold growth inhibition as IAE50 (for 50% inhibition);
#' passing "50" or "90" as the `leveltext` in such cases will produce the
#' desired appearance. If left unspecified, labels will simply be the string
#' representations of the paramter `levels`
#' * `metrics`: A named numeric or named character vector specifying additional
#' metrics for the combination; they will be added to the table containing the
#' calculated IAE values. Numeric values will be rounded to three significant
#' figures; string values will be included as is. Names will be parsed by
#' default, so the string "CI\[90\]" will be displayed with the "90" as a
#' subscript
#' * `surftheme`: A `ggplot2` `theme` object specifying any additional theme
#' adjustments to add to all response surface plots
#' * `curvetheme`: A `ggplot2` `theme` object specifying any additional theme
#' adjustments to add to both potentiation curve plots
#' * `layout`: The specific layout to be used; determines which report elements
#' are included. Can be "simple", "standard" (the default), or "dense"
#' * `orientation`: The expected orientantion of the rendered page; can be
#' "portrait" (the default) or "landscape"
#'
#'
#' @return A graphical object containing all plots and tables, arranged
#' according to the desired format. The resulting object is optimized for a
#' single page, either portrait or landscape as specified in `control`
#' @export
#'
#' @examples
#' surface <- synergisticExample
#' analysis <- runBraidAnalysis(measure~concA+concB, surface,
#' defaults=c(0,1), getCIs=FALSE)
#'
#' report <- makeBraidReport(analysis,c("A Drug","B Drug"),
#' levels=c(0.5, 0.9),limits=c(5,5))
#' print(report)
#'
#' control <- list(abbs=c("A","B"),units=c("\u00B5M"),leveltext=c("50","90"),
#' xscale=scale_x_log10(breaks=c(0.1,0.5,2,10),
#' labels=as.character),
#' fillscale=scale_fill_viridis_c(option="A"),
#' colorscale=scale_color_brewer(palette="Set1"),
#' title="Example Analysis")
#' nextReport <- makeBraidReport(analysis,c("A Drug","B Drug"),
#' levels=c(0.5, 0.9),limits=c(5,5),
#' control=control)
#' print(nextReport)
makeBraidReport <- function(analysis,compounds,levels,limits,control=list()) {
newerggplot <- utils::packageVersion("ggplot2")>="3.5.0"
if (!inherits(analysis,"braidAnalysis")) {
stop("Parameter 'analysis' must be an object of class 'braidAnalysis'.")
}
# Control options:
#
# xscale, yscale, fillscale, colorscale
# xname, yname, effectname, title
# contourcolor
# abbs, units
# irregular, swidth, sheight, npoints
# leveltext, metrics
# layout, orientation
defaultControls <- list(
xscale = "log10",
yscale = "log10",
contourcolor = "black",
npoints = 50,
layout = "standard",
orientation = "portrait"
)
defaultControls[names(control)] <- control
control <- defaultControls
if (!is.null(control$base_size)) {
base_size <- control$base_size
} else if (control$layout=="standard") {
base_size <- 8
} else if (control$layout=="dense") {
base_size <- 7
} else if (control$layout=="simple") {
base_size <- 9
} else { stop(paste0("Unrecognized layout options '",control$layout,"'.")) }
braidPar <- stats::coef(analysis$braidFit)
braidDir <- ifelse (braidPar[["Ef"]]>=braidPar[["E0"]],1,-1)
# build additive parameter vector
additivePar <- braidPar
if (is.null(analysis$hillFit1) || is.null(analysis$hillFit2)) {
additiveDir <- braidDir
} else {
hfit1 <- analysis$hillFit1
hfit2 <- analysis$hillFit2
hpar1 <- stats::coef(analysis$hillFit1)
hpar2 <- stats::coef(analysis$hillFit2)
hDir1 <- ifelse(hpar1[["Ef"]]>=hpar1[["E0"]],1,-1)
hDir2 <- ifelse(hpar2[["Ef"]]>=hpar2[["E0"]],1,-1)
if (stats::var(hfit1$fitted.values)==0 | stats::var(hfit1$act)==0) { cor1 <- 0 }
else { cor1 <- stats::cor(hfit1$act,hfit1$fitted.values) }
if (stats::var(hfit2$fitted.values)==0 | stats::var(hfit2$act)==0) { cor2 <- 0 }
else { cor2 <- stats::cor(hfit2$act,hfit2$fitted.values) }
if (cor1>=cor2) {
additiveDir <- hDir1
additivePar[["E0"]] <- hpar1[["E0"]]
} else {
additiveDir <- hDir2
additivePar[["E0"]] <- hpar2[["E0"]]
}
}
additivePar[["kappa"]] <- 0
if (!is.null(analysis$hillFit1)) {
hpar1 <- stats::coef(analysis$hillFit1)
if (sign(additiveDir*(hpar1[["Ef"]]-additivePar[["E0"]]))>=0) {
additivePar[c("IDMA","na","EfA")] <- hpar1[c("IDM","n","Ef")]
}
}
if (!is.null(analysis$hillFit2)) {
hpar2 <- stats::coef(analysis$hillFit2)
if (sign(additiveDir*(hpar2[["Ef"]]-additivePar[["E0"]]))>=0) {
additivePar[c("IDMB","nb","EfB")] <- hpar2[c("IDM","n","Ef")]
}
}
additivePar[["Ef"]] <- additiveDir*max(additiveDir*additivePar[c("EfA","EfB")])
mainDF <- data.frame(conc1=analysis$concs[,1],conc2=analysis$concs[,2],act=analysis$act)
mainDF$bfit <- evalBraidModel(mainDF$conc1,mainDF$conc2,braidPar)
mainDF$afit <- evalBraidModel(mainDF$conc1,mainDF$conc2,additivePar)
mainDF$berror <- mainDF$act-mainDF$bfit
mainDF$aerror <- mainDF$act-mainDF$afit
r2add <- stats::cov.wt(cbind(mainDF$act,mainDF$afit),wt=analysis$weights,cor=TRUE)$cor[1,2]
r2brd <- stats::cov.wt(cbind(mainDF$act,mainDF$bfit),wt=analysis$weights,cor=TRUE)$cor[1,2]
r2add <- signif(r2add^2,3)
r2brd <- signif(r2brd^2,3)
surflabels <- list()
surfscales <- list()
if (!is.null(control$units) && length(control$units)==1) {
control$units <- c(control$units,control$units)
}
if (is.null(control$xname)) {
if (is.null(control$abbs)) { surflabels$x <- substring(compounds[[1]],1,3) }
else ( surflabels$x <- control$abbs[[1]])
if (!is.null(control$units)) {
surflabels$x <- paste0(surflabels$x," (",control$units[[1]],")")
}
} else { surflabels$x <- control$xname }
if (is.null(control$yname)) {
if (is.null(control$abbs)) { surflabels$y <- substring(compounds[[2]],1,3) }
else ( surflabels$y <- control$abbs[[2]])
if (!is.null(control$units)) {
surflabels$y <- paste0(surflabels$y," (",control$units[[2]],")")
}
} else { surflabels$y <- control$yname }
if (is.null(control$effectname)) { surflabels$fill <- "Effect" }
else { surflabels$fill <- control$effectname }
if (is.character(control$xscale)) {
if (newerggplot) {
surfscales$x <- scale_x_continuous(transform=control$xscale)
} else {
surfscales$x <- scale_x_continuous(trans=control$xscale)
}
} else if (inherits(control$xscale,"ScaleContinuousPosition")) {
surfscales$x <- control$xscale
} else {
stop("Control parameter \"xscale\" must be a transform string or a ggplot2 position scale,")
}
if (is.character(control$yscale)) {
if (newerggplot) {
surfscales$y <- scale_y_continuous(transform=control$yscale)
} else {
surfscales$y <- scale_y_continuous(trans=control$yscale)
}
} else if (inherits(control$yscale,"ScaleContinuousPosition")) {
surfscales$y <- control$yscale
} else {
stop("Control parameter \"yscale\" must be a transform string or a ggplot2 position scale,")
}
if (is.null(control$fillscale)) {
surfscales$fill <- scale_fill_distiller(palette="RdYlBu",direction=-braidDir)
} else if (inherits(control$fillscale,"ScaleContinuous")) {
surfscales$fill <- control$fillscale
} else {
stop("Control parameter \"fillscale\" must be ggplot2 fill scale,")
}
valrange <- range(c(mainDF$act,mainDF$afit,mainDF$bfit))
errrange <- range(c(mainDF$aerror,mainDF$berror))
blankdf <- data.frame(conc1=stats::median(mainDF$conc1,na.rm=TRUE),
conc2=stats::median(mainDF$conc2,na.rm=TRUE),
act=valrange,bfit=valrange,afit=valrange,
berror=errrange,aerror=errrange)
if (is.null(control$irregular)) { irregular <- FALSE }
else { irregular <- control$irregular }
if (is.null(control$swidth)) { swidth <- NULL }
else { swidth <- control$swidth }
if (is.null(control$sheight)) { sheight <- NULL }
else { sheight <- control$sheight }
if (is.null(sheight)) {
if (is.null(swidth)) {
smooth_aes <- aes()
} else {
smooth_aes <- aes(width=swidth)
}
} else {
if (is.null(swidth)) {
smooth_aes <- aes(height=sheight)
} else {
smooth_aes <- aes(width=swidth,height=sheight)
}
}
crel <- is.finite(log(mainDF$conc1)) & is.finite(log(mainDF$conc2))
if (length(limits)==1) { limits <- c(limits,limits) }
limits <- limits[1:2]
limitdf <- data.frame(conc1=c(),conc2=c())
if (any(crel)) {
crng1 <- range(mainDF$conc1[crel])
crng2 <- range(mainDF$conc2[crel])
if (!any(limits<c(crng1[[1]],crng2[[1]])) &&
any(limits<=c(crng1[[2]],crng2[[2]]))) {
checkvec <- limits<=c(crng1[[2]],crng2[[2]])
if (all(checkvec)) {
limitdf <- data.frame(conc1=c(0,limits[[1]],limits[[1]]),
conc2=c(limits[[2]],limits[[2]],0))
} else if (checkvec[[1]]) {
limitdf <- data.frame(conc1=c(limits[[1]],limits[[1]]),
conc2=c(0,Inf))
} else {
limitdf <- data.frame(conc1=c(0,Inf),
conc2=c(limits[[2]],limits[[2]]))
}
}
}
plots <- list()
plots[[1]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="act"))+geom_blank(data=blankdf)
plots[[2]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="act",z="act"))+geom_blank(data=blankdf)
plots[[3]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="aerror"))+ geom_blank(data=blankdf)
plots[[4]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="afit",z="afit"))+geom_blank(data=blankdf)
plots[[5]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="berror"))+geom_blank(data=blankdf)
plots[[6]] <- ggplot(mainDF,aes_string(x="conc1",y="conc2",fill="bfit",z="bfit"))+geom_blank(data=blankdf)
if (irregular) {
plots[[1]] <- plots[[1]]+geom_braid_glass()
plots[[3]] <- plots[[3]]+geom_braid_glass()
plots[[5]] <- plots[[5]]+geom_braid_glass()
} else {
plots[[1]] <- plots[[1]]+geom_braid()
plots[[3]] <- plots[[3]]+geom_braid()
plots[[5]] <- plots[[5]]+geom_braid()
}
plots[[2]] <- plots[[2]]+
geom_braid_smooth(smooth_aes,npoints=control$npoints)+
geom_braid_contour(smooth_aes,npoints=control$npoints,breaks=levels,colour=control$contourcolor)+
annotate("contour",x=limitdf$conc1,y=limitdf$conc2,colour=control$contourcolor,linetype=2)
plots[[4]] <- plots[[4]]+
geom_braid_smooth(smooth_aes,npoints=control$npoints)+
geom_braid_contour(smooth_aes,npoints=control$npoints,breaks=levels,colour=control$contourcolor)+
annotate("contour",x=limitdf$conc1,y=limitdf$conc2,colour=control$contourcolor,linetype=2)
plots[[6]] <- plots[[6]]+
geom_braid_smooth(smooth_aes,npoints=control$npoints)+
geom_braid_contour(smooth_aes,npoints=control$npoints,breaks=levels,colour=control$contourcolor)+
annotate("contour",x=limitdf$conc1,y=limitdf$conc2,colour=control$contourcolor,linetype=2)
plots[[1]] <- plots[[1]]+
labs(x=surflabels$x,y=surflabels$y,fill=surflabels$fill,title="Measured Data")+
surfscales$x+surfscales$y+surfscales$fill
plots[[2]] <- plots[[2]]+
labs(x=surflabels$x,y=surflabels$y,fill=surflabels$fill,title="Smoothed Data")+
surfscales$x+surfscales$y+surfscales$fill
plots[[3]] <- plots[[3]]+
labs(x=surflabels$x,y=surflabels$y,title="Additive Error")+
surfscales$x+surfscales$y+scale_fill_gradient2("Error")
plots[[4]] <- plots[[4]]+
labs(x=surflabels$x,y=surflabels$y,fill=surflabels$fill,
title=paste0("Additive Fit (R\u00B2=",r2add,")"))+
surfscales$x+surfscales$y+surfscales$fill
plots[[5]] <- plots[[5]]+
labs(x=surflabels$x,y=surflabels$y,title="BRAID Error")+
surfscales$x+surfscales$y+scale_fill_gradient2("Error")
plots[[6]] <- plots[[6]]+
labs(x=surflabels$x,y=surflabels$y,fill=surflabels$fill,
title=paste0("BRAID Fit (R\u00B2=",r2brd,")"))+
surfscales$x+surfscales$y+surfscales$fill
for (i in 1:6) {
plots[[i]] <- plots[[i]]+
theme(text=element_text(size=base_size),
legend.key.size=grid::unit(0.8,"lines"))
if (!is.null(control$surftheme) && inherits(control$surftheme,"theme")) {
plots[[i]] <- plots[[i]]+control$surftheme
}
}
paramdf <- data.frame(name=c("ID['M,A']","ID['M,B']","n[a]","n[b]","kappa",
"E[0]","E['f,A']","E['f,B']","E[f]"),
value=signif(stats::coef(analysis$braidFit),3))
bmodel <- analysis$braidFit$model
if (!is.null(analysis$braidFit$ciMat)) {
paramdf$lo <- NA_real_
paramdf$hi <- NA_real_
paramdf$lo[bmodel] <- signif(pmin(analysis$braidFit$ciMat[,1],paramdf$value[bmodel]),3)
paramdf$hi[bmodel] <- signif(pmax(analysis$braidFit$ciMat[,2],paramdf$value[bmodel]),3)
paramdf$text <- ifelse(is.na(paramdf$lo),
as.character(paramdf$value),
paste0(as.character(paramdf$value)," (",
as.character(paramdf$lo)," \u2013 ",
as.character(paramdf$hi),")"))
} else {
paramdf$text <- as.character(paramdf$value)
}
if (!any((7:8)%in%bmodel)) {
paramdf <- paramdf[c(1:6,9),]
} else if (!(9%in%bmodel)) {
paramdf <- paramdf[1:8,]
}
paramvec <- paramdf$text
names(paramvec) <- paramdf$name
partab <- parameterTable(paramvec,parse=TRUE,title="Best-Fit BRAID",base_size=base_size)
if (is.null(control$leveltext)) {
leveltext <- as.character(levels)
} else {
leveltext <- control$leveltext
}
metricdf <- data.frame(name=paste0("IAE['",leveltext,"']"),
value=signif(estimateIAE(braidPar,levels,limits),3))
if (!is.null(analysis$braidFit$ciMat)) {
iaemat <- calcBraidConfInt(
analysis$braidFit,
function(par) estimateIAE(par,levels,limits)
)
metricdf$lo <- signif(iaemat[,1],3)
metricdf$hi <- signif(iaemat[,3],3)
} else {
metricdf$lo <- NA_real_
metricdf$hi <- NA_real_
}
metricdf$text <- ifelse(is.na(metricdf$lo),
as.character(metricdf$value),
paste0(as.character(metricdf$value)," (",
as.character(metricdf$lo)," \u2013 ",
as.character(metricdf$hi),")"))
if (!is.null(control$metrics)) {
if (is.null(names(control$metrics))) {
warning("Control value 'metrics' must be a named vector. Values ignored.")
} else if (is.numeric(control$metrics)) {
metricdf <- rbind(metricdf,data.frame(name=names(control$metrics),
value=as.numeric(control$metrics),
lo=NA_real_,hi=NA_real_,
text=as.character(control$metrics)))
} else if (is.character(control$metrics)) {
metricdf <- rbind(metricdf,data.frame(name=names(control$metrics),
value=NA_real_,
lo=NA_real_,hi=NA_real_,
text=control$metrics))
} else {
warning("Control value 'metrics' must either be numeric or a character vector. Values ignored.")
}
}
metricvec <- metricdf$text
names(metricvec) <- metricdf$name
metrictab <- parameterTable(metricvec,parse=TRUE,title="Summary Metrics",base_size=base_size)
spacetab1 <- gtable::gtable(widths=unit(10,"mm"),heights=unit(10,"mm"),name="spacer")
spacetab2 <- gtable::gtable(widths=unit(10,"mm"),heights=unit(5,"mm"),name="spacer")
if (control$layout != "simple") {
pls <- suppressWarnings(layer_scales(plots[[1]]))
plsx <- pls$x
xbreaks <- plsx$break_info()$major_source
xbreaks <- plsx$trans$inverse(xbreaks)
xvector <- plsx$range$range
xvector <- seq(xvector[[1]],xvector[[2]],length=control$npoints)
xvector <- plsx$trans$inverse(xvector)
xvector[xvector<0] <- 0
xvector <- unique(xvector)
plsy <- pls$y
ybreaks <- plsy$break_info()$major_source
ybreaks <- plsy$trans$inverse(ybreaks)
yvector <- plsy$range$range
yvector <- seq(yvector[[1]],yvector[[2]],length=control$npoints)
yvector <- plsy$trans$inverse(yvector)
yvector[yvector<0] <- 0
yvector <- unique(yvector)
xbreaks <- sort(unique(c(0,xbreaks[
xbreaks>=min(analysis$concs[,1],na.rm=TRUE) &
xbreaks<=max(analysis$concs[,1],na.rm=TRUE)
])))
if (length(xbreaks)==6) { xbreaks <- xbreaks[c(1,3,4,5)] }
else if (length(xbreaks)>6) { xbreaks <- xbreaks[seq(1,length(xbreaks),by=2)] }
ybreaks <- sort(unique(c(0,ybreaks[
ybreaks>=min(analysis$concs[,2],na.rm=TRUE) &
ybreaks<=max(analysis$concs[,2],na.rm=TRUE)
])))
if (length(ybreaks)==6) { ybreaks <- ybreaks[c(1,3,4,5)] }
else if (length(ybreaks)>6) { ybreaks <- ybreaks[seq(1,length(ybreaks),by=2)] }
iaedf1 <- data.frame(name=ybreaks,
value=signif(invertBraidModel(DB=ybreaks,
effect=levels[[1]],
bpar=braidPar),3))
if (!is.null(analysis$braidFit$ciMat)) {
iaeci1 <- calcBraidConfInt(
analysis$braidFit,
function(p) invertBraidModel(DB=ybreaks,effect=levels[[1]],bpar=p)
)
iaedf1$lo <- signif(iaeci1[,1],3)
iaedf1$hi <- signif(iaeci1[,3],3)
} else {
iaedf1$lo <- NA_real_
iaedf1$hi <- NA_real_
}
iaedf1$text <- ifelse(is.na(iaedf1$lo),
as.character(iaedf1$value),
paste0(as.character(iaedf1$value)," (",
as.character(iaedf1$lo)," \u2013 ",
as.character(iaedf1$hi),")"))
iaedf1 <- iaedf1[,c("name","text")]
iaetitle1 <- paste0("IC['",leveltext[[1]],"']~of~'",compounds[[1]],"'")
iaetab1 <- dataFrameTable(iaedf1,title=iaetitle1,names=c(surflabels$y,"Estimate"),parse="title",base_size=base_size)
iaedf2 <- data.frame(name=xbreaks,
value=signif(invertBraidModel(DA=xbreaks,
effect=levels[[1]],
bpar=braidPar),3))
if (!is.null(analysis$braidFit$ciMat)) {
iaeci2 <- calcBraidConfInt(
analysis$braidFit,
function(p) invertBraidModel(DA=xbreaks,effect=levels[[1]],bpar=p)
)
iaedf2$lo <- signif(iaeci2[,1],3)
iaedf2$hi <- signif(iaeci2[,3],3)
} else {
iaedf2$lo <- NA_real_
iaedf2$hi <- NA_real_
}
iaedf2$text <- ifelse(is.na(iaedf2$lo),
as.character(iaedf2$value),
paste0(as.character(iaedf2$value)," (",
as.character(iaedf2$lo)," \u2013 ",
as.character(iaedf2$hi),")"))
iaedf2 <- iaedf2[,c("name","text")]
iaetitle2 <- paste0("IC['",leveltext[[1]],"']~of~'",compounds[[2]],"'")
iaetab2 <- dataFrameTable(iaedf2,title=iaetitle2,names=c(surflabels$x,"Estimate"),parse="title",base_size=base_size)
plotscales <- list()
plotscales$x1 <- surfscales$x
plotscales$x2 <- recastPositionScale(surfscales$y,"x")
plotscales$y <- recastPositionScale(surfscales$fill,"y")
if (is.null(control$colorscale)) {
plotscales$color <- scale_color_discrete()
} else {
plotscales$color <- control$colorscale
}
plotscales$fill <- recastFillScale(plotscales$color)
curvedf1 <- data.frame(
conc1 = rep(xvector,times=length(ybreaks)),
conc2 = rep(ybreaks,each=length(xvector))
)
curvedf1$act <- evalBraidModel(curvedf1$conc1,curvedf1$conc2,braidPar)
if (!is.null(analysis$braidFit$ciMat)) {
curveci1 <- calcBraidConfInt(
analysis$braidFit,
function(p) evalBraidModel(curvedf1$conc1,curvedf1$conc2,bpar=p)
)
curvedf1$act_lo <- curveci1[,1]
curvedf1$act_hi <- curveci1[,3]
}
curvedf1$conc2 <- factor(curvedf1$conc2)
curveplot1 <- ggplot(curvedf1,aes_string(x="conc1",y="act"))+
geom_hline(yintercept=levels,linetype=2,colour="black")
if (!is.null(analysis$braidFit$ciMat)) {
curveplot1 <- curveplot1 +
geom_ribbon(aes_string(fill="conc2",ymin="act_lo",ymax="act_hi"),alpha=0.3)
}
curveplot1 <- curveplot1 +
geom_line(aes_string(colour="conc2"))+
labs(x=surflabels$x,y=surflabels$fill,color=surflabels$y,
title=paste0("Potentiation of ",compounds[[1]]))+
plotscales$x1+plotscales$y+plotscales$color
if (!is.null(analysis$braidFit$ciMat)) {
curveplot1 <- curveplot1 +
labs(fill=surflabels$y)+
plotscales$fill
}
curveplot1 <- curveplot1 +
theme(text=element_text(size=base_size),
legend.key.size=grid::unit(0.8,"lines"))
curvedf2 <- data.frame(
conc1 = rep(xbreaks,each=length(yvector)),
conc2 = rep(yvector,times=length(xbreaks))
)
curvedf2$act <- evalBraidModel(curvedf2$conc1,curvedf2$conc2,braidPar)
if (!is.null(analysis$braidFit$ciMat)) {
curveci2 <- calcBraidConfInt(
analysis$braidFit,
function(p) evalBraidModel(curvedf2$conc1,curvedf2$conc2,bpar=p)
)
curvedf2$act_lo <- curveci2[,1]
curvedf2$act_hi <- curveci2[,3]
}
curvedf2$conc1 <- factor(curvedf2$conc1)
curveplot2 <- ggplot(curvedf2,aes_string(x="conc2",y="act"))+
geom_hline(yintercept=levels,linetype=2,colour="black")
if (!is.null(analysis$braidFit$ciMat)) {
curveplot2 <- curveplot2 +
geom_ribbon(aes_string(fill="conc1",ymin="act_lo",ymax="act_hi"),alpha=0.3)
}
curveplot2 <- curveplot2 +
geom_line(aes_string(colour="conc1"))+
labs(x=surflabels$y,y=surflabels$fill,color=surflabels$x,
title=paste0("Potentiation of ",compounds[[2]]))+
plotscales$x2+plotscales$y+plotscales$color
if (!is.null(analysis$braidFit$ciMat)) {
curveplot2 <- curveplot2 +
labs(fill=surflabels$x)+
plotscales$fill
}
curveplot2 <- curveplot2 +
theme(text=element_text(size=base_size),
legend.key.size=grid::unit(0.8,"lines"))
if (!is.null(control$curvetheme) && inherits(control$curvetheme,"theme")) {
curveplot1 <- curveplot1 + control$curvetheme
curveplot2 <- curveplot2 + control$curvetheme
}
}
if (is.null(control$title)) {
control$title = paste0(compounds[[1]]," vs. ",compounds[[2]])
}
titleGrob <- grid::textGrob(control$title,gp=grid::gpar(fontsize=14,fontface="bold"))
if (control$layout=="standard") {
if (control$orientation=="portrait") {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots[c(1,2,5,6)],ncol=2))
g2 <- rbind(partab,spacetab1,metrictab,size="max")
# g3 <- rbind(iaetab1,spacetab2,iaetab2,size="max")
# g4 <- cowplot::plot_grid(curveplot1,curveplot2,nrow=1)
g3 <- cowplot::plot_grid(iaetab1,iaetab2,curveplot1,curveplot2,ncol=2,rel_heights=c(2,3))
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(g1,g2,rel_widths=c(2,1)),
g3,
# cowplot::plot_grid(g3,g4,rel_widths=c(2,5)),
ncol=1,rel_heights=c(1,11,11)
)
} else {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots[c(1,2,5,6)],ncol=2))
g2 <- cowplot::plot_grid(partab,metrictab)
g3 <- cowplot::plot_grid(iaetab1,iaetab2)
g4 <- cowplot::plot_grid(curveplot1,curveplot2,ncol=1)
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(
cowplot::plot_grid(g1,g2,ncol=1,rel_heights=c(2,1)),
cowplot::plot_grid(g3,g4,ncol=1,rel_heights=c(2,5)),
nrow=1,rel_widths=c(5,3)
),
ncol=1,rel_heights=c(1,16)
)
}
} else if (control$layout=="dense") {
if (control$orientation=="portrait") {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots,ncol=2))
g2 <- rbind(partab,spacetab1,metrictab,size="max")
g3 <- rbind(iaetab1,spacetab2,iaetab2,size="max")
g4 <- cowplot::plot_grid(curveplot1,curveplot2,nrow=1)
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(g1,g2,rel_widths=c(2,1)),
cowplot::plot_grid(g3,g4,rel_widths=c(2,5)),
ncol=1,rel_heights=c(1,16,6)
)
} else {
for (i in 1:6) {
plots[[i]] <- plots[[i]]+theme(legend.position="bottom")
}
curveplot1 <- curveplot1+theme(legend.position="bottom")
curveplot2 <- curveplot2+theme(legend.position="bottom")
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots,ncol=2))
g2 <- rbind(partab,spacetab2,metrictab,spacetab2,
iaetab1,spacetab2,iaetab2,size="max")
g4 <- cowplot::plot_grid(curveplot1,curveplot2,ncol=1)
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(g1,g2,g4,nrow=1,rel_widths=c(6,3,4)),
ncol=1,rel_heights=c(1,16)
)
}
} else if (control$layout=="simple") {
if (control$orientation=="portrait") {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=c(plots[c(1,2,5,6)],list(partab),list(metrictab)),ncol=2))
cowplot::plot_grid(
titleGrob,
g1,
ncol=1,rel_heights=c(1,20)
)
} else {
g1 <- suppressWarnings(cowplot::plot_grid(plotlist=plots[c(1,2,5,6)],ncol=2))
g2 <- rbind(partab,spacetab1,metrictab,size="max")
cowplot::plot_grid(
titleGrob,
cowplot::plot_grid(g1,g2,rel_widths=c(5,2)),
ncol=1,rel_heights=c(1,16)
)
}
}
}
recastPositionScale <- function(oldScale,dimension="x") {
if (utils::packageVersion("ggplot2")>="3.5.0") {
if (dimension=="x") {
newScale <- scale_x_continuous(
name=oldScale$name,
breaks=oldScale$breaks,
minor_breaks=oldScale$minor_breaks,
n.breaks=oldScale$n.breaks,
labels=oldScale$labels,
limits=oldScale$limits,
expand=oldScale$expand,
oob=oldScale$oob,
na.value=oldScale$na.value,
transform=oldScale$trans
)
} else {
newScale <- scale_y_continuous(
name=oldScale$name,
breaks=oldScale$breaks,
minor_breaks=oldScale$minor_breaks,
n.breaks=oldScale$n.breaks,
labels=oldScale$labels,
limits=oldScale$limits,
expand=oldScale$expand,
oob=oldScale$oob,
na.value=oldScale$na.value,
transform=oldScale$trans
)
}
} else {
if (dimension=="x") {
newScale <- scale_x_continuous(
name=oldScale$name,
breaks=oldScale$breaks,
minor_breaks=oldScale$minor_breaks,
n.breaks=oldScale$n.breaks,
labels=oldScale$labels,
limits=oldScale$limits,
expand=oldScale$expand,
oob=oldScale$oob,
na.value=oldScale$na.value,
trans=oldScale$trans
)
} else {
newScale <- scale_y_continuous(
name=oldScale$name,
breaks=oldScale$breaks,
minor_breaks=oldScale$minor_breaks,
n.breaks=oldScale$n.breaks,
labels=oldScale$labels,
limits=oldScale$limits,
expand=oldScale$expand,
oob=oldScale$oob,
na.value=oldScale$na.value,
trans=oldScale$trans
)
}
}
}
recastFillScale <- function(colorscale) {
discrete_scale(
"fill",
palette=colorscale$palette,
name=colorscale$name,
breaks=colorscale$breaks,
labels=colorscale$labels,
limits=colorscale$limits,
expand=colorscale$expand,
na.translate=colorscale$na.translate,
na.value=colorscale$na.value,
drop=colorscale$drop,
guide=colorscale$guide,
position=colorscale$position
)
}
#' BRAID Surface Analysis
#'
#' Performs a convenient pre-built set of BRAID and dose-response analysis
#' tasks
#'
#' @inheritParams braidrm::findBestBraid
#' @param ... Additional parameters to be passed to [braidrm::findBestBraid()]
#'
#' @return An object of class `braidAnalysis`, containing the following values:
#' * `concs`: a width-two array containing the two tested doses for each
#' measurement
#' * `act`: a numeric vector with as many values as `concs` has rows,
#' containing the measured values for each measurement
#' * `weights`: a numeric vector of weights, the same length as `act`,
#' specifying the weight given to each measurement in fitting. All weights are
#' 1 by default
#' * `braidFit`: a fit object of class `braidrm` containing the best-fit BRAID
#' surface according to the given constraints
#' * `hillFit1`: If the given data contains measurements of the first drug in
#' isolation, those measurements are fit using [basicdrm::findBestHillModel];
#' the results of this analysis are stored as an object of class `hillrm` as
#' `hillFit1`. If no such measurements are found, this will be `NULL`
#' * `hillFit2`: the corresponding fit for measurements of the second drug
#' alone, if they are included; `NULL` otherwise
#'
#' @export
#' @examples
#' surface <- synergisticExample
#'
#' analysis <- runBraidAnalysis(measure~concA+concB, surface, defaults=c(0,1))
#'
#' names(analysis)
runBraidAnalysis <- function(formula,data,defaults,weights=NULL,start=NULL,
direction=0,lower=NULL,upper=NULL,useBIC=TRUE,...) UseMethod("runBraidAnalysis")
#' @export
#' @rdname runBraidAnalysis
runBraidAnalysis.formula <- function(formula,data,defaults,weights=NULL,start=NULL,
direction=0,lower=NULL,upper=NULL,useBIC=TRUE,...) {
mf <- stats::model.frame(formula=formula, data=data)
concs <- stats::model.matrix(attr(mf, "terms"), data=mf)
tms <- attr(concs,"assign")
for (i in seq(length(tms),1,by=-1)) {
if (tms[i]==0) { concs <- concs[,-i] }
}
act <- stats::model.response(mf)
weights <- eval(substitute(weights),data)
bfit <- runBraidAnalysis.default(concs,act,defaults,weights,start,
direction,lower,upper,useBIC,...)
bfit$call <- match.call()
return(bfit)
}
#' @export
#' @rdname runBraidAnalysis
runBraidAnalysis.default <- function(formula,data,defaults,weights=NULL,start=NULL,
direction=0,lower=NULL,upper=NULL,useBIC=TRUE,...) {
concs <- formula
act <- data
if (is.null(weights)) { weights <- rep(1,length(act)) }
else if (length(weights)==1) { weights <- rep(weights,length(act)) }
if (any(is.finite(log(concs[,1])) & concs[,2]==0)) {
conc_r1 <- concs[concs[,2]==0,1]
act_r1 <- act[concs[,2]==0]
weights_r1 <- weights[concs[,2]==0]
if (is.null(start)) { start_r1 <- NULL }
else { start_r1 <- c(start[c(1,3)],defaults) }
if (is.null(lower)) { lower_r1 <- NULL }
else { upper_r1 <- lower[c(1,3,6,7)] }
if (is.null(lower)) { upper_r1 <- NULL }
else { upper_r1 <- upper[c(1,3,6,7)] }
hfit1 <- basicdrm::findBestHillModel(conc_r1,act_r1,defaults,weights_r1,start_r1,
direction,lower_r1,upper_r1,useBIC)
} else {
hfit1 <- NULL
}
if (any(is.finite(log(concs[,2])) & concs[,1]==0)) {
conc_r2 <- concs[concs[,1]==0,2]
act_r2 <- act[concs[,1]==0]
weights_r2 <- weights[concs[,1]==0]
if (is.null(start)) { start_r2 <- NULL }
else { start_r2 <- c(start[c(2,4)],defaults) }
if (is.null(lower)) { lower_r2 <- NULL }
else { upper_r2 <- lower[c(2,4,6,8)] }
if (is.null(lower)) { upper_r2 <- NULL }
else { upper_r2 <- upper[c(2,4,6,8)] }
hfit2 <- basicdrm::findBestHillModel(conc_r2,act_r2,defaults,weights_r2,start_r2,
direction,lower_r2,upper_r2,useBIC)
} else {
hfit2 <- NULL
}
bfit <- findBestBraid(concs,act,defaults,weights=weights,start=start,
direction=direction,lower=lower,upper=upper,useBIC=useBIC,...)
structure(
list(concs=concs,act=act,weights=weights,
braidFit=bfit,hillFit1=hfit1,hillFit2=hfit2),
class="braidAnalysis"
)
}
#' Basic BRAID Analysis Conversion
#'
#' @param bfit A BRAID fit object of class `braidrm`
#'
#' @details
#' While we strongly recommend using the [runBraidAnalysis()] function for a
#' more complete treatment of a combination; there may be circumstances in
#' which is necessary or preferable to use an existing BRAID fit object (of
#' class `braidrm`). Thsi function takes such a fit and wraps it in a minimal
#' `braidAnalysis` object which can then be passed to [makeBraidReport()]
#'
#' @return A BRAID analysis object of class `braidanalysis` containin only the
#' results from the given BRAID fit
#' @export
#' @examples
#' surface <- antagonisticExample
#' fit <- braidrm(measure ~ concA + concB, surface, model="kappa2")
#'
#' analysis <- basicBraidAnalysis(fit)
#'
#' names(analysis)
basicBraidAnalysis <- function(bfit) {
structure(
list(concs=bfit$concs,act=bfit$act,weights=bfit$weights,
braidFit=bfit,hillFit1=NULL,hillFit2=NULL),
class="braidAnalysis"
)
}
Try the braidReports package in your browser
Any scripts or data that you put into this service are public.
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.