,
Zhi Zhao [ctb],
Tero Aittokallio [ctb]| Title: | Statistical Framework for in Vivo Drug Combination Studies |
|---|---|
| Description: | A framework for evaluating drug combination effects in preclinical in vivo studies. 'SynergyLMM' provides functions to analyze longitudinal tumor growth experiments using linear mixed-effects models, perform time-dependent analyses of synergy and antagonism, evaluate model diagnostics and performance, and assess both post-hoc and a priori statistical power. The calculation of drug combination synergy follows the statistical framework provided by Demidenko and Miller (2019, <doi:10.1371/journal.pone.0224137>). The implementation and analysis of linear mixed-effect models is based on the methods described by Pinheiro and Bates (2000, <doi:10.1007/b98882>), and Gałecki and Burzykowski (2013, <doi:10.1007/978-1-4614-3900-4>). |
| Authors: | Rafael Romero-Becerra [aut, cre]
,
Zhi Zhao [ctb],
Tero Aittokallio [ctb] |
| Maintainer: | Rafael Romero-Becerra <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 1.0.1.9000 |
| Built: | 2025-02-11 10:25:23 UTC |
| Source: | https://github.com/rafromb/synergylmm |
A priori power calculation for a hypothetical two-drugs combination study of synergy using linear-mixed models with varying drug combination effect and/or experimental variability.
APrioriPwr( npg = 5, time = c(0, 3, 5, 10), grwrControl = 0.08, grwrA = 0.07, grwrB = 0.06, grwrComb = 0.03, sd_ranef = 0.01, sgma = 0.1, sd_eval = NULL, sgma_eval = NULL, grwrComb_eval = NULL, method = "Bliss", ... )APrioriPwr( npg = 5, time = c(0, 3, 5, 10), grwrControl = 0.08, grwrA = 0.07, grwrB = 0.06, grwrComb = 0.03, sd_ranef = 0.01, sgma = 0.1, sd_eval = NULL, sgma_eval = NULL, grwrComb_eval = NULL, method = "Bliss", ... )
npg |
Number of subjects per group. |
time |
Vector with the times at which the tumor volume measurements have been performed. |
grwrControl |
Coefficient for Control treatment group tumor growth rate. |
grwrA |
Coefficient for Drug A treatment group tumor growth rate. |
grwrB |
Coefficient for Drug B treatment group tumor growth rate. |
grwrComb |
Coefficient for Combination (Drug A + Drug B) treatment group tumor growth rate. |
sd_ranef |
Random effects standard deviation (between-subject variance) for the model. |
sgma |
Residuals standard deviation (within-subject variance) for the model. |
sd_eval |
A vector with values representing the standard deviations of random effects, through which the power for synergy calculation will be evaluated. |
sgma_eval |
A vector with values representing the standard deviations of the residuals, through which the power for synergy calculation will be evaluated. |
grwrComb_eval |
A vector with values representing the coefficients for Combination treatment group tumor growth rate, through which the power for synergy calculation will be evaluated. |
method |
String indicating the method for synergy calculation. Possible methods are "Bliss" and "HSA", corresponding to Bliss and highest single agent, respectively. |
... |
Additional parameters to be passed to nlmeU::Pwr.lme method. |
APrioriPwr allows for total customization of an hypothetical drug combination study and allows the user
to define several experimental parameters, such as the sample size, time of measurements, or drug effect,
for the power evaluation of synergy for Bliss and HSA reference models. The power calculation is
based on F-tests of the fixed effects of the model as previously described (Helms, R. W. (1992),
Verbeke and Molenberghs (2009), Gałecki and Burzykowski (2013)).
The focus the power analysis with APrioriPwr is on the drug combination effect and the variability in the
experiment. For other a priori power analysis see also PwrSampleSize() and PwrTime().
npg, time, grwrControl, grwrA, grwrB, grwrComb, sd_ranef and sgma are parameters referring to
the initial exemplary data set and corresponding fitted model. These values can be obtained from a fitted model, using lmmModel_estimates(),
or be defined by the user.
sd_eval, sgma_eval, and grwrComb_eval are the different values that will be modified in the initial
exemplary data set to fit the corresponding model and calculate the power.
The functions returns several plots:
A plot representing the hypothetical data, with the regression lines for each
treatment group according to grwrControl, grwrA, grwrB and grwrComb values. The values
assigned to sd_ranef and sgma are also shown.
A plot showing the values of the power calculation depending on the values assigned to
sd_eval and sgma_eval. The power result corresponding to the values assigned to sd_ranef and
sgma is shown with a red dot.
A plot showing the values of the power calculation depending on the values assigned to
grwrComb_eval. The vertical dashed line indicates the value of grwrComb. The horizontal
line indicates the power of 0.80.
The statistical power for the fitted model for the initial data set according to the values given by
npg, time, grwrControl, grwrA, grwrB, grwrComb, sd_ranef and sgma is also shown in the console.
Helms, R. W. (1992). Intentionally incomplete longitudinal designs: I. Methodology and comparison of some full span designs. Statistics in Medicine, 11(14–15), 1889–1913. https://doi.org/10.1002/sim.4780111411
Verbeke, G. & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer New York. https://doi.org/10.1007/978-1-4419-0300-6
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
PostHocPwr,PwrSampleSize(), PwrTime().
APrioriPwr( sd_eval = seq(0.01, 0.2, 0.01), sgma_eval = seq(0.01, 0.2, 0.01), grwrComb_eval = seq(0.01, 0.1, 0.005) )APrioriPwr( sd_eval = seq(0.01, 0.2, 0.01), sgma_eval = seq(0.01, 0.2, 0.01), grwrComb_eval = seq(0.01, 0.1, 0.005) )
CookDistance allows the user to identify those subjects with a greater influence in the estimation of the
(tumor growth rate) for the treatment group, based in the calculation of Cook's distances.
CookDistance(model, cook_thr = NA, label_angle = 0, verbose = TRUE)CookDistance(model, cook_thr = NA, label_angle = 0, verbose = TRUE)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
cook_thr |
Numeric value indicating the threshold for the Cook's distance. If not specified, the threshold is set to the 90% percentile of the Cook's distance values. |
label_angle |
Numeric value indicating the angle for the label of subjects with a Cook's distance greater than |
verbose |
Logical indicating if the subjects with a Cook's distance greater than |
The identification of the subjects with a greater influence in each estimated representing the tumor growth is based on the calculation of Cook's distances, as
described in Gałecki and Burzykowsk (2013). To compute the Cook's distance for the estimates (i.e., the contribution to each subject to the coefficient of its treatment group),
first a matrix containing the leave-one-subject-out estimates or is calculated. Then, the Cook's distances are calculated according to:
where is the estimate of the parameter vector obtained by fitting the model to the data with the -th subject excluded. The denominator of
the expression is equal to the number of the fixed-effects coefficients, which, under the assumption that the design matrix is of full rank, is equivalent to the rank of the design matrix.
A plot of the Cook's distance value for each subject, indicating those subjects
whose Cook's distance is greater than cook_thr.
If saved to a variable, the function returns a vector with the Cook's distances for each subject.
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
#' # Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Calulate Cook's distances for each subject CookDistance(model = lmm) # Change the Cook's distance threshold CookDistance(model = lmm, cook_thr = 0.15)#' # Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Calulate Cook's distances for each subject CookDistance(model = lmm) # Change the Cook's distance threshold CookDistance(model = lmm, cook_thr = 0.15)
getRTV is a helper function used inside lmmModel() to obtain a dataframe with a column RTV corresponding
to the relative tumor volume to time time_start, and a column logRTV with the logarithm of RTV.
getRTV(data, time_start)getRTV(data, time_start)
data |
Data frame with the tumor growth data. The input data frame columns have to have the following names:
|
time_start |
Numeric value indicating the time at which the treatment started. |
The function returns the original data frame of tumor growth data, with 3 additional columns, corresponding to:
RTV: Relative tumor volume to the tumor volume at start_time.
logRTV: Logarithm of RTV column.
TV0: Tumor volume at start_time.
# Load example dataset data("grwth_data") # Change column names colnames(grwth_data) <- c("SampleID", "Time", "Treatment", "TV") # Calculate relative tumor volume getRTV(data = grwth_data, time_start = 0)# Load example dataset data("grwth_data") # Change column names colnames(grwth_data) <- c("SampleID", "Time", "Treatment", "TV") # Calculate relative tumor volume getRTV(data = grwth_data, time_start = 0)
A long format data frame with example tumor growth data, generated with simulateTumorGrowth(),
representing a simulated in vivo two-drugs combination experiment, involving 8 subjects per group, with 30 days of follow-up, taking measurement
every 3 days.
data(grwth_data)data(grwth_data)
A data frame with 352 rows and 4 columns:
Subject IDs.
Time for each measurement in an arbitrary unit, for example, days.
Treatment group for each subject (Control, DrugA, DrugB, and Combination).
Measurement representing the tumor volume in an arbitrary unit, for example, mm^3.
lmmModel() fits a linear mixed effect model from a tumor growth dataset. The input data frame must be in long format and include at least the following columns: column with the sample ids,
column with the time at which each measurement has been done, a column indicating the treatment group, and a column with the tumor measurement (e.g. tumor volume).
lmmModel( data, sample_id = "SampleID", time = "Time", treatment = "Treatment", tumor_vol = "TV", trt_control = "Control", drug_a = "Drug_A", drug_b = "Drug_B", drug_c = NA, combination = "Combination", time_start = NULL, time_end = NULL, min_observations = 1, show_plot = TRUE, tum_vol_0 = "ignore", ... )lmmModel( data, sample_id = "SampleID", time = "Time", treatment = "Treatment", tumor_vol = "TV", trt_control = "Control", drug_a = "Drug_A", drug_b = "Drug_B", drug_c = NA, combination = "Combination", time_start = NULL, time_end = NULL, min_observations = 1, show_plot = TRUE, tum_vol_0 = "ignore", ... )
data |
A data frame with the tumor growth data, in long format. It must contain at least the following columns: mice IDs, time of follow-up (numeric number), treatment and tumor volume (numeric number). |
sample_id |
String indicating the name of the column with the mice IDs. |
time |
String indicating the name of the column with the time of follow-up. |
treatment |
String indicating the name of the column with the treatment corresponding to each sample. |
tumor_vol |
String indicating the name of the column with the tumor volume (or any other measurement representing the tumor growth). |
trt_control |
String indicating the name assigned to the 'Control' group. |
drug_a |
String indicating the name assigned to the 'Drug A' group. |
drug_b |
String indicating the name assigned to the 'Drug B' group. |
drug_c |
String indicating the name assigned to the 'Drug C' group (if present). |
combination |
String indicating the name assigned to the Combination ('Drug A' + 'Drug B', or 'Drug A' + 'Drug B' + 'Drug C') group. |
time_start |
Numeric value indicating the time point at which the treatment started. If not
specified, the minimum value in the |
time_end |
Numeric value indicating the last time point to be included in the analysis. If not
specified, the maximum value in the |
min_observations |
Minimum number of observation for each sample to be included in the analysis. |
show_plot |
Logical indicating if a plot for the log of the relative tumor volume (RTV) vs Time for each sample, and the model calculated marginal slope for each treatment, should be produced. |
tum_vol_0 |
Select the behavior of the function regarding measurements in which the tumor measurement is 0, and therefore the logarithmic transformation is not possible. Possible options are 'ignore' (default), to ignore these measurements, or 'transform', to add 1 unit to all measurements before the log transformation. |
... |
Additional arguments to be passed to |
lmmModel() relies in the assumption that tumor growth follows an exponential kinetics. Any departure from this assumption can be tested using the diagnostics functions ranefDiagnostics(),
residDiagnostics(), and ObsvsPred().
The model formula for the linear mixed-effect fitted model is:
where denotes the value of the logarithm of the relative tumor volume measured for subject at time . The coefficient
represents the fixed effectsat time for each treatment , where in the case of
two-drugs combination experiments, or in the case of three-drugs combination experiments, and indicates the
tumor-specific growth rate for each treatment group.
Term corresponds to the subject-specific random slope that takes into account the longitudinal nature of the data, where is the random effect for subject .
is the residual random term for subject at time .
The implementation of the linear mixed model in lmmModel() is done using nlme::lme, which also allows for the
specification of within-group correlations structures and/or unequal variances. These, and additional parameters,
can be passed to the nlme::lme function through the ... argument for fitting the model (see examples below).
An object of class "lme" (see nlme::lme for details) representing the linear mixed-effects model fit. If show_plot = TRUE, the plot
of the tumor growth data obtained with plot_lmmModel() is also shown.
Pinheiro JC, Bates DM (2000). Mixed-Effects Models in S and S-PLUS. Springer, New York. doi:10.1007/b98882 doi:10.1007/b98882.
Pinheiro J, Bates D, R Core Team (2024). nlme: Linear and Nonlinear Mixed Effects Models. R package version 3.1-166, https://CRAN.R-project.org/package=nlme.
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
data("grwth_data") # Most simple model lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Changing the last time point of follow-up lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", time_end = 21 ) # Adding additional parameters for model fitting lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", # Adding variance function to represent a different variance per subject weights = nlme::varIdent(form = ~1|SampleID), # Specifiying control values for lme Fit (useful when convergence problems appear) control = nlme::lmeControl(maxIter = 1000, msMaxIter = 1000, niterEM = 100, msMaxEval = 1000) ) # Fit a model specifying a different variance per Time lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", # Adding variance function to represent a different variance per Time weights = nlme::varIdent(form = ~1|Time) )data("grwth_data") # Most simple model lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Changing the last time point of follow-up lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", time_end = 21 ) # Adding additional parameters for model fitting lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", # Adding variance function to represent a different variance per subject weights = nlme::varIdent(form = ~1|SampleID), # Specifiying control values for lme Fit (useful when convergence problems appear) control = nlme::lmeControl(maxIter = 1000, msMaxIter = 1000, niterEM = 100, msMaxEval = 1000) ) # Fit a model specifying a different variance per Time lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", # Adding variance function to represent a different variance per Time weights = nlme::varIdent(form = ~1|Time) )
lmmModel_estimates allows the user to easily extract some of the interesting model estimates for further use in other functions,
such as for power calculation.
lmmModel_estimates(model)lmmModel_estimates(model)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
The model estimates provided by lmmModel_estimates include:
Fixed effect coefficients: , , , ,
which represent the estimated specific growth rates for the Control, Drug A, Drug B and Combination groups, respectively.
These are shown in columns control, drug_a, drug_b, and combination, respectively.
Standard deviation of the random effects (between-subject variance). Column sd_ranef.
Standard deviation of the residuals (within-subject variance). Column sd_resid.
A data frame with the estimated values for the coefficients of the tumor growth for each treatment,
the standard deviation of the random effects, and the standard deviation of the residuals of the model.
These values can be useful for the power analysis of the model using APrioriPwr().
data("grwth_data") # Fit example model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Get the estimates lmmModel_estimates(lmm)data("grwth_data") # Fit example model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Get the estimates lmmModel_estimates(lmm)
lmmSynergy allows for the calculation of synergy using 3 different references models: Bliss independence, highest single agent and
response additivity. The calculation of synergy is based on hypothesis testing on the coefficient estimates from the model fitted by
lmmModel().
lmmSynergy( model, method = "Bliss", min_time = 0, robust = FALSE, type = "CR2", ra_nsim = 1000, show_plot = TRUE, ... )lmmSynergy( model, method = "Bliss", min_time = 0, robust = FALSE, type = "CR2", ra_nsim = 1000, show_plot = TRUE, ... )
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
method |
String indicating the method for synergy calculation. Possible methods are "Bliss", "HSA" and "RA", corresponding to Bliss, highest single agent and response additivity, respectively. |
min_time |
Minimun time for which to start calculating synergy. |
robust |
If TRUE, uncertainty is estimated using sandwich-based robust estimators of the variance-covariance matrix of the regression coefficient estimates provided by clubSandwich::vcovCR.lme. |
type |
Character string specifying which small-sample adjustment should be used, with available options "CR0", "CR1", "CR1p", "CR1S", "CR2", or "CR3".
See "Details" section of |
ra_nsim |
Number of random sampling to calculate the synergy for Response Additivity model. |
show_plot |
Logical indicating if a plot with the results of the synergy calculation should be generated. |
... |
Additional arguments to be passed to |
lmmSynergy uses the statistical description provided by Demidenko and Miller (2019) for the calculation of synergy. It is based on hypothesis testing
on the coefficients estimates from the model fitted by lmmModel(): , , , ,
which represent the estimated specific growth rates for the Control, Drug A, Drug B and Combination groups, respectively.
Bliss Indepence Model
For Bliss model, lmmSynergy test the following null hypothesis:
Two-drugs combination experiment:
Three-drugs combination experiment:
Highes Single Agent (HSA)
Two-drugs combination experiment:
For the HSA model, lmmSynergy test the following null hypothesis:
Three-drugs combination experiment:
For the HSA model, lmmSynergy test the following null hypothesis:
Response Additivity (RA)
For the RA model, lmmSynergy test the following null hypothesis:
Two-drugs combination experiment:
Three-drugs combination experiment:
For Bliss and HSA models, lmmSynergy uses marginaleffects::hypotheses() to conduct hypothesis tests on the estimated coefficients of the model.
In the case of the RA model, the null hypothesis is tested comparing the area under the curve (i.e. cumulative effect from the beginning of a treatment to
a time point of interest) obtained from each side of the equation for the null hypothesis, based on ra_sim random samplings from the
distribution of the coefficients.
Combination Index and Synergy Score
The results obtained by lmmSynergy include the synergy score (SS) and combination index (CI) for the model, for each time point, together with their confidence interval,
and the corresponding p-value. The values of SS and CI provided by lmmSynergy follow previous definitions of these metrics so they have the same interpretation:
The SS has been defined as the excess response due to drug interaction compared to the reference model (Ianevski et al. (2017), Ianevski, Giri, and Aittokallio (2022), Mao and Guo (2023)).
Following this definition, a , , and , represent synergistic, additive and antagonistic effects, respectively.
According to the common definition of the CI, a , , and represent synergistic, additive and antagonistic effects, respectively (Yadav et al. (2015), Demidenko and Miller (2019),
Mao and Guo (2023)), and provides information about the observed drug combination effect versus the expected additive effect provided by the reference synergy model.
A drug combination effect larger than the expected () would indicate synergism, a drug combination effect equal to the expected () would indicate additivity,
and a lower drug combination effect than the expected () would indicate antagonism.
As mentioned above, the results include the synergy results for each day. This means that lmmSynergy refits the model using the data from time_start defined in lmmModel() until
each time point, providing the synergy results for each of these models and for that specific time point.
Uncertainty estimation using robust estimators
If robust = TRUE, lmmSynergy deals with possible model misspecifications, allowing for cluster-robust variance estimation using clubSandwich::vcovCR.lme.
When using robust = TRUE, setting type = "CR2" is recommended. See more details in clubSandwich::vcovCR().
Note: When a variance structure has been specified in the model it is recommended to use always robust = TRUE to get a better estimation.
The function returns a list with two elements:
Constrasts: List with the outputs of the linear test for the synergy null hypothesis obtained by marginaleffects::hypotheses() for each time.
See marginaleffects::hypotheses() for more details.
Synergy: Data frame with the synergy results, indicating the model of synergy ("Bliss", "HSA" or "RA"), the metric (combination index and synergy score),
the value of the metric estimate (with upper and lower confidence interval bounds) and the p-value, for each time.
If show_plot = TRUE, a plot with the synergy results obtained with plot_lmmSynergy() is also shown.
Demidenko, Eugene, and Todd W. Miller. 2019. Statistical Determination of Synergy Based on Bliss Definition of Drugs Independence. PLoS ONE 14 (November). https://doi.org/10.1371/journal.pone.0224137.
Yadav, Bhagwan, Krister Wennerberg, Tero Aittokallio, and Jing Tang. 2015. Searching for Drug Synergy in Complex Dose–Response Landscapes Using an Interaction Potency Model. Computational and Structural Biotechnology Journal 13: 504–13. https://doi.org/10.1016/j.csbj.2015.09.001.
Ianevski, Aleksandr, Liye He, Tero Aittokallio, and Jing Tang. 2017. SynergyFinder: A Web Application for Analyzing Drug Combination Dose–Response Matrix Data. Bioinformatics 33 (August): 2413–15. https://doi.org/10.1093/bioinformatics/btx162.
Ianevski, Aleksandr, Anil K Giri, and Tero Aittokallio. 2022. SynergyFinder 3.0: An Interactive Analysis and Consensus Interpretation of Multi-Drug Synergies Across Multiple Samples. Nucleic Acids Research 50 (July): W739–43. https://doi.org/10.1093/nar/gkac382.
Mao, Binchen, and Sheng Guo. 2023. Statistical Assessment of Drug Synergy from in Vivo Combination Studies Using Mouse Tumor Models. Cancer Research Communications 3 (October): 2146–57. https://doi.org/10.1158/2767-9764.CRC-23-0243.
Vincent Arel-Bundock, Noah Greifer, and Andrew Heiss. Forthcoming. How to Interpret Statistical Models Using marginaleffects in R and Python. Journal of Statistical Software. https://marginaleffects.com
# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Most simple use with default values syn <- lmmSynergy(lmm) # Accessing to synergy results data frame syn$Synergy # Selecting different reference models: ## Bliss lmmSynergy(lmm, method = "Bliss") ## HSA lmmSynergy(lmm, method = "HSA") ## RA lmmSynergy(lmm, method = "RA", ra_sim = 1000) # Only calculate synergy from Time 12 onwards lmmSynergy(lmm, min_time = 12) # Using robust standard errors lmmSynergy(lmm, method = "Bliss", robust = TRUE, type = "CR2")# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Most simple use with default values syn <- lmmSynergy(lmm) # Accessing to synergy results data frame syn$Synergy # Selecting different reference models: ## Bliss lmmSynergy(lmm, method = "Bliss") ## HSA lmmSynergy(lmm, method = "HSA") ## RA lmmSynergy(lmm, method = "RA", ra_sim = 1000) # Only calculate synergy from Time 12 onwards lmmSynergy(lmm, min_time = 12) # Using robust standard errors lmmSynergy(lmm, method = "Bliss", robust = TRUE, type = "CR2")
logLikSubjectDisplacements allows the user to evaluate the log-likelihood displacement for each subject,
indicating the influence of every subject to the model.
logLikSubjectDisplacements( model, disp_thrh = NA, label_angle = 0, var_name = NULL, verbose = TRUE, ... )logLikSubjectDisplacements( model, disp_thrh = NA, label_angle = 0, var_name = NULL, verbose = TRUE, ... )
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
disp_thrh |
Numeric value indicating the threshold of log-likelihood displacement. If not specified, the threshold is set to the 90% percentile of the log-likelihood displacement values. |
label_angle |
Numeric value indicating the angle for the label of subjects with a log-likelihood displacement greater than |
var_name |
Name of the variable for the weights of the model in the case that a variance structure has been specified using |
verbose |
Logical indicating if subjects with a log-likelihood displacement greater than |
... |
Extra arguments, if any, for lattice::panel.xyplot. |
The evaluation of the log-likelihood displacement is based in the analysis proposed in Verbeke and Molenberghs (2009) and Gałecki and Burzykowski (2013).
First, a list of models fitted to leave-one-subject-out datasets are obtained. Then, for each model, the maximum likelihood estimate obtained by fitting the
model to all data and the maximum likelihood estimate obtained by fitting the model to the data with the -th subject removed are obtained and used for the
log-likelihood displacement calculation. The likelihood displacement, , is defined as twice the difference between the log-likelihood computed at a
maximum and displaced values of estimated parameters (Verbeke and Molenberghs (2009), Gałecki and Burzykowsk (2013)):
where is the maximum-likelihood estimate of obtained by fitting the model to all data, while is
the maximum-likelihood estimate obtained by fitting the model to the data with the -subject excluded.
Returns a plot of the log-likelihood displacement values for each subject, indicating those subjects
whose contribution is greater than disp_thrh.
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
Verbeke, G. & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer New York. https://doi.org/10.1007/978-1-4419-0300-6
# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Obtain log-likelihood displacement for each subject logLikSubjectDisplacements(model = lmm) # Modifying the threshold for log-likelihood displacement logLikSubjectDisplacements(model = lmm, disp_thrh = 1) # Calculating the log-likelihood contribution in a model with a variance structure specified lmm_var <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", weights = nlme::varIdent(form = ~ 1|SampleID) ) # Calculate the log-likelihood contribution logLikSubjectDisplacements(model = lmm, var_name = "SampleID")# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Obtain log-likelihood displacement for each subject logLikSubjectDisplacements(model = lmm) # Modifying the threshold for log-likelihood displacement logLikSubjectDisplacements(model = lmm, disp_thrh = 1) # Calculating the log-likelihood contribution in a model with a variance structure specified lmm_var <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", weights = nlme::varIdent(form = ~ 1|SampleID) ) # Calculate the log-likelihood contribution logLikSubjectDisplacements(model = lmm, var_name = "SampleID")
ObsvsPred allows the user to have a straight forward idea about how the model is fitting the data, providing
plots of the predicted regression lines versus the actual data points.
ObsvsPred(model, nrow = 4, ncol = 5, ...)ObsvsPred(model, nrow = 4, ncol = 5, ...)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
nrow |
Number of rows of the layout to organize the observed vs predicted plots. |
ncol |
Number of columns of the layout to organize the observed vs predicted plots. |
... |
Additional arguments to be passed to |
The function provides visual and quantitative information about the performance of the model:
A layout of the observed and predicted values of (relative tumor volume) vs Time for each SampleID (i.e. subject),
with the actual measurements, the regression line for each SampleID, and the marginal, treatment-specific,
regression line for each treatment group.
Performance metrics of the model obtain calculated using performance::model_performance(). The maximum likelihood-based Akaike's Information Criterion (AIC),
small sample AIC (AICc), and Bayesian Information Criterion, and the Nakagawa's r-squared
root mean squared error (RMSE) of the residuals, and the standard deviation of the residuals (sigma) are provided.
Performance metrics of the model obtain calculated using performance::model_performance() and a layout of plots of the observed vs predicted values for each SampleID.
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
Lüdecke et al., (2021). performance: An R Package for Assessment, Comparison and Testing of Statistical Models. Journal of Open Source Software, 6(60), 3139. https://doi.org/10.21105/joss.03139
Sakamoto, Y., M. Ishiguro, and G. Kitagawa. 1984. Akaike Information Criterion Statistics. Mathematics and Its Applications. Reidel.
Nakagawa, Shinichi, and Holger Schielzeth. 2013. A General and Simple Method for Obtaining r2 from Generalized Linear Mixed-effects Models. Methods in Ecology and Evolution 4 (February): 133–42. https://doi.org/10.1111/j.2041-210x.2012.00261.x.
Johnson, Paul C. D. 2014. Extension of Nakagawa & Schielzeth’s r 2 GLMM to Random Slopes Models. Methods in Ecology and Evolution 5 (September): 944–46. https://doi.org/10.1111/2041-210X.12225.
Nakagawa, Shinichi, Paul C. D. Johnson, and Holger Schielzeth. 2017. The Coefficient of Determination r2 and Intra-Class Correlation Coefficient from Generalized Linear Mixed-Effects Models Revisited and Expanded. Journal of The Royal Society Interface 14 (September): 20170213. https://doi.org/10.1098/rsif.2017.0213.
# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Obtain Observed vs Predicted plots, and model performance metrics ObsvsPred(model = lmm, nrow = 4, ncol = 8)# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Obtain Observed vs Predicted plots, and model performance metrics ObsvsPred(model = lmm, nrow = 4, ncol = 8)
Vizualization of tumor growth data and linear mixed model fitted regression line for the fixed effects. This functions returns a ggplot2 plot, allowing for further personalization.
plot_lmmModel( model, trt_control = "Control", drug_a = "Drug_A", drug_b = "Drug_B", drug_c = NA, combination = "Combination" )plot_lmmModel( model, trt_control = "Control", drug_a = "Drug_A", drug_b = "Drug_B", drug_c = NA, combination = "Combination" )
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
trt_control |
String indicating the name assigned to the 'Control' group. |
drug_a |
String indicating the name assigned to the 'Drug A' group. |
drug_b |
String indicating the name assigned to the 'Drug B' group. |
drug_c |
String indicating the name assigned to the 'Drug C' group (if present). |
combination |
String indicating the name assigned to the Combination ('Drug A' + 'Drug B', or 'Drug A' + 'Drug B' + 'Drug C') group. |
A ggplot2 plot (see ggplot2::ggplot() for more details) showing the tumor growth data represented as log(relative tumor volume) versus time since treatment initiation.
The regression lines corresponding to the fixed effects for each treatment group are also plotted.
data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Default plot plot_lmmModel(lmm, trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Adding ggplot2 elements plot_lmmModel(lmm, trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) + ggplot2::labs(title = "Example Plot") + ggplot2::theme(legend.position = "top")data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Default plot plot_lmmModel(lmm, trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Adding ggplot2 elements plot_lmmModel(lmm, trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) + ggplot2::labs(title = "Example Plot") + ggplot2::theme(legend.position = "top")
Visualization of synergy results obtained by lmmSynergy(). This functions returns a ggplot2 plot, allowing for
further personalization.
plot_lmmSynergy(syn_data)plot_lmmSynergy(syn_data)
syn_data |
Object obtained by |
plot_lmmSynergy produces a ggplot2 plot with the results of the synergy calculation. Each dot represents the estimated combination index
or synergy score, and the gray lines represent the 95% confidence intervals, for each day. Each dot is colored based on the , with
purple colors indicating a , and green colors indicating a .
A list with ggplot2 plots (see ggplot2::ggplot() for more details) with the combination index (CI) and synergy score (SS)
estimates, confidence intervals and p-values for the synergy calculation using linear mixed models.
data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Obtain synergy results lmmSyn <- lmmSynergy(lmm) # Plot synergy results plot_lmmSynergy(lmmSyn) # Accessing to the combination index plot plot_lmmSynergy(lmmSyn)$CI # Accessing to only synergy score plot plot_lmmSynergy(lmmSyn)$SS # Accessing to the grid of both plots side by side plot_lmmSynergy(lmmSyn)$CI_SS # Adding ggplot2 elements plot_lmmSynergy(lmmSyn)$CI + ggplot2::labs(title = "Synergy Calculation for Bliss") + ggplot2::theme(legend.position = "top")data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Obtain synergy results lmmSyn <- lmmSynergy(lmm) # Plot synergy results plot_lmmSynergy(lmmSyn) # Accessing to the combination index plot plot_lmmSynergy(lmmSyn)$CI # Accessing to only synergy score plot plot_lmmSynergy(lmmSyn)$SS # Accessing to the grid of both plots side by side plot_lmmSynergy(lmmSyn)$CI_SS # Adding ggplot2 elements plot_lmmSynergy(lmmSyn)$CI + ggplot2::labs(title = "Synergy Calculation for Bliss") + ggplot2::theme(legend.position = "top")
Visualization of observed vs predicted values by a fitted linear mixed model of tumor growth data.
plot_ObsvsPred(model, nrow = 4, ncol = 5)plot_ObsvsPred(model, nrow = 4, ncol = 5)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
nrow |
Number of rows of the layout to organize the observed vs predicted plots. |
ncol |
Number of columns of the layout to organize the observed vs predicted plots. |
A layout (arranged in nrow rows and ncol columns) of the observed and predicted values of (relative tumor volume) vs Time for each SampleID (i.e. subject),
with the actual measurements, the regression line for each SampleID, and the marginal, treatment-specific,
regression line for each treatment group.
#' data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Obtain the plots plot_ObsvsPred(lmm, nrow = 4, ncol = 8)#' data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Obtain the plots plot_ObsvsPred(lmm, nrow = 4, ncol = 8)
Visualization of random effects diagnostics for a fitted linear mixed model of tumor growth data.
plot_ranefDiagnostics(model)plot_ranefDiagnostics(model)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
A list with different plots for evaluating the normality and homoscedasticity of the random effects, including:
A normal Q-Q plot of the random effects of the model.
A normal Q-Q plot of the residuals by sample.
Boxplots of the "raw" residuals (observed - fitted) by sample.
Scatter plots of the normalized residuals (standardized residuals pre-multiplied by the inverse square-root factor of the estimated error correlation matrix, see nlme::residuals.lme)
vs fitted values by sample. Observations with absolute standardized (normalized) residuals greater than the quantile of the standard normal distribution
are identified in the plots labelled with the time point corresponding to the observation.
data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Generate plots plot_ranefDiagnostics(lmm) # Access to specific plots plot_ranefDiagnostics(lmm)[[1]] plot_ranefDiagnostics(lmm)[[2]]data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Generate plots plot_ranefDiagnostics(lmm) # Access to specific plots plot_ranefDiagnostics(lmm)[[1]] plot_ranefDiagnostics(lmm)[[2]]
Visualization of residuals diagnostics for a fitted linear mixed model of tumor growth data.
plot_residDiagnostics(model)plot_residDiagnostics(model)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
A list with different plots for evaluating the normality and homoscedasticity of the normalized residuals (standardized residuals pre-multiplied by the inverse square-root factor of the estimated error correlation matrix, see nlme::residuals.lme), including:
A normal Q-Q plot of the normalized residuals of the model.
A normal Q-Q plot of the normalized residuals of the model by Time.
A normal Q-Q plot of the normalized residuals of the model by Treatment.
A dotplot of normalized residuals vs fitted values.
A dotplot of the normalized residuals by Time and Treatment.
data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Generate plots plot_residDiagnostics(lmm) # Access to specific plots plot_residDiagnostics(lmm)[[1]] plot_residDiagnostics(lmm)[[2]]data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Generate plots plot_residDiagnostics(lmm) # Access to specific plots plot_residDiagnostics(lmm)[[1]] plot_residDiagnostics(lmm)[[2]]
Generic function to generate different types of plots based on SynergyLMM outputs.
plot_SynergyLMM(object, plot_type = "lmmModel", ...)plot_SynergyLMM(object, plot_type = "lmmModel", ...)
object |
An object generated by |
plot_type |
String indicating the type of plot to generate. Possible options include:
|
... |
Additional arguments passed to the specific plot function. |
Different output plots are produced depending on the plot_type selected.
plot_lmmModel(), plot_lmmSynergy(), plot_ObsvsPred(), plot_ranefDiagnostics(), plot_residDiagnostics().
data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Plot lmmModel plot_SynergyLMM(lmm, plot_type = "lmmModel", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Plot ObsvsPred plot_SynergyLMM(lmm, plot_type = "ObsvsPred") # Plot ranefDiagnostics plot_SynergyLMM(lmm, plot_type = "ranefDiagnostics") # Plot residDiagnostics plot_SynergyLMM(lmm, plot_type = "residDiagnostics") # Plot lmmSynergy lmmSyn <- lmmSynergy(lmm) plot_SynergyLMM(lmmSyn, plot_type = "lmmSynergy")data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination", show_plot = FALSE ) # Plot lmmModel plot_SynergyLMM(lmm, plot_type = "lmmModel", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Plot ObsvsPred plot_SynergyLMM(lmm, plot_type = "ObsvsPred") # Plot ranefDiagnostics plot_SynergyLMM(lmm, plot_type = "ranefDiagnostics") # Plot residDiagnostics plot_SynergyLMM(lmm, plot_type = "residDiagnostics") # Plot lmmSynergy lmmSyn <- lmmSynergy(lmm) plot_SynergyLMM(lmmSyn, plot_type = "lmmSynergy")
PostHocPwr allows for the post hoc power analysis of the synergy hypothesis testing for Bliss and HSA refence
models for a given tumor growth data fitted model.
PostHocPwr(model, nsim = 1000, method = "Bliss", pvalue = 0.05, time = NA, ...)PostHocPwr(model, nsim = 1000, method = "Bliss", pvalue = 0.05, time = NA, ...)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
nsim |
Number of simulations to perform. |
method |
String indicating the method for synergy calculation. Possible methods are "Bliss" and "HSA", corresponding to Bliss and highest single agent, respectively. |
pvalue |
Threshold for the p-value of synergy calculation to be considered statistically significant. |
time |
Time point for which to calculate the statistical power. If not specified, the last time point is used by default. |
... |
Additional parameters to be passed to nlmeU::simulateY: |
The post hoc power calculation relies on simulation of the dependent variable, using nlmeU::simulateY.
For a given fitted model of the tumor growth data, nsim simulations of the dependent variable ()
are done, based on the marginal distribution implied by the fitted model.
The model is then fitted to the new values of the dependant variable.
For each simulation, the new estimates from each model are then used for the synergy hypothesis testing as explained in lmmSynergy, and the p-values stored.
The power is returned as the proportion of simulations resulting in a significant synergy hypothesis testing
(p-value < pvalue).
When time is specified, PostHocPwr refits the model using the data from the time_start time point defined in lmmModel() until time, and report the
statistical power for that model. If time is not specified, the model fitted using all data points is used for statistical power calculation.
Returns a numeric value of the power for the synergy calculation for the model using the method specified in method.
The power is expressed as the proportion of simulations that provides a p-value below the threshold specified in pvalue.
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
APrioriPwr(), PwrSampleSize(), PwrTime().
#' data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) PostHocPwr(lmm, nsim = 50) # 50 simulations for shorter computing time # Using a seed to obtain reproducible results PostHocPwr(lmm, seed = 123, nsim = 50) # Calculating the power for an specific day PostHocPwr(lmm, nsim = 50, time = 6)#' data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) PostHocPwr(lmm, nsim = 50) # 50 simulations for shorter computing time # Using a seed to obtain reproducible results PostHocPwr(lmm, seed = 123, nsim = 50) # Calculating the power for an specific day PostHocPwr(lmm, nsim = 50, time = 6)
This function is generic; method functions can be written to handle specific classes of objects.
Pwr(object, ...)Pwr(object, ...)
object |
an object containing the results returned by a model fitting function (e.g., |
... |
some methods for this generic function may require additional arguments. |
numeric scalar value.
Andrzej Galecki and Tomasz Burzykowski
data("grwth_data") lmm <- lmmModel(data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination") Pwr(lmm, L = c("Time:TreatmentControl" = 1, "Time:TreatmentDrugA" = -1, "Time:TreatmentDrugB" = -1, "Time:TreatmentCombination" = 1))data("grwth_data") lmm <- lmmModel(data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination") Pwr(lmm, L = c("Time:TreatmentControl" = 1, "Time:TreatmentDrugA" = -1, "Time:TreatmentDrugB" = -1, "Time:TreatmentCombination" = 1))
This is method for Pwr() generic function. It is a modified version from the
one described by Galecki and Burzykowski implemented in nlmeU package (nlmeU::Pwr.lme).
## S3 method for class 'lme' Pwr( object, ..., type = c("sequential", "marginal"), Terms, L, verbose = FALSE, sigma, ddf = numeric(0), alpha = 0.05, altB = NULL, tol = 1e-10 )## S3 method for class 'lme' Pwr( object, ..., type = c("sequential", "marginal"), Terms, L, verbose = FALSE, sigma, ddf = numeric(0), alpha = 0.05, altB = NULL, tol = 1e-10 )
object |
an object containing |
... |
some additional arguments may be required. |
type |
an optional character string specifying the type of sum of squares to be used in F-tests
needed for power calculations. Syntax is the same as for |
Terms |
an optional integer or character vector specifying which terms
in the model should be jointly tested to be zero using a Wald F-test. See
|
L |
an optional numeric vector or array specifying linear combinations
of the coefficients in the model that should be tested to be zero. See
|
verbose |
an optional logical value. See |
sigma |
numeric scalar value. |
ddf |
numeric scalar value. Argument can be used to redefine default number of denominator degrees of freedom |
alpha |
numeric scalar value. By default 0.05. |
altB |
matrix/vector containing alternative values for beta parameters |
tol |
numeric scalar value. |
a data frame inheriting from class Pwr.lme
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
nlme::anova.lme, nlmeU::Pwr.lme
A priori power calculation for a hypothetical two-drugs combination study of synergy evaluation using linear-mixed models depending on the sample size per group.
PwrSampleSize( npg = c(5, 8, 10), time = c(0, 3, 5, 10), grwrControl = 0.08, grwrA = 0.07, grwrB = 0.06, grwrComb = 0.03, sd_ranef = 0.01, sgma = 0.1, method = "Bliss", ... )PwrSampleSize( npg = c(5, 8, 10), time = c(0, 3, 5, 10), grwrControl = 0.08, grwrA = 0.07, grwrB = 0.06, grwrComb = 0.03, sd_ranef = 0.01, sgma = 0.1, method = "Bliss", ... )
npg |
A vector with the sample size (number of subjects) per group to calculate the power of the synergy analysis. |
time |
Vector with the times at which the tumor volume measurements have been performed. |
grwrControl |
Coefficient for Control treatment group tumor growth rate. |
grwrA |
Coefficient for Drug A treatment group tumor growth rate. |
grwrB |
Coefficient for Drug B treatment group tumor growth rate. |
grwrComb |
Coefficient for Combination (Drug A + Drug B) treatment group tumor growth rate. |
sd_ranef |
Random effects standard deviation for the model. |
sgma |
Residuals standard deviation for the model. |
method |
String indicating the method for synergy calculation. Possible methods are "Bliss" and "HSA", corresponding to Bliss and highest single agent, respectively. |
... |
Additional parameters to be passed to nlmeU::Pwr.lme method. |
PwrSampleSize allows the user to define an hypothetical drug combination study, customizing several
experimental parameters, such as the sample size, time of measurements, or drug effect,
for the power evaluation of synergy for Bliss and HSA reference models. The power calculation is
based on F-tests of the fixed effects of the model as previously described (Helms, R. W. (1992),
Verbeke and Molenberghs (2009), Gałecki and Burzykowski (2013)).
The focus the power analysis with PwrSampleSize is on the sample size per group. The function allows
for the evaluation of how the statistical power changes when the sample size per group varies while the
other parameters are kept constant. For other a priori power analysis see also APrioriPwr() and PwrTime().
time, grwrControl, grwrA, grwrB, grwrComb, sd_ranef and sgma are parameters referring to
the initial exemplary data set and corresponding fitted model. These values can be obtained from a fitted model, using lmmModel_estimates(),
or be defined by the user.
npg is a vector indicating the different sample sizes for which the statistical power is going to be evaluated, keeping the
rest of parameters fixed.
The functions returns two plots:
A plot representing the hypothetical data, with the regression lines for each
treatment group according to grwrControl, grwrA, grwrB and grwrComb values. The values
assigned to sd_ranef and sgma are also shown.
A plot showing the values of the power calculation depending on the values assigned to
npg.
The function also returns the data frame with the power for the analysis for each sample size
specified in npg.
Helms, R. W. (1992). Intentionally incomplete longitudinal designs: I. Methodology and comparison of some full span designs. Statistics in Medicine, 11(14–15), 1889–1913. https://doi.org/10.1002/sim.4780111411
Verbeke, G. & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer New York. https://doi.org/10.1007/978-1-4419-0300-6
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
PostHocPwr, APrioriPwr(), PwrTime().
PwrSampleSize(npg = 1:20)PwrSampleSize(npg = 1:20)
A priori power calculation for a hypothetical two-drugs combination study of synergy depending on the time of follow-up or the frequency of measurements.
PwrTime( npg = 5, time = list(seq(0, 9, 3), seq(0, 21, 3), seq(0, 30, 3)), type = "max", grwrControl = 0.08, grwrA = 0.07, grwrB = 0.06, grwrComb = 0.03, sd_ranef = 0.01, sgma = 0.1, method = "Bliss", ... )PwrTime( npg = 5, time = list(seq(0, 9, 3), seq(0, 21, 3), seq(0, 30, 3)), type = "max", grwrControl = 0.08, grwrA = 0.07, grwrB = 0.06, grwrComb = 0.03, sd_ranef = 0.01, sgma = 0.1, method = "Bliss", ... )
npg |
Number of mouse per group. |
time |
A list in which each element is a vector with the times at which the tumor volume measurements have been performed.
If |
type |
String indicating whether to calculate the power depending on the time of follow-up ("max"), or the frequency of measurements ("freq"). |
grwrControl |
Coefficient for Control treatment group tumor growth rate. |
grwrA |
Coefficient for Drug A treatment group tumor growth rate. |
grwrB |
Coefficient for Drug B treatment group tumor growth rate. |
grwrComb |
Coefficient for Combination (Drug A + Drug B) treatment group tumor growth rate. |
sd_ranef |
Random effects standard deviation for the model. |
sgma |
Residuals standard deviation for the model. |
method |
String indicating the method for synergy calculation. Possible methods are "Bliss" and "HSA", corresponding to Bliss and highest single agent, respectively. |
... |
Additional parameters to be passed to nlmeU::Pwr.lme method. |
PwrTime allows the user to define an hypothetical drug combination study, customizing several
experimental parameters, such as the sample size, time of measurements, or drug effect,
for the power evaluation of synergy for Bliss and HSA reference models. The power calculation is
based on F-tests of the fixed effects of the model as previously described (Helms, R. W. (1992),
Verbeke and Molenberghs (2009), Gałecki and Burzykowski (2013)).
The focus the power analysis with PwrTime is on the time at which the measurements are done. The function allows
for the evaluation of how the statistical power changes when the time of follow-up varies while the frequency
of measurements is keep constant. It also allows to how the statistical power changes when the time of follow-up is
kept constant, but the frequency of measurements varies.
For other a priori power analysis see also APrioriPwr() and PwrSampleSize().
npg, grwrControl, grwrA, grwrB, grwrComb, sd_ranef and sgma are parameters referring to
the initial exemplary data set and corresponding fitted model. These values can be obtained from a fitted model, using lmmModel_estimates(),
or be defined by the user.
time is a list in which each element is a vector with the times at which the tumor volume measurements have been performed, and for
which the statistical power is going to be evaluated, keeping the rest of parameters fixed.
The functions returns two plots:
A plot representing the hypothetical data, with the regression lines for each
treatment group according to grwrControl, grwrA, grwrB and grwrComb values. The values
assigned to sd_ranef and sgma are also shown.
A plot showing the values of the power calculation depending on the values assigned to
Time. If type is set to "max", the plot shows how the power varies depending on the maximum time of follow-up.
If type is set to "freq", the plot shows how the power varies depending on how frequently the measurements have
been performed.
The function also returns the data frame with the power for the analysis for each value specified in Time.
Helms, R. W. (1992). Intentionally incomplete longitudinal designs: I. Methodology and comparison of some full span designs. Statistics in Medicine, 11(14–15), 1889–1913. https://doi.org/10.1002/sim.4780111411
Verbeke, G. & Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. Springer New York. https://doi.org/10.1007/978-1-4419-0300-6
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
PostHocPwr, APrioriPwr(), PwrSampleSize().
# Power analysis maintaining the frequency of measurements # and varying the time of follow-up ('type = "max"') PwrTime(time = list(seq(0, 9, 3), seq(0, 12, 3), seq(0, 15, 3), seq(0, 21, 3), seq(0, 30, 3)), type = "max") # Power analysis maintaining the time of follow-up # and varying the frequency of measurements ('type = "freq"') PwrTime(time = list(seq(0, 10, 1), seq(0, 10, 2), seq(0, 10, 5), seq(0, 10, 10)), type = "freq")# Power analysis maintaining the frequency of measurements # and varying the time of follow-up ('type = "max"') PwrTime(time = list(seq(0, 9, 3), seq(0, 12, 3), seq(0, 15, 3), seq(0, 21, 3), seq(0, 30, 3)), type = "max") # Power analysis maintaining the time of follow-up # and varying the frequency of measurements ('type = "freq"') PwrTime(time = list(seq(0, 10, 1), seq(0, 10, 2), seq(0, 10, 5), seq(0, 10, 10)), type = "freq")
ranefDiagnostics provides several plots as well as statistical test for the examination
of the normality of the random effects of the input model.
ranefDiagnostics(model, verbose = TRUE)ranefDiagnostics(model, verbose = TRUE)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
verbose |
Logical indicating if the normality and homoscedasticity tests results should be printed to the console. |
One of the assumptions of the model obtained with lmmModel() (as in any other linear mixed model) is that
the random effects are normally distributed:
For the evaluation of this assumption, ranefDiagnostics provides Q-Q plots of random effects,
together with statistical assessment of their normality using Shapiro-Wilk, D'Agostini and Anderson-Darling normality tests.
Additionally, Q-Q plots of the normalized residuals
(standardized residuals pre-multiplied by the inverse square-root factor of the estimated error correlation matrix, see nlme::residuals.lme)
by sample are provided to allow for the identification of subjects
which could be notably different from the others and be affecting the adequacy of the model.
Additionally, boxplots of the "raw" residuals (observed - fitted) by sample and scatter plots of the normalized residuals versus fitted values by sample
are provided to give information about variability of the residuals by subject and possible outlier observations. Observations with absolute standardized (normalized)
residuals greater than the quantile of the standard normal distribution
are identified in the scatter plots labelled with the time point corresponding to the observation.
A list with different elements for the diagnostics of the random effects are produced:
plots: Different plots for evaluating the normality and homoscedasticity of the random effects are produced.
Normality: List with the results from 3 different test of the normality of the random effects: Shapiro - Wilk normality test,
D'Agostino normality test and Anderson - Darling normality test.
Levene.test: results from Levene homoscedasticity test (car::leveneTest()) of the normalized residuals by SampleID (i.e., by subject).
Fligner.test: results from Fligner homoscedasticity test (stats::fligner.test()) of the normalized residuals by SampleID (i.e., by subject).
Pinheiro JC, Bates DM (2000). Mixed-Effects Models in S and S-PLUS. Springer, New York. doi:10.1007/b98882.
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Run random effects diagnostics ranef_diag <- ranefDiagnostics(lmm) #Access to individual plots ranef_diag$Plots[1] ranef_diag$Plots[2] # Access to normality tests ranef_diag$Normality ranef_diag$Normality$Shapiro.test # Access to homoscedasticity tests of residuals by subject ranef_diag$Levene.test ranef_diag$Fligner.test# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Run random effects diagnostics ranef_diag <- ranefDiagnostics(lmm) #Access to individual plots ranef_diag$Plots[1] ranef_diag$Plots[2] # Access to normality tests ranef_diag$Normality ranef_diag$Normality$Shapiro.test # Access to homoscedasticity tests of residuals by subject ranef_diag$Levene.test ranef_diag$Fligner.test
residDiagnostics provides several plots as well as statistical test for the examination
of the normality and homoscedasticity of the residuals of the input model.
residDiagnostics(model, pvalue = 0.05, verbose = TRUE)residDiagnostics(model, pvalue = 0.05, verbose = TRUE)
model |
An object of class "lme" representing the linear mixed-effects model fitted by |
pvalue |
Threshold for the p-value of outlier observations based on their normalized residuals. |
verbose |
Logical indicating if the normality and homoscedasticity tests results, and the list of potential outlier observations should be printed to the console. |
One of the assumption of the model fit by lmmModel() is that the residuals are normally distributed.
For the evaluation of this assumption, residDiagnostics provides Q-Q plots of the normalized residuals
(standardized residuals pre-multiplied by the inverse square-root factor of the estimated error correlation matrix, see nlme::residuals.lme),
together with statistical assessment of their
normality using Shapiro-Wilk, D'Agostini and Anderson-Darling normality tests. Additionally, Q-Q plots of the normalized residuals by time point and
treatment group are provided to be able to detect time points or treatment groups which could be notably different from the others and be
affecting the adequacy of the model.
Scatter plots of the normalized residuals versus fitted values and normalized residuals per time and per treatment are also provided to give information about variability of the residuals and possible outlier observations. These plots are accompanied by Levene and Fligner-Killend homogeneity of variance test results.
Observations with absolute standardized (normalized) residuals greater than the quantile of the standard normal distribution
are identified and reported as potential outlier observations.
A list with different elements for the diagnostics of the residuals are produced:
plots: Different plots for evaluating the normality and homocedasticity of the residuals.
outliers: Data frame with the identified outliers based on the Pearson residuals and the value of pval. The column resid.p contains the
value of the Pearson residuals for each observation.
Normality: List with the results from 3 different test of the normality of the normalized residuals of the model: Shapiro - Wilk normality test,
D'Agostino normality test and Anderson - Darling normality test.
Levene.test: List with the Levene homoscedasticity test results of the normalized residuals by Time and Treatment.
Fligner.test: List with the Fligner-Killeen homoscedasticity test results of the normalized residuals by Time and Treatment.
Pinheiro JC, Bates DM (2000). Mixed-Effects Models in S and S-PLUS. Springer, New York. doi:10.1007/b98882.
Andrzej Galecki & Tomasz Burzykowski (2013) Linear Mixed-Effects Models Using R: A Step-by-Step Approach First Edition. Springer, New York. ISBN 978-1-4614-3899-1
# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Residuals diagnostics resid_diag <- residDiagnostics(model = lmm, pvalue = 0.05) # Access outliers data frame resid_diag$Outliers # Access individual plots resid_diag$Plots[1] resid_diag$Plots[2] # Access results of normality tests resid_diag$Normality resid_diag$Normality$Shapiro.test # Access to homoscedasticity test results resid_diag$Levene.test resid_diag$Fligner.test# Load the example data data(grwth_data) # Fit the model lmm <- lmmModel( data = grwth_data, sample_id = "subject", time = "Time", treatment = "Treatment", tumor_vol = "TumorVolume", trt_control = "Control", drug_a = "DrugA", drug_b = "DrugB", combination = "Combination" ) # Residuals diagnostics resid_diag <- residDiagnostics(model = lmm, pvalue = 0.05) # Access outliers data frame resid_diag$Outliers # Access individual plots resid_diag$Plots[1] resid_diag$Plots[2] # Access results of normality tests resid_diag$Normality resid_diag$Normality$Shapiro.test # Access to homoscedasticity test results resid_diag$Levene.test resid_diag$Fligner.test
Helper function to simulate tumor growth data for a two-drug combination experiment.
simulateTumorGrowth( npg = 5, timepoints = c(0, 3, 5, 10), initial_volume = 100, grwrControl = 0.08, grwrA = 0.07, grwrB = 0.06, grwrComb = 0.04, sd = 0.1 )simulateTumorGrowth( npg = 5, timepoints = c(0, 3, 5, 10), initial_volume = 100, grwrControl = 0.08, grwrA = 0.07, grwrB = 0.06, grwrComb = 0.04, sd = 0.1 )
npg |
Number of samples per group. |
timepoints |
Vector with the time points at which the tumor volume measurements have been performed. |
initial_volume |
Initial volume for the tumor growth. |
grwrControl |
Coefficient for Control treatment group tumor growth rate. |
grwrA |
Coefficient for Drug A treatment group tumor growth rate. |
grwrB |
Coefficient for Drug B treatment group tumor growth rate. |
grwrComb |
Coefficient for Combination (Drug A + Drug B) treatment group tumor growth rate. |
sd |
Variability for the tumor growth. |
The function simulates the tumor growth following exponential kinetics, given by
where is given by initial_volume, is given by timepoints
and are the coefficients given by grwrControl, grwrA, grwrB, and grwrComb.
The variability is simulated using the sd argument to add random noise from a normal distribution ,
with = sd.
A data frame with tumor growth data in long format.
# This code generates the 'grwth_data' example dataset: set.seed(123) grwth_data <- simulateTumorGrowth(npg = 8, timepoints = seq(0,30,3), initial_volume = 200, grwrControl = 0.08, grwrA = 0.07, grwrB = 0.065, grwrComb = 0.04, sd = 0.2)# This code generates the 'grwth_data' example dataset: set.seed(123) grwth_data <- simulateTumorGrowth(npg = 8, timepoints = seq(0,30,3), initial_volume = 200, grwrControl = 0.08, grwrA = 0.07, grwrB = 0.065, grwrComb = 0.04, sd = 0.2)