Optimization algorithms

We advise the user to employ the method = "nlslm" argument which is set as the default in monotherapy curve estimation. It is based on minpack.lm::nlsLM function with an underlying Levenberg-Marquardt algorithm for non-linear least squares estimation. This algorithm is known to be more robust than method = "nls" and its Gauss-Newton algorithm. In cases with nice sigmoid-shaped data, both methods should however lead to similar results.

method = "optim" is a simple sum-of-squared-residuals minimization driven by a default Nelder-Mead algorithm from optim minimizer. It is typically slower than non-linear least squares based estimation and can lead to a significant increase in computational time for larger datasets and bootstrapped statistics. In nice cases, Nelder-Mead algorithm and non-linear least squares can lead to rather similar estimates but this is not always the case as these algorithms are based on different techniques.

In general, we advise that in automated batch processing whenever method = "nlslm" does not converge fast enough and/or emits a warning, user should implement a fallback to method = "optim" and re-do the estimation. If none of these suggestions work, it might be useful to fiddle around and slightly perturb starting values for the algorithms as well. By default, these are obtained from the initialMarginal function.

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")
  })

Note as well that additional arguments to fitMarginals passed via ... ellipsis argument will be passed on to the respective solver function, i.e. minpack.lm::nlsLM, nls or optim.

Custom marginal fit

While BIGL package provides several routines to fit 4-parameter log-logistic dose-response models, some users may prefer to use their own optimizers to estimate the relevant parameters. It is rather easy to integrate this into the workflow by constructing a custom MarginalFit object. It is in practice a simple list with

• coef: named vector with coefficient estimates
• sigma: standard deviation of residuals
• df: degrees of freedom from monotherapy curve estimates
• model: model of the marginal estimation which allows imposing linear constraints on parameters. If no constraints are necessary, it can be left out or assigned the output of constructFormula function with no inputs.
• shared_asymptote: whether estimation is constrained to share the asymptote. During the estimation, this is deduced from model object.
• method: method used in dose-response curve estimation which will be re-used in bootstrapping
• transforms: power and biological transformation functions (and their inverses) used in monotherapy curve estimation. This should be a list in a format described above. If transforms is unspecified or NULL, no transformations will be used in statistical bootstrapping unless the user asks for it explicitly via one of the arguments to fitSurface.

Other elements in the MarginalFit are currently unused for evaluating synergy and can be disregarded. These elements, however, might be necessary to ensure proper working of available methods for the MarginalFit object.

As an example, the following code generates a custom MarginalFit object that can be passed further to estimate a response surface under the null hypothesis.

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")

Note that during bootstrapping this would use minpack.lm::nlsLM function to re-estimate parameters from data following the null. A custom optimizer for bootstrapping is currently not implemented.

(Generalized) Loewe model

If transformation functions were estimated using fitMarginals, these will be automatically recycled from the marginalFit object when doing calculations for the response surface fit. Alternatively, transformation functions can be passed by a separate argument. Since the marginalFit object was estimated without the shared asymptote constraint, the following will compute the response surface based on the generalized Loewe model.

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

Null model: Generalized Loewe Additivity
Variance assumption used: "equal"
Mean occupancy rate: 0.6732745

Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
Transformations: Yes

                 Drug A Drug B
Slope             2.132  1.201
Maximal response  0.088  0.091
log10(EC50)       1.008  0.977

Common baseline at: 0.132



No test statistics were computed.

The occcupancy matrix used in the expected response calculation for the Loewe models can be accessed with rs$occupancy.

For off-axis data and a fixed dose combination, the Z-score for that dose combination is defined to be the standardized difference between the observed effect and the effect predicted by a generalized Loewe model. If the observed effect differs significantly from the prediction, it might be due to the presence of synergy or antagonism. If multiple observations refer to the same combination of doses, then a mean is taken over these multiple standardized differences.

The following plot illustrates the isobologram of the chosen null model. Coloring and contour lines within the plot should help the user distinguish areas and dose combinations that generate similar response according to the null model. Note that the isobologram is plotted by default on a logarithmically scaled grid of doses.

isobologram(rs)

The plot below illustrates the above considerations in a 3-dimensional setting. In this plot, points refer to the observed effects whereas the surface is the model-predicted response. The surface is colored according to the median Z-scores where blue coloring indicates possible synergistic effects (red coloring would indicate possible antagonism).

plot(rs, legend = FALSE, main = "")
view3d(0, -75)
rglwidget()

Highest Single Agent

For the Highest Single Agent null model to work properly, it is expected that both marginal curves are either decreasing or increasing. Equivalent summary and plot methods are also available for this type of null model.

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

Null model: Highest Single Agent with differing maximal response
Variance assumption used: "equal"

Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
Transformations: Yes

                 Drug A Drug B
Slope             2.132  1.201
Maximal response  0.088  0.091
log10(EC50)       1.008  0.977

Common baseline at: 0.132

Exact meanR test (H0 = no synergy/antagonism):
    F(49,117) = 16.432 (p-value < 2e-16)

Evidence for effects in data: Syn
Points with significant deviations from the null: 
     d1     d2          R p-value call
3   2.0   0.12  -4.161608 0.00281  Syn
10  2.0   0.49  -4.095304 0.00358  Syn
17  2.0   2.00  -4.375919 0.00113  Syn
18  7.8   2.00  -8.192123  <2e-16  Syn
24  2.0   7.80  -4.197198 0.00238  Syn
25  7.8   7.80 -22.097201  <2e-16  Syn
32  7.8  31.00  -7.656737  <2e-16  Syn
40 31.0 130.00  -3.440117 0.03841  Syn
47 31.0 500.00  -3.668181 0.01791  Syn

Overall maxR summary:
 Call Syn Ant Total
  Syn   9   0    49

Bliss Independence

Also for the Bliss independence null model to work properly, it is expected that both marginal curves are either decreasing or increasing. Equivalent summary and plot methods are also available for this type of null model.

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

Null model: Bliss independence with differing maximal response
Variance assumption used: "equal"

Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
Transformations: Yes

                 Drug A Drug B
Slope             2.132  1.201
Maximal response  0.088  0.091
log10(EC50)       1.008  0.977

Common baseline at: 0.132

Exact meanR test (H0 = no synergy/antagonism):
    F(49,117) = 5.3573 (p-value = 0)

Evidence for effects in data: Syn
Points with significant deviations from the null: 
    d1    d2          R p-value call
3  2.0  0.12  -3.501361 0.03219  Syn
18 7.8  2.00  -3.514212 0.03067  Syn
24 2.0  7.80  -5.029746   1e-04  Syn
25 7.8  7.80 -10.130484  <2e-16  Syn
31 2.0 31.00  -3.522727 0.02967  Syn
32 7.8 31.00  -3.656988 0.01887  Syn

Overall maxR summary:
 Call Syn Ant Total
  Syn   6   0    49

Alternative Loewe Generalization

Also for the Alternative Loewe Generalization null model to work properly, it is expected that both marginal curves are either decreasing or increasing. Equivalent summary and plot methods are also available for this type of null model.

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

Null model: Alternative generalization of Loewe Additivity
Variance assumption used: "equal"

Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
Transformations: Yes

                 Drug A Drug B
Slope             2.132  1.201
Maximal response  0.088  0.091
log10(EC50)       1.008  0.977

Common baseline at: 0.132

Exact meanR test (H0 = no synergy/antagonism):
    F(49,117) = 4.6588 (p-value = 0)

Evidence for effects in data: Syn
Points with significant deviations from the null: 
    d1    d2         R p-value call
3  2.0  0.12 -3.890538 0.00861  Syn
25 7.8  7.80 -9.296311  <2e-16  Syn
32 7.8 31.00 -4.444400 0.00096  Syn

Overall maxR summary:
 Call Syn Ant Total
  Syn   3   0    49

meanR

Here we will use the previously computed CP covariance matrix to speed up the process.

• normal errors

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

• non-normal errors

The previous piece of code assumes normal errors. If we drop this assumption, we can use bootstrap methods to resample from the observed errors. Other parameters for bootstrapping, such as additional distribution for errors, wild bootstrapping to account for heteroskedasticity, are also available. See help(fitSurface).

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

Both tests use the same calculated F-statistic but compare it to different null distributions. In this particular case, both tests lead to identical results.

maxR

The meanR statistic can be complemented by the maxR statistic for each of available dose combinations. We will do this once again by assuming both normal and non-normal errors similar to the computation of the meanR statistic.

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)

Here is the summary of maxR statistics. It lists the total number of dose combinations listed as synergistic or antagonistic for Experiment 1 given the above calculations.

By using the outsidePoints function, we can obtain a quick summary indicating which dose combinations in Experiment 1 appear to deviate significantly from the null model according to the maxR statistic.

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

Synergistic effects of drug combinations can be depicted in a bi-dimensional contour plot where the x-axis and y-axis represent doses of Compound 1 and Compound 2 respectively and each point is colored based on the p-value and sign of the respective maxR statistic.

contour(maxR_B,
         ## colorPalette = c("blue", "black", "black", "red"),
        main = paste0(" Experiment ", i, " contour plot for maxR"),
        scientific = TRUE, digits = 3, cutoff = cutoff)

Previously, we had colored the 3-dimensional predicted response surface plot based on its Z-score, i.e. deviation of the predicted versus the observed effect. We can also easily color it based on the computed maxR statistic to account for additional statistical variation.

plot(maxR_B, color = "maxR", legend = FALSE, main = "")
view3d(0, -75)
rglwidget()