Neuenschwander, Branson, and Gsponer (2008) (NBG) introduced a derivative of the CRM for dose-escalation clinical trials using the model:

\[ \text{logit} p_i = \alpha + \exp{(\beta)} \log{(x_i / d^*)}, \]

where \(p_i\) is the probability of toxicity at the \(i\)th dose, \(x_i\), and \(d^*\) is a reference dose. Here \(\alpha\) and \(\beta\) are model parameters on which the authors place a bivariate normal prior. This model is very similar to the two-parameter logistic CRM, implemented with stan_crm(model = 'logistic2'). However, a notable difference is that the dose, \(x_i\), enters the model as a covariate. This dispenses with the toxicity skeleton that is used in the CRM.

The authors introduce their approach in a paper that argues for greater parameterisation in dose-escalation models so that they might more accurately estimate the entire dose-toxicity curve than simplistic one-parameter models.

Usage

Let’s run some examples.

To access the Stan implementations of the NBG model, we must load trialr:

library(trialr)

For illustration, let us reproduce the model fit in the lower right pane in Figure 1 of Neuenschwander, Branson, and Gsponer (2008). The authors fit their model to a partially-completed historic trial investigating 15 doses, using the highest dose as the reference dose:

dose <- c(1, 2.5, 5, 10, 15, 20, 25, 30, 40, 50, 75, 100, 150, 200, 250)
d_star <- 250

The original investigators sought a dose associated with 30% toxicity:

target <- 0.30

and the analysis concerns which dose to give after observing the following outcomes:

outcomes <- '1NNN 2NNNN 3NNNN 4NNNN 7TT'

We see that the original investigators escalated through doses 1 to 3, but after having observed no toxicity, seem to have thrown caution to the wind and escalated straight to dose 7. Unfortunately, that move seems to have be imprudent because two toxicities were seen in two patients.

Neuenschwander et al. introduce normal priors for the parameters \(\alpha\) and \(\beta\) which we specify when we call stan_nbg to fit the model:

fit <- stan_nbg(outcome_str = outcomes, real_doses = dose, d_star = d_star,
                target = target, alpha_mean = 2.15, alpha_sd = 0.84,
                beta_mean = 0.52, beta_sd = 0.8, seed = 2020, refresh = 0)

Note that presently, the implementation in trialr uses independent normal priors on \(\alpha\) and \(\beta\), and not the bivariate normal that the authors used. This will hopefully be addressed in a future release of trialr. Nevertheless, we see that the small difference of prior apparently makes little difference to the posterior because the mean estimates of the probability of toxicity are very close to that shown in Figure 1 of NBG’s manuscript:

fit
#>    Patient Dose Toxicity Weight
#> 1        1    1        0      1
#> 2        2    1        0      1
#> 3        3    1        0      1
#> 4        4    2        0      1
#> 5        5    2        0      1
#> 6        6    2        0      1
#> 7        7    2        0      1
#> 8        8    3        0      1
#> 9        9    3        0      1
#> 10      10    3        0      1
#> 11      11    3        0      1
#> 12      12    4        0      1
#> 13      13    4        0      1
#> 14      14    4        0      1
#> 15      15    4        0      1
#> 16      16    7        1      1
#> 17      17    7        1      1
#> 
#>    Dose N Tox ProbTox MedianProbTox ProbMTD
#> 1     1 3   0  0.0126       0.00563 0.00025
#> 2     2 4   0  0.0325       0.01967 0.00050
#> 3     3 4   0  0.0675       0.04934 0.02250
#> 4     4 4   0  0.1385       0.11863 0.10450
#> 5     5 0   0  0.2064       0.18937 0.16775
#> 6     6 0   0  0.2692       0.25641 0.17100
#> 7     7 2   2  0.3264       0.31878 0.14425
#> 8     8 0   0  0.3780       0.37409 0.16850
#> 9     9 0   0  0.4661       0.46884 0.12475
#> 10   10 0   0  0.5371       0.54698 0.07125
#> 11   11 0   0  0.6616       0.67700 0.02000
#> 12   12 0   0  0.7390       0.75753 0.00375
#> 13   13 0   0  0.8259       0.84691 0.00075
#> 14   14 0   0  0.8717       0.89297 0.00025
#> 15   15 0   0  0.8993       0.91986 0.00000
#> 
#> The model targets a toxicity level of 0.3.
#> The dose with estimated toxicity probability closest to target is 7.
#> The dose most likely to be the MTD is 6.
#> Model entropy: 2.06

We see that the design advocates selecting dose 7 for the next patients. This can be verified in the manuscript.

To illustrate the researchers’ motivation for developing their method, the original one-parameter CRM model using the original investigators’ skeleton was:

skeleton = c(0.01, 0.015, 0.02, 0.025, 0.03, 0.04, 0.05, 0.10, 0.17, 0.30)
fit2 <- stan_crm(outcomes, skeleton = skeleton, target = target,
                 model = 'empiric', beta_sd = 1.34, seed = 2020, refresh = 0)
fit2
#>    Patient Dose Toxicity Weight
#> 1        1    1        0      1
#> 2        2    1        0      1
#> 3        3    1        0      1
#> 4        4    2        0      1
#> 5        5    2        0      1
#> 6        6    2        0      1
#> 7        7    2        0      1
#> 8        8    3        0      1
#> 9        9    3        0      1
#> 10      10    3        0      1
#> 11      11    3        0      1
#> 12      12    4        0      1
#> 13      13    4        0      1
#> 14      14    4        0      1
#> 15      15    4        0      1
#> 16      16    7        1      1
#> 17      17    7        1      1
#> 
#>    Dose Skeleton N Tox ProbTox MedianProbTox ProbMTD
#> 1     1    0.010 3   0  0.0754        0.0592 0.00725
#> 2     2    0.015 4   0  0.0924        0.0759 0.00600
#> 3     3    0.020 4   0  0.1071        0.0906 0.01125
#> 4     4    0.025 4   0  0.1201        0.1038 0.01150
#> 5     5    0.030 0   0  0.1321        0.1161 0.02425
#> 6     6    0.040 0   0  0.1538        0.1386 0.02650
#> 7     7    0.050 2   2  0.1733        0.1589 0.10250
#> 8     8    0.100 0   0  0.2533        0.2432 0.27850
#> 9     9    0.170 0   0  0.3419        0.3369 0.36500
#> 10   10    0.300 0   0  0.4763        0.4775 0.16725
#> 
#> The model targets a toxicity level of 0.3.
#> The dose with estimated toxicity probability closest to target is 9.
#> The dose most likely to be the MTD is 9.
#> Model entropy: 1.61

Note that the above fit uses just the lowest ten doses, to be consistent with Table 1 in Neuenschwander, Branson, and Gsponer (2008). Incredibly, this design advocates escalating two more doses to dose 9, despite the outcomes 7TT in the previous cohort.

To be fair to the CRM, this scenario has become a didactic example for how not to choose a skeleton. The prior probabilities of toxicity are far too close together. Choosing a perhaps more sensible skeleton, the model advocates more defensible behaviour:

skeleton = c(0.03, 0.06, 0.12, 0.20, 0.30, 0.40, 0.50, 0.59, 0.67, 0.74)
# Obtained using dfcrm::getprior(0.05, 0.3, 5, 10)
fit3 <- stan_crm(outcomes, skeleton = skeleton, target = target,
                 model = 'empiric', beta_sd = 1.34, seed = 2020, refresh = 0)
fit3
#>    Patient Dose Toxicity Weight
#> 1        1    1        0      1
#> 2        2    1        0      1
#> 3        3    1        0      1
#> 4        4    2        0      1
#> 5        5    2        0      1
#> 6        6    2        0      1
#> 7        7    2        0      1
#> 8        8    3        0      1
#> 9        9    3        0      1
#> 10      10    3        0      1
#> 11      11    3        0      1
#> 12      12    4        0      1
#> 13      13    4        0      1
#> 14      14    4        0      1
#> 15      15    4        0      1
#> 16      16    7        1      1
#> 17      17    7        1      1
#> 
#>    Dose Skeleton N Tox ProbTox MedianProbTox ProbMTD
#> 1     1     0.03 3   0  0.0177       0.00744 0.00025
#> 2     2     0.06 4   0  0.0344       0.01959 0.00250
#> 3     3     0.12 4   0  0.0700       0.05163 0.01625
#> 4     4     0.20 4   0  0.1229       0.10544 0.08625
#> 5     5     0.30 0   0  0.1977       0.18583 0.19575
#> 6     6     0.40 0   0  0.2820       0.27782 0.26925
#> 7     7     0.50 2   2  0.3759       0.37951 0.24000
#> 8     8     0.59 0   0  0.4689       0.47830 0.13325
#> 9     9     0.67 0   0  0.5584       0.57133 0.04650
#> 10   10     0.74 0   0  0.6421       0.65647 0.01000
#> 
#> The model targets a toxicity level of 0.3.
#> The dose with estimated toxicity probability closest to target is 6.
#> The dose most likely to be the MTD is 6.
#> Model entropy: 1.77

With this particular skeleton that places the prior guess of MTD at dose 5 and spaces out the prior probabilities of toxicity, de-escalation to dose 6 is suggested.

trialr and the escalation package

escalation is an R package that provides a grammar for specifying dose-finding clinical trials. For instance, it is common for trialists to say something like ‘I want to use this published design… but I want it to stop once \(n\) patients have been treated at the recommended dose’ or ‘…but I want to prevent dose skipping’ or ‘…but I want to select dose using a more risk-averse metric than merely closest-to-target’.

trialr and escalation work together to achieve these goals. trialr provides model-fitting capabilities to escalation, including the NBG method described here. escalation then provides additional classes to achieve all of the above custom behaviours, and more.

escalation also provides methods for running simulations and calculating dose-paths. Simulations are regularly used to appraise the operating characteristics of adaptive clinical trial designs. Dose-paths are a tool for analysing and visualising all possible future trial behaviours. Both are provided for a wide array of dose-finding designs, with or without custom behaviours like those identified above. There are many examples in the escalation vignettes at https://cran.r-project.org/package=escalation.

References

Neuenschwander, Beat, Michael Branson, and Thomas Gsponer. 2008. “Critical Aspects of the Bayesian Approach to Phase I Cancer Trials.” Statistics in Medicine 27 (13): 2420–39. https://doi.org/10.1002/sim.3230.