oCRM demonstration
Test data
Based on a real (non-Roche) study, using a three category toxicity scale (0=“No event”, 1=“sub-toxic AE”, 2=“DLT”).
The doses used were 1, 2, 4, 8, 16 and 30 units. (Doses of 40 and 50 were also available, but were not used.)
The prior distribution of model parameters was
\[ \begin{bmatrix} \alpha_1 \\ \alpha_2 \\ log(\beta) \end{bmatrix} \sim N \left( \begin{bmatrix} 4 \\ 3 \\ 0 \end{bmatrix}, \begin{bmatrix} 3 & 0 & 0\\ 0 & 4 & 0 \\ 0 & 0 & 1 \end{bmatrix} \right) \]
| Subject | Cohort | Dose | Tox status |
|---|---|---|---|
| 1 | 1 | 1 | 0 |
| 2 | 1 | 1 | 0 |
| 3 | 2 | 2 | 0 |
| 4 | 2 | 2 | 0 |
| 5 | 3 | 4 | 0 |
| 6 | 3 | 4 | 0 |
| 7 | 4 | 8 | 0 |
| 8 | 4 | 8 | 0 |
| 9 | 5 | 16 | 1 |
| 10 | 5 | 16 | 2 |
| 11 | 6 | 30 | 2 |
| 12 | 6 | 30 | 1 |
| 13 | 7 | 30 | 1 |
| 14 | 7 | 30 | 1 |
| 15 | 8 | 30 | 1 |
The Ordinal CRM model
Define the oCRM model objects (almost) as specified in the design document
alpha[j] ~ dnorm(meanAlpha[j], precAlpha[j]) T(, alpha[j-1]) is valid JAGS syntax, but not valid R syntax. Note the use of the %_% operator to overcome this. For details, see the online doc for h_jags_write_model().
TODO: Need to convert the vector of integer categories
yto a matrix of indicator valuesDLT[,]whereDLT[i, j] = (y[i] >= j). I’m not sure how best to do this within thecrmPackecosystem. See the validation section below for an example implementation.
.DataOrdinal <- setClass(
Class = "DataOrdinal",
contains = "GeneralData",
slots = c(
params = "ModelParamsNormal",
x = "numeric",
y = "integer",
doseGrid = "numeric",
nGrid = "integer",
xLevel = "integer",
yCategories = "integer",
placebo = "logical"
),
prototype = prototype(
x = numeric(),
y = integer(),
doseGrid = numeric(),
nGrid = 0L,
xLevel = integer(),
yCategories = c("No DLT" = 0L, "DLT" = 1L),
placebo = FALSE
)
# ,
# validity = v_data_cat_tox
)
.Data <- setClass(
Class = "Data",
contains = "DataOrdinal"
# ,
# validity = v_data
)
.LogisticLogNormOrd <- setClass(
Class = "LogisticLogNormalOrd",
contains = "ModelLogNormal"
)
LogisticLogNormOrd <- function(meanAlpha, varAlpha, meanBeta, varBeta, ref_dose = 0) {
model_ln <- ModelParamsNormal(
mean = c(meanAlpha, meanBeta),
cov = diag(c(varAlpha, varBeta))
)
.LogisticLogNormOrd(
params = model_ln,
ref_dose = crmPack:::positive_number(ref_dose),
datamodel = function() {
for (i in 1:15) {
x_rel[i] <- log(x[i] / ref_dose)
for (j in 1:length(meanAlpha)) {
logit(p[i, j]) <- alpha[j] + beta * x_rel[i]
DLT[i, j] ~ dbern(p[i, j])
}
}
},
priormodel = function() {
alpha[1] ~ dnorm(meanAlpha[1], precAlpha[1, 1])
for (j in 2:length(meanAlpha)) {
alpha[j] ~ dnorm(meanAlpha[j], precAlpha[j, j]) %_% T( , alpha[j-1])
}
beta ~ dnorm(meanBeta, precBeta)
},
modelspecs = function(from_prior) {
ms <- list(
meanAlpha = meanAlpha,
precAlpha = solve(diag(varAlpha)),
meanBeta = meanBeta,
precBeta = 1/varBeta
)
if (!from_prior) {
ms$ref_dose <- ref_dose
}
ms
},
init = function(yCategories) {
list(alpha = seq(from=5, to=3, length.out=length(yCategories)-1), beta = 0)
},
datanames = c("nObs", "yCategories", "y", "x"),
sample = c("alpha", "beta")
)
}Check consequential code modifications
Modify the ParamsLogNormal class to allow for dimensions other than 2
model_ln <- ModelParamsNormal(
mean = c(c(3, 4), 0),
cov = diag(c(c(4, 3), 1))
)
model_lnAn object of class "ModelParamsNormal"
Slot "mean":
[1] 3 4 0
Slot "cov":
[,1] [,2] [,3]
[1,] 4 0 0
[2,] 0 3 0
[3,] 0 0 1
Slot "prec":
[,1] [,2] [,3]
[1,] 0.25 0.0000000 0
[2,] 0.00 0.3333333 0
[3,] 0.00 0.0000000 1Change the assertion in h_jags_get_model_inits to allow init to be a closure in .LogisticLogNormOrd.
standardCRM <- ModelLogNormal(
mean = c(3, 0),
cov = matrix(c(3, 0, 0, 1), ncol = 2),
ref_dose = 1
)
validObject(standardCRM)[1] TRUEoCRM <- LogisticLogNormOrd(
meanAlpha = c(3, 4),
varAlpha = c(4, 3),
meanBeta = 0,
varBeta = 1,
ref_dose = 1
)
validObject(oCRM)[1] TRUEAlso:
- Changed the assertion in
h_jags_add_dummyfromassert_class(object, "Data")toassert_class(object, "DataOrdinal"). - Loops over
yCategoriesshould I think be from2tolength(yCategories)and not any other variation. - Temporarily removed validity checking for simplicity
Evaluate the prior
mcmcSummary and summaryTable are utility functions defined earlier in this document.
emptyData <- new(
"DataOrdinal",
doseGrid=c(1, 2, 4, 8, 16, 30, 40, 50),
yCategories = c("No event" = 0L, "sub-toxic AE" = 1L, "DLT" = 2L)
)
model <- LogisticLogNormOrd(
meanAlpha = c(-2, -1),
meanBeta = 0,
varAlpha = c(8, 4),
varBeta = 2,
ref_dose = 25
)
prior <- mcmc(emptyData, model, McmcOptions(samples = 25000))
priorSummary <- prior@data %>% mcmcSummary()
priorSummary %>%
ggplot(aes(x=Dose)) +
geom_ribbon(aes(ymin=Q10, ymax=Q90), alpha=0.3, fill="steelblue") +
geom_line(aes(y=Median)) +
theme_light() +
facet_wrap(vars(Category), labeller=label_both) +
labs(title = "Mean prior probability of toxicity category > x, with 80% ci")
priorSummary %>% summaryTable(emptyData, "Prior probability of toxicity category > x")| Dose | Median | 10th centile | 90th centile |
|---|---|---|---|
| sub-toxic AE | |||
| 1 | 0.00 | 0.00 | 0.03 |
| 2 | 0.00 | 0.00 | 0.04 |
| 4 | 0.01 | 0.00 | 0.04 |
| 8 | 0.02 | 0.00 | 0.05 |
| 16 | 0.04 | 0.00 | 0.05 |
| 30 | 0.07 | 0.06 | 0.16 |
| 40 | 0.09 | 0.06 | 0.52 |
| 50 | 0.11 | 0.06 | 0.82 |
| DLT | |||
| 1 | 0.00 | 0.00 | 0.04 |
| 2 | 0.01 | 0.00 | 0.05 |
| 4 | 0.01 | 0.00 | 0.05 |
| 8 | 0.02 | 0.00 | 0.06 |
| 16 | 0.05 | 0.00 | 0.06 |
| 30 | 0.08 | 0.07 | 0.19 |
| 40 | 0.11 | 0.07 | 0.58 |
| 50 | 0.13 | 0.08 | 0.85 |
Evaluate the posterior
postData <- new(
"DataOrdinal",
x = df$Dose,
y = as.integer(df$ToxStatus),
doseGrid=c(1, 2, 4, 8, 16, 30, 40, 50),
yCategories = c("No event" = 0L, "sub-toxic AE" = 1L, "DLT" = 2L)
)
post <- mcmc(postData, model, McmcOptions(samples = 25000))
postSummary <- post@data %>% mcmcSummary()
postSummary %>%
ggplot(aes(x=Dose)) +
geom_ribbon(aes(ymin=Q10, ymax=Q90), alpha=0.3, fill="steelblue") +
geom_line(aes(y=Median)) +
theme_light() +
facet_wrap(vars(Category), labeller=label_both) +
labs(title = "Mean posterior probability of toxicity category > x, with 80% ci")
postSummary %>% summaryTable(emptyData, "Posterior probability of toxicity category > x")| Dose | Median | 10th centile | 90th centile |
|---|---|---|---|
| sub-toxic AE | |||
| 1 | 0.00 | 0.00 | 0.04 |
| 2 | 0.00 | 0.00 | 0.04 |
| 4 | 0.01 | 0.00 | 0.05 |
| 8 | 0.02 | 0.00 | 0.05 |
| 16 | 0.04 | 0.00 | 0.06 |
| 30 | 0.07 | 0.06 | 0.17 |
| 40 | 0.10 | 0.07 | 0.55 |
| 50 | 0.12 | 0.07 | 0.83 |
| DLT | |||
| 1 | 0.45 | 0.00 | 0.93 |
| 2 | 0.62 | 0.00 | 0.93 |
| 4 | 0.77 | 0.00 | 0.94 |
| 8 | 0.87 | 0.02 | 0.95 |
| 16 | 0.93 | 0.58 | 0.95 |
| 30 | 0.96 | 0.96 | 0.99 |
| 40 | 0.97 | 0.96 | 1.00 |
| 50 | 0.98 | 0.96 | 1.00 |
Validation
To confirm the accuracy of the results above, perform the same analysis from first principles.
Prior
priorModelString <- "
model {
#Independent univariate parameters for clarity
alpha[1] ~ dnorm(meanAlpha1, 1/(sdAlpha1*sdAlpha1))
alpha[2] ~ dnorm(meanAlpha2, 1/(sdAlpha2*sdAlpha2)) T(, alpha[1])
#Common slope.
gamma ~ dnorm(meanLogBeta, 1/(sdLogBeta*sdLogBeta))
beta <- exp(gamma)
}
Inits {
list(alpha=c(5, 3), gamma=0)
}
#monitor# alpha[1], alpha[2], beta
#data# meanAlpha1, meanAlpha2, meanLogBeta, sdAlpha1, sdAlpha2, sdLogBeta
"
meanAlpha1 <- 5
meanAlpha2 <- 3
meanLogBeta <- log(1)
sdAlpha1 <- 4
sdAlpha2 <- 4
sdLogBeta <- 3
priorValidation <- as_tibble(run.jags(priorModelString)$mcmc[[1]]) %>%
rename(Alpha1=`alpha[1]`, Alpha2=`alpha[2]`, Beta=beta) %>%
expand(nesting(Alpha1, Alpha2, Beta), Dose=emptyData@doseGrid) %>%
mutate(
Z1=exp(Alpha1 + exp(Beta)*log(Dose/model@ref_dose)),
Z2=exp(Alpha2 + exp(Beta)*log(Dose/model@ref_dose)),
Prob1=Z1 / (1 + Z1),
Prob2=Z2 / (1 + Z2)
) %>%
pivot_longer(
starts_with("Prob"),
names_to="Category",
values_to="Prob",
names_prefix="Prob"
) %>%
group_by(Dose, Category) %>%
summarise(
Median=median(Prob, na.rm = TRUE),
Q10=quantile(Prob, probs=0.1, na.rm = TRUE),
Q90=quantile(Prob, probs=0.9, na.rm = TRUE),
.groups="drop"
)Note: Transposing BUGS inits into R format
Compiling rjags model...
Calling the simulation using the rjags method...
Note: the model did not require adaptation
Burning in the model for 4000 iterations...
Running the model for 10000 iterations...
Simulation complete
Calculating summary statistics...Warning: Convergence cannot be assessed with only 1 chainFinished running the simulationpriorValidation %>%
ggplot(aes(x=Dose)) +
geom_ribbon(aes(ymin=Q10, ymax=Q90), alpha=0.3, fill="steelblue") +
geom_line(aes(y=Median)) +
theme_light() +
facet_wrap(vars(Category), labeller=label_both) +
labs(title = "Mean prior probability of toxicity category > x, with 80% ci")
priorValidation %>% summaryTable(emptyData, "Prior probability of toxicity category > x")| Dose | Median | 10th centile | 90th centile |
|---|---|---|---|
| sub-toxic AE | |||
| 1 | 0.00 | 0.00 | 0.99 |
| 2 | 0.01 | 0.00 | 0.99 |
| 4 | 0.05 | 0.00 | 1.00 |
| 8 | 0.20 | 0.00 | 1.00 |
| 16 | 0.63 | 0.00 | 1.00 |
| 30 | 1.00 | 0.64 | 1.00 |
| 40 | 1.00 | 0.74 | 1.00 |
| 50 | 1.00 | 0.80 | 1.00 |
| DLT | |||
| 1 | 0.00 | 0.00 | 0.37 |
| 2 | 0.00 | 0.00 | 0.58 |
| 4 | 0.00 | 0.00 | 0.77 |
| 8 | 0.00 | 0.00 | 0.89 |
| 16 | 0.03 | 0.00 | 0.96 |
| 30 | 0.86 | 0.03 | 1.00 |
| 40 | 0.92 | 0.05 | 1.00 |
| 50 | 0.95 | 0.07 | 1.00 |
Posterior
postModelString <- "
data {
for (i in 1:length(d)) {
for (j in 1:2) {
DLT[i, j] <- r[i] >= j
}
}
}
model {
#Independent univariate parameters for clarity
alpha[1] ~ dnorm(meanAlpha1, 1/(sdAlpha1*sdAlpha1))
alpha[2] ~ dnorm(meanAlpha2, 1/(sdAlpha2*sdAlpha2))
#Common slope.
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
"
meanAlpha1 <- 4
meanAlpha2 <- 3
meanLogBeta <- log(1)
sdAlpha1 <- 3
sdAlpha2 <- 4
sdLogBeta <- 1
d <- df$Dose
r <- df$ToxStatus
dRef <- model@ref_dose
postValidation <- as_tibble(run.jags(priorModelString)$mcmc[[1]]) %>%
rename(Alpha1=`alpha[1]`, Alpha2=`alpha[2]`, Beta=beta) %>%
expand(nesting(Alpha1, Alpha2, Beta), Dose=emptyData@doseGrid) %>%
mutate(
Z1=exp(Alpha1 + exp(Beta)*log(Dose/model@ref_dose)),
Z2=exp(Alpha2 + exp(Beta)*log(Dose/model@ref_dose)),
Prob1=Z1 / (1 + Z1),
Prob2=Z2 / (1 + Z2)
) %>%
pivot_longer(
starts_with("Prob"),
names_to="Category",
values_to="Prob",
names_prefix="Prob"
) %>%
group_by(Dose, Category) %>%
summarise(
Median=median(Prob, na.rm = TRUE),
Q10=quantile(Prob, probs=0.1, na.rm = TRUE),
Q90=quantile(Prob, probs=0.9, na.rm = TRUE),
.groups="drop"
)Note: Transposing BUGS inits into R format
Compiling rjags model...
Calling the simulation using the rjags method...
Note: the model did not require adaptation
Burning in the model for 4000 iterations...
Running the model for 10000 iterations...
Simulation complete
Calculating summary statistics...Warning: Convergence cannot be assessed with only 1 chainFinished running the simulationpostValidation %>%
ggplot(aes(x=Dose)) +
geom_ribbon(aes(ymin=Q10, ymax=Q90), alpha=0.3, fill="steelblue") +
geom_line(aes(y=Median)) +
theme_light() +
facet_wrap(vars(Category), labeller=label_both) +
labs(title = "Mean posterior probability of toxicity category > x, with 80% ci")
postValidation %>% summaryTable(emptyData, "Posterior probability of toxicity category > x")| Dose | Median | 10th centile | 90th centile |
|---|---|---|---|
| sub-toxic AE | |||
| 1 | 0.00 | 0.00 | 0.77 |
| 2 | 0.02 | 0.00 | 0.91 |
| 4 | 0.09 | 0.00 | 0.97 |
| 8 | 0.41 | 0.00 | 0.99 |
| 16 | 0.85 | 0.00 | 1.00 |
| 30 | 0.99 | 0.73 | 1.00 |
| 40 | 1.00 | 0.86 | 1.00 |
| 50 | 1.00 | 0.91 | 1.00 |
| DLT | |||
| 1 | 0.00 | 0.00 | 0.08 |
| 2 | 0.00 | 0.00 | 0.21 |
| 4 | 0.00 | 0.00 | 0.47 |
| 8 | 0.02 | 0.00 | 0.77 |
| 16 | 0.12 | 0.00 | 0.93 |
| 30 | 0.80 | 0.05 | 1.00 |
| 40 | 0.92 | 0.10 | 1.00 |
| 50 | 0.96 | 0.15 | 1.00 |
Environment
R version 4.2.1 (2022-06-23)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.5 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] runjags_2.2.1-7 broom_1.0.1 checkmate_2.1.0
[4] crmPack_2.0.0.9133 truncnorm_1.0-8 kableExtra_1.3.4
[7] knitr_1.40 forcats_0.5.2 stringr_1.4.1
[10] dplyr_1.0.10 purrr_0.3.5 readr_2.1.3
[13] tidyr_1.2.1 tibble_3.1.8 ggplot2_3.3.6
[16] tidyverse_1.3.2
loaded via a namespace (and not attached):
[1] httr_1.4.4 jsonlite_1.8.2 viridisLite_0.4.1
[4] modelr_0.1.9 assertthat_0.2.1 highr_0.9
[7] googlesheets4_1.0.1 cellranger_1.1.0 yaml_2.3.6
[10] lattice_0.20-45 pillar_1.8.1 backports_1.4.1
[13] glue_1.6.2 digest_0.6.30 rvest_1.0.3
[16] colorspace_2.0-3 htmltools_0.5.3 pkgconfig_2.0.3
[19] haven_2.5.1 GenSA_1.1.7 mvtnorm_1.1-3
[22] scales_1.2.1 webshot_0.5.4 svglite_2.1.0
[25] rjags_4-13 tzdb_0.3.0 googledrive_2.0.0
[28] farver_2.1.1 generics_0.1.3 ellipsis_0.3.2
[31] withr_2.5.0 cli_3.4.1 magrittr_2.0.3
[34] crayon_1.5.2 readxl_1.4.1 evaluate_0.17
[37] fs_1.5.2 fansi_1.0.3 parallelly_1.32.1
[40] xml2_1.3.3 tools_4.2.1 hms_1.1.2
[43] gargle_1.2.1 formatR_1.12 lifecycle_1.0.3
[46] munsell_0.5.0 reprex_2.0.2 lambda.r_1.2.4
[49] compiler_4.2.1 systemfonts_1.0.4 rlang_1.0.6
[52] futile.logger_1.4.3 grid_4.2.1 rstudioapi_0.14
[55] htmlwidgets_1.5.4 labeling_0.4.2 rmarkdown_2.20.1
[58] gtable_0.3.1 DBI_1.1.3 R6_2.5.1
[61] gridExtra_2.3 lubridate_1.8.0 fastmap_1.1.0
[64] utf8_1.2.2 futile.options_1.0.1 stringi_1.7.8
[67] parallel_4.2.1 vctrs_0.4.2 coda_0.19-4
[70] dbplyr_2.2.1 tidyselect_1.2.0 xfun_0.37.1