In Roche/crmPack: Object-Oriented Implementation of CRM Designs

knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
library(kableExtra)
library(runjags)

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~0~(d) = 1 by definition.]

where d~ref~ is a reference dose. The αs are constrained such that α~1~ > α~2~ > ... > α~K~.

The model is a constrained parallel lines regression on the logit scale with appropriate priors imposed on the α~k~s 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

doseList <- c(5, 15, 45, 70, 100, 220, 300, 600, 1000, 1800, 4000, 10000, 16000)

observedData <- tibble(
  Dose = c(
    rep(5, 1),
    rep(15, 4),
    rep(45, 5),
    rep(70, 5),
    rep(100, 5),
    rep(220, 8),
    rep(300, 6),
    rep(600, 15),
    rep(1000, 8),
    rep(1800, 9),
    rep(4000, 10),
    rep(10000, 14),
    rep(16000, 2)
  ),
  Status = c(
    rep(0, 1),
    rep(0, 4),
    rep(0, 5),
    rep(0, 5),
    rep(0, 5),
    rep(0, 7), 1,
    rep(0, 5), 1,
    rep(0, 10), 1, 1, 1, 1, 1,
    rep(0, 5), 1, 1, 1,
    rep(0, 8), 1,
    rep(0, 8), 1, 1,
    rep(0, 10), 1, 1, 1, 1,
    rep(1, 1),
    rep(2, 1)
  )
)

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)

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")

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"))

post$graph

Recommended dose

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

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

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~1~s [that is, on p(DLT or subDLT | dose)] and that the toxicity bands are defined as

tibble(
  Band = c("Underdosing", "Target toxicity", "Excess toxicity", "Unacceptable toxicity"),
  Lower = c(0, 0.15, 0.25, 0.4),
  Upper = c(0.15, 0.25, 0.4, 1)
) %>%
  kable(table.attr = "style='width:30%;'") %>%
  add_header_above(c(" " = 1, "Bounds" = 2)) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"))

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))

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?