synergylmm

A comprehensive statistical framework for designing and analyzing in vivo drug combination experiments.

https://github.com/rafromb/synergylmm

Repository

A comprehensive statistical framework for designing and analyzing in vivo drug combination experiments.

README.Rmd

---
output: github_document
---



```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  out.width = "100%"
)
```

# SynergyLMM  




[![CRAN status](https://www.r-pkg.org/badges/version/SynergyLMM)](https://CRAN.R-project.org/package=SynergyLMM)
[![Lifecycle: stable](https://img.shields.io/badge/lifecycle-stable-brightgreen.svg)](https://lifecycle.r-lib.org/articles/stages.html#stable)


A comprehensive statistical framework for designing and analyzing _in vivo_ drug combination experiments. 

## Table of Contents

- [Installation](#installation)
- [1. Fit Model](#1-fit-model)
- [1.1 Gompertz Growth Model](#11-gompertz-growth-model)
- [1.2 Specify Residual Variance Structure](#12-specify-residual-variance-structure)
- [1.3 Adjust Optimization Controls](#13-adjust-optimization-controls-to-solve-convergence-problems)
- [2. Synergy Analysis](#2-synergy-analysis)
- [3. Model Diagnosis](#3-model-diagnosis)
- [4. Power Analysis](#4-power-analysis)

## Installation

You can install the development version of SynergyLMM from [GitHub](https://github.com/) with:

``` r
# install.packages("pak")
pak::pak("RafRomB/SynergyLMM")
```

Or you can install the [CRAN-released version](https://cran.r-project.org/package=SynergyLMM) with:

``` r
install.packages("SynergyLMM")
```

## Wep-App

You can also use **SynergyLMM** directly in your browser at: https://synergylmm.uiocloud.no/ 

## Example Use of SynergyLMM

This is a basic example which shows how to use SynergyLMM to analyze synergy in a 2-drug combination _in vivo_ experiment.

```{r load_SynergyLMM}
library(SynergyLMM)
```

We start by loading the data (in long format). We will use the example data provided in the package:

```{r load_data}
data("grwth_data")
```

### 1. Fit Model

The first step is fitting the model from our data:

```{r fit_model, fig.width=8}
# Most simple model
lmm <- lmmModel(
  data = grwth_data,
  grwth_model = "exp",
  sample_id = "subject",
  time = "Time",
  treatment = "Treatment",
  tumor_vol = "TumorVolume",
  trt_control = "Control",
  drug_a = "DrugA",
  drug_b = "DrugB",
  combination = "Combination"
)
```

#### 1.1 Gompertz Growth Model
The previous example fits the model using an exponential growth. Users can also use a Gompertz growth model using `grwth_model = "gompertz"`:


```{r fit_model_gomp, fig.width=8}
lmm_gomp <- lmmModel(
  data = grwth_data,
  grwth_model = "gompertz", # Specify Gompertz growth model
  start_values = "selfStart",
  sample_id = "subject",
  time = "Time",
  treatment = "Treatment",
  tumor_vol = "TumorVolume",
  trt_control = "Control",
  drug_a = "DrugA",
  drug_b = "DrugB",
  combination = "Combination"
)
```

When using the Gompertz growth model, setting `start_values = "selfStart"` derives the initial values from a call to `stats::nls`, which can facilitate successful model fitting.

#### 1.2 Specify Residual Variance Structure

To address heteroscedasticity, users can define a variance structure using the argument `weights` and functions from the [`nlme`](https://cran.r-project.org/web/packages/nlme/index.html) R package. For example, a treatment-specific variance structure can be specify using `weights = nlme::varIdent(form = ~1|Treatment)`:

```{r varIdent}
lmmVar <- lmmModel(
  data = grwth_data,
  grwth_model = "exp",
  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|Treatment),
  show_plot = FALSE
)
```

Other options are:

- Time-specific variance structure: `weights = nlme::varIdent(form = ~1|Time)`
- Subject-specific variance structure: `weights = nlme::varIdent(form = ~1|SampleID)`

Combinations of variance functions can also be specified using `nlme::varComb()`. For example, to model residual variance as varying independently across both individual samples and time points, users can specify:

- `weights = nlme::varComb(nlme::varIdent(form = ~1|SampleID), nlme::varIdent(form = ~1|Time))`

#### 1.3 Adjust Optimization Controls to Solve Convergence Problems

Convergence issues can appear when using complex variance or correlation structures, particularly with non-linear models like Gompertz.

User can customize optimization settings—e.g., number of iterations, convergence tolerance, or optimizer method—via the `control` argument in `lmmModel()` function. 

For example, to increase the maximum number of iterations for the optimization algorithm: `control = nlme::lmeControl(maxIter = 1000)`

```{r lmeControl}
lmmCont <- lmmModel(
  data = grwth_data,
  grwth_model = "exp",
  sample_id = "subject",
  time = "Time",
  treatment = "Treatment",
  tumor_vol = "TumorVolume",
  trt_control = "Control",
  drug_a = "DrugA",
  drug_b = "DrugB",
  combination = "Combination",
  control = nlme::lmeControl(maxIter = 1000),
  show_plot = FALSE
)
```

- For the exponential model use: `control = nlme::lmeControl()`
- For the Gompertz model use: `control = nlme::nlmeControl()`

#### 1.4 Model Estimates

We can obtain the model estimates and their standard deviation (sd) using:

```{r model_estimates}
lmmModel_estimates(lmm)
```

Also for the Gompertz model:

```{r}
lmmModel_estimates(lmm_gomp)
```

### 2. Synergy Analysis

Synergy calculation can be done using `lmmSynergy`. The reference model to assess the drug combination effect can be changed to Bliss independence, highest single agent, or response additivity, by setting `method = "Bliss"`, `method = "HSA"`, and `method = "RA"`, respectively.

**Bliss independence model**

```{r bliss_syn, fig.width=10}
lmmSynergy(lmm, method = "Bliss")
```

**Highest Single Agent model**

```{r hsa_syn, fig.width=10}
lmmSynergy(lmm, method = "HSA")
```

**Response Additivity**

```{r ra_syn, fig.width=10}
set.seed(123)
lmmSynergy(lmm, method = "RA", ra_nsim = 1000)
```

**Using Robust Estimates**

```{r robustSE, fig.width=10}
lmmSynergy(lmm, method = "Bliss", robust = TRUE)
```
- _Note: Sandwich-based robust estimates are only available for the exponential model._


### 3. Model Diagnostics

We can perform the model diagnostics using the following functions:

**Random Effects**

```{r ranef_diag, fig.width=8, fig.height=10}
ranefDiagnostics(lmm)
```
**Residuals Diagnostics**

```{r resid_diagnostics, fig.width=8, fig.height=12}
residDiagnostics(lmm)
```

**Observed versus Predicted Values**

```{r obs_vs_pred, fig.height=8, fig.width=10}
ObsvsPred(lmm)
```

**Influential Diagnostics**

_Cook's distances based on the change of fitted values_

```{r CooksD_fitted}
CookDistance(lmm)
```


_Cook's distances based on the change of fixed effects_

```{r CooksD_fixef}
CookDistance(lmm, type = "fixef")
```

_log likelihood displacements_

- _Note: The calculation of log likelihood displacements is only available for the exponential growth model._

```{r logLikDisp}
logLikSubjectDisplacements(lmm)
```


### 4. Power Analysis

- _Note: The_ post-hoc _power analysis and the_ a priori _power analysis are only available for the exponential growth model._

**Post-Hoc Power Analysis**

```{r posthoc_pwr}
set.seed(123)
PostHocPwr(lmm, method = "Bliss", time = 30)
```

**A Priori Power Analysis**

We will estimate the effect of sample size on statistical power based on the estimates from the model:

```{r pwrsamplesize, fig.width=12}
# Vector with different sample sizes per group
npg <- 3:15

# Obtain model estimates
(lmmestim <- lmmModel_estimates(lmm))

# Obtain time points

(timepoints <- unique(lmm$dt1$Time))

# Calculate power depending on sample size per group

PwrSampleSize(
  npg = npg,
  time = timepoints,
  grwrControl = round(lmmestim$Control,3),
  grwrA = round(lmmestim$DrugA,3),
  grwrB = round(lmmestim$DrugB,3),
  grwrComb = round(lmmestim$Combination,3),
  sd_ranef = round(lmmestim$sd_ranef,3),
  sgma = round(lmmestim$sd_resid,3),
  method = "Bliss"
)
```

Owner

GitHub Events

Issues and Pull Requests

Last synced: 10 months ago

Packages