analysis.R
In BIGL: Biochemically Intuitive Generalized Loewe Model

## ----init, message = FALSE----------------------------------------------------
library(BIGL)
library(knitr)
library(ggplot2)
set.seed(12345)

if (!requireNamespace("rmarkdown", quietly = TRUE) || !rmarkdown::pandoc_available("1.14")) {
  warning(call. = FALSE, "These vignettes assume rmarkdown and pandoc
          version 1.14. These were not found. Older versions will not work.")
  knitr::knit_exit()
}

## ----settings-----------------------------------------------------------------
nExp <- 4             # Dataset has 11 experiments, we consider only 4
cutoff <- 0.95        # Cutoff for p-values to use in plot.maxR() function

## ----data---------------------------------------------------------------------
data("directAntivirals", package = "BIGL")
head(directAntivirals)

## -----------------------------------------------------------------------------
subsetData <- function(data, i) {
  ## Subset data to a single experiment and, optionally, select the necessary
  ## columns only
  subset(data, experiment == i)[, c("effect", "d1", "d2")]
}

## ----subset, out.width="100%"-------------------------------------------------
i <- 1
data <- subsetData(directAntivirals, i)

## ----transformations----------------------------------------------------------
## Define forward and reverse transform functions
transforms <- list(
  "BiolT" = function(y, args) with(args, N0*exp(y*time.hours)),
  "InvBiolT" = function(T, args) with(args, 1/time.hours*log(T/N0)),
  "PowerT" = function(y, args) with(args, log(y)),
  "InvPowerT" = function(T, args) with(args, exp(T)),
  "compositeArgs" = list(N0 = 1,
                         time.hours = 72)
)

## ----autotransform, eval=FALSE------------------------------------------------
#  transforms_auto <- getTransformations(data)
#  fitMarginals(data, transforms = transforms_auto)
#  
#  ## In the case of 1-parameter Box-Cox transformation, it is easy
#  ## to retrieve the power parameter by evaluating the function at 0.
#  ## If parameter is 0, then it is a log-transformation.
#  with(transforms_auto, -1 / PowerT(0, compositeArgs))

## ----marginalFit--------------------------------------------------------------
## Fitting marginal models
marginalFit <- fitMarginals(data, transforms = transforms, method = "nls", 
    names = c("Drug A", "Drug B"))
summary(marginalFit)

## ----marginalPlot, fig.align="center", fig.height = 4, fig.width = 6----------
## Plotting marginal models
plot(marginalFit) + ggtitle(paste("Direct-acting antivirals - Experiment" , i))

## ----marginalFitC, eval = FALSE-----------------------------------------------
#  ## Parameter ordering: h1, h2, b, m1, m2, e1, e2
#  ## Constraint 1: m1 = m2. Constraint 2: b = 0.1
#  constraints <- list("matrix" = rbind(c(0, 0, 0, -1, 1, 0, 0),
#                                       c(0, 0, 1, 0, 0, 0, 0)),
#                      "vector" = c(0, 0.1))
#  
#  ## Parameter estimates will now satisfy equality:
#  ##   constraints$matrix %*% pars == constraints$vector
#  fitMarginals(data, transforms = transforms,
#               constraints = constraints)

## ----marginalFitFixed, eval = FALSE-------------------------------------------
#  ## Set baseline at 0.1 and maximal responses at 0.
#  fitMarginals(data, transforms = transforms,
#               fixed = c("m1" = 0, "m2" = 0, "b" = 0.1))

## ----fallback, eval = FALSE---------------------------------------------------
#  nlslmFit <- tryCatch({
#    fitMarginals(data, transforms = transforms,
#                 method = "nlslm")
#  }, warning = function(w) w, error = function(e) e)
#  
#  if (inherits(nlslmFit, c("warning", "error")))
#    optimFit <- tryCatch({
#      fitMarginals(data, transforms = transforms,
#                   method = "optim")
#    })

## ----eval=FALSE---------------------------------------------------------------
#  marginalFit <- list("coef" = c("h1" = 1, "h2" = 2, "b" = 0,
#                                 "m1" = 1.2, "m2" = 1, "e1" = 0.5, "e2" = 0.5),
#                      "sigma" = 0.1,
#                      "df" = 123,
#                      "model" = constructFormula(),
#                      "shared_asymptote" = FALSE,
#                      "method" = "nlslm",
#                      "transforms" = transforms)
#  class(marginalFit) <- append(class(marginalFit), "MarginalFit")

## ----analysis, message=FALSE, comment = NA------------------------------------
rs <- fitSurface(data, marginalFit,
                 null_model = "loewe",
                 B.CP = 50, statistic = "none", parallel = FALSE)
summary(rs)

## ----image, warning=FALSE, comment = NA, fig.width = 6, fig.height = 4, fig.align = "center"----
isobologram(rs)

## ----plot3d, warning=FALSE, fig.align="center", fig.height=7, fig.width=7-----
plot(rs, legend = FALSE, main = "")

## ----analysis_hsa, message=FALSE, comment = NA--------------------------------
rsh <- fitSurface(data, marginalFit,
                  null_model = "hsa",
                  B.CP = 50, statistic = "both", parallel = FALSE)
summary(rsh)

## ----analysis_bliss, message=FALSE, comment = NA------------------------------
rsb <- fitSurface(data, marginalFit, 
                  null_model = "bliss",
                  B.CP = 50, statistic = "both", parallel = FALSE)
summary(rsb)

## ----analysis_loewe2, message=FALSE, comment = NA-----------------------------
rsl2 <- fitSurface(data, marginalFit, 
                  null_model = "loewe2",
                  B.CP = 50, statistic = "both", parallel = FALSE)
summary(rsl2)

## ----plot_2d_cross_section, message=FALSE, comment = NA, fig.width = 8, fig.height = 6----
nullModels <- c("loewe", "loewe2", "bliss", "hsa")
rs_list <- Map(fitSurface, null_model = nullModels, MoreArgs = list(
        data = data, fitResult = marginalFit, 
        B.CP = 50, statistic = "none", parallel = FALSE)
)

synergy_plot_bycomp(rs_list, ylab = "Response", plotBy = "Drug A", color = TRUE)

## ----meanrnorm, message = FALSE-----------------------------------------------
meanR_N <- fitSurface(data, marginalFit,
                      statistic = "meanR", CP = rs$CP, B.B = NULL,
                      parallel = FALSE)

## ----meanrnonnorm, message = FALSE--------------------------------------------
meanR_B <- fitSurface(data, marginalFit,
                      statistic = "meanR", CP = rs$CP, B.B = 20,
                      parallel = FALSE)

## ----meanresults, echo=FALSE--------------------------------------------------
MeanR_both <- rbind("Normal errors" = c(meanR_N$meanR$FStat, meanR_N$meanR$p.value),
                    "Bootstrapped errors" = c(meanR_B$meanR$FStat, meanR_B$meanR$p.value))
colnames(MeanR_both) <- c("F-statistic", "p-value")
kable(MeanR_both)

## ----maxboth, message = FALSE-------------------------------------------------
maxR_N <- fitSurface(data, marginalFit,
                     statistic = "maxR", CP = rs$CP, B.B = NULL,
                     parallel = FALSE)
maxR_B <- fitSurface(data, marginalFit,
                     statistic = "maxR", CP = rs$CP, B.B = 20,
                     parallel = FALSE)
maxR_both <- rbind(summary(maxR_N$maxR)$totals,
                   summary(maxR_B$maxR)$totals)

## ----printmax, echo = FALSE---------------------------------------------------
rownames(maxR_both) <- c("Normal errors", "Bootstrapped errors")
kable(maxR_both)

## ----maxoutside, results="asis"-----------------------------------------------
outPts <- outsidePoints(maxR_B$maxR$Ymean)
kable(outPts, caption = paste0("Non-additive points for Experiment ", i))

## ----maxcontour, fig.align="center", fig.width=6, fig.height=5----------------
contour(maxR_B,
       colorPalette = c("blue", "white", "red"),
        main = paste0(" Experiment ", i, " contour plot for maxR"),
        scientific = TRUE, digits = 3, cutoff = cutoff
)

## ----plot3dmax, warning=FALSE, fig.height=7, fig.width=7----------------------
plot(maxR_B, color = "maxR", legend = FALSE, main = "")

## ----summarySingleConfInt-----------------------------------------------------
summary(maxR_B$confInt)

## ----plotSingleConfInt, fig.height=5, fig.width=8-----------------------------
plotConfInt(maxR_B, color = "effect-size")

## ----contour_effectsize, warning=FALSE, fig.align="center", fig.width=6, fig.height=5, message=FALSE, comment = NA----
contour(
    maxR_B,
    colorPalette = c("Syn" = "blue", "None" = "white", "Ant" = "red"),
    main = paste0(" Experiment ", i, " contour plot for effect size"),
    color = "effect-size",
    scientific = TRUE, digits = 3, cutoff = cutoff
)

## ----plot3d_effectsize, warning=FALSE, fig.height=7, fig.width=7, message=FALSE, comment = NA----
plot(maxR_B, color = "effect-size", legend = FALSE, main = "", gradient = FALSE)

## ----heterogenanalysis, fig.width=6, fig.height=5-----------------------------
marginalFit <- fitMarginals(data, transforms = NULL)
summary(marginalFit)

resU <- fitSurface(data, marginalFit, method = "unequal", 
    statistic = "both", B.CP = 20, B.B = 20, parallel = FALSE)
summary(resU)

## ----modelVariancePlot, fig.width=6, fig.height=5-----------------------------
plotMeanVarFit(data)
plotMeanVarFit(data, log = "xy") #Clearer on the log-scale
plotMeanVarFit(data, trans = "log") #Thresholded at maximum observed variance

## ----modelVarianceSum, fig.width=6, fig.height=5------------------------------
resM <- fitSurface(data, marginalFit, method = "model", 
    statistic = "both", B.CP = 20, B.B = 20, parallel = FALSE)

## ----modelVarianceSumLogTransform, fig.width=6, fig.height=5, eval = FALSE----
#  resL <- fitSurface(data, marginalFit, method = "model", trans = "log",
#      statistic = "both", B.CP = 20, B.B = 20, parallel = FALSE)

## ----resM---------------------------------------------------------------------
summary(resM) 

## ----fullanalysis, message=FALSE----------------------------------------------
marginalFits <- list()
datasets <- list()
respSurfaces <- list()
maxR.summary <- list()
for (i in seq_len(nExp)) {
  ## Select experiment
  data <- subsetData(directAntivirals, i)
  ## Fit joint marginal model
  marginalFit <- fitMarginals(data, transforms = transforms,
                              method = "nlslm")
  ## Predict response surface based on generalized Loewe model
  respSurface <- fitSurface(data, marginalFit,
                            statistic = "maxR", B.CP = 20,
                            parallel = FALSE)

  datasets[[i]] <- data
  marginalFits[[i]] <- marginalFit
  respSurfaces[[i]] <- respSurface
  maxR.summary[[i]] <- summary(respSurface$maxR)$totals
}

## ----maxrfull, echo=FALSE-----------------------------------------------------
allMaxR <- do.call(rbind, maxR.summary)
rownames(allMaxR) <- paste("Experiment", 1:nrow(allMaxR))
kable(allMaxR, row.names = TRUE)

## ----tabs, echo = FALSE, results = "asis"-------------------------------------
i <- 4
genCaption <- function(k) paste("Non-additive points for Experiment", k)
outPts <- outsidePoints(respSurfaces[[i]]$maxR$Ymean)
print(kable(outPts, caption = genCaption(i)))

## ----fullcontour, echo=FALSE, fig.align = "center", fig.width = 6, fig.height = 5----
i <- 4
contour(respSurfaces[[i]],
        main = paste("Experiment", i),
        scientific = TRUE, digits = 3, cutoff = cutoff)