In the general introductory CRM vignette, we introduced the different flavours of the Continual Reassessment Method (CRM) implmented in trialr. In this vignette, we demonstrate some of trialr’s capabilities for analysing dose transition pathways (Yap et al. 2017) with the general CRM model.

## Introduction

A pathway in a dose-finding trial represents a sequence of doses that were delivered to patients and the outcomes that they experienced. trialr provides functions for calculating pathways, analysing model behaviour, and visualising results.

A syntax for succinctly representing pathways will be useful. Brock et al. (2017) introduced a method for describing dose pathways in efficacy and toxicity dose-finding trials using the characters E, T, N & B to represent patients that experienced efficacy only, toxicity only, neither or both. This syntax is used in the EffTox vignette. In the CRM setting, where patients experience only toxicity or no toxicity, we can restrict the character set to simply T and N.

Thus, a dose-finding trial that treated a cohort of three patients at dose-level 1, with none experiencing the dose-limiting toxicity (DLT) event, may be represented as 1NNN. If the next cohort of three was treated at dose-level 2 and the first of these three experienced DLT, we could describe the entire pathway as 1NNN 2TNN. Cohorts are assumed to be separated by spaces. This syntax can be used to describe the outcomes observed hitherto in a trial, the notional future outcomes, or some combination of the two.

We might analyse the model fits on a pathway to learn how model estimation has evolved in a trial. That is the topic of the next section. We might also like to investigate what a CRM model would recommend in each feasible future pathway, as described by Yap et al. (2017). Used in this way, the analysis of dose transition pathways is a valuable step for honing a trial design, or anticipating future trial directions. This is the topic of the latter section.

## Dose pathways with CRM in trialr

Suppose we are mid-way through a dose-finding trial using a CRM design. We have treated and evaluated nine patients at two dose-levels, with one toxicity seen at the second dose:

outcome_str <- '1NNN 2NTN 2NNN'

Seeking a dose associated with a risk of toxicity close to 25%, we commenced the trial with the following dose-toxicity skeleton:

skeleton <- c(0.05, 0.15, 0.25, 0.4, 0.6)
target <- 0.25

The crm_path_analysis function provides a convenient way of fitting the CRM model to each cohort in a pathway. This allows us to see how our model has evolved in its inferences.

library(trialr)

path <- crm_path_analysis(
outcome_str = outcome_str,
skeleton = skeleton, target = target, model = 'empiric',
beta_sd = 1, seed = 123, refresh = 0)

The code above fits four CRM models:

• using the prior beliefs only, before any patients have been treated (i.e. the pathway is the empty string "")
• after outcomes 1NNN
• after outcomes 1NNN 2NTN
• after outcomes 1NNN 2NTN 2NNN

The returned object has type dose_finding_paths, a simple list-like object containing the nodes on the dose pathway. These nodes contain the crm_fit objects that would have been returned by calling stan_crm on the prevailing outcomes. Nodes also contain references to their parent so that the sequence of a pathway is known. This is pertinent below in the analysis of future pathways where many different future outcomes are possible, creating rich tree-like structures.

For all dose_finding_paths objects, the nodes are keyed by the pathway:

names(path)
## [1] ""               "1NNN"           "1NNN 2NTN"      "1NNN 2NTN 2NNN"

We can convert the path object to a tibble for flexible analysis:

library(tibble)

df <- as_tibble(path)
df
## # A tibble: 4 x 8
##   .node .parent .depth outcomes    next_dose fit      parent_fit dose_index
##   <dbl>   <dbl>  <dbl> <chr>           <dbl> <named > <named li> <named li>
## 1     1      NA      0 ""                  2 <crm_fi… <NULL>     <dbl [5]>
## 2     2       1      1 "1NNN"              4 <crm_fi… <crm_fit>  <dbl [5]>
## 3     3       2      2 "1NNN 2NTN"         2 <crm_fi… <crm_fit>  <dbl [5]>
## 4     4       3      3 "1NNN 2NTN…         3 <crm_fi… <crm_fit>  <dbl [5]>

tibble was chosen over the base R data.frame class because of its flexibility. For instance, the fit column contains the crm_fit object associated with the outcomes. The dose_index column contains a vector of integers reflecting the dose-level under study. This is useful if we unnest the data to analyse statistics pertaining to specific doses. For instance, we might be interested in how our beliefs on the rate of toxicity at dose-level 2 have evolved. Using functions from dplyr, tidyr and purrr, we can extract the prob_tox vector from each fit, unnest, and filter to retain only rows that pertain to dose 2.

The behaviour of unnest changed in v1.0 of tidyr. As of version v1.0, the command to do this is:

library(tidyr)
library(purrr)
library(dplyr)

df %>%
mutate(prob_tox = fit %>% map('prob_tox')) %>%
select(outcomes, dose_index, prob_tox) %>%
unnest(cols = c(dose_index, prob_tox)) %>%
filter(dose_index == 2)
## # A tibble: 4 x 3
##   outcomes         dose_index prob_tox
##   <chr>                 <dbl>    <dbl>
## 1 ""                        2    0.226
## 2 "1NNN"                    2    0.128
## 3 "1NNN 2NTN"               2    0.232
## 4 "1NNN 2NTN 2NNN"          2    0.168

For older versions of tidyr, that command is:

library(tidyr)
library(purrr)
library(dplyr)

df %>%
mutate(prob_tox = fit %>% map('prob_tox')) %>%
select(outcomes, dose_index, prob_tox) %>%
unnest %>%
filter(dose_index == 2)

We see that at each update, the toxicity estimate moves in the direction we would expect in response to the outcomes observed.

Having the crm_fit object available allows us to analyse alternative algorithms for choosing dose. For instance, the careful_escalation function will avoid skipping doses in escalation and halt a trial when there is sufficient evidence that a particular dose is too toxic.

df %>%
mutate(
recommended_dose = fit %>% map_int('recommended_dose'),
careful_dose = fit %>% map_dbl(careful_escalation,
tox_threshold = target + 0.1,
certainty_threshold = 0.7)
) %>%
select(outcomes, recommended_dose, careful_dose)
## # A tibble: 4 x 3
##   outcomes         recommended_dose careful_dose
##   <chr>                       <int>        <dbl>
## 1 ""                              2            1
## 2 "1NNN"                          4            2
## 3 "1NNN 2NTN"                     2            2
## 4 "1NNN 2NTN 2NNN"                3            3

We will use dose selection functions like careful_escalation again below when calculating future dose transition pathways. To complete this example, let us observe how careful_escalation will eventually recommended the dose-level NA to signify that the trial should be stopped.

paths <- crm_path_analysis(
outcome_str = '1NNN 2NTN 2NNN 3TTT 1TTT 1TNT',
skeleton = skeleton, target = target, model = 'logistic',
a0 = 3, beta_mean = 0,beta_sd = 1,
seed = 123, refresh = 0)

df <- as_tibble(paths)

df %>%
mutate(
recommended_dose = fit %>% map_int('recommended_dose'),
careful_dose = fit %>% map_dbl(careful_escalation,
tox_threshold = target + 0.1,
certainty_threshold = 0.7)
) %>%
select(outcomes, recommended_dose, careful_dose)
## # A tibble: 7 x 3
##   outcomes                        recommended_dose careful_dose
##   <chr>                                      <int>        <dbl>
## 1 ""                                             1            1
## 2 "1NNN"                                         5            2
## 3 "1NNN 2NTN"                                    2            2
## 4 "1NNN 2NTN 2NNN"                               3            3
## 5 "1NNN 2NTN 2NNN 3TTT"                          1            1
## 6 "1NNN 2NTN 2NNN 3TTT 1TTT"                     1            1
## 7 "1NNN 2NTN 2NNN 3TTT 1TTT 1TNT"                1            1

After these inopportune outcomes, we see that the careful_escalation function now advocates stopping. It does this because there is a greater than 70% posterior probability that the risk of toxicity at the lowest dose-level exceeds the target toxicity rate plus 10%. It stops with reference to toxicity risk at the lowest dose by default. To scrutinise another dose-level, we could have used the reference_dose parameter.

## Dose transition pathways for future cohorts

Statisticians have demonstrated again and again the poor performance of the perennial 3+3 dose-finding design (O’Quigley, Pepe, and Fisher 1990; Iasonos et al. 2008; Le Tourneau, Lee, and Siu 2009). Despite this, the approach has remained stubbornly popular for decades (Rogatko et al. 2007; Chiuzan et al. 2017). Why is this?

It is surely due in part to 3+3’s simplicity. Dose selection decisions are governed by simple rules based on the outcomes of cohorts of three patients. Once familiar with these rules, perfect foresight on the dose selection pathways is possible. For instance 1NNN, will be followed by a cohort at dose 2. 1NNT will be followed by another cohort at dose 1. 1NNT 1NNN will be followed by a cohort at dose 2. 1NNT 1NNN 2NTT will result in the trial being stopped and dose 1 being declared the maximum tolerable dose (MTD). These simple rules can be generalised ad-nauseam to foresee every conceivable trial pathway.

Dose transition pathways (DTPs) are a general tool for furnishing relatively complex statistical dose-finding designs with the same level of foresight and transparency. Yap et al. (2017) introduced DTPs with the CRM as a tool to aid trial design and planning.

The trialr function crm_dtps calculates DTPs for an arbitrary sequence of future cohorts of your choosing. The function is essentially a sequence of calls to stan_crm, with some logic to record the path structure and avoid redundant invocations. Thus, each of the CRM model types supported by stan_crm is supported here.

To calculate paths for two future cohorts of two patients, we use the parameter cohort_sizes = c(2, 2):

paths <- crm_dtps(skeleton = skeleton,
target = target,
model = 'empiric',
cohort_sizes = c(2, 2),
next_dose = 2,
beta_sd = 1,
refresh = 0)

The parameters skeleton, target and model are mandatory because they are passed to stan_crm. Depending on the model type chosen, further parameters will be required. For instance, under the empiric model we must provide a value for the prior standard deviation of $$\beta$$ via the beta_sd parameter. Refer to the introductory CRM vignette for more details.

In the example above, the parameter next_dose determines the dose-level that will be given to the first cohort. If omitted, the first cohort is treated at the dose suggested by the model fit only to the prior information. The option of overriding this is provided because in such experimental medical settings, trialists may choose to start at a dose they firmly believe will be safe, even if their prior belief is that some higher doses will also be tolerable.

Further parameters can be passed ultimately to rstan::sampling via the ellipsis operator to tailor the MCMC sampling. Here, we provide refresh = 0 to suppress sampling messages. We might also have adjusted the number of cores to use for sampling, or the number of warmup samples, for instance. See the rstan documentation for full details.

As with the examples in the previous section, a dose_finding_paths object is returned that may be converted to a tibble for further analysis.

df <- as_tibble(paths)
df
## # A tibble: 13 x 8
##    .node .parent .depth outcomes next_dose fit       parent_fit  dose_index
##    <dbl>   <dbl>  <dbl> <chr>        <dbl> <named l> <named lis> <named li>
##  1     1      NA      0 ""               2 <crm_fit> <NULL>      <dbl [5]>
##  2     2       1      1 "NN"             4 <crm_fit> <crm_fit>   <dbl [5]>
##  3     3       2      2 "NN"             5 <crm_fit> <crm_fit>   <dbl [5]>
##  4     4       1      1 "NT"             1 <crm_fit> <crm_fit>   <dbl [5]>
##  5     5       4      2 "NN"             2 <crm_fit> <crm_fit>   <dbl [5]>
##  6     6       1      1 "TT"             1 <crm_fit> <crm_fit>   <dbl [5]>
##  7     7       6      2 "NN"             1 <crm_fit> <crm_fit>   <dbl [5]>
##  8     8       2      2 "NT"             3 <crm_fit> <crm_fit>   <dbl [5]>
##  9     9       4      2 "NT"             1 <crm_fit> <crm_fit>   <dbl [5]>
## 10    10       6      2 "NT"             1 <crm_fit> <crm_fit>   <dbl [5]>
## 11    11       2      2 "TT"             2 <crm_fit> <crm_fit>   <dbl [5]>
## 12    12       4      2 "TT"             1 <crm_fit> <crm_fit>   <dbl [5]>
## 13    13       6      2 "TT"             1 <crm_fit> <crm_fit>   <dbl [5]>

We see that our analysis of cohort_sizes = c(2, 2) has yielded 13 model fits. There is one node of depth 0. This is the trial starting point, where next_dose was provided manually. There are three nodes of depth 1, corresponding to paths 2NN, 2NT and 2TT. Finally there are nine nodes of depth 2.

We may find it helpful to see the data in a wide format:

spread_paths(df %>% select(-fit, -parent_fit, -dose_index))
## # A tibble: 9 x 6
##   outcomes0 next_dose0 outcomes1 next_dose1 outcomes2 next_dose2
##   <chr>          <dbl> <chr>          <dbl> <chr>          <dbl>
## 1 ""                 2 NN                 4 NN                 5
## 2 ""                 2 NN                 4 NT                 3
## 3 ""                 2 NN                 4 TT                 2
## 4 ""                 2 NT                 1 NN                 2
## 5 ""                 2 NT                 1 NT                 1
## 6 ""                 2 NT                 1 TT                 1
## 7 ""                 2 TT                 1 NN                 1
## 8 ""                 2 TT                 1 NT                 1
## 9 ""                 2 TT                 1 TT                 1

This arranges the paths in rows. Similar columns are output for each node in the path. The column names are suffixed so they can be distinguised from one another. The table above confirms that there are nine ways to traverse these first two cohort of two patients.

The pathways we have calculated thus far were calculated from a blank slate - no patients had yet been treated. Pathways may also be calculated for trials in progress. All that is required is that we specify the trial path already observed using the previous_outcomes parameter. In the following example, we calculate pathways for the next two cohorts of three patients, having observed outcomes 2NN 3TN.

paths2 <- crm_dtps(skeleton = skeleton,
target = target,
model = 'empiric',
cohort_sizes = c(3, 3),
previous_outcomes = '2NN 3TN',
next_dose = 2,
beta_sd = 1,
refresh = 0)
spread_paths(as_tibble(paths2) %>% select(-fit, -parent_fit, -dose_index))
## # A tibble: 16 x 6
##    outcomes0 next_dose0 outcomes1 next_dose1 outcomes2 next_dose2
##    <chr>          <dbl> <chr>          <dbl> <chr>          <dbl>
##  1 ""                 2 NNN                3 NNN                4
##  2 ""                 2 NNN                3 NNT                3
##  3 ""                 2 NNN                3 NTT                2
##  4 ""                 2 NNN                3 TTT                1
##  5 ""                 2 NNT                2 NNN                3
##  6 ""                 2 NNT                2 NNT                2
##  7 ""                 2 NNT                2 NTT                1
##  8 ""                 2 NNT                2 TTT                1
##  9 ""                 2 NTT                1 NNN                2
## 10 ""                 2 NTT                1 NNT                1
## 11 ""                 2 NTT                1 NTT                1
## 12 ""                 2 NTT                1 TTT                1
## 13 ""                 2 TTT                1 NNN                1
## 14 ""                 2 TTT                1 NNT                1
## 15 ""                 2 TTT                1 NTT                1
## 16 ""                 2 TTT                1 TTT                1

When calculaing pathways, the dose given for the next cohort is that with posterior mean probability of toxicity closest to the target. This is the default behaviour for the CRM model. However, we might wish to tailor the algorithm used to select doses. For instance, we might wish to avoid the skipping of doses in escalation and incorporate a mechanism that advises stopping the trial if excess toxicity is seen. As described above, these behaviours are provided by the careful_escalation function. We provide a custom dose selection function via the user_dose_func parameter:

paths3 <- crm_dtps(
skeleton = skeleton,
target = target,
model = 'empiric',
cohort_sizes = c(3, 3),
previous_outcomes = '2NN 3TN',
next_dose = 2,
beta_sd = 1,
user_dose_func = function(x) {
careful_escalation(x, tox_threshold = target + 0.1,
certainty_threshold = 0.7)
},
seed = 123, refresh = 0)

df3 <- as_tibble(paths3)
spread_paths(df3 %>% select(-fit, -parent_fit, -dose_index))
## # A tibble: 16 x 6
##    outcomes0 next_dose0 outcomes1 next_dose1 outcomes2 next_dose2
##    <chr>          <dbl> <chr>          <dbl> <chr>          <dbl>
##  1 ""                 2 NNN                3 NNN                4
##  2 ""                 2 NNN                3 NNT                3
##  3 ""                 2 NNN                3 NTT                2
##  4 ""                 2 NNN                3 TTT                2
##  5 ""                 2 NNT                2 NNN                3
##  6 ""                 2 NNT                2 NNT                2
##  7 ""                 2 NNT                2 NTT                1
##  8 ""                 2 NNT                2 TTT                1
##  9 ""                 2 NTT                1 NNN                2
## 10 ""                 2 NTT                1 NNT                1
## 11 ""                 2 NTT                1 NTT                1
## 12 ""                 2 NTT                1 TTT                1
## 13 ""                 2 TTT                1 NNN                1
## 14 ""                 2 TTT                1 NNT                1
## 15 ""                 2 TTT                1 NTT                1
## 16 ""                 2 TTT                1 TTT               NA

We see that the custom dose function yields a single difference compared to paths2. After observing 2TTT 1TTT, the custom function stops because $$Prob(DLT_1 > 0.35 | X) > 0.7$$, i.e. it is likely that even dose-level 1 exceeds our target level of toxicity.

We use careful_escalation here for illustration. user_dose_func can be any function that takes a crm_fit as its single parameter and returns either an integer dose-level, or NA to signify that the trial should stop.

The analysis of DTPs becomes bewildering as the volume of information grows. A clear method of visualisation is extremely valuable. The DiagrammeR package creates graphs of nodes and edges. This is perfect for visualising tree-like structures like those created by crm_dtps.

To create a graph, DiagrammeR requires a data.frame of nodes and another of edges. The nodes are the shapes in a graph. The edges are the lines that join them. trialr has already given us all of the information we require to define these elements. All that remains is to choose a pleasing colour scheme.

# This section of code outputs to the Viewer pane in RStudio
if(Sys.getenv("RSTUDIO") == "1") {

library(DiagrammeR)

df3 %>%
transmute(id = .node,
type = NA,
label = case_when(
is.na(next_dose) ~ 'Stop',
TRUE ~ next_dose %>% as.character()),
shape = 'circle',
fillcolor = case_when(
next_dose == 1 ~ 'slategrey',
next_dose == 2 ~ 'skyblue1',
next_dose == 3 ~ 'royalblue1',
next_dose == 4 ~ 'orchid4',
next_dose == 5 ~ 'royalblue4',
is.na(next_dose) ~ 'red'
)
) -> ndf

df3 %>%
filter(!is.na(.parent)) %>%
select(from = .parent, to = .node, label = outcomes) %>%
mutate(rel = "leading_to") -> edf

graph <- create_graph(nodes_df = ndf, edges_df = edf)
render_graph(graph)
}

Compared to the tabular form above, the graph representation of DTPs is far superior. The central node shows that dose 2 will be given to the next cohort. Displaying the nodes that require stopping in a bold and symbolic colour like red immediately conveys useful information. The only path that will see the model advocate stopping in the next two cohorts is if every patient experiences toxicity.

Using graphs and DiagrammeR to visualise paths is potentially very flexible and powerful. Obviously different colours and shapes may be chosen. The opacity (or alpha) of the nodes and paths may be adjusted by their probability of occurrence so that paths that are less likely appear more feint. Opportunities abound.

# Summary

That concludes this vignette on the analysis of dose transitions with CRM in trialr. The functions described will help investigators design and conduct high-quality dose-finding clinical trials. Behind the scenes, this package leverages the computational power of rstan to fit its models. The objects returned work nicely with modern tidyverse packages and programming-styles to offer a flexible suite of tools for clinical trialists.

## Other CRM vignettes

There are many vignettes illuminating the CRM in trialr:

# trialr

trialr is available at https://github.com/brockk/trialr and https://CRAN.R-project.org/package=trialr

# References

Brock, Kristian, Lucinda Billingham, Mhairi Copland, Shamyla Siddique, Mirjana Sirovica, and Christina Yap. 2017. “Implementing the EffTox Dose-Finding Design in the Matchpoint Trial.” BMC Medical Research Methodology 17 (1): 112. https://doi.org/10.1186/s12874-017-0381-x.

Chiuzan, Cody, Jonathan Shtaynberger, Gulam A. Manji, Jimmy K. Duong, Gary K. Schwartz, Anastasia Ivanova, and Shing M. Lee. 2017. “Dose-Finding Designs for Trials of Molecularly Targeted Agents and Immunotherapies.” Journal of Biopharmaceutical Statistics 27 (3): 477–94. https://doi.org/10.1080/10543406.2017.1289952.

Iasonos, Alexia, Andrew S Wilton, Elyn R Riedel, Venkatraman E Seshan, and David R Spriggs. 2008. “A Comprehensive Comparison of the Continual Reassessment Method to the Standard 3 + 3 Dose Escalation Scheme in Phase I Dose-Finding Studies.” Clinical Trials (London, England) 5 (5): 465–77. https://doi.org/10.1177/1740774508096474.

Le Tourneau, Christophe, J. Jack Lee, and Lillian L. Siu. 2009. “Dose Escalation Methods in Phase I Cancer Clinical Trials.” Journal of the National Cancer Institute 101 (10): 708–20. https://doi.org/10.1093/jnci/djp079.

O’Quigley, J, M Pepe, and L Fisher. 1990. “Continual Reassessment Method: A Practical Design for Phase 1 Clinical Trials in Cancer.” Biometrics 46 (1): 33–48. https://doi.org/10.2307/2531628.

Rogatko, André, David Schoeneck, William Jonas, Mourad Tighiouart, Fadlo R. Khuri, and Alan Porter. 2007. “Translation of Innovative Designs into Phase I Trials.” Journal of Clinical Oncology 25 (31): 4982–6. https://doi.org/10.1200/JCO.2007.12.1012.

Yap, Christina, Lucinda J. Billingham, Ying Kuen Cheung, Charlie Craddock, and John O’Quigley. 2017. “Dose Transition Pathways: The Missing Link Between Complex Dose-Finding Designs and Simple Decision-Making.” Clinical Cancer Research 23 (24): 7440–7. https://doi.org/10.1158/1078-0432.CCR-17-0582.