Show Packages
source(here::here("scripts", "load_packages.R"))
library(tidyverse)
library(rethinking)
library(cmdstanr)
library(posterior)Bayesian ORR Modeling — Categorical Dose Specification
1 Overview
This document presents the Bayesian objective response rate (ORR) analysis under the categorical dose specification, in which treatment exposure is modeled using discrete dose thresholds rather than assuming a continuous relationship. This approach allows for the possibility that response gains occur in a stepwise manner—for example, after reaching clinically meaningful exposure levels—rather than increasing uniformly with each additional dose.
All models are fit in R using rethinking::ulam(). The model structure was developed and evaluated in a separate generative modeling workflow prior to analysis of the observed data. In that stage, simulated datasets were used to assess prior behavior, confirm parameter recovery, and evaluate whether the categorical formulation could recover threshold-based exposure–response patterns under controlled conditions. The analyses presented here therefore represent the application of a pre-specified and previously validated modeling framework to the clinical cohort.
The primary adjusted ORR model estimates the association between dose category and response while accounting for key clinical covariates, including age, immunosuppression status, tumor location, stage, treatment agent, and ECOG performance status.
The categorical dose specification differs from the linear and log-dose formulations by relaxing assumptions about smooth dose–response relationships. Instead, it allows the data to inform whether response probability changes meaningfully at specific exposure levels (e.g., 0–1 doses, 2 doses, ≥3 doses). This makes it particularly useful for detecting threshold effects or diminishing returns that may not be well captured by continuous dose models.
To evaluate robustness, we also include prior-sensitivity analyses and dose-intensity sensitivity analyses derived from treatment timing and pre-response exposure. These complementary analyses assess whether observed threshold effects are stable across reasonable prior choices and alternative representations of treatment delivery.
Taken together, this document complements the primary manuscript by presenting the ORR estimates from the categorical dose model alongside the modeling assumptions, diagnostic checks, and sensitivity analyses that support their interpretation.
2 Load and Wrangle Modeling Dataset
We begin by loading the processed modeling dataset (df_model_raw), which is generated by the data-wrangling script:
dose_intensity_and_pre_response_exposure_sensitivity.R
Show R Code
load_latest_rds <- function(folder, pattern) {
file <-
list.files(folder, pattern = pattern, full.names = TRUE) |>
tibble(file = _) |>
mutate(
timestamp = stringr::str_extract(
basename(file),
"\\d{4}-\\d{2}-\\d{2}-\\d{2}-\\d{2}-\\d{2}"
),
timestamp = lubridate::ymd_hms(timestamp)
) |>
arrange(desc(timestamp)) |>
slice(1) |>
pull(file)
readRDS(file)
}
df_model_raw <- load_latest_rds(
here::here("files", "Dose Intensity Sensitivity"),
"^Dose Intensity Sensitivity Model Data.*\\.rds$"
)Show R Code
df_analysis <- df_model_raw |>
mutate(
response = case_when(
sys_resp_path %in% c("Partial response", "Major response", "Complete response") ~ 1,
sys_resp_path == "Non-responder" ~ 0,
(is.na(sys_resp_path) | sys_resp_path == "Not assessed") &
sys_resp_clin %in% c("Partial response", "Complete response") ~ 1,
TRUE ~ 0
),
dose_minus_1 = num_doses - 1,
log_dose = log(num_doses),
dose_group = if_else(num_doses <= 2, "2 doses", "3+ doses"),
dose_2 = as.integer(num_doses >= 2),
dose_3 = as.integer(num_doses >= 3),
dose_4 = as.integer(num_doses >= 4),
dose_5plus = as.integer(num_doses >= 5),
age_centered = age_at_treatment - 76,
immunosuppressed = factor(
immune_status == "Immunosuppressed",
levels = c(FALSE, TRUE)
),
location_condensed = case_when(
str_detect(primary_tumor_location, "Head|Neck|Face|Scalp|Ear") ~ "Head/Neck",
str_detect(primary_tumor_location, "Conjunctiva|Eye") ~ "Conjunctiva",
str_detect(primary_tumor_location, "Unknown") ~ "Unknown Primary",
TRUE ~ "Other"
),
location_condensed = factor(
location_condensed,
levels = c("Head/Neck", "Conjunctiva", "Other", "Unknown Primary")
),
stage_condensed = case_when(
clinical_stage %in% c("I", "II") ~ "Stage I-II",
clinical_stage == "III" ~ "Stage III",
clinical_stage == "IV" ~ "Stage IV",
TRUE ~ "Unstaged"
),
stage_condensed = factor(
stage_condensed,
levels = c("Stage I-II", "Stage III", "Stage IV", "Unstaged")
),
ecog_num = suppressWarnings(as.numeric(ecog)),
ecog_condensed = case_when(
ecog_num %in% c(0, 1) ~ "ECOG 0-1",
ecog_num %in% c(2, 3) ~ "ECOG 2-3",
TRUE ~ NA_character_
),
ecog_condensed = factor(
ecog_condensed,
levels = c("ECOG 0-1", "ECOG 2-3")
),
agent = case_when(
str_detect(sys_med, "Cemiplimab") ~ "Cemiplimab",
str_detect(sys_med, "Pembrolizumab") ~ "Pembrolizumab",
str_detect(sys_med, "Ipilimumab.*Nivolumab|Nivolumab.*Ipilimumab") ~ "Ipi+Nivo",
TRUE ~ NA_character_
),
agent = factor(
agent,
levels = c("Cemiplimab", "Pembrolizumab", "Ipi+Nivo")
),
agent_schedule = factor(
agent_schedule,
levels = c("cemiplimab", "pembro_q3", "pembro_q6", "pembro_mixed", "Ipi-nivo")
),
surgery = as.integer(has_surg)
) |>
filter(!is.na(num_doses), num_doses >= 1)Show R Code
df_analysis <- df_analysis |>
select(
patient,
MRN,
response,
num_doses,
log_dose,
dose_group,
dose_2,
dose_3,
dose_4,
dose_5plus,
age_centered,
sex,
immunosuppressed,
location_condensed,
stage_condensed,
agent,
agent_schedule,
ecog_condensed,
expected_interval_days,
likely_schedule,
num_doses_recorded,
num_doses_pre_response,
dose_intensity_final,
dose_intensity_pre_response,
dose_intensity_pre_response_uncapped,
mean_interval_days,
interval_deviation_days,
extra_doses_after_response,
dose_intensity_source,
surgery,
response_cutoff_date
)Show R Code
df_analysis <- df_analysis |>
mutate(
location_condensed = relevel(location_condensed, ref = "Head/Neck"),
stage_condensed = relevel(stage_condensed, ref = "Stage I-II"),
agent = relevel(agent, ref = "Cemiplimab"),
ecog_condensed = relevel(ecog_condensed, ref = "ECOG 0-1")
)Show R Code
df_analysis <- df_analysis |>
mutate(
ecog23 = as.integer(ecog_condensed == "ECOG 2-3"),
loc_conj = as.integer(location_condensed == "Conjunctiva"),
loc_other = as.integer(location_condensed == "Other"),
loc_unknown = as.integer(location_condensed == "Unknown Primary"),
stage_III = as.integer(stage_condensed == "Stage III"),
stage_IV = as.integer(stage_condensed == "Stage IV"),
stage_unstaged = as.integer(stage_condensed == "Unstaged"),
agent_pembro_q3 = if_else(agent_schedule == "pembro_q3", 1L, 0L, missing = 0L),
agent_pembro_q6 = if_else(agent_schedule == "pembro_q6", 1L, 0L, missing = 0L),
agent_pembro_mixed = if_else(agent_schedule == "pembro_mixed", 1L, 0L, missing = 0L),
agent_ipinivo = if_else(agent_schedule == "Ipi-nivo", 1L, 0L, missing = 0L)
) |>
filter(!is.na(agent_schedule))Show R Code
dat_cat <- list(
N = nrow(df_analysis),
y = df_analysis$response,
dose_2 = df_analysis$dose_2,
dose_3 = df_analysis$dose_3,
dose_4 = df_analysis$dose_4,
dose_5plus = df_analysis$dose_5plus,
age = df_analysis$age_centered,
imm = as.integer(df_analysis$immunosuppressed) - 1,
ecog23 = df_analysis$ecog23,
loc_conj = df_analysis$loc_conj,
loc_other = df_analysis$loc_other,
loc_unknown = df_analysis$loc_unknown,
stage_III = df_analysis$stage_III,
stage_IV = df_analysis$stage_IV,
stage_unstaged = df_analysis$stage_unstaged,
agent_pembro_q3 = df_analysis$agent_pembro_q3,
agent_pembro_q6 = df_analysis$agent_pembro_q6,
agent_pembro_mixed = df_analysis$agent_pembro_mixed,
agent_ipinivo = df_analysis$agent_ipinivo,
alpha_mu = qlogis(0.40)
)2.1 Fit Categorical Threshold Model Using Skeptical Priors
This section fits the Bayesian ORR model under the categorical threshold dose specification using skeptical priors on the dose terms. Rather than assuming that response probability changes smoothly with each additional dose, this formulation treats dose exposure as a set of discrete categories, allowing the model to detect possible threshold effects. In practical terms, this asks whether response odds differ meaningfully once patients reach specific exposure levels, such as 2 doses, 3 doses, 4 doses, or 5 or more doses, relative to the reference group.
The use of skeptical priors is important in this setting. Because categorical dose models are more flexible than linear dose models, they can more easily over-interpret random variation as apparent threshold effects. To guard against that, each dose-category coefficient is assigned a prior centered at zero with relatively tight dispersion. This expresses prior skepticism that any single threshold, on its own, produces a large shift in response probability unless the observed data provide convincing evidence. In other words, the model is allowed to detect threshold behavior, but it is discouraged from inferring abrupt dose-response jumps without support from the cohort.
As in the other ORR models, the categorical dose effects are estimated while adjusting for clinically relevant covariates, including age, immunosuppression, ECOG performance status, tumor location, disease stage, and treatment regimen. The intercept is anchored to a clinically plausible baseline response probability through alpha_mu, defined earlier in the workflow.
This model serves two purposes in the broader analysis. First, it provides a direct test of whether the exposure-response relationship may be nonlinear or threshold-based rather than approximately linear. Second, when paired with the weakly informative categorical model, it functions as a prior-sensitivity analysis, helping assess whether any apparent threshold pattern is robust or instead dependent on less regularizing prior assumptions.
# Fit Bayesian logistic regression model using a categorical threshold
# dose specification with skeptical priors.
#
# Outcome:
# y = binary response indicator
#
# Link:
# logit(p), where p is the probability of objective response
#
# Dose representation:
# Dose is modeled as a set of indicator variables rather than as a
# continuous term. This allows the model to estimate separate shifts
# in log-odds of response for specific exposure categories.
#
# Reference category:
# The omitted reference group is the lowest-dose category defined
# during data preparation (for example, 0-1 doses, depending on how
# dat_cat was constructed). Each beta_dose coefficient is therefore
# interpreted relative to that reference group.
m_cat_real_skeptical <- ulam(
alist(
# Likelihood:
# Each observed response outcome follows a Bernoulli distribution
# with patient-specific response probability p.
y ~ bernoulli(p),
# Linear predictor:
# Response probability is modeled on the log-odds scale as a function
# of dose category plus adjustment covariates.
logit(p) <- alpha +
beta_dose2 * dose_2 +
beta_dose3 * dose_3 +
beta_dose4 * dose_4 +
beta_dose5plus * dose_5plus +
beta_age * age +
beta_imm * imm +
beta_ecog * ecog23 +
beta_conj * loc_conj +
beta_other * loc_other +
beta_unknown * loc_unknown +
beta_stageIII * stage_III +
beta_stageIV * stage_IV +
beta_stageUn * stage_unstaged +
beta_pembro_q3 * agent_pembro_q3 +
beta_pembro_q6 * agent_pembro_q6 +
beta_pembro_mixed * agent_pembro_mixed +
beta_ipinivo * agent_ipinivo,
# Intercept prior:
# alpha_mu is defined earlier from the expected baseline response rate.
# This anchors the model to a clinically plausible starting probability
# before covariate adjustments are applied.
alpha ~ normal(alpha_mu, 1.5),
# Skeptical priors for categorical dose effects:
# These priors are centered at zero and fairly tight, expressing the
# belief that large threshold-specific effects should not be assumed
# unless clearly supported by the observed data.
#
# Interpretation:
# Each coefficient represents the change in log-odds of response for
# that dose category relative to the reference dose group.
beta_dose2 ~ normal(0, 0.2),
beta_dose3 ~ normal(0, 0.2),
beta_dose4 ~ normal(0, 0.2),
beta_dose5plus ~ normal(0, 0.2),
# Age effect:
# Centered at no effect, with narrow scale because very large changes
# in response odds per one-year increase in age are not expected.
beta_age ~ normal(0, 0.02),
# Immunosuppression:
# Prior centered below zero based on the expectation that
# immunosuppressed patients may have lower response probability.
beta_imm ~ normal(-0.8, 0.2),
# ECOG 2-3 versus reference:
# Mildly negative prior reflecting the possibility of worse response
# in patients with poorer performance status.
beta_ecog ~ normal(-0.3, 0.3),
# Tumor location priors:
# Conjunctival disease receives a stronger negative prior based on
# clinical expectations, whereas other and unknown locations are
# weakly regularized around no effect.
beta_conj ~ normal(-2.3, 0.5),
beta_other ~ normal(0, 0.3),
beta_unknown ~ normal(0, 0.3),
# Stage effects:
# Weakly regularizing priors centered at zero, allowing the data to
# determine whether stage shifts response probability.
beta_stageIII ~ normal(0, 0.3),
beta_stageIV ~ normal(0, 0.3),
beta_stageUn ~ normal(0, 0.3),
# Treatment regimen effects:
# Pembrolizumab schedule indicators are centered at no effect.
# Ipilimumab+nivolumab is given a mildly positive prior reflecting a
# plausible expectation of greater activity, while retaining
# substantial uncertainty.
beta_pembro_q3 ~ normal(0, 0.3),
beta_pembro_q6 ~ normal(0, 0.3),
beta_pembro_mixed ~ normal(0, 0.3),
beta_ipinivo ~ normal(0.4, 0.3)
),
# Analysis dataset prepared earlier for the categorical dose model
data = dat_cat,
# MCMC settings
chains = 4,
cores = 4,
iter = 2000
)2.2 Posterior predictive check
Show R Code
post_cat_skeptical <- extract.samples(m_cat_real_skeptical)Show R Code
# ------------------------------------------------------------------------------
# Manual posterior predictive check from fitted probabilities
# ------------------------------------------------------------------------------
manual_ppc_orr <- function(model, analysis_df, response_col = "response",
n_ppc_draws = 1000, seed = 123) {
set.seed(seed)
# Posterior fitted probabilities: usually draws x patients
p_hat <- link(model)
p_hat <- as.matrix(p_hat)
n_obs <- nrow(analysis_df)
# Defensive orientation check
if (ncol(p_hat) == n_obs) {
p_hat_checked <- p_hat
} else if (nrow(p_hat) == n_obs) {
p_hat_checked <- t(p_hat)
} else {
stop(
"Dimensions of link(model) do not match analysis_df. ",
"Expected one dimension to equal nrow(analysis_df)."
)
}
n_ppc_draws <- min(n_ppc_draws, nrow(p_hat_checked))
draw_ids <- sample(
seq_len(nrow(p_hat_checked)),
size = n_ppc_draws,
replace = FALSE
)
p_hat_ppc <- p_hat_checked[draw_ids, , drop = FALSE]
# Simulate binary outcomes from posterior fitted probabilities
yrep <- matrix(
rbinom(
n = length(p_hat_ppc),
size = 1,
prob = as.vector(p_hat_ppc)
),
nrow = nrow(p_hat_ppc),
ncol = ncol(p_hat_ppc)
)
observed_orr <- mean(analysis_df[[response_col]])
ppc_orr <- tibble(
draw = seq_len(nrow(yrep)),
simulated_orr = rowMeans(yrep)
)
ppc_summary <- ppc_orr |>
summarise(
observed_orr = observed_orr,
mean_simulated_orr = mean(simulated_orr),
median_simulated_orr = median(simulated_orr),
lower_89 = quantile(simulated_orr, 0.055),
upper_89 = quantile(simulated_orr, 0.945),
p_sim_le_observed = mean(simulated_orr <= observed_orr),
mean_fitted_probability = mean(p_hat_checked)
)
list(
p_hat = p_hat_checked,
yrep = yrep,
ppc_orr = ppc_orr,
summary = ppc_summary,
observed_orr = observed_orr
)
}Show R Code
plot_ppc_orr <- function(ppc_obj, title, subtitle = "Red dashed line = observed ORR") {
ggplot(ppc_obj$ppc_orr, aes(x = simulated_orr)) +
geom_histogram(
bins = 50,
fill = "steelblue",
alpha = 0.7,
color = "white"
) +
geom_vline(
xintercept = ppc_obj$observed_orr,
color = "red",
linewidth = 1.2,
linetype = "dashed"
) +
scale_x_continuous(
limits = c(0, 1),
labels = scales::percent_format(accuracy = 1)
) +
labs(
title = title,
subtitle = subtitle,
x = "Simulated cohort ORR",
y = "Count"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, face = "bold")
)
}Show R Code
p_cat_skeptical <- link(m_cat_real_skeptical)Show R Code
ppc_cat_skeptical <- manual_ppc_orr(
model = m_cat_real_skeptical,
analysis_df = df_analysis,
response_col = "response",
n_ppc_draws = 1000,
seed = 123
)
ppc_cat_skeptical$summary# A tibble: 1 × 7
observed_orr mean_simulated_orr median_simulated_orr lower_89 upper_89
<dbl> <dbl> <dbl> <dbl> <dbl>
1 0.646 0.641 0.640 0.566 0.714
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>Show R Code
plot_ppc_orr(
ppc_cat_skeptical,
title = "Posterior Predictive Check: Categorical Dose Model",
subtitle = "Skeptical prior; red dashed line = observed ORR"
)
2.3 Diagnostics
Show R Code
precis(m_cat_real_skeptical, depth = 2) mean sd 5.5% 94.5% rhat
alpha 0.9212924488 0.31370557 0.43564481 1.41603073 1.0001627
beta_dose2 0.0944282435 0.18632759 -0.20437697 0.39056651 1.0043255
beta_dose3 -0.0909574730 0.17662294 -0.37041050 0.19083989 1.0039848
beta_dose4 0.0763885370 0.17414012 -0.20222503 0.34949458 1.0007888
beta_dose5plus 0.1413137766 0.17746788 -0.14330677 0.42722212 1.0023395
beta_age 0.0003794618 0.01172043 -0.01845963 0.01895054 1.0026269
beta_imm -0.9039015698 0.17914548 -1.19332077 -0.61570449 1.0005856
beta_ecog -0.3105355846 0.26472746 -0.72107078 0.11286032 1.0020048
beta_conj -2.3053962997 0.44435803 -3.02299258 -1.60778603 1.0004844
beta_other -0.1374801433 0.26598282 -0.56070650 0.28213232 1.0041119
beta_unknown 0.0517497510 0.29422699 -0.41377543 0.51174581 1.0019665
beta_stageIII -0.0048443800 0.23586808 -0.38541381 0.35928007 0.9993818
beta_stageIV -0.3119673194 0.23315006 -0.67617208 0.05871795 0.9997218
beta_stageUn 0.2591618317 0.26681187 -0.16885673 0.69339521 1.0002587
beta_pembro_q3 -0.0586417628 0.22952634 -0.42999443 0.31502882 1.0013708
beta_pembro_q6 -0.0520201632 0.28289212 -0.50856005 0.40245863 1.0011811
beta_pembro_mixed 0.2536983917 0.27507505 -0.18736594 0.69330863 1.0010207
beta_ipinivo 0.4207126401 0.29037120 -0.04331428 0.88800205 1.0014758
ess_bulk
alpha 4397.255
beta_dose2 6457.648
beta_dose3 8556.934
beta_dose4 7091.864
beta_dose5plus 8993.712
beta_age 8179.245
beta_imm 8923.631
beta_ecog 7218.684
beta_conj 8472.089
beta_other 9532.843
beta_unknown 9024.214
beta_stageIII 6116.512
beta_stageIV 5760.632
beta_stageUn 7292.570
beta_pembro_q3 7755.704
beta_pembro_q6 8304.800
beta_pembro_mixed 9527.818
beta_ipinivo 8337.355Show R Code
plot(precis(m_cat_real_skeptical))
2.4 Baseline Probability
Show R Code
baseline_prob_cat_skeptical <- plogis(post_cat_skeptical$alpha)
mean(baseline_prob_cat_skeptical)[1] 0.7111326Show R Code
quantile(baseline_prob_cat_skeptical, c(0.055, 0.945)) 5.5% 94.5%
0.6072208 0.8047154 2.5 Posterior predicted dose-response curve
Show R Code
dose_number <- 1:15
new_data_cat <- tibble(
dose_number = dose_number,
dose_2 = as.integer(dose_number >= 2),
dose_3 = as.integer(dose_number >= 3),
dose_4 = as.integer(dose_number >= 4),
dose_5plus = as.integer(dose_number >= 5),
age = 0,
imm = 0,
ecog23 = 0,
loc_conj = 0,
loc_other = 0,
loc_unknown = 0,
stage_III = 0,
stage_IV = 0,
stage_unstaged = 0,
agent_pembro_q3 = 0,
agent_pembro_q6 = 0,
agent_pembro_mixed = 0,
agent_ipinivo = 0
)
p_cat_skeptical_new <- link(m_cat_real_skeptical, data = new_data_cat)
mean_p_cat_skeptical <- apply(p_cat_skeptical_new, 2, mean)
ci_cat_skeptical <- apply(p_cat_skeptical_new, 2, PI, 0.89)Show R Code
plot_df_cat_skeptical <- tibble(
dose_number = dose_number,
mean_p = mean_p_cat_skeptical,
lower = ci_cat_skeptical[1, ],
upper = ci_cat_skeptical[2, ]
)Show R Code
plot_predicted_response_cat_skeptical <-
ggplot(plot_df_cat_skeptical, aes(x = dose_number, y = mean_p)) +
geom_ribbon(
aes(ymin = lower, ymax = upper),
fill = "steelblue",
alpha = 0.22
) +
geom_step(linewidth = 1.2, color = "steelblue4") +
geom_point(
size = 2.8,
shape = 21,
stroke = 0.9,
fill = "white",
color = "steelblue4"
) +
scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, by = 0.1),
labels = scales::percent_format(accuracy = 1),
expand = expansion(mult = c(0.01, 0.02))
) +
scale_x_continuous(
breaks = 1:15,
expand = expansion(mult = c(0.01, 0.02))
) +
labs(
title = "Predicted Response Probability by Number of Doses",
subtitle = "Posterior mean predictions from Bayesian logistic regression\nCategorical threshold model with skeptical prior; shaded region = 89% CrI",
x = "Number of doses",
y = "Predicted response probability"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
plot_predicted_response_cat_skeptical
2.6 Extract posterior threshold effects
Show R Code
posterior_cat_skeptical <- tibble(
beta_dose2 = post_cat_skeptical$beta_dose2,
beta_dose3 = post_cat_skeptical$beta_dose3,
beta_dose4 = post_cat_skeptical$beta_dose4,
beta_dose5plus = post_cat_skeptical$beta_dose5plus
)Show R Code
threshold_summary_cat_skeptical <- tibble(
threshold = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4"),
mean_beta = c(
mean(post_cat_skeptical$beta_dose2),
mean(post_cat_skeptical$beta_dose3),
mean(post_cat_skeptical$beta_dose4),
mean(post_cat_skeptical$beta_dose5plus)
),
median_beta = c(
median(post_cat_skeptical$beta_dose2),
median(post_cat_skeptical$beta_dose3),
median(post_cat_skeptical$beta_dose4),
median(post_cat_skeptical$beta_dose5plus)
),
ci_lower_89 = c(
quantile(post_cat_skeptical$beta_dose2, 0.055),
quantile(post_cat_skeptical$beta_dose3, 0.055),
quantile(post_cat_skeptical$beta_dose4, 0.055),
quantile(post_cat_skeptical$beta_dose5plus, 0.055)
),
ci_upper_89 = c(
quantile(post_cat_skeptical$beta_dose2, 0.945),
quantile(post_cat_skeptical$beta_dose3, 0.945),
quantile(post_cat_skeptical$beta_dose4, 0.945),
quantile(post_cat_skeptical$beta_dose5plus, 0.945)
),
p_gt_0 = c(
mean(post_cat_skeptical$beta_dose2 > 0),
mean(post_cat_skeptical$beta_dose3 > 0),
mean(post_cat_skeptical$beta_dose4 > 0),
mean(post_cat_skeptical$beta_dose5plus > 0)
)
)
threshold_summary_cat_skeptical# A tibble: 4 × 6
threshold mean_beta median_beta ci_lower_89 ci_upper_89 p_gt_0
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Dose 2 vs 1 0.0944 0.0915 -0.204 0.391 0.701
2 Dose 3 vs 2 -0.0910 -0.0898 -0.370 0.191 0.306
3 Dose 4 vs 3 0.0764 0.0744 -0.202 0.349 0.669
4 Dose 5+ vs 4 0.141 0.140 -0.143 0.427 0.788Threshold-effect plot
Show R Code
plot_cat_skeptical_df <-
bind_rows(
tibble(beta = post_cat_skeptical$beta_dose2, threshold = "Dose 2 vs 1"),
tibble(beta = post_cat_skeptical$beta_dose3, threshold = "Dose 3 vs 2"),
tibble(beta = post_cat_skeptical$beta_dose4, threshold = "Dose 4 vs 3"),
tibble(beta = post_cat_skeptical$beta_dose5plus, threshold = "Dose 5+ vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4")
)
)
threshold_summary_plot_cat_skeptical <-
plot_cat_skeptical_df |>
group_by(threshold) |>
summarise(
mean_beta = mean(beta),
lower = quantile(beta, 0.055),
upper = quantile(beta, 0.945),
p_gt_0 = mean(beta > 0),
.groups = "drop"
)
plot_cat_skeptical <-
ggplot(threshold_summary_plot_cat_skeptical,
aes(x = threshold, y = mean_beta)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_pointrange(
aes(ymin = lower, ymax = upper),
linewidth = 0.9
) +
geom_text(
aes(
y = upper + 0.05,
label = paste0("P(>0) = ", round(p_gt_0 * 100, 0), "%")
),
size = 3.5
) +
labs(
title = "Posterior Threshold Effects: Categorical Dose Model",
subtitle = "Skeptical prior; points = posterior means, bars = 89% CrI",
x = NULL,
y = "Log-odds increment at threshold"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.text.x = element_text(face = "bold")
)
plot_cat_skeptical
2.7 Posterior distributions: skeptical categorical model
Show R Code
posterior_cat_skeptical_long <- bind_rows(
tibble(beta = post_cat_skeptical$beta_dose2, threshold = "Dose 2 vs 1"),
tibble(beta = post_cat_skeptical$beta_dose3, threshold = "Dose 3 vs 2"),
tibble(beta = post_cat_skeptical$beta_dose4, threshold = "Dose 4 vs 3"),
tibble(beta = post_cat_skeptical$beta_dose5plus, threshold = "Dose 5+ vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4")
)
)
cat_skeptical_summary <- posterior_cat_skeptical_long |>
group_by(threshold) |>
summarise(
mean_beta = mean(beta),
lower = quantile(beta, 0.055),
upper = quantile(beta, 0.945),
p_gt_0 = mean(beta > 0),
.groups = "drop"
)
plot_cat_density_skeptical <-
ggplot(posterior_cat_skeptical_long, aes(x = beta)) +
geom_density(fill = "steelblue", alpha = 0.35, color = "steelblue4", linewidth = 0.9) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 0.9) +
geom_vline(
data = cat_skeptical_summary,
aes(xintercept = lower),
linetype = "dotted",
color = "grey40",
linewidth = 0.7
) +
geom_vline(
data = cat_skeptical_summary,
aes(xintercept = upper),
linetype = "dotted",
color = "grey40",
linewidth = 0.7
) +
geom_text(
data = cat_skeptical_summary,
aes(
x = Inf,
y = Inf,
label = paste0(
"89% CrI: [", round(lower, 2), ", ", round(upper, 2), "]\n",
"P(β > 0) = ", round(p_gt_0 * 100, 0), "%"
)
),
hjust = 1.05,
vjust = 1.5,
size = 3.3,
inherit.aes = FALSE
) +
facet_wrap(~ threshold, scales = "free_y") +
labs(
title = "Posterior Distributions: Categorical Threshold Effects",
subtitle = "Skeptical prior",
x = "Threshold-specific log-odds effect",
y = "Density"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
strip.text = element_text(face = "bold")
)
plot_cat_density_skeptical
Show R Code
fmt_pct <- function(x) sprintf("%.0f%%", 100 * x)
adjust_k <- 1.1
cols_cat <- c(
"Dose 2 vs 1" = "#1F78B4",
"Dose 3 vs 2" = "#1B9E77",
"Dose 4 vs 3" = "#D95F02",
"Dose ≥5 vs 4" = "#7B2CBF"
)Show R Code
# ---- Combine posteriors + ordering ----
posterior_cat_combined_skeptical <-
bind_rows(
tibble(beta = post_cat_skeptical$beta_dose2, threshold = "Dose 2 vs 1"),
tibble(beta = post_cat_skeptical$beta_dose3, threshold = "Dose 3 vs 2"),
tibble(beta = post_cat_skeptical$beta_dose4, threshold = "Dose 4 vs 3"),
tibble(beta = post_cat_skeptical$beta_dose5plus, threshold = "Dose ≥5 vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose ≥5 vs 4")
)
)
# ---- Colors ----
cols_cat <- c(
"Dose 2 vs 1" = "#1F78B4",
"Dose 3 vs 2" = "#1B9E77",
"Dose 4 vs 3" = "#D95F02",
"Dose ≥5 vs 4" = "#7B2CBF"
)
# ---- Summaries for annotation ----
summ_cat_skeptical <-
posterior_cat_combined_skeptical |>
group_by(threshold) |>
summarise(
p_gt0 = mean(beta > 0),
OR_md = exp(median(beta)),
.groups = "drop"
)
fmt_pct <- function(x) sprintf("%.1f%%", 100 * x)
ann_text_skeptical <- paste0(
"Posterior Pr(β > 0):\n",
"Dose 2 vs 1: ", fmt_pct(summ_cat_skeptical$p_gt0[summ_cat_skeptical$threshold == "Dose 2 vs 1"]), "\n",
"Dose 3 vs 2: ", fmt_pct(summ_cat_skeptical$p_gt0[summ_cat_skeptical$threshold == "Dose 3 vs 2"]), "\n",
"Dose 4 vs 3: ", fmt_pct(summ_cat_skeptical$p_gt0[summ_cat_skeptical$threshold == "Dose 4 vs 3"]), "\n",
"Dose ≥5 vs 4: ", fmt_pct(summ_cat_skeptical$p_gt0[summ_cat_skeptical$threshold == "Dose ≥5 vs 4"]), "\n\n",
"Median OR (incremental):\n",
"Dose 2 vs 1: ×", sprintf("%.2f", summ_cat_skeptical$OR_md[summ_cat_skeptical$threshold == "Dose 2 vs 1"]), "\n",
"Dose 3 vs 2: ×", sprintf("%.2f", summ_cat_skeptical$OR_md[summ_cat_skeptical$threshold == "Dose 3 vs 2"]), "\n",
"Dose 4 vs 3: ×", sprintf("%.2f", summ_cat_skeptical$OR_md[summ_cat_skeptical$threshold == "Dose 4 vs 3"]), "\n",
"Dose ≥5 vs 4: ×", sprintf("%.2f", summ_cat_skeptical$OR_md[summ_cat_skeptical$threshold == "Dose ≥5 vs 4"])
)
# ---- Compute ymax ----
adjust_k <- 1.1
dens_list_skeptical <- posterior_cat_combined_skeptical |>
group_by(threshold) |>
summarise(
dens = list(stats::density(beta, adjust = adjust_k)),
.groups = "drop"
)
ymax_skeptical <- max(vapply(dens_list_skeptical$dens, function(d) max(d$y), numeric(1)))
plot_categorical_skeptical <-
ggplot(posterior_cat_combined_skeptical, aes(x = beta)) +
geom_density(aes(fill = threshold), alpha = 0.18, color = NA, adjust = 1.1) +
geom_density(aes(color = threshold), fill = NA, linewidth = 1.15, adjust = 1.1) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
scale_fill_manual(values = cols_cat) +
scale_color_manual(values = cols_cat) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
scale_y_continuous(expand = expansion(mult = c(0.02, 0.02))) +
annotate(
"label",
x = 1.25,
y = ymax_skeptical,
label = ann_text_skeptical,
hjust = 0, vjust = 1,
size = 3.1,
label.size = 0.25,
fill = "white", alpha = 0.92
) +
labs(
title = "Posterior Distributions: Categorical Dose Threshold Effects",
subtitle = "Skeptical Prior",
x = "β (log-odds; incremental effect of crossing each threshold)",
y = "Density",
color = "Increment",
fill = "Increment"
) +
coord_cartesian(xlim = c(NA, 2.1), clip = "off") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, face = "bold"),
legend.position = "bottom",
plot.margin = margin(t = 10, r = 40, b = 10, l = 10)
)
plot_categorical_skeptical
2.8 Fit model with weakly informative prior
m_cat_real_weak <- ulam(
alist(
y ~ bernoulli(p),
logit(p) <- alpha +
beta_dose2 * dose_2 +
beta_dose3 * dose_3 +
beta_dose4 * dose_4 +
beta_dose5plus * dose_5plus +
beta_age * age +
beta_imm * imm +
beta_ecog * ecog23 +
beta_conj * loc_conj +
beta_other * loc_other +
beta_unknown * loc_unknown +
beta_stageIII * stage_III +
beta_stageIV * stage_IV +
beta_stageUn * stage_unstaged +
beta_pembro_q3 * agent_pembro_q3 +
beta_pembro_q6 * agent_pembro_q6 +
beta_pembro_mixed * agent_pembro_mixed +
beta_ipinivo * agent_ipinivo,
alpha ~ normal(alpha_mu, 1.5),
beta_dose2 ~ normal(0.25, 0.25),
beta_dose3 ~ normal(0.15, 0.25),
beta_dose4 ~ normal(0.10, 0.25),
beta_dose5plus ~ normal(0.05, 0.25),
beta_age ~ normal(0, 0.02),
beta_imm ~ normal(-0.8, 0.2),
beta_ecog ~ normal(-0.3, 0.3),
beta_conj ~ normal(-2.3, 0.5),
beta_other ~ normal(0, 0.3),
beta_unknown ~ normal(0, 0.3),
beta_stageIII ~ normal(0, 0.3),
beta_stageIV ~ normal(0, 0.3),
beta_stageUn ~ normal(0, 0.3),
beta_pembro_q3 ~ normal(0, 0.3),
beta_pembro_q6 ~ normal(0, 0.3),
beta_pembro_mixed ~ normal(0, 0.3),
beta_ipinivo ~ normal(0.4, 0.3)
),
data = dat_cat,
chains = 4,
cores = 4,
iter = 2000
)2.9 Posterior predictive check of model with weak prior
Show R Code
post_cat_weak <- extract.samples(m_cat_real_weak)Show R Code
p_cat_weak <- link(m_cat_real_weak)Show R Code
ppc_cat_weak <- manual_ppc_orr(
model = m_cat_real_weak,
analysis_df = df_analysis,
response_col = "response",
n_ppc_draws = 1000,
seed = 123
)
ppc_cat_weak$summary# A tibble: 1 × 7
observed_orr mean_simulated_orr median_simulated_orr lower_89 upper_89
<dbl> <dbl> <dbl> <dbl> <dbl>
1 0.646 0.643 0.646 0.566 0.714
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>Show R Code
plot_ppc_orr(
ppc_cat_weak,
title = "Posterior Predictive Check: Categorical Dose Model",
subtitle = "Weakly informative prior; red dashed line = observed ORR"
)
2.10 Posterior parameter interpretation
Show R Code
precis(m_cat_real_weak, depth = 2) mean sd 5.5% 94.5% rhat
alpha 0.641885324 0.33965084 0.09541816 1.18023375 1.0016023
beta_dose2 0.354882797 0.23326749 -0.02315859 0.73342757 0.9996452
beta_dose3 -0.068141054 0.20929002 -0.40637631 0.27619348 1.0017430
beta_dose4 0.136321391 0.21097370 -0.20023085 0.47234156 1.0000074
beta_dose5plus 0.200572739 0.21483659 -0.13802426 0.55200393 0.9997660
beta_age 0.001351118 0.01136991 -0.01647778 0.01955368 1.0007983
beta_imm -0.908659444 0.17490247 -1.18172778 -0.62867420 1.0034165
beta_ecog -0.309262227 0.26145953 -0.72028226 0.11497642 1.0002346
beta_conj -2.311931021 0.42858545 -2.98553693 -1.63249015 1.0005781
beta_other -0.138734676 0.26767774 -0.56162062 0.28952495 1.0011760
beta_unknown 0.054923300 0.27945309 -0.38764670 0.49321553 0.9996675
beta_stageIII -0.012655996 0.23918689 -0.39525832 0.37956691 1.0005987
beta_stageIV -0.322060196 0.23655158 -0.69376606 0.06086378 1.0012935
beta_stageUn 0.251976652 0.26458921 -0.16710190 0.66668997 1.0010447
beta_pembro_q3 -0.072870688 0.23443051 -0.45553127 0.29923232 1.0003603
beta_pembro_q6 -0.051053263 0.28229687 -0.49700621 0.39166873 1.0018507
beta_pembro_mixed 0.238899247 0.27750463 -0.21030381 0.67471893 1.0004910
beta_ipinivo 0.426045087 0.30622075 -0.07120204 0.91587557 1.0022985
ess_bulk
alpha 3809.084
beta_dose2 4966.564
beta_dose3 5977.923
beta_dose4 6254.582
beta_dose5plus 7527.020
beta_age 8984.921
beta_imm 7140.958
beta_ecog 7307.466
beta_conj 5322.043
beta_other 8216.679
beta_unknown 8001.107
beta_stageIII 5579.720
beta_stageIV 4924.334
beta_stageUn 6609.002
beta_pembro_q3 7280.697
beta_pembro_q6 7157.381
beta_pembro_mixed 6738.184
beta_ipinivo 8079.685Show R Code
plot(precis(m_cat_real_weak))
Posterior predicted dose-response curve
Show R Code
p_cat_weak_new <- link(m_cat_real_weak, data = new_data_cat)
mean_p_cat_weak <- apply(p_cat_weak_new, 2, mean)
ci_cat_weak <- apply(p_cat_weak_new, 2, PI, 0.89)Show R Code
plot_df_cat_weak <- tibble(
dose_number = dose_number,
mean_p = mean_p_cat_weak,
lower = ci_cat_weak[1, ],
upper = ci_cat_weak[2, ]
)Show R Code
plot_predicted_response_cat_weak <-
ggplot(plot_df_cat_weak, aes(x = dose_number, y = mean_p)) +
geom_ribbon(
aes(ymin = lower, ymax = upper),
fill = "steelblue",
alpha = 0.22
) +
geom_step(linewidth = 1.2, color = "steelblue4") +
geom_point(
size = 2.8,
shape = 21,
stroke = 0.9,
fill = "white",
color = "steelblue4"
) +
scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, by = 0.1),
labels = scales::percent_format(accuracy = 1),
expand = expansion(mult = c(0.01, 0.02))
) +
scale_x_continuous(
breaks = 1:15,
expand = expansion(mult = c(0.01, 0.02))
) +
labs(
title = "Predicted Response Probability by Number of Doses",
subtitle = "Posterior mean predictions from Bayesian logistic regression\nCategorical threshold model with weakly informative prior; shaded region = 89% CrI",
x = "Number of doses",
y = "Predicted response probability"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
plot_predicted_response_cat_weak
2.11 Extract and summarize posterior threshold effects
Show R Code
posterior_cat_weak <- tibble(
beta_dose2 = post_cat_weak$beta_dose2,
beta_dose3 = post_cat_weak$beta_dose3,
beta_dose4 = post_cat_weak$beta_dose4,
beta_dose5plus = post_cat_weak$beta_dose5plus
)Show R Code
threshold_summary_cat_weak <- tibble(
threshold = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4"),
mean_beta = c(
mean(post_cat_weak$beta_dose2),
mean(post_cat_weak$beta_dose3),
mean(post_cat_weak$beta_dose4),
mean(post_cat_weak$beta_dose5plus)
),
median_beta = c(
median(post_cat_weak$beta_dose2),
median(post_cat_weak$beta_dose3),
median(post_cat_weak$beta_dose4),
median(post_cat_weak$beta_dose5plus)
),
ci_lower_89 = c(
quantile(post_cat_weak$beta_dose2, 0.055),
quantile(post_cat_weak$beta_dose3, 0.055),
quantile(post_cat_weak$beta_dose4, 0.055),
quantile(post_cat_weak$beta_dose5plus, 0.055)
),
ci_upper_89 = c(
quantile(post_cat_weak$beta_dose2, 0.945),
quantile(post_cat_weak$beta_dose3, 0.945),
quantile(post_cat_weak$beta_dose4, 0.945),
quantile(post_cat_weak$beta_dose5plus, 0.945)
),
p_gt_0 = c(
mean(post_cat_weak$beta_dose2 > 0),
mean(post_cat_weak$beta_dose3 > 0),
mean(post_cat_weak$beta_dose4 > 0),
mean(post_cat_weak$beta_dose5plus > 0)
)
)
threshold_summary_cat_weak# A tibble: 4 × 6
threshold mean_beta median_beta ci_lower_89 ci_upper_89 p_gt_0
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Dose 2 vs 1 0.355 0.353 -0.0232 0.733 0.936
2 Dose 3 vs 2 -0.0681 -0.0654 -0.406 0.276 0.368
3 Dose 4 vs 3 0.136 0.137 -0.200 0.472 0.742
4 Dose 5+ vs 4 0.201 0.199 -0.138 0.552 0.819Show R Code
plot_cat_weak_df <-
bind_rows(
tibble(beta = post_cat_weak$beta_dose2, threshold = "Dose 2 vs 1"),
tibble(beta = post_cat_weak$beta_dose3, threshold = "Dose 3 vs 2"),
tibble(beta = post_cat_weak$beta_dose4, threshold = "Dose 4 vs 3"),
tibble(beta = post_cat_weak$beta_dose5plus, threshold = "Dose 5+ vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4")
)
)
threshold_summary_plot_cat_weak <-
plot_cat_weak_df |>
group_by(threshold) |>
summarise(
mean_beta = mean(beta),
lower = quantile(beta, 0.055),
upper = quantile(beta, 0.945),
p_gt_0 = mean(beta > 0),
.groups = "drop"
)
plot_cat_weak <-
ggplot(threshold_summary_plot_cat_weak,
aes(x = threshold, y = mean_beta)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_pointrange(
aes(ymin = lower, ymax = upper),
linewidth = 0.9
) +
geom_text(
aes(
y = upper + 0.05,
label = paste0("P(>0) = ", round(p_gt_0 * 100, 0), "%")
),
size = 3.5
) +
labs(
title = "Posterior Threshold Effects: Categorical Dose Model",
subtitle = "Weakly informative prior; points = posterior means, bars = 89% CrI",
x = NULL,
y = "Log-odds increment at threshold"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.text.x = element_text(face = "bold")
)
plot_cat_weak
2.12 Compare categorical threshold effects across priors
Show R Code
prior_sensitivity_cat <- tibble(
threshold = rep(c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4"), 2),
prior = rep(c("Skeptical", "Weakly informative"), each = 4),
mean_beta = c(
mean(post_cat_skeptical$beta_dose2),
mean(post_cat_skeptical$beta_dose3),
mean(post_cat_skeptical$beta_dose4),
mean(post_cat_skeptical$beta_dose5plus),
mean(post_cat_weak$beta_dose2),
mean(post_cat_weak$beta_dose3),
mean(post_cat_weak$beta_dose4),
mean(post_cat_weak$beta_dose5plus)
),
median_beta = c(
median(post_cat_skeptical$beta_dose2),
median(post_cat_skeptical$beta_dose3),
median(post_cat_skeptical$beta_dose4),
median(post_cat_skeptical$beta_dose5plus),
median(post_cat_weak$beta_dose2),
median(post_cat_weak$beta_dose3),
median(post_cat_weak$beta_dose4),
median(post_cat_weak$beta_dose5plus)
),
ci_lower_89 = c(
quantile(post_cat_skeptical$beta_dose2, 0.055),
quantile(post_cat_skeptical$beta_dose3, 0.055),
quantile(post_cat_skeptical$beta_dose4, 0.055),
quantile(post_cat_skeptical$beta_dose5plus, 0.055),
quantile(post_cat_weak$beta_dose2, 0.055),
quantile(post_cat_weak$beta_dose3, 0.055),
quantile(post_cat_weak$beta_dose4, 0.055),
quantile(post_cat_weak$beta_dose5plus, 0.055)
),
ci_upper_89 = c(
quantile(post_cat_skeptical$beta_dose2, 0.945),
quantile(post_cat_skeptical$beta_dose3, 0.945),
quantile(post_cat_skeptical$beta_dose4, 0.945),
quantile(post_cat_skeptical$beta_dose5plus, 0.945),
quantile(post_cat_weak$beta_dose2, 0.945),
quantile(post_cat_weak$beta_dose3, 0.945),
quantile(post_cat_weak$beta_dose4, 0.945),
quantile(post_cat_weak$beta_dose5plus, 0.945)
),
p_gt_0 = c(
mean(post_cat_skeptical$beta_dose2 > 0),
mean(post_cat_skeptical$beta_dose3 > 0),
mean(post_cat_skeptical$beta_dose4 > 0),
mean(post_cat_skeptical$beta_dose5plus > 0),
mean(post_cat_weak$beta_dose2 > 0),
mean(post_cat_weak$beta_dose3 > 0),
mean(post_cat_weak$beta_dose4 > 0),
mean(post_cat_weak$beta_dose5plus > 0)
)
)
prior_sensitivity_cat# A tibble: 8 × 7
threshold prior mean_beta median_beta ci_lower_89 ci_upper_89 p_gt_0
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Dose 2 vs 1 Skeptical 0.0944 0.0915 -0.204 0.391 0.701
2 Dose 3 vs 2 Skeptical -0.0910 -0.0898 -0.370 0.191 0.306
3 Dose 4 vs 3 Skeptical 0.0764 0.0744 -0.202 0.349 0.669
4 Dose 5+ vs 4 Skeptical 0.141 0.140 -0.143 0.427 0.788
5 Dose 2 vs 1 Weakly info… 0.355 0.353 -0.0232 0.733 0.936
6 Dose 3 vs 2 Weakly info… -0.0681 -0.0654 -0.406 0.276 0.368
7 Dose 4 vs 3 Weakly info… 0.136 0.137 -0.200 0.472 0.742
8 Dose 5+ vs 4 Weakly info… 0.201 0.199 -0.138 0.552 0.819Show R Code
posterior_cat_weak_long <- bind_rows(
tibble(beta = post_cat_weak$beta_dose2, threshold = "Dose 2 vs 1"),
tibble(beta = post_cat_weak$beta_dose3, threshold = "Dose 3 vs 2"),
tibble(beta = post_cat_weak$beta_dose4, threshold = "Dose 4 vs 3"),
tibble(beta = post_cat_weak$beta_dose5plus, threshold = "Dose 5+ vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4")
)
)
cat_weak_summary <- posterior_cat_weak_long |>
group_by(threshold) |>
summarise(
mean_beta = mean(beta),
lower = quantile(beta, 0.055),
upper = quantile(beta, 0.945),
p_gt_0 = mean(beta > 0),
.groups = "drop"
)
plot_cat_density_weak <-
ggplot(posterior_cat_weak_long, aes(x = beta)) +
geom_density(fill = "steelblue", alpha = 0.35, color = "steelblue4", linewidth = 0.9) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 0.9) +
geom_vline(
data = cat_weak_summary,
aes(xintercept = lower),
linetype = "dotted",
color = "grey40",
linewidth = 0.7
) +
geom_vline(
data = cat_weak_summary,
aes(xintercept = upper),
linetype = "dotted",
color = "grey40",
linewidth = 0.7
) +
geom_text(
data = cat_weak_summary,
aes(
x = Inf,
y = Inf,
label = paste0(
"89% CrI: [", round(lower, 2), ", ", round(upper, 2), "]\n",
"P(β > 0) = ", round(p_gt_0 * 100, 0), "%"
)
),
hjust = 1.05,
vjust = 1.5,
size = 3.3,
inherit.aes = FALSE
) +
facet_wrap(~ threshold, scales = "free_y") +
labs(
title = "Posterior Distributions: Categorical Threshold Effects",
subtitle = "Weakly informative prior",
x = "Threshold-specific log-odds effect",
y = "Density"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
strip.text = element_text(face = "bold")
)
plot_cat_density_weak
Show R Code
posterior_cat_combined_weak <-
bind_rows(
tibble(beta = post_cat_weak$beta_dose2, threshold = "Dose 2 vs 1"),
tibble(beta = post_cat_weak$beta_dose3, threshold = "Dose 3 vs 2"),
tibble(beta = post_cat_weak$beta_dose4, threshold = "Dose 4 vs 3"),
tibble(beta = post_cat_weak$beta_dose5plus, threshold = "Dose ≥5 vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose ≥5 vs 4")
)
)
summ_cat_weak <-
posterior_cat_combined_weak |>
group_by(threshold) |>
summarise(
p_gt0 = mean(beta > 0),
OR_md = exp(median(beta)),
.groups = "drop"
)
ann_text_weak <- paste0(
"Posterior Pr(β > 0):\n",
"Dose 2 vs 1: ", fmt_pct(summ_cat_weak$p_gt0[summ_cat_weak$threshold == "Dose 2 vs 1"]), "\n",
"Dose 3 vs 2: ", fmt_pct(summ_cat_weak$p_gt0[summ_cat_weak$threshold == "Dose 3 vs 2"]), "\n",
"Dose 4 vs 3: ", fmt_pct(summ_cat_weak$p_gt0[summ_cat_weak$threshold == "Dose 4 vs 3"]), "\n",
"Dose ≥5 vs 4: ", fmt_pct(summ_cat_weak$p_gt0[summ_cat_weak$threshold == "Dose ≥5 vs 4"]), "\n\n",
"Median OR (incremental):\n",
"Dose 2 vs 1: ×", sprintf("%.2f", summ_cat_weak$OR_md[summ_cat_weak$threshold == "Dose 2 vs 1"]), "\n",
"Dose 3 vs 2: ×", sprintf("%.2f", summ_cat_weak$OR_md[summ_cat_weak$threshold == "Dose 3 vs 2"]), "\n",
"Dose 4 vs 3: ×", sprintf("%.2f", summ_cat_weak$OR_md[summ_cat_weak$threshold == "Dose 4 vs 3"]), "\n",
"Dose ≥5 vs 4: ×", sprintf("%.2f", summ_cat_weak$OR_md[summ_cat_weak$threshold == "Dose ≥5 vs 4"])
)
dens_list_weak <- posterior_cat_combined_weak |>
group_by(threshold) |>
summarise(
dens = list(stats::density(beta, adjust = adjust_k)),
.groups = "drop"
)
ymax_weak <- max(vapply(dens_list_weak$dens, function(d) max(d$y), numeric(1)))
plot_categorical_weak <-
ggplot(posterior_cat_combined_weak, aes(x = beta)) +
geom_density(aes(fill = threshold), alpha = 0.18, color = NA, adjust = 1.1) +
geom_density(aes(color = threshold), fill = NA, linewidth = 1.15, adjust = 1.1) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
scale_fill_manual(values = cols_cat) +
scale_color_manual(values = cols_cat) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
scale_y_continuous(expand = expansion(mult = c(0.02, 0.02))) +
annotate(
"label",
x = 1.25,
y = ymax_weak,
label = ann_text_weak,
hjust = 0, vjust = 1,
size = 3.1,
label.size = 0.25,
fill = "white", alpha = 0.92
) +
labs(
title = "Posterior Distributions: Categorical Dose Threshold Effects",
subtitle = "Weakly Informative Prior",
x = "β (log-odds; incremental effect of crossing each threshold)",
y = "Density",
color = "Increment",
fill = "Increment"
) +
coord_cartesian(xlim = c(NA, 2.1), clip = "off") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, face = "bold"),
legend.position = "bottom",
plot.margin = margin(t = 10, r = 40, b = 10, l = 10)
)
plot_categorical_weak
3 Dose-Intensity–Adjusted Categorical Model
This section extends the categorical threshold dose model to incorporate dose-intensity–adjusted exposure, using the same conceptual framework introduced in the linear and log-dose sensitivity analysis. Rather than counting all administered doses, exposure is redefined to include only those doses delivered prior to the clinical response timepoint, thereby reducing the risk of reverse causation (i.e., patients receiving more doses simply because they have already responded).
As in the primary categorical model, treatment exposure is represented using discrete dose categories. However, these categories are now constructed from pre-response dose counts, yielding a time-aware representation of exposure that more closely reflects the causal question of interest: whether greater treatment intensity before response is associated with a higher probability of response.
The model specification differs from the log-dose version only in how dose is parameterized. Instead of a single continuous dose term, the model includes a set of indicator variables corresponding to clinically interpretable exposure thresholds:
- 2 doses
- 3 doses
- 4 doses
- ≥5 doses
Each coefficient therefore represents the change in log-odds of response for that exposure category relative to the reference group (typically the lowest exposure category).
This formulation preserves the flexibility of the categorical approach—allowing for nonlinear or threshold effects—while incorporating a more temporally aligned definition of treatment exposure. As such, it serves as a key sensitivity analysis to evaluate whether previously observed dose–response patterns are robust to a more causally coherent definition of dose intensity.
Show R Code
df_analysis <- df_analysis |>
mutate(
dose_intensity_pre_response_uncapped_centered =
dose_intensity_pre_response_uncapped -
mean(dose_intensity_pre_response_uncapped, na.rm = TRUE)
)
df_analysis |>
summarise(
n_non_missing = sum(!is.na(dose_intensity_pre_response_uncapped)),
mean_raw = mean(dose_intensity_pre_response_uncapped, na.rm = TRUE),
sd_raw = sd(dose_intensity_pre_response_uncapped, na.rm = TRUE)
)# A tibble: 1 × 3
n_non_missing mean_raw sd_raw
<int> <dbl> <dbl>
1 189 0.977 0.113Show R Code
df_analysis_cat_di <- df_analysis |>
filter(
!is.na(response),
!is.na(dose_2),
!is.na(dose_3),
!is.na(dose_4),
!is.na(dose_5plus),
!is.na(age_centered),
!is.na(immunosuppressed),
!is.na(ecog23),
!is.na(loc_conj),
!is.na(loc_other),
!is.na(loc_unknown),
!is.na(stage_III),
!is.na(stage_IV),
!is.na(stage_unstaged),
!is.na(agent_pembro_q3),
!is.na(agent_pembro_q6),
!is.na(agent_pembro_mixed),
!is.na(agent_ipinivo),
!is.na(dose_intensity_pre_response_uncapped_centered)
)Show R Code
df_analysis_cat_di |>
summarise(
n = n(),
mean_response = mean(response),
mean_di = mean(dose_intensity_pre_response_uncapped_centered),
sd_di = sd(dose_intensity_pre_response_uncapped_centered),
min_di = min(dose_intensity_pre_response_uncapped_centered),
max_di = max(dose_intensity_pre_response_uncapped_centered)
)# A tibble: 1 × 6
n mean_response mean_di sd_di min_di max_di
<int> <dbl> <dbl> <dbl> <dbl> <dbl>
1 189 0.646 -3.81e-17 0.113 -0.624 0.3563.1 Build the data list
Show R Code
dat_cat_di_uncapped <- list(
N = nrow(df_analysis_cat_di),
y = df_analysis_cat_di$response,
d2 = df_analysis_cat_di$dose_2,
d3 = df_analysis_cat_di$dose_3,
d4 = df_analysis_cat_di$dose_4,
d5 = df_analysis_cat_di$dose_5plus,
di = df_analysis_cat_di$dose_intensity_pre_response_uncapped_centered,
age = df_analysis_cat_di$age_centered,
imm = as.integer(df_analysis_cat_di$immunosuppressed) - 1,
ecog = df_analysis_cat_di$ecog23,
lc = df_analysis_cat_di$loc_conj,
lo = df_analysis_cat_di$loc_other,
lu = df_analysis_cat_di$loc_unknown,
s3 = df_analysis_cat_di$stage_III,
s4 = df_analysis_cat_di$stage_IV,
su = df_analysis_cat_di$stage_unstaged,
aq3 = df_analysis_cat_di$agent_pembro_q3,
aq6 = df_analysis_cat_di$agent_pembro_q6,
amx = df_analysis_cat_di$agent_pembro_mixed,
ani = df_analysis_cat_di$agent_ipinivo,
alpha_mu = qlogis(0.40)
)Show R Code
tibble(
variable = names(dat_cat_di_uncapped),
class = purrr::map_chr(dat_cat_di_uncapped, ~ paste(class(.x), collapse = ",")),
n = purrr::map_int(dat_cat_di_uncapped, length),
n_na = purrr::map_int(dat_cat_di_uncapped, ~ sum(is.na(.))),
n_inf = purrr::map_int(dat_cat_di_uncapped, ~ sum(is.infinite(.)))
)# A tibble: 21 × 5
variable class n n_na n_inf
<chr> <chr> <int> <int> <int>
1 N integer 1 0 0
2 y numeric 189 0 0
3 d2 integer 189 0 0
4 d3 integer 189 0 0
5 d4 integer 189 0 0
6 d5 integer 189 0 0
7 di numeric 189 0 0
8 age numeric 189 0 0
9 imm numeric 189 0 0
10 ecog integer 189 0 0
# ℹ 11 more rows3.2 Fit the weak-prior categorized ORR model with uncapped pre-response dose intensity
m_cat_real_weak_di_uncapped <- ulam(
alist(
y ~ bernoulli(p),
logit(p) <- alpha +
beta_d2 * d2 +
beta_d3 * d3 +
beta_d4 * d4 +
beta_d5 * d5 +
beta_di * di +
beta_age * age +
beta_imm * imm +
beta_ecog * ecog +
beta_lc * lc +
beta_lo * lo +
beta_lu * lu +
beta_s3 * s3 +
beta_s4 * s4 +
beta_su * su +
beta_aq3 * aq3 +
beta_aq6 * aq6 +
beta_amx * amx +
beta_ani * ani,
alpha ~ normal(alpha_mu, 1.5),
beta_d2 ~ normal(0.25, 0.25),
beta_d3 ~ normal(0.15, 0.25),
beta_d4 ~ normal(0.10, 0.25),
beta_d5 ~ normal(0.05, 0.25),
beta_di ~ normal(0.10, 0.12),
beta_age ~ normal(0, 0.02),
beta_imm ~ normal(-0.8, 0.2),
beta_ecog ~ normal(-0.3, 0.3),
beta_lc ~ normal(-2.3, 0.5),
beta_lo ~ normal(0, 0.3),
beta_lu ~ normal(0, 0.3),
beta_s3 ~ normal(0, 0.3),
beta_s4 ~ normal(0, 0.3),
beta_su ~ normal(0, 0.3),
beta_aq3 ~ normal(0, 0.3),
beta_aq6 ~ normal(0, 0.3),
beta_amx ~ normal(0, 0.3),
beta_ani ~ normal(0.4, 0.3)
),
data = dat_cat_di_uncapped,
chains = 4,
cores = 4,
iter = 2000
)3.3 Diagnostics
Show R Code
precis(m_cat_real_weak_di_uncapped, depth = 2) mean sd 5.5% 94.5% rhat ess_bulk
alpha 0.644222293 0.33514124 0.10702123 1.18910292 1.0004800 3408.720
beta_d2 0.350652913 0.23109996 -0.01614075 0.72176095 1.0004992 4660.824
beta_d3 -0.068545069 0.20697914 -0.39610787 0.25842264 0.9998906 7039.331
beta_d4 0.138391629 0.21608310 -0.21271147 0.48759719 1.0014870 7396.213
beta_d5 0.192718833 0.21598461 -0.15038819 0.53938171 0.9994993 6822.119
beta_di 0.104702802 0.11692087 -0.08064457 0.29013206 1.0005062 7436.582
beta_age 0.001370511 0.01190506 -0.01781793 0.02028470 1.0025485 7357.254
beta_imm -0.907738829 0.17712012 -1.20108110 -0.62707837 1.0005153 7101.471
beta_ecog -0.308090361 0.25071616 -0.70477556 0.09047891 1.0015742 7310.262
beta_lc -2.312182211 0.42327369 -2.98292518 -1.64873189 1.0003462 8580.399
beta_lo -0.139449634 0.26747533 -0.56289565 0.30001823 0.9998579 7467.994
beta_lu 0.055973731 0.28936293 -0.40747342 0.51135862 1.0018404 7189.880
beta_s3 -0.010185488 0.23266907 -0.38248924 0.35197766 1.0011849 5371.190
beta_s4 -0.316896805 0.23447813 -0.68129440 0.05473524 1.0009808 5907.081
beta_su 0.253377327 0.26713188 -0.16919944 0.68271941 1.0011745 8693.717
beta_aq3 -0.073253705 0.23244125 -0.44577364 0.28800665 1.0019000 7016.036
beta_aq6 -0.052635363 0.29056824 -0.51784350 0.41205756 1.0008659 8676.296
beta_amx 0.240208777 0.28150827 -0.20704940 0.69579290 0.9998606 7966.446
beta_ani 0.424686321 0.30676839 -0.05258281 0.90750624 1.0003964 7418.775Show R Code
trankplot(m_cat_real_weak_di_uncapped)
Show R Code
plot(precis(m_cat_real_weak_di_uncapped))
3.4 Posterior predictive check
Show R Code
post_cat_weak_di_uncapped <- extract.samples(m_cat_real_weak_di_uncapped)Show R Code
p_cat_weak_di_uncapped <- link(m_cat_real_weak_di_uncapped)Show R Code
ppc_cat_di_uncapped <- manual_ppc_orr(
model = m_cat_real_weak_di_uncapped,
analysis_df = df_analysis_cat_di,
response_col = "response",
n_ppc_draws = 1000,
seed = 123
)
ppc_cat_di_uncapped$summary# A tibble: 1 × 7
observed_orr mean_simulated_orr median_simulated_orr lower_89 upper_89
<dbl> <dbl> <dbl> <dbl> <dbl>
1 0.646 0.642 0.646 0.566 0.714
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>Show R Code
plot_ppc_orr(
ppc_cat_di_uncapped,
title = "Posterior Predictive Check: Categorical Dose Model + Dose Intensity",
subtitle = "Weak prior with uncapped pre-response dose intensity; red dashed line = observed ORR"
)
3.5 Extract posterior threshold and dose-intensity effects
Show R Code
post_cat_weak_di_uncapped <- extract.samples(m_cat_real_weak_di_uncapped)
posterior_dose2_cat_weak_di <- post_cat_weak_di_uncapped$beta_d2
posterior_dose3_cat_weak_di <- post_cat_weak_di_uncapped$beta_d3
posterior_dose4_cat_weak_di <- post_cat_weak_di_uncapped$beta_d4
posterior_dose5_cat_weak_di <- post_cat_weak_di_uncapped$beta_d5
posterior_di_cat_weak <- post_cat_weak_di_uncapped$beta_di
posterior_cat_weak_di_uncapped <- tibble(
beta_dose2 = posterior_dose2_cat_weak_di,
beta_dose3 = posterior_dose3_cat_weak_di,
beta_dose4 = posterior_dose4_cat_weak_di,
beta_dose5plus = posterior_dose5_cat_weak_di,
beta_di = posterior_di_cat_weak
)Show R Code
p_dose2_gt_0_cat_weak_di <- mean(posterior_dose2_cat_weak_di > 0)
p_dose3_gt_0_cat_weak_di <- mean(posterior_dose3_cat_weak_di > 0)
p_dose4_gt_0_cat_weak_di <- mean(posterior_dose4_cat_weak_di > 0)
p_dose5_gt_0_cat_weak_di <- mean(posterior_dose5_cat_weak_di > 0)
p_di_gt_0_cat_weak <- mean(posterior_di_cat_weak > 0)
p_di_gt_005_cat_weak <- mean(posterior_di_cat_weak > 0.05)
p_di_gt_010_cat_weak <- mean(posterior_di_cat_weak > 0.10)Threshold terms after adjustment for uncapped pre-response dose intensity P(β_dose2 > 0) = 93% P(β_dose3 > 0) = 38% P(β_dose4 > 0) = 74% P(β_dose5plus > 0) = 82%
Uncapped pre-response dose-intensity term P(β_DI > 0) = 81% P(β_DI > 0.05) = 68% P(β_DI > 0.10) = 51%
3.6 Compare threshold effects with and without dose intensity
Show R Code
delta_beta_dose2_di_adjustment <- posterior_dose2_cat_weak_di - post_cat_weak$beta_dose2
delta_beta_dose3_di_adjustment <- posterior_dose3_cat_weak_di - post_cat_weak$beta_dose3
delta_beta_dose4_di_adjustment <- posterior_dose4_cat_weak_di - post_cat_weak$beta_dose4
delta_beta_dose5_di_adjustment <- posterior_dose5_cat_weak_di - post_cat_weak$beta_dose5plus
delta_cat_di_comparison <- tibble(
threshold = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4"),
mean_delta = c(
mean(delta_beta_dose2_di_adjustment),
mean(delta_beta_dose3_di_adjustment),
mean(delta_beta_dose4_di_adjustment),
mean(delta_beta_dose5_di_adjustment)
),
median_delta = c(
median(delta_beta_dose2_di_adjustment),
median(delta_beta_dose3_di_adjustment),
median(delta_beta_dose4_di_adjustment),
median(delta_beta_dose5_di_adjustment)
),
ci_lower_89 = c(
quantile(delta_beta_dose2_di_adjustment, 0.055),
quantile(delta_beta_dose3_di_adjustment, 0.055),
quantile(delta_beta_dose4_di_adjustment, 0.055),
quantile(delta_beta_dose5_di_adjustment, 0.055)
),
ci_upper_89 = c(
quantile(delta_beta_dose2_di_adjustment, 0.945),
quantile(delta_beta_dose3_di_adjustment, 0.945),
quantile(delta_beta_dose4_di_adjustment, 0.945),
quantile(delta_beta_dose5_di_adjustment, 0.945)
),
p_delta_gt_0 = c(
mean(delta_beta_dose2_di_adjustment > 0),
mean(delta_beta_dose3_di_adjustment > 0),
mean(delta_beta_dose4_di_adjustment > 0),
mean(delta_beta_dose5_di_adjustment > 0)
)
)
delta_cat_di_comparison# A tibble: 4 × 6
threshold mean_delta median_delta ci_lower_89 ci_upper_89 p_delta_gt_0
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Dose 2 vs 1 -0.00423 -0.00434 -0.523 0.508 0.496
2 Dose 3 vs 2 -0.000404 -0.000907 -0.471 0.472 0.499
3 Dose 4 vs 3 0.00207 -0.00190 -0.480 0.491 0.498
4 Dose 5+ vs 4 -0.00785 -0.00442 -0.481 0.479 0.494Show R Code
cat_di_comparison <- tibble(
threshold = rep(c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4"), 2),
model = rep(
c(
"Weak prior categorical model",
"Weak prior categorical model + uncapped pre-response dose intensity"
),
each = 4
),
mean_beta = c(
mean(post_cat_weak$beta_dose2),
mean(post_cat_weak$beta_dose3),
mean(post_cat_weak$beta_dose4),
mean(post_cat_weak$beta_dose5plus),
mean(posterior_dose2_cat_weak_di),
mean(posterior_dose3_cat_weak_di),
mean(posterior_dose4_cat_weak_di),
mean(posterior_dose5_cat_weak_di)
),
median_beta = c(
median(post_cat_weak$beta_dose2),
median(post_cat_weak$beta_dose3),
median(post_cat_weak$beta_dose4),
median(post_cat_weak$beta_dose5plus),
median(posterior_dose2_cat_weak_di),
median(posterior_dose3_cat_weak_di),
median(posterior_dose4_cat_weak_di),
median(posterior_dose5_cat_weak_di)
),
ci_lower_89 = c(
quantile(post_cat_weak$beta_dose2, 0.055),
quantile(post_cat_weak$beta_dose3, 0.055),
quantile(post_cat_weak$beta_dose4, 0.055),
quantile(post_cat_weak$beta_dose5plus, 0.055),
quantile(posterior_dose2_cat_weak_di, 0.055),
quantile(posterior_dose3_cat_weak_di, 0.055),
quantile(posterior_dose4_cat_weak_di, 0.055),
quantile(posterior_dose5_cat_weak_di, 0.055)
),
ci_upper_89 = c(
quantile(post_cat_weak$beta_dose2, 0.945),
quantile(post_cat_weak$beta_dose3, 0.945),
quantile(post_cat_weak$beta_dose4, 0.945),
quantile(post_cat_weak$beta_dose5plus, 0.945),
quantile(posterior_dose2_cat_weak_di, 0.945),
quantile(posterior_dose3_cat_weak_di, 0.945),
quantile(posterior_dose4_cat_weak_di, 0.945),
quantile(posterior_dose5_cat_weak_di, 0.945)
),
p_gt_0 = c(
mean(post_cat_weak$beta_dose2 > 0),
mean(post_cat_weak$beta_dose3 > 0),
mean(post_cat_weak$beta_dose4 > 0),
mean(post_cat_weak$beta_dose5plus > 0),
mean(posterior_dose2_cat_weak_di > 0),
mean(posterior_dose3_cat_weak_di > 0),
mean(posterior_dose4_cat_weak_di > 0),
mean(posterior_dose5_cat_weak_di > 0)
)
)
cat_di_comparison# A tibble: 8 × 7
threshold model mean_beta median_beta ci_lower_89 ci_upper_89 p_gt_0
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Dose 2 vs 1 Weak prior … 0.355 0.353 -0.0232 0.733 0.936
2 Dose 3 vs 2 Weak prior … -0.0681 -0.0654 -0.406 0.276 0.368
3 Dose 4 vs 3 Weak prior … 0.136 0.137 -0.200 0.472 0.742
4 Dose 5+ vs 4 Weak prior … 0.201 0.199 -0.138 0.552 0.819
5 Dose 2 vs 1 Weak prior … 0.351 0.351 -0.0161 0.722 0.934
6 Dose 3 vs 2 Weak prior … -0.0685 -0.0712 -0.396 0.258 0.377
7 Dose 4 vs 3 Weak prior … 0.138 0.138 -0.213 0.488 0.744
8 Dose 5+ vs 4 Weak prior … 0.193 0.191 -0.150 0.539 0.8183.7 Posterior predicted response curve from the model including dose intensity
For the dose-response plot, hold di = 0, which corresponds to the mean centered uncapped pre-response dose intensity in the cohort.
Show R Code
dose_number <- 1:15
new_data_cat_di <- tibble(
d2 = as.integer(dose_number >= 2),
d3 = as.integer(dose_number >= 3),
d4 = as.integer(dose_number >= 4),
d5 = as.integer(dose_number >= 5),
di = 0,
age = 0,
imm = 0,
ecog = 0,
lc = 0,
lo = 0,
lu = 0,
s3 = 0,
s4 = 0,
su = 0,
aq3 = 0,
aq6 = 0,
amx = 0,
ani = 0
)
p_cat_weak_di_new <- link(
m_cat_real_weak_di_uncapped,
data = new_data_cat_di
)
mean_p_cat_weak_di <- apply(p_cat_weak_di_new, 2, mean)
ci_cat_weak_di <- apply(p_cat_weak_di_new, 2, PI, 0.89)Show R Code
plot_df_cat_weak_di <- tibble(
dose_number = dose_number,
mean_p = mean_p_cat_weak_di,
lower = ci_cat_weak_di[1, ],
upper = ci_cat_weak_di[2, ]
)Show R Code
plot_predicted_response_cat_weak_di_uncapped <-
ggplot(plot_df_cat_weak_di, aes(x = dose_number, y = mean_p)) +
geom_ribbon(
aes(ymin = lower, ymax = upper),
fill = "steelblue",
alpha = 0.22
) +
geom_step(
linewidth = 1.2,
color = "steelblue4"
) +
geom_point(
size = 2.8,
shape = 21,
stroke = 0.9,
fill = "white",
color = "steelblue4"
) +
scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, by = 0.1),
labels = scales::percent_format(accuracy = 1),
expand = expansion(mult = c(0.01, 0.02))
) +
scale_x_continuous(
breaks = 1:15,
expand = expansion(mult = c(0.01, 0.02))
) +
labs(
title = "Predicted Response Probability by Number of Doses",
subtitle = "Weak prior categorical model adjusted for uncapped pre-response dose intensity\nDose intensity held at cohort mean; shaded region = 89% CrI",
x = "Number of doses",
y = "Predicted response probability"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
plot_predicted_response_cat_weak_di_uncapped
Show R Code
posterior_cat_weak_di_long <- bind_rows(
tibble(beta = posterior_dose2_cat_weak_di, threshold = "Dose 2 vs 1"),
tibble(beta = posterior_dose3_cat_weak_di, threshold = "Dose 3 vs 2"),
tibble(beta = posterior_dose4_cat_weak_di, threshold = "Dose 4 vs 3"),
tibble(beta = posterior_dose5_cat_weak_di, threshold = "Dose 5+ vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4")
)
)
cat_weak_di_summary <- posterior_cat_weak_di_long |>
group_by(threshold) |>
summarise(
mean_beta = mean(beta),
lower = quantile(beta, 0.055),
upper = quantile(beta, 0.945),
p_gt_0 = mean(beta > 0),
.groups = "drop"
)
plot_cat_density_weak_di <-
ggplot(posterior_cat_weak_di_long, aes(x = beta)) +
geom_density(fill = "steelblue", alpha = 0.35, color = "steelblue4", linewidth = 0.9) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 0.9) +
geom_vline(
data = cat_weak_di_summary,
aes(xintercept = lower),
linetype = "dotted",
color = "grey40",
linewidth = 0.7
) +
geom_vline(
data = cat_weak_di_summary,
aes(xintercept = upper),
linetype = "dotted",
color = "grey40",
linewidth = 0.7
) +
geom_text(
data = cat_weak_di_summary,
aes(
x = Inf,
y = Inf,
label = paste0(
"89% CrI: [", round(lower, 2), ", ", round(upper, 2), "]\n",
"P(β > 0) = ", round(p_gt_0 * 100, 0), "%"
)
),
hjust = 1.05,
vjust = 1.5,
size = 3.3,
inherit.aes = FALSE
) +
facet_wrap(~ threshold, scales = "free_y") +
labs(
title = "Posterior Distributions: Categorical Threshold Effects",
subtitle = "Weak prior model adjusted for uncapped pre-response dose intensity",
x = "Threshold-specific log-odds effect",
y = "Density"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
strip.text = element_text(face = "bold")
)
plot_cat_density_weak_di
Show R Code
posterior_cat_combined_weak_di <-
bind_rows(
tibble(beta = posterior_dose2_cat_weak_di, threshold = "Dose 2 vs 1"),
tibble(beta = posterior_dose3_cat_weak_di, threshold = "Dose 3 vs 2"),
tibble(beta = posterior_dose4_cat_weak_di, threshold = "Dose 4 vs 3"),
tibble(beta = posterior_dose5_cat_weak_di, threshold = "Dose ≥5 vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose ≥5 vs 4")
)
)
summ_cat_weak_di <-
posterior_cat_combined_weak_di |>
group_by(threshold) |>
summarise(
p_gt0 = mean(beta > 0),
OR_md = exp(median(beta)),
.groups = "drop"
)
ann_text_weak_di <- paste0(
"Posterior Pr(β > 0):\n",
"Dose 2 vs 1: ", fmt_pct(summ_cat_weak_di$p_gt0[summ_cat_weak_di$threshold == "Dose 2 vs 1"]), "\n",
"Dose 3 vs 2: ", fmt_pct(summ_cat_weak_di$p_gt0[summ_cat_weak_di$threshold == "Dose 3 vs 2"]), "\n",
"Dose 4 vs 3: ", fmt_pct(summ_cat_weak_di$p_gt0[summ_cat_weak_di$threshold == "Dose 4 vs 3"]), "\n",
"Dose ≥5 vs 4: ", fmt_pct(summ_cat_weak_di$p_gt0[summ_cat_weak_di$threshold == "Dose ≥5 vs 4"]), "\n\n",
"Median OR (incremental):\n",
"Dose 2 vs 1: ×", sprintf("%.2f", summ_cat_weak_di$OR_md[summ_cat_weak_di$threshold == "Dose 2 vs 1"]), "\n",
"Dose 3 vs 2: ×", sprintf("%.2f", summ_cat_weak_di$OR_md[summ_cat_weak_di$threshold == "Dose 3 vs 2"]), "\n",
"Dose 4 vs 3: ×", sprintf("%.2f", summ_cat_weak_di$OR_md[summ_cat_weak_di$threshold == "Dose 4 vs 3"]), "\n",
"Dose ≥5 vs 4: ×", sprintf("%.2f", summ_cat_weak_di$OR_md[summ_cat_weak_di$threshold == "Dose ≥5 vs 4"])
)
dens_list_weak_di <- posterior_cat_combined_weak_di |>
group_by(threshold) |>
summarise(
dens = list(stats::density(beta, adjust = adjust_k)),
.groups = "drop"
)
ymax_weak_di <- max(vapply(dens_list_weak_di$dens, function(d) max(d$y), numeric(1)))
plot_categorical_weak_di <-
ggplot(posterior_cat_combined_weak_di, aes(x = beta)) +
geom_density(aes(fill = threshold), alpha = 0.18, color = NA, adjust = 1.1) +
geom_density(aes(color = threshold), fill = NA, linewidth = 1.15, adjust = 1.1) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
scale_fill_manual(values = cols_cat) +
scale_color_manual(values = cols_cat) +
scale_x_continuous(expand = expansion(mult = c(0.02, 0.02))) +
scale_y_continuous(expand = expansion(mult = c(0.02, 0.02))) +
annotate(
"label",
x = 1.25,
y = ymax_weak_di,
label = ann_text_weak_di,
hjust = 0, vjust = 1,
size = 3.1,
label.size = 0.25,
fill = "white", alpha = 0.92
) +
labs(
title = "Posterior Distributions: Categorical Dose Threshold Effects",
subtitle = "Weak Prior Model Adjusted for Uncapped Pre-Response Dose Intensity",
x = "β (log-odds; incremental effect of crossing each threshold)",
y = "Density",
color = "Increment",
fill = "Increment"
) +
coord_cartesian(xlim = c(NA, 2.1), clip = "off") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, face = "bold"),
legend.position = "bottom",
plot.margin = margin(t = 10, r = 40, b = 10, l = 10)
)
plot_categorical_weak_di
3.8 Posterior threshold-effect plot after dose-intensity adjustment
Show R Code
threshold_summary_plot_cat_weak_di <- tibble(
threshold = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4"),
mean_beta = c(
mean(posterior_dose2_cat_weak_di),
mean(posterior_dose3_cat_weak_di),
mean(posterior_dose4_cat_weak_di),
mean(posterior_dose5_cat_weak_di)
),
lower = c(
quantile(posterior_dose2_cat_weak_di, 0.055),
quantile(posterior_dose3_cat_weak_di, 0.055),
quantile(posterior_dose4_cat_weak_di, 0.055),
quantile(posterior_dose5_cat_weak_di, 0.055)
),
upper = c(
quantile(posterior_dose2_cat_weak_di, 0.945),
quantile(posterior_dose3_cat_weak_di, 0.945),
quantile(posterior_dose4_cat_weak_di, 0.945),
quantile(posterior_dose5_cat_weak_di, 0.945)
),
p_gt_0 = c(
mean(posterior_dose2_cat_weak_di > 0),
mean(posterior_dose3_cat_weak_di > 0),
mean(posterior_dose4_cat_weak_di > 0),
mean(posterior_dose5_cat_weak_di > 0)
)
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4")
)
)
plot_cat_weak_di <-
ggplot(threshold_summary_plot_cat_weak_di,
aes(x = threshold, y = mean_beta)) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_pointrange(
aes(ymin = lower, ymax = upper),
linewidth = 0.9
) +
geom_text(
aes(
y = upper + 0.05,
label = paste0("P(>0) = ", round(p_gt_0 * 100, 0), "%")
),
size = 3.5
) +
labs(
title = "Posterior Threshold Effects: Categorical Dose Model",
subtitle = "Weak prior model adjusted for uncapped pre-response dose intensity; points = posterior means, bars = 89% CrI",
x = NULL,
y = "Log-odds increment at threshold"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.text.x = element_text(face = "bold")
)
plot_cat_weak_di
3.9 Companion plot for the dose-intensity coefficient itself
Show R Code
dens_di <- density(posterior_di_cat_weak)
dens_df_di <- data.frame(x = dens_di$x, y = dens_di$y)
ymax_di <- max(dens_df_di$y)
ci_lower_di <- quantile(posterior_di_cat_weak, 0.055)
ci_upper_di <- quantile(posterior_di_cat_weak, 0.945)
plot_di_cat_weak <-
ggplot() +
geom_line(
data = dens_df_di,
aes(x = x, y = y),
linewidth = 0.9,
color = "#08306B"
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
geom_vline(
xintercept = c(ci_lower_di, ci_upper_di),
linetype = "dotted",
color = "grey40",
linewidth = 0.8
) +
annotate(
"text",
x = mean(posterior_di_cat_weak),
y = ymax_di * 0.80,
label = paste0(
"89% CrI: [",
round(ci_lower_di, 2),
", ",
round(ci_upper_di, 2),
"]"
),
size = 3.4
) +
labs(
title = "Posterior Distribution: Uncapped Pre-Response Dose Intensity Effect",
subtitle = "Weak prior categorical model",
x = expression(beta[DI]),
y = "Density"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, face = "bold")
)
plot_di_cat_weak
3.10 Diagnostics
Show R Code
df_analysis_cat_di |>
summarise(
corr_numdoses_di = cor(
num_doses,
dose_intensity_pre_response_uncapped_centered,
use = "complete.obs"
),
corr_dose2_di = cor(
dose_2,
dose_intensity_pre_response_uncapped_centered,
use = "complete.obs"
),
corr_dose3_di = cor(
dose_3,
dose_intensity_pre_response_uncapped_centered,
use = "complete.obs"
),
corr_dose4_di = cor(
dose_4,
dose_intensity_pre_response_uncapped_centered,
use = "complete.obs"
),
corr_dose5_di = cor(
dose_5plus,
dose_intensity_pre_response_uncapped_centered,
use = "complete.obs"
)
)# A tibble: 1 × 5
corr_numdoses_di corr_dose2_di corr_dose3_di corr_dose4_di corr_dose5_di
<dbl> <dbl> <dbl> <dbl> <dbl>
1 -0.0704 -0.0506 -0.156 -0.0780 -0.0630Show R Code
ggplot(
df_analysis_cat_di,
aes(x = num_doses, y = dose_intensity_pre_response_uncapped_centered)
) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = FALSE, color = "red") +
labs(
title = "Relationship Between Number of Doses and Pre-Response Uncapped Dose Intensity",
x = "Number of doses",
y = "Centered uncapped pre-response dose intensity"
) +
theme_minimal()
Show R Code
df_joint_cat <- tibble(
beta_dose2 = posterior_dose2_cat_weak_di,
beta_di = posterior_di_cat_weak
)
joint_posterior_plot_cat <-
ggplot(df_joint_cat, aes(x = beta_dose2, y = beta_di)) +
geom_point(alpha = 0.15, size = 1) +
stat_density_2d(
aes(fill = after_stat(level)),
geom = "polygon",
alpha = 0.3
) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
labs(
title = "Joint Posterior of Early Dose Threshold Effect and Dose Intensity",
subtitle = "Weak prior categorical model with uncapped pre-response dose intensity",
x = expression(beta[dose2]),
y = expression(beta[DI])
) +
theme_classic() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(face = "bold", hjust = 0.5)
)
joint_posterior_plot_cat
Show R Code
delta_cat_df <- bind_rows(
tibble(delta_beta = delta_beta_dose2_di_adjustment, threshold = "Dose 2 vs 1"),
tibble(delta_beta = delta_beta_dose3_di_adjustment, threshold = "Dose 3 vs 2"),
tibble(delta_beta = delta_beta_dose4_di_adjustment, threshold = "Dose 4 vs 3"),
tibble(delta_beta = delta_beta_dose5_di_adjustment, threshold = "Dose 5+ vs 4")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4")
)
)
ggplot(delta_cat_df, aes(x = delta_beta)) +
geom_density(fill = "steelblue", alpha = 0.4) +
geom_vline(xintercept = 0, linetype = "dashed") +
facet_wrap(~ threshold, scales = "free_y") +
labs(
title = "Difference in Posterior Threshold Effects After Adding Dose Intensity",
subtitle = "Adjusted model minus unadjusted model",
x = expression(Delta * beta),
y = "Density"
) +
theme_minimal() +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(face = "bold", hjust = 0.5)
)
Dose-intensity adjustment was evaluated by adding the uncapped pre-response dose-intensity variable to the weak-prior categorical ORR model and comparing the posterior threshold-effect estimates before and after adjustment. This comparison shows how much inference on the threshold-specific dose effects changes after accounting for treatment timing intensity prior to response.
Show R Code
delta_cat_di_comparison# A tibble: 4 × 6
threshold mean_delta median_delta ci_lower_89 ci_upper_89 p_delta_gt_0
<chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Dose 2 vs 1 -0.00423 -0.00434 -0.523 0.508 0.496
2 Dose 3 vs 2 -0.000404 -0.000907 -0.471 0.472 0.499
3 Dose 4 vs 3 0.00207 -0.00190 -0.480 0.491 0.498
4 Dose 5+ vs 4 -0.00785 -0.00442 -0.481 0.479 0.494Show R Code
posterior_cat_compare_long <- bind_rows(
tibble(beta = post_cat_skeptical$beta_dose2, threshold = "Dose 2 vs 1", model = "Skeptical"),
tibble(beta = post_cat_skeptical$beta_dose3, threshold = "Dose 3 vs 2", model = "Skeptical"),
tibble(beta = post_cat_skeptical$beta_dose4, threshold = "Dose 4 vs 3", model = "Skeptical"),
tibble(beta = post_cat_skeptical$beta_dose5plus, threshold = "Dose 5+ vs 4", model = "Skeptical"),
tibble(beta = post_cat_weak$beta_dose2, threshold = "Dose 2 vs 1", model = "Weak"),
tibble(beta = post_cat_weak$beta_dose3, threshold = "Dose 3 vs 2", model = "Weak"),
tibble(beta = post_cat_weak$beta_dose4, threshold = "Dose 4 vs 3", model = "Weak"),
tibble(beta = post_cat_weak$beta_dose5plus, threshold = "Dose 5+ vs 4", model = "Weak"),
tibble(beta = posterior_dose2_cat_weak_di, threshold = "Dose 2 vs 1", model = "Weak + DI"),
tibble(beta = posterior_dose3_cat_weak_di, threshold = "Dose 3 vs 2", model = "Weak + DI"),
tibble(beta = posterior_dose4_cat_weak_di, threshold = "Dose 4 vs 3", model = "Weak + DI"),
tibble(beta = posterior_dose5_cat_weak_di, threshold = "Dose 5+ vs 4", model = "Weak + DI")
) |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose 5+ vs 4")
),
model = factor(model, levels = c("Skeptical", "Weak", "Weak + DI"))
)
ggplot(posterior_cat_compare_long, aes(x = beta, color = model, fill = model)) +
geom_density(alpha = 0.18, linewidth = 0.9) +
geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
facet_wrap(~ threshold, scales = "free_y") +
labs(
title = "Posterior Distributions Across Categorical Model Specifications",
subtitle = "Comparison of skeptical, weak, and weak + dose-intensity-adjusted models",
x = "Threshold-specific log-odds effect",
y = "Density"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
strip.text = element_text(face = "bold")
)
Show R Code
# ------------------------------------------------------------------------------
# Sensitivity analysis: categorical ORR model using pre-response dose count
# ------------------------------------------------------------------------------
df_analysis_cat_pre_response <- df_analysis |>
mutate(
# Number of doses administered before response assessment / response documentation
dose_pre_response = num_doses_pre_response,
dose_pre_response_minus_1 = dose_pre_response - 1,
log_dose_pre_response = log(dose_pre_response),
dose_pre_response_group = if_else(dose_pre_response <= 2, "2 doses", "3+ doses"),
# Cumulative categorical thresholds based on pre-response dose count
dose_pre_response_2 = as.integer(dose_pre_response >= 2),
dose_pre_response_3 = as.integer(dose_pre_response >= 3),
dose_pre_response_4 = as.integer(dose_pre_response >= 4),
dose_pre_response_5plus = as.integer(dose_pre_response >= 5)
) |>
filter(
!is.na(response),
!is.na(dose_pre_response),
dose_pre_response >= 1,
!is.na(dose_pre_response_2),
!is.na(dose_pre_response_3),
!is.na(dose_pre_response_4),
!is.na(dose_pre_response_5plus),
!is.na(age_centered),
!is.na(immunosuppressed),
!is.na(ecog23),
!is.na(loc_conj),
!is.na(loc_other),
!is.na(loc_unknown),
!is.na(stage_III),
!is.na(stage_IV),
!is.na(stage_unstaged),
!is.na(agent_pembro_q3),
!is.na(agent_pembro_q6),
!is.na(agent_pembro_mixed),
!is.na(agent_ipinivo)
)Show R Code
df_analysis_cat_pre_response |>
summarise(
n = n(),
observed_orr = mean(response),
min_pre_response_doses = min(dose_pre_response, na.rm = TRUE),
max_pre_response_doses = max(dose_pre_response, na.rm = TRUE),
mean_total_doses = mean(num_doses, na.rm = TRUE),
mean_pre_response_doses = mean(dose_pre_response, na.rm = TRUE),
n_total_differs_from_pre_response = sum(num_doses != dose_pre_response, na.rm = TRUE)
)# A tibble: 1 × 7
n observed_orr min_pre_response_doses max_pre_response_doses
<int> <dbl> <int> <int>
1 189 0.646 1 20
# ℹ 3 more variables: mean_total_doses <dbl>, mean_pre_response_doses <dbl>,
# n_total_differs_from_pre_response <int>Show R Code
dat_cat_pre_response <- list(
N = nrow(df_analysis_cat_pre_response),
y = df_analysis_cat_pre_response$response,
# Pre-response categorical thresholds
dose_2 = df_analysis_cat_pre_response$dose_pre_response_2,
dose_3 = df_analysis_cat_pre_response$dose_pre_response_3,
dose_4 = df_analysis_cat_pre_response$dose_pre_response_4,
dose_5plus = df_analysis_cat_pre_response$dose_pre_response_5plus,
age = df_analysis_cat_pre_response$age_centered,
imm = as.integer(df_analysis_cat_pre_response$immunosuppressed) - 1,
ecog23 = df_analysis_cat_pre_response$ecog23,
loc_conj = df_analysis_cat_pre_response$loc_conj,
loc_other = df_analysis_cat_pre_response$loc_other,
loc_unknown = df_analysis_cat_pre_response$loc_unknown,
stage_III = df_analysis_cat_pre_response$stage_III,
stage_IV = df_analysis_cat_pre_response$stage_IV,
stage_unstaged = df_analysis_cat_pre_response$stage_unstaged,
agent_pembro_q3 = df_analysis_cat_pre_response$agent_pembro_q3,
agent_pembro_q6 = df_analysis_cat_pre_response$agent_pembro_q6,
agent_pembro_mixed = df_analysis_cat_pre_response$agent_pembro_mixed,
agent_ipinivo = df_analysis_cat_pre_response$agent_ipinivo,
alpha_mu = qlogis(0.40)
)m_cat_real_weak_pre_response <- ulam(
alist(
y ~ bernoulli(p),
logit(p) <- alpha +
beta_dose2 * dose_2 +
beta_dose3 * dose_3 +
beta_dose4 * dose_4 +
beta_dose5plus * dose_5plus +
beta_age * age +
beta_imm * imm +
beta_ecog * ecog23 +
beta_conj * loc_conj +
beta_other * loc_other +
beta_unknown * loc_unknown +
beta_stageIII * stage_III +
beta_stageIV * stage_IV +
beta_stageUn * stage_unstaged +
beta_pembro_q3 * agent_pembro_q3 +
beta_pembro_q6 * agent_pembro_q6 +
beta_pembro_mixed * agent_pembro_mixed +
beta_ipinivo * agent_ipinivo,
alpha ~ normal(alpha_mu, 1.5),
# Same weakly informative categorical priors as total-dose model
beta_dose2 ~ normal(0.25, 0.25),
beta_dose3 ~ normal(0.15, 0.25),
beta_dose4 ~ normal(0.10, 0.25),
beta_dose5plus ~ normal(0.05, 0.25),
beta_age ~ normal(0, 0.02),
beta_imm ~ normal(-0.8, 0.2),
beta_ecog ~ normal(-0.3, 0.3),
beta_conj ~ normal(-2.3, 0.5),
beta_other ~ normal(0, 0.3),
beta_unknown ~ normal(0, 0.3),
beta_stageIII ~ normal(0, 0.3),
beta_stageIV ~ normal(0, 0.3),
beta_stageUn ~ normal(0, 0.3),
beta_pembro_q3 ~ normal(0, 0.3),
beta_pembro_q6 ~ normal(0, 0.3),
beta_pembro_mixed ~ normal(0, 0.3),
beta_ipinivo ~ normal(0.4, 0.3)
),
data = dat_cat_pre_response,
chains = 4,
cores = 4,
iter = 2000
)Show R Code
precis(m_cat_real_weak_pre_response, depth = 2) mean sd 5.5% 94.5% rhat
alpha 0.5615149963 0.33231272 0.04489284 1.10696011 1.0005308
beta_dose2 0.3885779504 0.22583689 0.03642770 0.74713004 1.0001001
beta_dose3 0.0048120898 0.20887528 -0.32520766 0.33588297 0.9998988
beta_dose4 0.1983962836 0.21469035 -0.13999484 0.54028789 1.0007690
beta_dose5plus 0.2330430935 0.21657112 -0.11138654 0.58194435 1.0011561
beta_age 0.0015693797 0.01192352 -0.01764487 0.02059912 1.0010167
beta_imm -0.8996524063 0.18187148 -1.18857816 -0.60533469 1.0004533
beta_ecog -0.3092374895 0.26053882 -0.71684104 0.10839448 1.0019321
beta_conj -2.2943287518 0.42437909 -2.99530125 -1.63051428 0.9999577
beta_other -0.1455416818 0.26778788 -0.58061376 0.28032629 1.0020528
beta_unknown 0.0541900769 0.29341716 -0.42070476 0.52255488 1.0007876
beta_stageIII -0.0008699984 0.24309116 -0.38639217 0.39146155 0.9993231
beta_stageIV -0.3114039816 0.23957836 -0.69942099 0.06464211 1.0002509
beta_stageUn 0.2452288892 0.27662837 -0.20239055 0.69544968 1.0023207
beta_pembro_q3 -0.0574590745 0.23443632 -0.42997859 0.32114561 1.0013645
beta_pembro_q6 -0.0560191788 0.29212711 -0.52301830 0.41807913 0.9998776
beta_pembro_mixed 0.2239854781 0.28163072 -0.21704470 0.67009799 1.0009072
beta_ipinivo 0.4228140394 0.29744511 -0.04520905 0.89182765 1.0015257
ess_bulk
alpha 3682.649
beta_dose2 6176.627
beta_dose3 6997.907
beta_dose4 6262.687
beta_dose5plus 7586.298
beta_age 9323.948
beta_imm 8730.188
beta_ecog 7510.286
beta_conj 9841.129
beta_other 8377.742
beta_unknown 7762.264
beta_stageIII 5217.635
beta_stageIV 5477.764
beta_stageUn 8265.998
beta_pembro_q3 7378.094
beta_pembro_q6 7457.922
beta_pembro_mixed 9702.041
beta_ipinivo 9017.915Show R Code
plot(precis(m_cat_real_weak_pre_response))
Show R Code
post_cat_weak_pre_response <- extract.samples(m_cat_real_weak_pre_response)
posterior_cat_weak_pre_response <- tibble(
beta_dose2 = post_cat_weak_pre_response$beta_dose2,
beta_dose3 = post_cat_weak_pre_response$beta_dose3,
beta_dose4 = post_cat_weak_pre_response$beta_dose4,
beta_dose5plus = post_cat_weak_pre_response$beta_dose5plus
)
save_files(
save_object = posterior_cat_weak_pre_response,
filename = "Posterior Categorical Threshold Effects Weak Prior PreResponse Dose",
subD = "files/Modeling/Posterior Categorical Threshold Effects Weak Prior PreResponse Dose"
)File saved to: /Users/davidmiller/Partners HealthCare Dropbox/David Miller/mLab/Projects/FIRST and Neoadjuvant Outcomes in CSCC at MGH-MEEI/files/Modeling/Posterior Categorical Threshold Effects Weak Prior PreResponse Dose/Posterior Categorical Threshold Effects Weak Prior PreResponse Dose-2026-05-08-05-13-52.rds Show R Code
threshold_summary_cat_weak_pre_response <- tibble(
threshold = c(
"Pre-response dose 2 vs 1",
"Pre-response dose 3 vs 2",
"Pre-response dose 4 vs 3",
"Pre-response dose ≥5 vs 4"
),
mean_beta = c(
mean(post_cat_weak_pre_response$beta_dose2),
mean(post_cat_weak_pre_response$beta_dose3),
mean(post_cat_weak_pre_response$beta_dose4),
mean(post_cat_weak_pre_response$beta_dose5plus)
),
median_beta = c(
median(post_cat_weak_pre_response$beta_dose2),
median(post_cat_weak_pre_response$beta_dose3),
median(post_cat_weak_pre_response$beta_dose4),
median(post_cat_weak_pre_response$beta_dose5plus)
),
ci_lower_89 = c(
quantile(post_cat_weak_pre_response$beta_dose2, 0.055),
quantile(post_cat_weak_pre_response$beta_dose3, 0.055),
quantile(post_cat_weak_pre_response$beta_dose4, 0.055),
quantile(post_cat_weak_pre_response$beta_dose5plus, 0.055)
),
ci_upper_89 = c(
quantile(post_cat_weak_pre_response$beta_dose2, 0.945),
quantile(post_cat_weak_pre_response$beta_dose3, 0.945),
quantile(post_cat_weak_pre_response$beta_dose4, 0.945),
quantile(post_cat_weak_pre_response$beta_dose5plus, 0.945)
),
p_gt_0 = c(
mean(post_cat_weak_pre_response$beta_dose2 > 0),
mean(post_cat_weak_pre_response$beta_dose3 > 0),
mean(post_cat_weak_pre_response$beta_dose4 > 0),
mean(post_cat_weak_pre_response$beta_dose5plus > 0)
),
median_or = exp(median_beta)
)
threshold_summary_cat_weak_pre_response# A tibble: 4 × 7
threshold mean_beta median_beta ci_lower_89 ci_upper_89 p_gt_0 median_or
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Pre-response d… 0.389 0.383 0.0364 0.747 0.961 1.47
2 Pre-response d… 0.00481 0.00321 -0.325 0.336 0.506 1.00
3 Pre-response d… 0.198 0.200 -0.140 0.540 0.820 1.22
4 Pre-response d… 0.233 0.234 -0.111 0.582 0.856 1.26Show R Code
# ------------------------------------------------------------------------------
# Manual posterior predictive check from fitted probabilities
# ------------------------------------------------------------------------------
set.seed(123)
p_hat_cat_pre_response <- link(m_cat_real_weak_pre_response)
stopifnot(ncol(p_hat_cat_pre_response) == nrow(df_analysis_cat_pre_response))
n_ppc_draws <- 1000
draw_ids <- sample(
seq_len(nrow(p_hat_cat_pre_response)),
size = n_ppc_draws,
replace = FALSE
)
p_hat_ppc_cat_pre_response <-
p_hat_cat_pre_response[draw_ids, , drop = FALSE]
yrep_cat_pre_response <- matrix(
rbinom(
n = length(p_hat_ppc_cat_pre_response),
size = 1,
prob = as.vector(p_hat_ppc_cat_pre_response)
),
nrow = nrow(p_hat_ppc_cat_pre_response),
ncol = ncol(p_hat_ppc_cat_pre_response)
)
observed_orr_cat_pre_response <-
mean(df_analysis_cat_pre_response$response)
ppc_cat_pre_response_orr <- tibble(
draw = seq_len(nrow(yrep_cat_pre_response)),
simulated_orr = rowMeans(yrep_cat_pre_response)
)
ppc_cat_pre_response_orr |>
summarise(
observed_orr = observed_orr_cat_pre_response,
mean_simulated_orr = mean(simulated_orr),
median_simulated_orr = median(simulated_orr),
lower_89 = quantile(simulated_orr, 0.055),
upper_89 = quantile(simulated_orr, 0.945),
p_sim_le_observed = mean(simulated_orr <= observed_orr_cat_pre_response),
mean_fitted_probability = mean(p_hat_cat_pre_response)
)# A tibble: 1 × 7
observed_orr mean_simulated_orr median_simulated_orr lower_89 upper_89
<dbl> <dbl> <dbl> <dbl> <dbl>
1 0.646 0.642 0.640 0.571 0.714
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>Show R Code
ggplot(ppc_cat_pre_response_orr, aes(x = simulated_orr)) +
geom_histogram(
bins = 50,
fill = "steelblue",
alpha = 0.7,
color = "white"
) +
geom_vline(
xintercept = observed_orr_cat_pre_response,
color = "red",
linewidth = 1.2,
linetype = "dashed"
) +
scale_x_continuous(
limits = c(0, 1),
labels = scales::percent_format(accuracy = 1)
) +
labs(
title = "Posterior Predictive Check: Pre-Response Categorical Dose Model",
subtitle = "Red dashed line = observed ORR",
x = "Simulated cohort ORR",
y = "Count"
) +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5, face = "bold")
)
Show R Code
max_pre_response_dose_for_plot <- min(
15,
max(df_analysis_cat_pre_response$dose_pre_response, na.rm = TRUE)
)
dose_number_pre_response <- 1:max_pre_response_dose_for_plot
new_data_cat_pre_response <- tibble(
dose_number = dose_number_pre_response,
dose_2 = as.integer(dose_number_pre_response >= 2),
dose_3 = as.integer(dose_number_pre_response >= 3),
dose_4 = as.integer(dose_number_pre_response >= 4),
dose_5plus = as.integer(dose_number_pre_response >= 5),
age = 0,
imm = 0,
ecog23 = 0,
loc_conj = 0,
loc_other = 0,
loc_unknown = 0,
stage_III = 0,
stage_IV = 0,
stage_unstaged = 0,
agent_pembro_q3 = 0,
agent_pembro_q6 = 0,
agent_pembro_mixed = 0,
agent_ipinivo = 0
)
p_cat_weak_pre_response_new <- link(
m_cat_real_weak_pre_response,
data = new_data_cat_pre_response
)
mean_p_cat_weak_pre_response <- apply(
p_cat_weak_pre_response_new,
2,
mean
)
ci_cat_weak_pre_response <- apply(
p_cat_weak_pre_response_new,
2,
PI,
0.89
)
plot_df_cat_weak_pre_response <- tibble(
dose_number = dose_number_pre_response,
mean_p = mean_p_cat_weak_pre_response,
lower = ci_cat_weak_pre_response[1, ],
upper = ci_cat_weak_pre_response[2, ]
)Show R Code
plot_predicted_response_cat_weak_pre_response <-
ggplot(plot_df_cat_weak_pre_response, aes(x = dose_number, y = mean_p)) +
geom_ribbon(
aes(ymin = lower, ymax = upper),
fill = "steelblue",
alpha = 0.22
) +
geom_step(linewidth = 1.2, color = "steelblue4") +
geom_point(
size = 2.8,
shape = 21,
stroke = 0.9,
fill = "white",
color = "steelblue4"
) +
scale_y_continuous(
limits = c(0, 1),
breaks = seq(0, 1, by = 0.1),
labels = scales::percent_format(accuracy = 1),
expand = expansion(mult = c(0.01, 0.02))
) +
scale_x_continuous(
breaks = seq(
min(plot_df_cat_weak_pre_response$dose_number),
max(plot_df_cat_weak_pre_response$dose_number),
by = 1
),
expand = expansion(mult = c(0.01, 0.02))
) +
labs(
title = "Predicted Response Probability by Pre-Response Dose Count",
subtitle = "Weak-prior categorical model using doses before response assessment; shaded region = 89% CrI",
x = "Number of pre-response doses",
y = "Predicted response probability"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.title = element_text(face = "bold"),
panel.grid.minor = element_blank()
)
plot_predicted_response_cat_weak_pre_response
Show R Code
cat_total_vs_pre_response_comparison <- tibble(
threshold = rep(
c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose ≥5 vs 4"),
2
),
exposure_definition = rep(
c("Total cumulative dose", "Pre-response dose"),
each = 4
),
mean_beta = c(
mean(post_cat_weak$beta_dose2),
mean(post_cat_weak$beta_dose3),
mean(post_cat_weak$beta_dose4),
mean(post_cat_weak$beta_dose5plus),
mean(post_cat_weak_pre_response$beta_dose2),
mean(post_cat_weak_pre_response$beta_dose3),
mean(post_cat_weak_pre_response$beta_dose4),
mean(post_cat_weak_pre_response$beta_dose5plus)
),
median_beta = c(
median(post_cat_weak$beta_dose2),
median(post_cat_weak$beta_dose3),
median(post_cat_weak$beta_dose4),
median(post_cat_weak$beta_dose5plus),
median(post_cat_weak_pre_response$beta_dose2),
median(post_cat_weak_pre_response$beta_dose3),
median(post_cat_weak_pre_response$beta_dose4),
median(post_cat_weak_pre_response$beta_dose5plus)
),
ci_lower_89 = c(
quantile(post_cat_weak$beta_dose2, 0.055),
quantile(post_cat_weak$beta_dose3, 0.055),
quantile(post_cat_weak$beta_dose4, 0.055),
quantile(post_cat_weak$beta_dose5plus, 0.055),
quantile(post_cat_weak_pre_response$beta_dose2, 0.055),
quantile(post_cat_weak_pre_response$beta_dose3, 0.055),
quantile(post_cat_weak_pre_response$beta_dose4, 0.055),
quantile(post_cat_weak_pre_response$beta_dose5plus, 0.055)
),
ci_upper_89 = c(
quantile(post_cat_weak$beta_dose2, 0.945),
quantile(post_cat_weak$beta_dose3, 0.945),
quantile(post_cat_weak$beta_dose4, 0.945),
quantile(post_cat_weak$beta_dose5plus, 0.945),
quantile(post_cat_weak_pre_response$beta_dose2, 0.945),
quantile(post_cat_weak_pre_response$beta_dose3, 0.945),
quantile(post_cat_weak_pre_response$beta_dose4, 0.945),
quantile(post_cat_weak_pre_response$beta_dose5plus, 0.945)
),
p_gt_0 = c(
mean(post_cat_weak$beta_dose2 > 0),
mean(post_cat_weak$beta_dose3 > 0),
mean(post_cat_weak$beta_dose4 > 0),
mean(post_cat_weak$beta_dose5plus > 0),
mean(post_cat_weak_pre_response$beta_dose2 > 0),
mean(post_cat_weak_pre_response$beta_dose3 > 0),
mean(post_cat_weak_pre_response$beta_dose4 > 0),
mean(post_cat_weak_pre_response$beta_dose5plus > 0)
),
median_or = exp(median_beta)
)
cat_total_vs_pre_response_comparison# A tibble: 8 × 8
threshold exposure_definition mean_beta median_beta ci_lower_89 ci_upper_89
<chr> <chr> <dbl> <dbl> <dbl> <dbl>
1 Dose 2 vs 1 Total cumulative d… 0.355 0.353 -0.0232 0.733
2 Dose 3 vs 2 Total cumulative d… -0.0681 -0.0654 -0.406 0.276
3 Dose 4 vs 3 Total cumulative d… 0.136 0.137 -0.200 0.472
4 Dose ≥5 vs 4 Total cumulative d… 0.201 0.199 -0.138 0.552
5 Dose 2 vs 1 Pre-response dose 0.389 0.383 0.0364 0.747
6 Dose 3 vs 2 Pre-response dose 0.00481 0.00321 -0.325 0.336
7 Dose 4 vs 3 Pre-response dose 0.198 0.200 -0.140 0.540
8 Dose ≥5 vs 4 Pre-response dose 0.233 0.234 -0.111 0.582
# ℹ 2 more variables: p_gt_0 <dbl>, median_or <dbl>Show R Code
cat_total_vs_pre_response_comparison |>
mutate(
threshold = factor(
threshold,
levels = c("Dose 2 vs 1", "Dose 3 vs 2", "Dose 4 vs 3", "Dose ≥5 vs 4")
)
) |>
ggplot(
aes(
x = threshold,
y = median_beta,
ymin = ci_lower_89,
ymax = ci_upper_89,
color = exposure_definition
)
) +
geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
geom_pointrange(
position = position_dodge(width = 0.45),
linewidth = 0.8
) +
labs(
title = "Categorical Threshold Effects by Dose Definition",
subtitle = "Total cumulative dose vs pre-response dose; points = posterior medians, bars = 89% CrI",
x = NULL,
y = "Threshold-specific log-odds effect",
color = "Exposure definition"
) +
theme_minimal(base_size = 14) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
axis.text.x = element_text(face = "bold")
)
Show R Code
# ------------------------------------------------------------------------------
# Posterior distributions: pre-response categorical threshold effects
# ------------------------------------------------------------------------------
posterior_cat_pre_response_long <-
bind_rows(
tibble(
beta = post_cat_weak_pre_response$beta_dose2,
threshold = "Dose 2 vs 1"
),
tibble(
beta = post_cat_weak_pre_response$beta_dose3,
threshold = "Dose 3 vs 2"
),
tibble(
beta = post_cat_weak_pre_response$beta_dose4,
threshold = "Dose 4 vs 3"
),
tibble(
beta = post_cat_weak_pre_response$beta_dose5plus,
threshold = "Dose ≥5 vs 4"
)
) |>
mutate(
threshold = factor(
threshold,
levels = c(
"Dose 2 vs 1",
"Dose 3 vs 2",
"Dose 4 vs 3",
"Dose ≥5 vs 4"
)
)
)
cat_pre_response_summary <-
posterior_cat_pre_response_long |>
group_by(threshold) |>
summarise(
median_beta = median(beta),
lower = quantile(beta, 0.055),
upper = quantile(beta, 0.945),
p_gt0 = mean(beta > 0),
median_or = exp(median_beta),
.groups = "drop"
)
cat_pre_response_summary# A tibble: 4 × 6
threshold median_beta lower upper p_gt0 median_or
<fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 Dose 2 vs 1 0.383 0.0364 0.747 0.961 1.47
2 Dose 3 vs 2 0.00321 -0.325 0.336 0.506 1.00
3 Dose 4 vs 3 0.200 -0.140 0.540 0.820 1.22
4 Dose ≥5 vs 4 0.234 -0.111 0.582 0.856 1.26Show R Code
plot_cat_density_weak_pre_response <-
ggplot(posterior_cat_pre_response_long, aes(x = beta)) +
geom_density(
fill = "steelblue",
alpha = 0.35,
color = "steelblue4",
linewidth = 0.9
) +
geom_vline(
xintercept = 0,
linetype = "dashed",
color = "red",
linewidth = 0.9
) +
geom_vline(
data = cat_pre_response_summary,
aes(xintercept = lower),
linetype = "dotted",
color = "grey40",
linewidth = 0.7
) +
geom_vline(
data = cat_pre_response_summary,
aes(xintercept = upper),
linetype = "dotted",
color = "grey40",
linewidth = 0.7
) +
geom_text(
data = cat_pre_response_summary,
aes(
x = Inf,
y = Inf,
label = paste0(
"89% CrI: [",
round(lower, 2),
", ",
round(upper, 2),
"]\n",
"P(β > 0) = ",
round(p_gt0 * 100, 0),
"%\n",
"Median OR = ",
round(median_or, 2)
)
),
hjust = 1.05,
vjust = 1.4,
size = 3.2,
inherit.aes = FALSE
) +
facet_wrap(~ threshold, scales = "free_y") +
labs(
title = "Posterior Distributions: Pre-Response Categorical Dose Threshold Effects",
subtitle = "Weakly informative prior; vertical dotted lines = 89% credible interval",
x = "Threshold-specific β coefficient",
y = "Density"
) +
theme_minimal(base_size = 13) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
strip.text = element_text(face = "bold")
)
plot_cat_density_weak_pre_response
4 Conclusion
Across these analyses, the categorical dose specification provides a useful complement to the linear and log-dose formulations by asking a different clinical question: not whether response increases smoothly with each additional dose, but whether meaningful gains appear after crossing specific exposure thresholds. In that sense, the model is less about estimating a single slope and more about evaluating whether the treatment–response relationship may be organized around clinically relevant milestones in exposure.
Several broad themes emerge from this workflow. First, the posterior predictive checks and model diagnostics support the adequacy of the categorical framework as a reasonable representation of the observed ORR data. Second, comparison of skeptical and weakly informative priors helps clarify which threshold-specific signals are robust to prior choice and which remain more uncertain. This is particularly important in a threshold model, where additional flexibility can improve interpretability but also increase susceptibility to over-reading sparse data patterns. Third, the dose-intensity sensitivity analysis adds an important temporal refinement by restricting attention to treatment intensity accrued before the response assessment window, thereby addressing the concern that higher cumulative exposure may partly reflect the consequence rather than the cause of response.
Taken together, these analyses suggest that any apparent dose–response pattern should be interpreted as graded and uncertain rather than as evidence for a single sharp biological cutoff. The categorical formulation is valuable because it can reveal where the posterior support for benefit appears strongest, but the threshold terms should not be mistaken for independently proven mechanistic breakpoints. Instead, they are best understood as model-based summaries of where clinically meaningful changes in response probability may be occurring within the observed cohort.
More broadly, the convergence or divergence of results across the skeptical, weakly informative, and dose-intensity–adjusted models is itself informative. When a threshold effect persists across specifications, confidence increases that the signal reflects a stable feature of the data rather than an artifact of a particular prior or exposure definition. When the posterior shifts materially after prior relaxation or dose-intensity adjustment, that finding is also scientifically useful, because it identifies which inferences remain sensitive to modeling assumptions and therefore should be interpreted cautiously.
In this way, the categorical dose analysis contributes not just a set of ORR estimates, but a more nuanced understanding of exposure–response structure. It helps distinguish between three possibilities: a broadly monotonic association, a threshold-like pattern, or substantial uncertainty that precludes strong claims about either. That perspective is especially useful for a clinical audience, where decisions are rarely based on whether a coefficient is simply above or below zero, but rather on whether the totality of evidence supports a credible and interpretable treatment pattern.
Overall, this document supports the main manuscript by showing that the categorical dose framework can be fit coherently, checked transparently, and stress-tested through prior and dose-intensity sensitivity analyses. The resulting inferences should be viewed as complementary to the continuous-dose models rather than competitive with them. Together, these approaches provide a richer picture of how treatment exposure may relate to response, while keeping uncertainty, model dependence, and causal timing considerations in clear view.