The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.
We first load the necessary packages and set some pre-defined values needed to replicate the analysis.
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.")
::knit_exit()
knitr }
<- 4 # Dataset has 11 experiments, we consider only 4
nExp <- 0.95 # Cutoff for p-values to use in plot.maxR() function cutoff
The data for the analysis must come in a data-frame with required columns d1
, d2
and effect
for doses of two compounds and observed cell counts respectively. The effect
column may represent also a type of normalized data and subsequent transformation functions should be adjusted.
We will use sample data included in the package - directAntivirals
.
data("directAntivirals", package = "BIGL")
head(directAntivirals)
## experiment cpd1 cpd2 effect d1 d2
## 1 1 cpd1_A cpd2_A 509.64 500 500
## 2 1 cpd1_A cpd2_A 589.67 500 500
## 3 1 cpd1_A cpd2_A 524.71 500 130
## 4 1 cpd1_A cpd2_A 575.45 500 130
## 5 1 cpd1_A cpd2_A 634.53 500 31
## 6 1 cpd1_A cpd2_A 702.35 500 31
This data consists of 11 experiments that can be processed separately. For initial illustration purposes we choose just one experiment and retain only the columns of interest. We define a simple function to do just that.
<- function(data, i) {
subsetData ## Subset data to a single experiment and, optionally, select the necessary
## columns only
subset(data, experiment == i)[, c("effect", "d1", "d2")]
}
Now let us only pick Experiment 4
to illustrate the functionality of the package.
<- 4
i <- subsetData(directAntivirals, i) data
Dose-response data for Experiment 4
will be used for a large share of the analysis presented here, therefore this subset is stored in a dataframe called data
. Later, we will run the analysis for some other experiments as well.
If raw data is measured in cell counts, data transformation might be of interest to improve accuracy and interpretation of the model. Of course, this will depend on the model specification. For example, if a generalized Loewe model is assumed on the growth rate of the cell count, the appropriate conversion should be made from the observed cell counts. The formula used would be \[y = N_0\exp\left(kt\right)\] where \(k\) is a growth rate, \(t\) is time (fixed) and \(y\) is the observed cell count. If such a transformation is specified, it is referred to as the biological transformation.
In certain cases, variance-stabilizing transformations (Box-Cox) can also be useful. We refer to these transformations as power transformations. In many cases, a simple logarithmic transformation can be sufficient but, if desired, a helper function optim.boxcox
is available to automate the selection of Box-Cox transformation parameters.
In addition to specifying biological and power transformations, users are also asked to specify their inverses. These are later used in the bootstrapping procedure and plotting methods.
As an example, we might define a transforms
list that will be passed to the fitting functions. It contains both biological growth rate and power transformations along with their inverses.
## Define forward and reverse transform functions
<- list(
transforms "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)
)
compositeArgs
contains the initial cell counts (N0
) and incubation time (time.hours
). In certain cases, the getTransformations
wrapper function can be employed to automatically obtain a prepared list with biological growth rate and power transformations based on results from optim.boxcox
. Its output will also contain the inverses of these transforms.
<- getTransformations(data)
transforms_auto 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))
Once dose-response dataframe is correctly set up, we may proceed onto synergy analysis. We will use transforms
as defined above with a logarithmic transformation. If not desired, transforms
can be set to NULL
and would be ignored.
Synergy analysis is quite modular and is divided into 3 parts:
The first step of the fitting procedure will consist in treating marginal data only, i.e. those observations within the experiment where one of the compounds is dosed at zero. For each compound the corresponding marginal doses are modelled using a 4-parameter logistic model.
The marginal models will be estimated together using non-linear least squares estimation procedure. Estimation of both marginal models needs to be simultaneous since it is assumed they share a common baseline that also needs to be estimated. The fitMarginals
function and other marginal estimation routines will automatically extract marginal data from the dose-response data frame.
Before proceeding onto the estimation, we get a rough guess of the parameters to use as starting values in optimization and then we fit the model. marginalFit
, returned by the fitMarginals
routine, is an object of class MarginalFit
which is essentially a list containing the main information about the marginal models, in particular the estimated coefficients.
The optional names
argument allows to specify the names of the compounds to be shown on the plots and in the summary. If not defined, the defaults (“Compound 1” and “Compound 2”) are used.
## Fitting marginal models
<- fitMarginals(data, transforms = transforms, method = "nls",
marginalFit names = c("Drug A", "Drug B"))
## Warning: The `transformations` argument of `fitMarginals()` will be deprecated as of
## BIGL newer versions.
## ℹ Please transform response before all analysis.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
summary(marginalFit)
## Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
## Transformations: Yes
##
## Drug A Drug B
## Slope 1.253 0.843
## Maximal response 0.076 0.082
## log10(EC50) 0.915 -1.041
##
## Common baseline at: 0.115
marginalFit
object retains the data that was supplied and the transformation functions used in the fitting procedure. It also has a plot
method which allows for a quick visualization of the fitting results.
## Plotting marginal models
plot(marginalFit) + ggtitle(paste("Direct-acting antivirals - Experiment" , i))
Note as well that the fitMarginals
function allows specifying linear constraints on parameters. This provides an easy way for the user to impose asymptote equality, specific baseline value and other linear constraints that might be useful. See help(constructFormula)
for more details.
## Parameter ordering: h1, h2, b, m1, m2, e1, e2
## Constraint 1: m1 = m2. Constraint 2: b = 0.1
<- list("matrix" = rbind(c(0, 0, 0, -1, 1, 0, 0),
constraints 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)
The fitMarginals
function allows an alternative user-friendly way to specify one or more fixed-value constraints using a named vector passed to the function via fixed
argument.
## Set baseline at 0.1 and maximal responses at 0.
fitMarginals(data, transforms = transforms,
fixed = c("m1" = 0, "m2" = 0, "b" = 0.1))
By default, no constraints are set, thus asymptotes are not shared and so a generalized Loewe model will be estimated.
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.
<- tryCatch({
nlslmFit fitMarginals(data, transforms = transforms,
method = "nlslm")
warning = function(w) w, error = function(e) e)
},
if (inherits(nlslmFit, c("warning", "error")))
<- tryCatch({
optimFit 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
.
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 estimatessigma
: standard deviation of residualsdf
: degrees of freedom from monotherapy curve estimatesmodel
: 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 bootstrappingtransforms
: 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.
<- list("coef" = c("h1" = 1, "h2" = 2, "b" = 0,
customMarginalFit "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(customMarginalFit) <- append(class(customMarginalFit), "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.
Five types of null models are available for calculating expected response surfaces.
shared_asymptote = FALSE
, in the marginal fitting procedure and null_model = "loewe"
in response calculation.shared_asymptote = TRUE
in the marginal fitting procedure and null_model = "loewe"
in response calculation.null_model = "hsa"
irrespective of the value of shared_asymptote
.null_model = "bliss"
. In the situations when maximal responses are constrained to be equal, the classical Bliss independence approach is used, when they are not equal, the Bliss independence calculation is performed on responses rescaled to the maximum range (i.e. absolute difference between baseline and maximal response).null_model = "loewe2"
. If the asymptotes are constrained to be equal, this reduces to the classical Loewe. Note that if shared_asymptote = TRUE
constraints are used, this also reduces to classical Loewe model.Three methods are available to control for errors
control = "FWER"
control = "FCR"
control = "dFCR"
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.
<- fitSurface(data, marginalFit,
rs null_model = "loewe",
B.CP = 50, statistic = "none", parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
summary(rs)
Null model: Generalized Loewe Additivity
Variance assumption used: "equal"
Mean occupancy rate: 0.5948603
Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
Transformations: Yes
Drug A Drug B
Slope 1.253 0.843
Maximal response 0.076 0.082
log10(EC50) 0.915 -1.041
Common baseline at: 0.115
No test statistics were computed.
CONFIDENCE INTERVALS
Overall effect
Estimated mean departure from null response surface with 95% confidence interval:
-0.2208 [-0.2594, -0.1757]
Evidence for effects in data: Syn
Significant pointwise effects
estimate lower upper call
0.12_0.004 -0.2404 -0.4094 -0.0713 Syn
0.12_1 -0.2675 -0.4724 -0.0625 Syn
0.49_0.016 -0.3341 -0.5100 -0.1581 Syn
0.49_0.063 -0.3613 -0.5420 -0.1805 Syn
0.49_0.25 -0.2129 -0.3929 -0.0328 Syn
0.49_1 -0.3840 -0.5892 -0.1788 Syn
130_0.25 -0.1904 -0.3637 -0.0171 Syn
2_0.004 -0.4308 -0.6034 -0.2581 Syn
2_0.016 -0.5848 -0.7535 -0.4162 Syn
2_0.063 -0.7147 -0.8899 -0.5394 Syn
2_0.25 -0.5367 -0.7123 -0.3611 Syn
2_1 -0.4389 -0.6452 -0.2326 Syn
31_0.016 -0.3208 -0.4970 -0.1447 Syn
31_0.063 -0.3865 -0.5621 -0.2109 Syn
31_0.25 -0.4510 -0.6350 -0.2671 Syn
31_1 -0.4359 -0.6663 -0.2054 Syn
7.8_0.004 -0.2392 -0.4212 -0.0571 Syn
7.8_0.016 -0.6429 -0.8180 -0.4677 Syn
7.8_0.063 -0.8496 -1.0195 -0.6797 Syn
7.8_0.25 -0.8124 -0.9824 -0.6424 Syn
7.8_1 -0.6211 -0.8335 -0.4087 Syn
Pointwise 95% confidence intervals summary:
Syn Ant Total
21 0 49
The occupancy 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 = "")
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.
<- fitSurface(data, marginalFit,
rsh null_model = "hsa",
B.CP = 50, statistic = "both", parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
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 1.253 0.843
Maximal response 0.076 0.082
log10(EC50) 0.915 -1.041
Common baseline at: 0.115
Exact meanR test (H0 = no synergy/antagonism):
F(49,117) = 24.8702 (p-value < 2e-16)
Evidence for effects in data: Syn
Points with significant deviations from the null:
d1 d2 absR p-value call
0.12_0.25 0.12 0.25000 3.984755 0.00571 Ant
0.49_0.016 0.49 0.01600 3.654198 0.01882 Syn
2_0.004 2.00 0.00400 7.030019 <2e-16 Syn
2_0.016 2.00 0.01600 12.372232 <2e-16 Syn
2_0.063 2.00 0.06300 11.845195 <2e-16 Syn
2_0.25 2.00 0.25000 6.820106 <2e-16 Syn
31_0.063 31.00 0.06300 4.048320 0.00433 Syn
31_0.25 31.00 0.25000 5.350282 1e-05 Syn
31_1 31.00 1.00000 5.224748 1e-05 Syn
7.8_0.00024 7.80 0.00024 6.175151 <2e-16 Ant
7.8_0.016 7.80 0.01600 8.091246 <2e-16 Syn
7.8_0.063 7.80 0.06300 14.996961 <2e-16 Syn
7.8_0.25 7.80 0.25000 15.296929 <2e-16 Syn
7.8_1 7.80 1.00000 7.654020 <2e-16 Syn
Overall maxR summary:
Call Syn Ant Total
Syn 12 2 49
CONFIDENCE INTERVALS
Overall effect
Estimated mean departure from null response surface with 95% confidence interval:
-0.2612 [-0.3024, -0.2024]
Evidence for effects in data: Syn
Significant pointwise effects
estimate lower upper call
0.12_0.004 -0.2593 -0.4216 -0.0969 Syn
0.12_1 -0.2693 -0.4610 -0.0775 Syn
0.49_0.00024 -0.1711 -0.3225 -0.0196 Syn
0.49_0.001 -0.1987 -0.3491 -0.0484 Syn
0.49_0.004 -0.2269 -0.3893 -0.0645 Syn
0.49_0.016 -0.4097 -0.5848 -0.2346 Syn
0.49_0.063 -0.4127 -0.5837 -0.2418 Syn
0.49_0.25 -0.2350 -0.4079 -0.0622 Syn
0.49_1 -0.3913 -0.5831 -0.1996 Syn
2_0.004 -0.5442 -0.7267 -0.3617 Syn
2_0.016 -0.8954 -1.0633 -0.7276 Syn
2_0.063 -0.9172 -1.0881 -0.7462 Syn
2_0.25 -0.6241 -0.7969 -0.4513 Syn
2_1 -0.4681 -0.6599 -0.2763 Syn
31_0.016 -0.3227 -0.5051 -0.1403 Syn
31_0.063 -0.3929 -0.5753 -0.2104 Syn
31_0.25 -0.4675 -0.6500 -0.2851 Syn
31_1 -0.4603 -0.6428 -0.2779 Syn
7.8_0.004 -0.2689 -0.4490 -0.0887 Syn
7.8_0.016 -0.7481 -0.9282 -0.5679 Syn
7.8_0.063 -1.1442 -1.3244 -0.9641 Syn
7.8_0.25 -1.1104 -1.2832 -0.9376 Syn
7.8_1 -0.7267 -0.9185 -0.5350 Syn
Pointwise 95% confidence intervals summary:
Syn Ant Total
23 0 49
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.
<- fitSurface(data, marginalFit,
rsb null_model = "bliss",
B.CP = 50, statistic = "both", parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
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 1.253 0.843
Maximal response 0.076 0.082
log10(EC50) 0.915 -1.041
Common baseline at: 0.115
Exact meanR test (H0 = no synergy/antagonism):
F(49,117) = 8.9099 (p-value < 2e-16)
Evidence for effects in data: Syn
Points with significant deviations from the null:
d1 d2 absR p-value call
0.49_1 0.49 1.00000 3.578273 0.02423 Syn
2_0.004 2.00 0.00400 3.568116 0.02524 Syn
2_0.016 2.00 0.01600 5.417975 2e-05 Syn
2_0.063 2.00 0.06300 7.533176 <2e-16 Syn
2_0.25 2.00 0.25000 5.089923 9e-05 Syn
31_0.016 31.00 0.01600 3.799524 0.01139 Syn
7.8_0.00024 7.80 0.00024 3.610220 0.02203 Ant
7.8_0.016 7.80 0.01600 5.931782 <2e-16 Syn
7.8_0.063 7.80 0.06300 8.103392 <2e-16 Syn
7.8_0.25 7.80 0.25000 7.244128 <2e-16 Syn
7.8_1 7.80 1.00000 4.009338 0.00526 Syn
Overall maxR summary:
Call Syn Ant Total
Syn 10 1 49
CONFIDENCE INTERVALS
Overall effect
Estimated mean departure from null response surface with 95% confidence interval:
-0.1734 [-0.2166, -0.1107]
Evidence for effects in data: Syn
Significant pointwise effects
estimate lower upper call
0.12_0.004 -0.2460 -0.3966 -0.0954 Syn
0.12_1 -0.2656 -0.4367 -0.0944 Syn
0.49_0.00024 -0.1558 -0.3010 -0.0106 Syn
0.49_0.001 -0.1488 -0.2958 -0.0019 Syn
0.49_0.016 -0.3423 -0.5027 -0.1818 Syn
0.49_0.063 -0.3610 -0.5230 -0.1990 Syn
0.49_0.25 -0.2018 -0.3630 -0.0406 Syn
0.49_1 -0.3702 -0.5416 -0.1989 Syn
2_0.004 -0.4096 -0.5827 -0.2365 Syn
2_0.016 -0.5498 -0.7213 -0.3783 Syn
2_0.063 -0.6520 -0.8166 -0.4874 Syn
2_0.25 -0.4539 -0.6150 -0.2927 Syn
2_1 -0.3600 -0.5299 -0.1901 Syn
31_0.016 -0.2525 -0.4075 -0.0974 Syn
31_0.063 -0.2346 -0.3846 -0.0846 Syn
31_0.25 -0.2055 -0.3543 -0.0566 Syn
7.8_0.004 -0.1875 -0.3492 -0.0259 Syn
7.8_0.016 -0.5204 -0.6784 -0.3624 Syn
7.8_0.063 -0.6312 -0.7824 -0.4799 Syn
7.8_0.25 -0.5444 -0.6944 -0.3944 Syn
7.8_1 -0.3674 -0.5232 -0.2115 Syn
Pointwise 95% confidence intervals summary:
Syn Ant Total
21 0 49
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.
<- fitSurface(data, marginalFit,
rsl2 null_model = "loewe2",
B.CP = 50, statistic = "both", parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
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 1.253 0.843
Maximal response 0.076 0.082
log10(EC50) 0.915 -1.041
Common baseline at: 0.115
Exact meanR test (H0 = no synergy/antagonism):
F(49,117) = 17.0629 (p-value < 2e-16)
Evidence for effects in data: Syn
Points with significant deviations from the null:
d1 d2 absR p-value call
0.49_0.016 0.49 0.01600 3.729228 0.01492 Syn
2_0.004 2.00 0.00400 5.769174 <2e-16 Syn
2_0.016 2.00 0.01600 7.898791 <2e-16 Syn
2_0.063 2.00 0.06300 9.581066 <2e-16 Syn
2_0.25 2.00 0.25000 6.165876 <2e-16 Syn
2_1 2.00 1.00000 3.865374 0.00945 Syn
31_0.063 31.00 0.06300 4.399434 0.0012 Syn
31_0.25 31.00 0.25000 5.523963 <2e-16 Syn
31_1 31.00 1.00000 4.938002 0.00013 Syn
7.8_0.00024 7.80 0.00024 4.589089 0.00053 Ant
7.8_0.016 7.80 0.01600 8.052364 <2e-16 Syn
7.8_0.063 7.80 0.06300 11.876320 <2e-16 Syn
7.8_0.25 7.80 0.25000 11.257939 <2e-16 Syn
7.8_1 7.80 1.00000 7.335702 <2e-16 Syn
Overall maxR summary:
Call Syn Ant Total
Syn 13 1 49
CONFIDENCE INTERVALS
Overall effect
Estimated mean departure from null response surface with 95% confidence interval:
-0.224 [-0.2669, -0.1913]
Evidence for effects in data: Syn
Significant pointwise effects
estimate lower upper call
0.12_0.004 -0.2416 -0.4030 -0.0802 Syn
0.12_1 -0.2679 -0.4499 -0.0860 Syn
0.49_0.00024 -0.1581 -0.3117 -0.0045 Syn
0.49_0.016 -0.3391 -0.5056 -0.1727 Syn
0.49_0.063 -0.3662 -0.5490 -0.1834 Syn
0.49_0.25 -0.2163 -0.4066 -0.0261 Syn
0.49_1 -0.3860 -0.5680 -0.2039 Syn
130_0.25 -0.1756 -0.3471 -0.0041 Syn
2_0.004 -0.4378 -0.6046 -0.2710 Syn
2_0.016 -0.5983 -0.7587 -0.4379 Syn
2_0.063 -0.7308 -0.9054 -0.5562 Syn
2_0.25 -0.5495 -0.7303 -0.3686 Syn
2_1 -0.4464 -0.6289 -0.2639 Syn
31_0.016 -0.3227 -0.5022 -0.1432 Syn
31_0.063 -0.3929 -0.5715 -0.2142 Syn
31_0.25 -0.4675 -0.6452 -0.2899 Syn
31_1 -0.4603 -0.6451 -0.2756 Syn
7.8_0.004 -0.2434 -0.4269 -0.0600 Syn
7.8_0.016 -0.6563 -0.8308 -0.4819 Syn
7.8_0.063 -0.8775 -1.0414 -0.7137 Syn
7.8_0.25 -0.8453 -1.0064 -0.6842 Syn
7.8_1 -0.6454 -0.8323 -0.4584 Syn
Pointwise 95% confidence intervals summary:
Syn Ant Total
22 0 49
Α 2-dimensional predicted response surface plot can be generated for a group of null models. In the plot, the points refer to the observed effects whereas the colored lines are the model-predicted responses for the different null models.
The panels correspond to concentration levels of one compound and the x-axis shows the concentration levels of the second compound. The user has the option to define which compound will be shown in the panels and in the x-axis.
<- c("loewe", "loewe2", "bliss", "hsa")
nullModels <- Map(fitSurface, null_model = nullModels, MoreArgs = list(
rs_list data = data, fitResult = marginalFit,
B.CP = 50, statistic = "none", parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
)
synergy_plot_bycomp(rs_list, ylab = "Response", plotBy = "Drug A", color = TRUE)
Presence of synergistic or antagonistic effects can be formalized by means of statistical tests. Two types of tests are considered here and are discussed in more details in the methodology vignette as well as the accompanying paper.
meanR
test evaluates how the predicted response surface based on a specified null model differs from the observed one. If the null hypothesis is rejected, this test suggests that at least some dose combinations may exhibit synergistic or antagonistic behaviour. The meanR
test is not designed to pinpoint which combinations produce these effects nor what type of deviating effect is present.
maxR
test allows to evaluate presence of synergistic/antagonistic effects for each dose combination and as such provides a point-by-point classification.
Both of the above test statistics have a well specified null distribution under a set of assumptions, namely normality of Z-scores. If this assumption is not satisfied, distribution of these statistics can be estimated using bootstrap. Normal approximation is significantly faster whereas bootstrapped distribution of critical values is likely to be more accurate in many practical cases.
Here we will use the previously computed CP
covariance matrix to speed up the process.
<- fitSurface(data, marginalFit,
meanR_N statistic = "meanR", CP = rs$CP, B.B = NULL,
parallel = FALSE)
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)
.
<- fitSurface(data, marginalFit,
meanR_B statistic = "meanR", CP = rs$CP, B.B = 20,
parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
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.
F-statistic | p-value | |
---|---|---|
Normal errors | 16.12679 | 0 |
Bootstrapped errors | 16.12679 | 0 |
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.
<- fitSurface(data, marginalFit,
maxR_N statistic = "maxR", CP = rs$CP, B.B = NULL,
parallel = FALSE)
<- fitSurface(data, marginalFit,
maxR_B statistic = "maxR", CP = rs$CP, B.B = 20,
parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
<- rbind(summary(maxR_N$maxR)$totals,
maxR_both 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 4 given the above calculations.
Call | Syn | Ant | Total | |
---|---|---|---|---|
Normal errors | Syn | 12 | 1 | 49 |
Bootstrapped errors | Syn | 10 | 1 | 49 |
By using the outsidePoints
function, we can obtain a quick summary indicating which dose combinations in Experiment 4 appear to deviate significantly from the null model according to the maxR
statistic.
<- outsidePoints(maxR_B$maxR$Ymean)
outPts kable(outPts, caption = paste0("Non-additive points for Experiment ", i))
d1 | d2 | absR | p-value | call | |
---|---|---|---|---|---|
2_0.004 | 2.0 | 0.00400 | 6.133101 | 0.00 | Syn |
2_0.016 | 2.0 | 0.01600 | 8.011602 | 0.00 | Syn |
2_0.063 | 2.0 | 0.06300 | 10.042789 | 0.00 | Syn |
2_0.25 | 2.0 | 0.25000 | 6.602912 | 0.00 | Syn |
31_0.063 | 31.0 | 0.06300 | 4.077214 | 0.02 | Syn |
31_0.25 | 31.0 | 0.25000 | 4.393387 | 0.00 | Syn |
7.8_0.00024 | 7.8 | 0.00024 | 4.436426 | 0.00 | Ant |
7.8_0.016 | 7.8 | 0.01600 | 8.157618 | 0.00 | Syn |
7.8_0.063 | 7.8 | 0.06300 | 11.914330 | 0.00 | Syn |
7.8_0.25 | 7.8 | 0.25000 | 10.903613 | 0.00 | Syn |
7.8_1 | 7.8 | 1.00000 | 6.158407 | 0.00 | Syn |
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", "white", "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 = "")
The BIGL package also yields effect sizes and corresponding confidence intervals with respect to any response surface. The overall effect size and confidence interval is output in the summary of the ResponseSurface
, but can also be called directly:
summary(maxR_B$confInt)
## Overall effect
## Estimated mean departure from null response surface with 95% confidence interval:
## -0.2208 [-0.2606, -0.1899]
## Evidence for effects in data: Syn
##
## Significant pointwise effects
## estimate lower upper call
## 0.12_0.004 -0.2404 -0.4031 -0.0776 Syn
## 0.12_1 -0.2675 -0.4649 -0.0701 Syn
## 0.49_0.00024 -0.1573 -0.3115 -0.0030 Syn
## 0.49_0.016 -0.3341 -0.5035 -0.1646 Syn
## 0.49_0.063 -0.3613 -0.5353 -0.1872 Syn
## 0.49_0.25 -0.2129 -0.3862 -0.0395 Syn
## 0.49_1 -0.3840 -0.5816 -0.1864 Syn
## 130_0.25 -0.1904 -0.3573 -0.0236 Syn
## 2_0.004 -0.4308 -0.5970 -0.2645 Syn
## 2_0.016 -0.5848 -0.7472 -0.4224 Syn
## 2_0.063 -0.7147 -0.8834 -0.5459 Syn
## 2_0.25 -0.5367 -0.7057 -0.3676 Syn
## 2_1 -0.4389 -0.6376 -0.2402 Syn
## 31_0.016 -0.3208 -0.4905 -0.1512 Syn
## 31_0.063 -0.3865 -0.5556 -0.2174 Syn
## 31_0.25 -0.4510 -0.6282 -0.2739 Syn
## 31_1 -0.4359 -0.6578 -0.2139 Syn
## 7.8_0.004 -0.2392 -0.4145 -0.0638 Syn
## 7.8_0.016 -0.6429 -0.8115 -0.4742 Syn
## 7.8_0.063 -0.8496 -1.0132 -0.6860 Syn
## 7.8_0.25 -0.8124 -0.9761 -0.6487 Syn
## 7.8_1 -0.6211 -0.8256 -0.4166 Syn
##
## Pointwise 95% confidence intervals summary:
## Syn Ant Total
## 22 0 49
In addition, a contour plot can be made with pointwise confidence intervals. Contour plot colouring can be defined according to the effect sizes or according to maxR results.
plotConfInt(maxR_B, color = "effect-size")
You can also customize the coloring of the contour plot and 3-dimensional predicted response surface plot based on effect sizes:
contour(
maxR_B,colorPalette = c("Syn" = "blue", "None" = "white", "Ant" = "red"),
main = paste0(" Experiment ", i, " contour plot for effect size"),
colorBy = "effect-size",
scientific = TRUE, digits = 3, cutoff = cutoff
)
plot(maxR_B, color = "effect-size", legend = FALSE, main = "", gradient = FALSE,
colorPalette = c("Ant" = "red", "None" = "white", "Syn" = "blue"),
colorPaletteNA = "white")
Starting from the package version 1.2.0
the variance can be estimated separately for on-axis (monotherapy) and off-axis points using method
argument to fitSurface
. The possible values for method
are:
"equal"
, equal variances assumed (as above, default),"unequal"
, variance is estimated separately for on-axis and off-axis points,"model"
, the variance is modelled as a function of the mean.Please see the methodology vignette for details. Below we show an example analysis in such case. Note that transformations are not possible if variances are not assumed equal.
<- fitMarginals(data, transforms = NULL)
marginalFit summary(marginalFit)
## Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
## Transformations: No
##
## Compound 1 Compound 2
## Slope 1.223 1.008
## Maximal response 214.310 462.322
## log10(EC50) 0.457 -1.639
##
## Common baseline at: 3979.458
<- fitSurface(data, marginalFit, method = "unequal",
resU statistic = "both", B.CP = 20, B.B = 20, parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
summary(resU)
## Null model: Generalized Loewe Additivity
## Variance assumption used: "unequal"
## Mean occupancy rate: 0.7222181
##
## Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
## Transformations: No
##
## Compound 1 Compound 2
## Slope 1.223 1.008
## Maximal response 214.310 462.322
## log10(EC50) 0.457 -1.639
##
## Common baseline at: 3979.458
##
## Bootstrapped meanR test (H0 = no synergy/antagonism):
## F(49,117) = 7.948 (p-value < 2e-16)
##
## Evidence for effects in data: Syn
## Points with significant deviations from the null:
## d1 d2 absR p-value call
## 0.12_0.004 0.12 0.00400 7.666810 <2e-16 Syn
## 0.49_0.00024 0.49 0.00024 4.896459 <2e-16 Syn
## 0.49_0.001 0.49 0.00100 4.378196 <2e-16 Syn
## 0.49_0.016 0.49 0.01600 5.902539 <2e-16 Syn
## 2_0.004 2.00 0.00400 7.032048 <2e-16 Syn
## 2_0.016 2.00 0.01600 7.318842 <2e-16 Syn
## 2_0.063 2.00 0.06300 5.169768 <2e-16 Syn
## 7.8_0.016 7.80 0.01600 4.265172 <2e-16 Syn
## 7.8_0.063 7.80 0.06300 4.378360 <2e-16 Syn
##
## Overall maxR summary:
## Call Syn Ant Total
## Syn 9 0 49
##
## CONFIDENCE INTERVALS
## Overall effect
## Estimated mean departure from null response surface with 95% confidence interval:
## -194.23 [-293.6487, 25.3107]
## Evidence for effects in data: None
##
## Significant pointwise effects
## estimate lower upper call
## 0.12_0.004 -744.1210 -957.2780 -530.9641 Syn
## 0.49_0.00024 -465.6427 -688.9717 -242.3138 Syn
## 0.49_0.001 -427.0063 -646.9899 -207.0227 Syn
## 0.49_0.004 -425.2081 -640.8506 -209.5656 Syn
## 0.49_0.016 -636.8597 -853.7796 -419.9397 Syn
## 0.49_0.063 -363.1038 -583.2415 -142.9662 Syn
## 2_0.004 -745.5425 -962.5064 -528.5787 Syn
## 2_0.016 -781.8747 -995.8941 -567.8553 Syn
## 2_0.063 -567.9238 -783.9720 -351.8756 Syn
## 2_0.25 -297.3783 -513.0082 -81.7483 Syn
## 7.8_0.004 -263.1519 -493.1373 -33.1664 Syn
## 7.8_0.016 -483.6000 -707.7305 -259.4696 Syn
## 7.8_0.063 -475.0748 -691.5275 -258.6220 Syn
## 7.8_0.25 -355.7178 -570.5703 -140.8652 Syn
## 7.8_1 -263.4871 -482.1530 -44.8211 Syn
##
## Pointwise 95% confidence intervals summary:
## Syn Ant Total
## 15 0 49
For the variance model, an exploratory plotting function is available to explore the relationship between the mean and the variance.
plotMeanVarFit(data)
plotMeanVarFit(data, log = "xy") #Clearer on the log-scale
plotMeanVarFit(data, trans = "log") #Thresholded at maximum observed variance
The linear fit seems fine in this case.
<- fitSurface(data, marginalFit, method = "model",
resM statistic = "both", B.CP = 20, B.B = 20, parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
If the log transformation yielded a better fit, then this could be achieved by using the following option.
<- fitSurface(data, marginalFit, method = "model", trans = "log",
resL statistic = "both", B.CP = 20, B.B = 20, parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR")
Negative variances were modelled, but variance model has the smallest observed variances as minimum so we can proceed.`
summary(resM)
## Null model: Generalized Loewe Additivity
## Variance assumption used: "model"
## Mean occupancy rate: 0.7222181
##
## Formula: effect ~ fn(h1, h2, b, m1, m2, e1, e2, d1, d2)
## Transformations: No
##
## Compound 1 Compound 2
## Slope 1.223 1.008
## Maximal response 214.310 462.322
## log10(EC50) 0.457 -1.639
##
## Common baseline at: 3979.458
##
## Bootstrapped meanR test (H0 = no synergy/antagonism):
## F(49,117) = 63.5616 (p-value < 2e-16)
##
## Evidence for effects in data: Syn
## Points with significant deviations from the null:
## d1 d2 absR p-value call
## 0.49_0.016 0.49 0.01600 4.002029 0.05 Syn
## 130_0.00024 130.00 0.00024 20.154443 <2e-16 Ant
## 130_0.001 130.00 0.00100 12.860045 <2e-16 Ant
## 130_0.016 130.00 0.01600 5.446635 <2e-16 Syn
## 130_0.063 130.00 0.06300 4.270034 0.05 Syn
## 130_0.25 130.00 0.25000 4.873464 <2e-16 Syn
## 130_1 130.00 1.00000 13.713871 <2e-16 Ant
## 2_0.004 2.00 0.00400 5.370226 <2e-16 Syn
## 2_0.016 2.00 0.01600 6.836464 <2e-16 Syn
## 2_0.063 2.00 0.06300 6.749047 <2e-16 Syn
## 2_1 2.00 1.00000 8.380021 <2e-16 Ant
## 31_0.004 31.00 0.00400 17.764933 <2e-16 Ant
## 31_0.016 31.00 0.01600 11.383025 <2e-16 Syn
## 31_0.063 31.00 0.06300 10.803941 <2e-16 Syn
## 31_1 31.00 1.00000 4.043937 0.05 Syn
## 500_0.00024 500.00 0.00024 3.973185 0.05 Ant
## 500_0.001 500.00 0.00100 4.881484 <2e-16 Syn
## 500_0.063 500.00 0.06300 6.560268 <2e-16 Syn
## 500_0.25 500.00 0.25000 10.550823 <2e-16 Syn
## 500_1 500.00 1.00000 8.463668 <2e-16 Ant
## 7.8_0.016 7.80 0.01600 6.851392 <2e-16 Syn
## 7.8_0.063 7.80 0.06300 27.710960 <2e-16 Syn
## 7.8_0.25 7.80 0.25000 15.181416 <2e-16 Syn
## 7.8_1 7.80 1.00000 10.426276 <2e-16 Syn
##
## Overall maxR summary:
## Call Syn Ant Total
## Syn 17 7 49
##
## CONFIDENCE INTERVALS
## Overall effect
## Estimated mean departure from null response surface with 95% confidence interval:
## -194.23 [-296.6101, -115.3013]
## Evidence for effects in data: Syn
##
## Significant pointwise effects
## estimate lower upper call
## 2_0.004 -745.5425 -1408.0272 -83.0579 Syn
## 2_0.016 -781.8747 -1296.3699 -267.3795 Syn
## 2_0.063 -567.9238 -916.9363 -218.9113 Syn
## 2_0.25 -297.3783 -543.5988 -51.1578 Syn
## 31_0.016 -139.2683 -253.9192 -24.6174 Syn
## 31_0.063 -160.0900 -282.7304 -37.4496 Syn
## 31_0.25 -183.3832 -356.3636 -10.4029 Syn
## 7.8_0.016 -483.6000 -755.1243 -212.0758 Syn
## 7.8_0.063 -475.0748 -579.7293 -370.4203 Syn
## 7.8_0.25 -355.7178 -521.7351 -189.7004 Syn
##
## Pointwise 95% confidence intervals summary:
## Syn Ant Total
## 10 0 49
In order to proceed with multiple experiments, we repeat the same procedure as previously. We collect all the necessary objects for which estimations do not have to be repeated to generate meanR
and maxR
statistics in a simple list.
<- list()
marginalFits <- list()
datasets <- list()
respSurfaces <- list()
maxR.summary for (i in seq_len(nExp)) {
## Select experiment
<- subsetData(directAntivirals, i)
data ## Fit joint marginal model
<- fitMarginals(data, transforms = transforms,
marginalFit method = "nlslm")
## Predict response surface based on generalized Loewe model
<- fitSurface(data, marginalFit,
respSurface statistic = "maxR", B.CP = 20,
parallel = FALSE,
wild_bootstrap = TRUE, wild_bootType = "normal",
control = "dFCR"
)
<- data
datasets[[i]] <- marginalFit
marginalFits[[i]] <- respSurface
respSurfaces[[i]] <- summary(respSurface$maxR)$totals
maxR.summary[[i]] }
We use the maxR
procedure with a chosen p-value cutoff of 0.95. If maxR
statistic falls outside the 95th percentile of its distribution (either bootstrapped or not), the respective off-axis dose combination is said to deviate significantly from the generalized Loewe model and the algorithm determines whether it deviates in a synergistic or antagonistic way.
Below is the summary of overall calls and number of deviating points for each experiment.
Call | Syn | Ant | Total | |
---|---|---|---|---|
Experiment 1 | Syn | 3 | 0 | 49 |
Experiment 2 | Syn | 3 | 0 | 49 |
Experiment 3 | Syn | 6 | 0 | 49 |
Experiment 4 | Syn | 14 | 1 | 49 |
Previous summarizing and visual analysis can be repeated on each of the newly defined experiments. For example, Experiment 4
indicates a total of 16 combinations that were called synergistic according to the maxR
test.
d1 | d2 | absR | p-value | call | |
---|---|---|---|---|---|
0.49_0.016 | 0.49 | 0.01600 | 3.461997 | 0.03624 | Syn |
0.49_0.063 | 0.49 | 0.06300 | 3.719475 | 0.01500 | Syn |
2_0.004 | 2.00 | 0.00400 | 5.558917 | 0.00000 | Syn |
2_0.016 | 2.00 | 0.01600 | 7.617450 | 0.00000 | Syn |
2_0.063 | 2.00 | 0.06300 | 9.780369 | 0.00000 | Syn |
2_0.25 | 2.00 | 0.25000 | 6.430839 | 0.00000 | Syn |
2_1 | 2.00 | 1.00000 | 3.565155 | 0.02538 | Syn |
31_0.016 | 31.00 | 0.01600 | 3.560541 | 0.02575 | Syn |
31_0.063 | 31.00 | 0.06300 | 4.447188 | 0.00102 | Syn |
31_0.25 | 31.00 | 0.25000 | 4.806716 | 0.00025 | Syn |
7.8_0.00024 | 7.80 | 0.00024 | 4.748063 | 0.00034 | Ant |
7.8_0.016 | 7.80 | 0.01600 | 7.879730 | 0.00000 | Syn |
7.8_0.063 | 7.80 | 0.06300 | 11.721510 | 0.00000 | Syn |
7.8_0.25 | 7.80 | 0.25000 | 10.852142 | 0.00000 | Syn |
7.8_1 | 7.80 | 1.00000 | 6.538138 | 0.00000 | Syn |
Consequently, above table for Experiment 4
can be illustrated in a contour plot.
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.