Design of ordinal CRM

The Ordinal CRM

The ordinal CRM, or oCRM, is an extension to the standard CRM in which the patient’s response is not binary (no DLT, DLT) but ordinal (for example: None, sub-DLT, DLT). The number of categories in the response scale is arbitrary. The responses are modelled using a standard parameterisation for ordered logistic regression:

Let p_k(d) be the probability that the response of a patient treated at dose d is in category k or higher, k=0, …, K; d=1, …, D.

Then logit(p_k(d)) = α_k + β log(d/d_ref), k=1, …, K. [p₀(d) = 1 by definition.]

where d_ref is a reference dose. The αs are constrained such that α₁ > α₂ > … > α_K.

The model is a constrained parallel lines regression on the logit scale with appropriate priors imposed on the α_ks and β.

Unfortunately, I know of no references for the oCRM, other than the FACTS user manual.

Example analysis from first prinicples

This example is motivated by a real Roche study in which the study team became concerned by the number of episodes of cytokine release syndrome (CRS) that did not strictly meet the definition of a DLT given in the protocol, yet which were serious enough to cause some disquiet.

In the example data below, the dose grid reflects reality, but the DLT and sub-DLT counts are entirely arbitrary.

Fictitious observed data

observedData %>%
  group_by(Dose, Status) %>%
  summarise(N = n(), .groups = "drop") %>%
  pivot_wider(
    names_from = Status,
    values_from = N,
    values_fill = 0
  ) %>%
  rename(NoDLT = `0`, SubDLT = `1`, DLT = `2`) %>%
  mutate(Treated = NoDLT + SubDLT + DLT, .before = 2)

## # A tibble: 13 × 5
##     Dose Treated NoDLT SubDLT   DLT
##    <dbl>   <int> <int>  <int> <int>
##  1     5       1     1      0     0
##  2    15       4     4      0     0
##  3    45       5     5      0     0
##  4    70       5     5      0     0
##  5   100       5     5      0     0
##  6   220       8     7      1     0
##  7   300       6     5      1     0
##  8   600      15    10      5     0
##  9  1000       8     5      3     0
## 10  1800       9     8      1     0
## 11  4000      10     8      2     0
## 12 10000      14    10      4     0
## 13 16000       2     0      1     1

Analysis

Define the model.

modelString <- "
data
{
  for (i in 1:length(d))
  {
    for (j in 1:2)
    {
      DLT[i, j] <- r[i] >= j
    }
  }
}
model
{
  #Independent univariate parameters for simplicity.
  #Recall JAGS uses precision, not variance.
  alpha[1] ~ dnorm(meanAlpha1, 1/(sdAlpha1*sdAlpha1))
  #                                                    Constrain the model
  alpha[2] ~ dnorm(meanAlpha2, 1/(sdAlpha2*sdAlpha2))  T(, alpha[1])
  #Common slope.  LogNormal distribution ensures slope is positive
  gamma ~ dnorm(meanLogBeta, 1/(sdLogBeta*sdLogBeta))
  beta <- exp(gamma)
  for (i in 1:length(d))
  {
    xhat[i] <- log(d[i] / dRef)
    for (j in 1:2)
    {
      z[i, j] <- alpha[j] + beta * xhat[i]
      p[i, j] <- exp(z[i, j]) / (1 + exp(z[i, j]))
      DLT[i, j] ~ dbern(p[i, j])
    }
  }
}
Inits
{
  list(alpha=c(5, 3), gamma=0)
}
#monitor# alpha[1], alpha[2], beta
#data# meanAlpha1, meanAlpha2, meanLogBeta, sdAlpha1, sdAlpha2, sdLogBeta, d, r, dRef
"

Define a reporting function.

analyse <- function(m, caption) {
  # Fit the model
  results <- as_tibble(run.jags(m)$mcmc[[1]]) %>%
    rename(alpha1 = `alpha[1]`, alpha2 = `alpha[2]`) # Fix awkwardly named columns

  # Derive p(DLT) [P2] and p(DLT or subDLT) [P1]
  results <- results %>%
    expand(nesting(alpha1, alpha2, beta), Dose = doseList) %>%
    mutate(
      XHat = log(Dose / dRef),
      z1 = alpha1 + beta * XHat,
      z2 = alpha2 + beta * XHat,
      P1 = exp(z1) / (1 + exp(z1)),
      P2 = exp(z2) / (1 + exp(z2))
    ) %>%
    # Old code.  Should update to pivot_longer
    gather(key = Type, value = Prob, P1, P2)

  # Summarise
  summary <- results %>%
    group_by(Dose, Type) %>%
    summarise(
      N = n(),
      Mean = mean(Prob, na.rm = TRUE),
      q50 = median(Prob, na.rm = TRUE),
      q05 = quantile(Prob, 0.05, na.rm = TRUE),
      q10 = quantile(Prob, 0.10, na.rm = TRUE),
      q20 = quantile(Prob, 0.20, na.rm = TRUE),
      q80 = quantile(Prob, 0.80, na.rm = TRUE),
      q90 = quantile(Prob, 0.90, na.rm = TRUE),
      q95 = quantile(Prob, 0.95, na.rm = TRUE),
      .groups = "drop"
    )

  table <- summary %>%
    # Old code.  Should update to pivot_longer
    gather(variable, value, -(Dose:Type)) %>%
    unite(temp, Type, variable) %>%
    spread(temp, value) %>%
    select(
      Dose, P1_N, P1_Mean, starts_with("P1_q"),
      P2_N, P2_Mean, starts_with("P2_q")
    )

  plot <- summary %>%
    ggplot() +
    geom_line(aes(x = Dose, y = Mean, colour = Type)) +
    geom_ribbon(aes(x = Dose, ymin = q05, ymax = q95, fill = Type), alpha = 0.1) +
    geom_ribbon(aes(x = Dose, ymin = q10, ymax = q90, fill = Type), alpha = 0.1) +
    geom_ribbon(aes(x = Dose, ymin = q20, ymax = q80, fill = Type), alpha = 0.1) +
    theme_light() +
    scale_colour_manual(name = " ", labels = c("Mean p(DLT or subDLT)", "Mean p(DLT)"), values = c("darkblue", "red")) +
    scale_fill_manual(name = " ", labels = c("Mean p(DLT or subDLT)", "Mean p(DLT)"), values = c("darkblue", "red")) +
    labs(
      y = "Probability",
      x = "Dose",
      title = caption,
      caption = "Shading shows the 90%, 80% and 60% central credible intervals"
    )
  rv <- NULL
  rv$summary <- summary
  rv$graph <- plot
  rv$tableWide <- table
  rv$samples <- results
  return(rv)
}

Define the hyperpriors and other model parameters. Values are entirely arbitrary.

dRef <- 450
meanAlpha1 <- 5
meanAlpha2 <- 3
meanLogBeta <- log(1)
sdAlpha1 <- 4
sdAlpha2 <- 4
sdLogBeta <- 3

# Hack the observed data
d <- observedData$Dose
r <- observedData$Status

Obtain the posterior

post <- analyse(modelString, caption = "Posterior distribution")

## Note: Transposing BUGS inits into R format

## Warning: Inits specified in the model file or using #inits# are ignored when a
## character string is given as the argument to inits

## Loading required namespace: rjags

## Warning: No initial values were provided - JAGS will use the same initial values
## for all chains

## Compiling rjags model...
## Calling the simulation using the rjags method...
## Adapting the model for 1000 iterations...
## Burning in the model for 4000 iterations...
## Running the model for 10000 iterations...
## Simulation complete
## Calculating summary statistics...
## Calculating the Gelman-Rubin statistic for 3 variables....
## Finished running the simulation

## Warning: attributes are not identical across measure variables;
## they will be dropped

post$tableWide %>%
  select(-P2_N) %>%
  kable(
    digits = c(0, 0, rep(2, 12)),
    col.names = c("Dose", "N", rep(c("Mean", "Q05", "Q10", "Q20", "Median", "Q80", "Q90", "Q95"), 2))
  ) %>%
  add_header_above(c(" " = 2, "p(DLT or subDLT)" = 8, "p(DLT)" = 8)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

		p(DLT or subDLT)								p(DLT)
Dose	N	Mean	Q05	Q10	Q20	Median	Q80	Q90	Q95	Mean	Q05	Q10	Q20	Median	Q80	Q90	Q95
5	10000	0.04	0.01	0.01	0.01	0.03	0.07	0.10	0.13	0.00	0	0.00	0.00	0	0	0.01	0.01
15	10000	0.06	0.01	0.02	0.02	0.05	0.09	0.12	0.15	0.00	0	0.00	0.00	0	0	0.01	0.01
45	10000	0.08	0.02	0.03	0.04	0.07	0.11	0.14	0.17	0.01	0	0.00	0.00	0	0	0.01	0.02
70	10000	0.09	0.03	0.04	0.05	0.08	0.13	0.16	0.18	0.01	0	0.00	0.00	0	0	0.01	0.02
100	10000	0.10	0.04	0.05	0.06	0.09	0.14	0.16	0.19	0.01	0	0.00	0.00	0	0	0.02	0.02
220	10000	0.13	0.06	0.07	0.09	0.12	0.16	0.19	0.21	0.01	0	0.00	0.00	0	0	0.02	0.02
300	10000	0.14	0.07	0.08	0.10	0.14	0.18	0.20	0.22	0.01	0	0.00	0.00	0	0	0.02	0.03
600	10000	0.17	0.10	0.12	0.13	0.17	0.21	0.23	0.25	0.01	0	0.00	0.00	0	0	0.03	0.03
1000	10000	0.20	0.13	0.14	0.16	0.20	0.24	0.26	0.28	0.01	0	0.00	0.00	0	0	0.03	0.04
1800	10000	0.24	0.16	0.18	0.20	0.24	0.28	0.30	0.32	0.02	0	0.00	0.01	0	0	0.04	0.05
4000	10000	0.30	0.20	0.22	0.24	0.30	0.35	0.38	0.41	0.02	0	0.00	0.01	0	0	0.05	0.06
10000	10000	0.38	0.22	0.25	0.29	0.38	0.46	0.51	0.54	0.04	0	0.01	0.01	0	0	0.07	0.09
16000	10000	0.42	0.24	0.27	0.32	0.42	0.52	0.57	0.61	0.04	0	0.01	0.01	0	0	0.09	0.11

post$graph

Recommended dose

The recommended dose (and stopping rules) can be based on any (combination) of the p_ks. For example, if the recommended dose is defined to be the highest dose with posterior mean p₁ less than 0.2, then the recommended dose is

post$summary %>%
  ungroup() %>%
  filter(
    Mean < 0.2,
    Type == "P1"
  ) %>%
  summarise(Recommended = max(Dose))

## # A tibble: 1 × 1
##   Recommended
##         <dbl>
## 1        1000

Integration with nCRM

Although an extension to the standard CRM model, there’s no reason why oCRM can’t be applied within an nCRM framework.

For example, suppose the recommended dose is determined on the basis of the p₁s [that is, on p(DLT or subDLT | dose)] and that the toxicity bands are defined as

	Bounds
Band	Lower	Upper
Underdosing	0.00	0.15
Target toxicity	0.15	0.25
Excess toxicity	0.25	0.40
Unacceptable toxicity	0.40	1.00

Then band membership can be calculated for each of the MCMC samples…

nCRM <- post$samples %>%
  filter(Type == "P1") %>%
  mutate(
    Under = Prob <= 0.15,
    Target = Prob > 0.15 & Prob <= 0.25,
    Excess = Prob > 0.25 & Prob <= 0.4,
    Unacceptable = Prob > 0.4,
    Check = Under + Target + Excess + Unacceptable
  )

… and summarised by dose

nCRM %>%
  select(-Check) %>%
  group_by(Dose) %>%
  summarise(
    Under = mean(Under),
    Target = mean(Target),
    Excess = mean(Excess),
    Unacceptable = mean(Unacceptable),
    .groups = "drop"
  ) %>%
  kable(
    table.attr = "style='width:40%;'",
    digits = c(0, 3, 3, 3, 3)
  ) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed")) %>%
  add_header_above(c(" " = 1, "Probability of band membership" = 4))

	Probability of band membership
Dose	Under	Target	Excess	Unacceptable
5	0.964	0.034	0.002	0.000
15	0.950	0.048	0.002	0.000
45	0.914	0.082	0.003	0.000
70	0.885	0.111	0.004	0.000
100	0.855	0.140	0.005	0.000
220	0.722	0.269	0.009	0.000
300	0.625	0.361	0.014	0.000
600	0.335	0.622	0.043	0.000
1000	0.133	0.732	0.135	0.000
1800	0.028	0.582	0.388	0.002
4000	0.005	0.236	0.691	0.068
10000	0.003	0.088	0.503	0.406
16000	0.002	0.062	0.374	0.563

Design discussion

Which Models should be extended to allow ordinal models?
Strictly speaking, the existing Model classes are special cases of the more general ordinal classes, with k=2 rather than k > 1 categories. However, I suspect refactoring the existing classes to be sub-classes of the corresponding ordinal Model class.
nextBest methods that handle the new Models will be needed.
New Stopping rules will be required.
Realistically, I expect the maximum number of categories required will be 3 or 4. Should we impose a maximum or allow an arbitrary number?
How best to provide arbitrary labels for DLT categories, if at all? The obvious options seem to me to be to provide alternate parameters in the constructor function(s). For example OrdinalLogisticLogNormal(k=3) to provide unnamed categories or default labels, or OrdinalLogisticLogNormal(dltGrades=c(0="None", 1="subDLT", 2="DLT")) to provide custom labels and ordering.
Other issues?

Environment

sessionInfo()

## R version 4.2.1 (2022-06-23)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur ... 10.16
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] runjags_2.2.1-7  kableExtra_1.3.4 forcats_0.5.2    stringr_1.4.1   
##  [5] dplyr_1.0.10     purrr_0.3.5      readr_2.1.3      tidyr_1.2.1     
##  [9] tibble_3.1.8     ggplot2_3.3.6    tidyverse_1.3.2 
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-45     svglite_2.1.0       lubridate_1.8.0    
##  [4] assertthat_0.2.1    digest_0.6.30       utf8_1.2.2         
##  [7] R6_2.5.1            cellranger_1.1.0    backports_1.4.1    
## [10] reprex_2.0.2        coda_0.19-4         evaluate_0.17      
## [13] highr_0.9           httr_1.4.4          pillar_1.8.1       
## [16] rlang_1.0.6         googlesheets4_1.0.1 readxl_1.4.1       
## [19] rstudioapi_0.14     jquerylib_0.1.4     rjags_4-13         
## [22] rmarkdown_2.17      labeling_0.4.2      webshot_0.5.4      
## [25] googledrive_2.0.0   munsell_0.5.0       broom_1.0.1        
## [28] compiler_4.2.1      modelr_0.1.9        xfun_0.34          
## [31] pkgconfig_2.0.3     systemfonts_1.0.4   htmltools_0.5.3    
## [34] tidyselect_1.2.0    viridisLite_0.4.1   fansi_1.0.3        
## [37] crayon_1.5.2        tzdb_0.3.0          dbplyr_2.2.1       
## [40] withr_2.5.0         grid_4.2.1          jsonlite_1.8.2     
## [43] gtable_0.3.1        lifecycle_1.0.3     DBI_1.1.3          
## [46] magrittr_2.0.3      scales_1.2.1        cli_3.4.1          
## [49] stringi_1.7.8       cachem_1.0.6        farver_2.1.1       
## [52] fs_1.5.2            xml2_1.3.3          bslib_0.4.0        
## [55] ellipsis_0.3.2      generics_0.1.3      vctrs_0.5.1        
## [58] tools_4.2.1         glue_1.6.2          hms_1.1.2          
## [61] parallel_4.2.1      fastmap_1.1.0       yaml_2.3.6         
## [64] colorspace_2.0-3    gargle_1.2.1        rvest_1.0.3        
## [67] knitr_1.40          haven_2.5.1         sass_0.4.2