Bayesian ORR Modeling — Linear Dose Specification

1 Overview

This page presents the Bayesian objective response rate (ORR) analysis for the linear dose specification. The goal is to estimate how treatment exposure relates to response in the observed study cohort under the assumption that each additional dose contributes a constant incremental effect on the log-odds of response.

All models are fit in R using rethinking::ulam(). The model structure shown here was developed in a separate generative workflow before being applied to the real data. In that preliminary stage, simulated datasets were used to examine prior behavior, evaluate parameter recovery, and assess whether the linear dose formulation could recover plausible exposure–response patterns under controlled conditions. The analyses below therefore represent the application of a pre-specified and previously stress-tested linear modeling framework to the clinical cohort.

The primary adjusted ORR model evaluates dose exposure while accounting for key clinical covariates, including age, immunosuppression status, tumor location, stage, treatment agent, and ECOG performance status.

The linear dose specification is the most direct and interpretable of the dose–response formulations considered in this project. Under this model, each additional dose is assumed to shift the log-odds of response by the same amount, regardless of whether that increase occurs early or late in the treatment course. This makes the linear formulation a useful reference model against which more flexible alternatives can be compared.

We also include prior-sensitivity analyses and dose-intensity sensitivity analyses derived from treatment timing and pre-response exposure. These analyses test whether the estimated linear association between number of doses and response is robust to reasonable changes in prior assumptions and to alternative representations of treatment delivery.

Taken together, this document is designed to complement the main paper by showing not only the ORR estimates from the linear dose model, but also the modeling choices, diagnostic checks, and sensitivity analyses that support their interpretation.

source(here::here("scripts", "load_packages.R"))

library(tidyverse)
library(rethinking)
library(cmdstanr)
library(posterior)

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

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

df_model_raw_2 <- readRDS(
  here::here("files", "Dose Intensity Sensitivity",
             "Dose Intensity Sensitivity Model Data-2026-04-29-08-52-22.rds")
)

extra <- df_model_raw |> anti_join(df_model_raw_2, by = join_by(MRN)
                              )

df_analysis <- df_model_raw |>
  select(-patient, -MRN) |>
  mutate(record_id = row_number())

df_analysis <- df_analysis |>
  mutate(
    # -----------------------------
    # Outcome
    # -----------------------------
    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 variables
    # -----------------------------
    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
    # -----------------------------
    age_centered = age_at_treatment - 76,

    # -----------------------------
    # Immunosuppression
    # -----------------------------
    immunosuppressed = factor(
      immune_status == "Immunosuppressed",
      levels = c(FALSE, TRUE)
    ),

    # -----------------------------
    # Location
    # -----------------------------
    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
    # -----------------------------
    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
    # -----------------------------
    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 variables
    # -----------------------------
    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")
    ),

    # Newer schedule-aware agent variable from yesterday
    agent_schedule = case_when(
      agent_schedule %in% c("cemiplimab", "Cemiplimab") ~ "cemiplimab",
      agent_schedule %in% c("pembro_q3", "Pembrolizumab_q3", "pembro q3") ~ "pembro_q3",
      agent_schedule %in% c("pembro_q6", "Pembrolizumab_q6", "pembro q6") ~ "pembro_q6",
      agent_schedule %in% c("pembro_mixed", "pembro mixed") ~ "pembro_mixed",
      agent_schedule %in% c("Ipi-nivo", "ipi-nivo", "Ipi+Nivo", "ipinivo") ~ "Ipi-nivo",
      TRUE ~ NA_character_
    ),
    agent_schedule = factor(
      agent_schedule,
      levels = c("cemiplimab", "pembro_q3", "pembro_q6", "pembro_mixed", "Ipi-nivo")
      ),

    # -----------------------------
    # Surgery
    # -----------------------------
    surgery = as.integer(has_surg)
  )

df_analysis <- df_analysis |> 
  filter(
    !is.na(num_doses),
    num_doses >= 1
  )

df_analysis <- df_analysis |>
  select(
    record_id,
    response,
    num_doses,
    dose_minus_1,
    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
  )

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_centered)),
    mean_centered = mean(dose_intensity_pre_response_uncapped_centered, na.rm = TRUE),
    sd_centered = sd(dose_intensity_pre_response_uncapped_centered, na.rm = TRUE),
    min_centered = min(dose_intensity_pre_response_uncapped_centered, na.rm = TRUE),
    max_centered = max(dose_intensity_pre_response_uncapped_centered, na.rm = TRUE)
  )

2.0.1 Variables

colnames(df_analysis)

 [1] "record_id"                                    
 [2] "response"                                     
 [3] "num_doses"                                    
 [4] "dose_minus_1"                                 
 [5] "log_dose"                                     
 [6] "dose_group"                                   
 [7] "dose_2"                                       
 [8] "dose_3"                                       
 [9] "dose_4"                                       
[10] "dose_5plus"                                   
[11] "age_centered"                                 
[12] "sex"                                          
[13] "immunosuppressed"                             
[14] "location_condensed"                           
[15] "stage_condensed"                              
[16] "agent"                                        
[17] "agent_schedule"                               
[18] "ecog_condensed"                               
[19] "expected_interval_days"                       
[20] "likely_schedule"                              
[21] "num_doses_recorded"                           
[22] "num_doses_pre_response"                       
[23] "dose_intensity_final"                         
[24] "dose_intensity_pre_response"                  
[25] "dose_intensity_pre_response_uncapped"         
[26] "mean_interval_days"                           
[27] "interval_deviation_days"                      
[28] "extra_doses_after_response"                   
[29] "dose_intensity_source"                        
[30] "surgery"                                      
[31] "response_cutoff_date"                         
[32] "dose_intensity_pre_response_uncapped_centered"

vars_main <- c(
  "response",
  "num_doses",
  "dose_minus_1",
  "log_dose",
  "age_centered",
  "immunosuppressed",
  "location_condensed",
  "stage_condensed",
  "agent",
  "agent_schedule",
  "ecog_condensed"
)

df_analysis |>
  select(all_of(vars_main)) |>
  glimpse()

Rows: 190
Columns: 11
$ response           <dbl> 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 1, …
$ num_doses          <int> 5, 16, 2, 5, 46, 2, 6, 8, 3, 4, 2, 2, 2, 4, 2, 2, 2…
$ dose_minus_1       <dbl> 4, 15, 1, 4, 45, 1, 5, 7, 2, 3, 1, 1, 1, 3, 1, 1, 2…
$ log_dose           <dbl> 1.6094379, 2.7725887, 0.6931472, 1.6094379, 3.82864…
$ age_centered       <dbl> 10, 11, -1, -43, 7, 4, 4, 17, 17, -5, -10, -9, 16, …
$ immunosuppressed   <fct> FALSE, TRUE, FALSE, FALSE, FALSE, FALSE, FALSE, TRU…
$ location_condensed <fct> Unknown Primary, Unknown Primary, Head/Neck, Head/N…
$ stage_condensed    <fct> Unstaged, Unstaged, Stage IV, Stage I-II, Stage III…
$ agent              <fct> Cemiplimab, Cemiplimab, Cemiplimab, Pembrolizumab, …
$ agent_schedule     <fct> cemiplimab, cemiplimab, cemiplimab, pembro_q3, cemi…
$ ecog_condensed     <fct> ECOG 0-1, ECOG 0-1, ECOG 0-1, ECOG 2-3, ECOG 0-1, E…

2.0.2 Agents Used

df_analysis <- df_analysis |>
  mutate(
    location_condensed = relevel(location_condensed, ref = "Head/Neck"),
    stage_condensed = relevel(stage_condensed, ref = "Stage I-II"),
    ecog_condensed = relevel(ecog_condensed, ref = "ECOG 0-1"),
    agent_schedule = relevel(agent_schedule, ref = "cemiplimab")
  )

df_analysis |>
  count(agent_schedule, sort = TRUE)

# A tibble: 6 × 2
  agent_schedule     n
  <fct>          <int>
1 cemiplimab       129
2 pembro_q3         45
3 pembro_mixed       9
4 pembro_q6          5
5 Ipi-nivo           1
6 <NA>               1

df_analysis <- df_analysis |>
  mutate(

    # ECOG dummy
    ecog23 = as.integer(ecog_condensed == "ECOG 2-3"),

    # Location dummies
    loc_conj = as.integer(location_condensed == "Conjunctiva"),
    loc_other = as.integer(location_condensed == "Other"),
    loc_unknown = as.integer(location_condensed == "Unknown Primary"),

    # Stage dummies
    stage_III = as.integer(stage_condensed == "Stage III"),
    stage_IV = as.integer(stage_condensed == "Stage IV"),
    stage_unstaged = as.integer(stage_condensed == "Unstaged"),

    # Schedule-aware agent dummies
    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))

df_analysis |>
  summarise(
    n_agent_schedule_na = sum(is.na(agent_schedule)),
    n_q3_na = sum(is.na(agent_pembro_q3)),
    n_q6_na = sum(is.na(agent_pembro_q6)),
    n_mixed_na = sum(is.na(agent_pembro_mixed)),
    n_ipi_na = sum(is.na(agent_ipinivo))
  )

# A tibble: 1 × 5
  n_agent_schedule_na n_q3_na n_q6_na n_mixed_na n_ipi_na
                <int>   <int>   <int>      <int>    <int>
1                   0       0       0          0        0

df_analysis |>
  count(agent_schedule, agent_pembro_q3, agent_pembro_q6, agent_pembro_mixed, agent_ipinivo)

# A tibble: 5 × 6
  agent_schedule agent_pembro_q3 agent_pembro_q6 agent_pembro_mixed
  <fct>                    <int>           <int>              <int>
1 cemiplimab                   0               0                  0
2 pembro_q3                    1               0                  0
3 pembro_q6                    0               1                  0
4 pembro_mixed                 0               0                  1
5 Ipi-nivo                     0               0                  0
# ℹ 2 more variables: agent_ipinivo <int>, n <int>

# ------------------------------------------------------------------------------
# Diagnose why the linear model includes fewer patients than df_analysis
# ------------------------------------------------------------------------------

model_vars <- c(
  "response",
  "dose_minus_1",
  "age_centered",
  "immunosuppressed",
  "ecog_condensed",
  "location_condensed",
  "stage_condensed",
  "agent_schedule"
)

# Count missing values in each variable needed for the model
df_analysis |>
  summarise(
    n_total = n(),
    across(
      all_of(model_vars),
      ~ sum(is.na(.x)),
      .names = "missing_{.col}"
    )
  )

# A tibble: 1 × 9
  n_total missing_response missing_dose_minus_1 missing_age_centered
    <int>            <int>                <int>                <int>
1     189                0                    0                    0
# ℹ 5 more variables: missing_immunosuppressed <int>,
#   missing_ecog_condensed <int>, missing_location_condensed <int>,
#   missing_stage_condensed <int>, missing_agent_schedule <int>

dat_linear <- list(

  N = nrow(df_analysis),
  y = df_analysis$response,

  dose = df_analysis$dose_minus_1,
  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)
)

---

3 Fit model using skeptical priors

m_linear_real_skeptical <- ulam(
  alist(
    y ~ bernoulli(p),

    logit(p) <- alpha +
      beta_dose * dose +
      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_dose ~ normal(0, 0.15),
    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_linear,
  chains = 4,
  cores = 4,
  iter = 2000
)

save_files(
  save_object = m_linear_real_skeptical,
  filename = "Model Linear Skeptical Prior",
  subD = "files/Modeling/Model Linear Skeptical Prior"
)

Step 2 — Posterior Predictive Check

This asks:

If the model were true, could it generate data that look like our data?

For your binary outcome, a simple check is the predicted response rate distribution.

Extract posterior predictions

post <- extract.samples(m_linear_real_skeptical)

p <- link(m_linear_real_skeptical)

link() gives posterior predicted probabilities for every patient.

# ------------------------------------------------------------------------------
# Manual posterior predictive check from fitted probabilities
# ------------------------------------------------------------------------------

set.seed(123)

p_hat_skeptical <- link(m_linear_real_skeptical)

stopifnot(ncol(p_hat_skeptical) == nrow(df_analysis))

n_ppc_draws <- 1000

draw_ids <- sample(
  seq_len(nrow(p_hat_skeptical)),
  size = n_ppc_draws,
  replace = FALSE
)

p_hat_ppc_skeptical <- p_hat_skeptical[draw_ids, , drop = FALSE]

yrep_skeptical <- matrix(
  rbinom(
    n = length(p_hat_ppc_skeptical),
    size = 1,
    prob = as.vector(p_hat_ppc_skeptical)
  ),
  nrow = nrow(p_hat_ppc_skeptical),
  ncol = ncol(p_hat_ppc_skeptical)
)

observed_orr_skeptical <- mean(df_analysis$response)

ppc_skeptical_orr <- tibble(
  draw = seq_len(nrow(yrep_skeptical)),
  simulated_orr = rowMeans(yrep_skeptical)
)

ppc_skeptical_orr |>
  summarise(
    observed_orr = observed_orr_skeptical,
    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_skeptical),
    mean_fitted_probability = mean(p_hat_skeptical)
  )

# 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.720
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>

This generates simulated response outcomes under the model.

Compare simulated vs observed response rate

ggplot(ppc_skeptical_orr, aes(x = simulated_orr)) +
  geom_histogram(
    bins = 50,
    fill = "steelblue",
    alpha = 0.7,
    color = "white"
  ) +
  geom_vline(
    xintercept = observed_orr_skeptical,
    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: Linear 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")
  )

## Posterior predictive check: rank distribution

3.1 Posterior parameter interpretation

## Model diagnostics

precis(m_linear_real_skeptical, depth = 2)

                          mean         sd        5.5%       94.5%      rhat
alpha              0.783717406 0.28210796  0.34726377  1.23507416 1.0004369
beta_dose          0.080135475 0.04294753  0.01292032  0.15114153 1.0021024
beta_age           0.001860411 0.01203287 -0.01713777  0.02106461 1.0009571
beta_imm          -0.908347149 0.17933848 -1.18748043 -0.61986475 1.0015891
beta_ecog         -0.306479570 0.26375136 -0.72138834  0.11565054 1.0029821
beta_conj         -2.294029784 0.42253180 -2.97195818 -1.62232577 1.0018507
beta_other        -0.148455789 0.28100103 -0.60698136  0.30482237 1.0013502
beta_unknown       0.054025007 0.28817600 -0.40469888  0.52466166 1.0016416
beta_stageIII     -0.008866530 0.23966303 -0.37789314  0.37348807 1.0004649
beta_stageIV      -0.319365735 0.23815099 -0.69922535  0.04842991 1.0027716
beta_stageUn       0.246070613 0.27435976 -0.19059212  0.68770620 1.0004301
beta_pembro_q3    -0.044829495 0.23677315 -0.42519951  0.33408142 1.0009871
beta_pembro_q6    -0.050229657 0.28510201 -0.50575405  0.40274564 1.0024390
beta_pembro_mixed  0.228135910 0.27356704 -0.21263900  0.66210886 0.9996578
beta_ipinivo       0.425415847 0.30627141 -0.06863530  0.90356926 0.9998916
                  ess_bulk
alpha             4564.860
beta_dose         5949.483
beta_age          6244.270
beta_imm          6340.662
beta_ecog         6061.341
beta_conj         6452.860
beta_other        5872.522
beta_unknown      6182.496
beta_stageIII     4399.533
beta_stageIV      5006.182
beta_stageUn      6249.769
beta_pembro_q3    6028.561
beta_pembro_q6    6546.577
beta_pembro_mixed 6943.551
beta_ipinivo      6201.655

plot(precis(m_linear_real_skeptical))

3.2 Extract paramaters from the model

post <- extract.samples(m_linear_real_skeptical)

3.3 Intercept

# Convert intercept to baseline probability
baseline_prob <- plogis(post$alpha)

# Summarize
mean(baseline_prob)

[1] 0.6833745

quantile(baseline_prob, c(0.055, 0.945))

     5.5%     94.5% 
0.5859539 0.7747054 

3.4 Step 4 — Posterior predicted dose–response curve

doses <- 0:14

# ------------------------------------------------------------------------------
# Create prediction dataset for the dose-response curve
# ------------------------------------------------------------------------------
# The goal is to estimate the model-adjusted expected response probability across
# increasing cumulative ICI doses while holding all other covariates fixed.
#
# In this model, dose is coded as dose_minus_1:
#   dose = 0 corresponds to 1 actual dose
#   dose = 1 corresponds to 2 actual doses
#   dose = 2 corresponds to 3 actual doses
#   etc.
#
# All other covariates are set to 0 to define the model reference population.
#
# Because age was centered as:
#   age_centered = age_at_treatment - 76
# age = 0 corresponds to age 76.
#
# For binary indicator variables, 0 means absence of that characteristic or
# membership in the omitted/reference category.
#
# Therefore, this prediction dataset represents the following reference-profile
# patient:
#   - 1 to 15 cumulative ICI doses, depending on the row
#   - age 76
#   - not immunosuppressed
#   - ECOG 0-1
#   - Head/Neck primary tumor location
#   - Stage I-II disease
#   - cemiplimab treatment schedule
#
# The resulting predictions should be interpreted as model-adjusted expected
# response probabilities at each dose level for this reference covariate profile,
# not as marginal predictions averaged over the observed cohort distribution.

new_data <- data.frame(
  dose = doses,

  # age = 0 means age 76, because age was centered at 76.
  age = 0,

  # imm = 0 means not immunosuppressed.
  imm = 0,

  # ecog23 = 0 means ECOG 0-1, the reference ECOG category.
  ecog23 = 0,

  # All location indicators set to 0 means Head/Neck,
  # the reference location category.
  loc_conj = 0,
  loc_other = 0,
  loc_unknown = 0,

  # All stage indicators set to 0 means Stage I-II,
  # the reference stage category.
  stage_III = 0,
  stage_IV = 0,
  stage_unstaged = 0,

  # All agent/schedule indicators set to 0 means cemiplimab,
  # the reference treatment category.
  agent_pembro_q3 = 0,
  agent_pembro_q6 = 0,
  agent_pembro_mixed = 0,
  agent_ipinivo = 0
)


# ------------------------------------------------------------------------------
# Generate posterior expected response probabilities for the reference population
# ------------------------------------------------------------------------------
# This line uses the fitted Bayesian logistic regression model to predict the
# expected probability of response for each row in `new_data`.
#
# `m_linear_real_skeptical` is the fitted model. It contains many posterior draws
# of the model parameters, such as:
#   - alpha
#   - beta_dose
#   - beta_age
#   - beta_imm
#   - etc.
#
# Each posterior draw represents one plausible version of the dose-response model
# after combining the observed data with the priors.
#
# `new_data` is the artificial prediction dataset created above. In that dataset,
# dose varies across the desired dose range, while all adjustment covariates are
# held fixed at their reference values.
#
# `link()` takes each posterior draw of the fitted model, plugs in the covariate
# values from `new_data`, and calculates the model-implied response probability.
#
# Intuitively, this asks:
#   "For a reference-profile patient, what response probability would the model
#    predict at each cumulative dose level?"
#
# Because the model is logistic, predictions are first calculated on the log-odds
# scale and then automatically transformed back to the probability scale.
#
# The result, `p`, is a matrix:
#   - each row is one posterior draw, or one plausible fitted model
#   - each column is one row of `new_data`, here one dose level
#
# So if there are 4,000 posterior draws and 15 dose levels, `p` will have
# approximately 4,000 rows and 15 columns.
#
# The values inside `p` are expected response probabilities, not simulated
# individual response outcomes. For example, a value of 0.72 means that, under
# that posterior draw and dose level, the expected probability of response is 72%.
#
# Later, we summarize `p` across posterior draws to obtain the posterior mean
# predicted response probability and the 89% credible interval at each dose level.

p <- link(m_linear_real_skeptical, data = new_data)

mean_p <- apply(p, 2, mean)
ci <- apply(p, 2, PI, 0.89)

plot_df_linear_skeptical<- tibble(
  dose_number = doses + 1,   # if doses = 0:15 and represents dose_minus_1
  mean_p = mean_p,
  lower = ci[1, ],
  upper = ci[2, ]
)

save_files(
  save_object = plot_df_linear_skeptical,
  filename = "Plot DF Predicted Response Linear Skeptical Prior",
  subD = "files/Modeling/Plot DF Predicted Response Linear Skeptical Prior"
)

plot_predicted_response_linear_skeptical <-
  ggplot(plot_df_linear_skeptical, aes(x = dose_number, y = mean_p)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    fill = "steelblue",
    alpha = 0.22
  ) +
  geom_line(
    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_linear_skeptical$dose_number),
    max(plot_df_linear_skeptical$dose_number),
    by = 1
  ),
  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\nLinear dose 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_linear_skeptical

save_files(
  save_object = plot_predicted_response_linear_skeptical,
  filename = "Predicted Response Probability Linear Model Skeptical Prior",
  subD = "files/Modeling/Predicted Response Probability Linear Model Skeptical Prior"
)

save_files(
  save_object = plot_predicted_response_linear_skeptical,
  filename = "Predicted Response Probability Linear Model Skeptical Prior",
  subD = "files/Modeling/Predicted Response Probability Linear Model Skeptical Prior",
  extension = ".png",
  width = 8,
  height = 6,
  dpi = 300,
  units = "in"
)

post_linear_skeptical <- extract.samples(m_linear_real_skeptical)

posterior_dose_linear_skeptical <- post_linear_skeptical$beta_dose

save_files(
  save_object = posterior_dose_linear_skeptical,
  filename = "Posterior Dose Linear Skeptical Prior",
  subD = "files/Modeling/Posterior Dose Linear Skeptical Prior"
)

# Probability statements

or_threshold_105 <- 1.05
or_threshold_110 <- 1.10

beta_threshold_105 <- log(or_threshold_105)
beta_threshold_110 <- log(or_threshold_110)

p_gt_0   <- mean(posterior_dose_linear_skeptical > 0)
p_or_105 <- mean(posterior_dose_linear_skeptical >= beta_threshold_105)
p_or_110 <- mean(posterior_dose_linear_skeptical >= beta_threshold_110)

P(β > 0) = 98% → Strong evidence for a positive dose effect

P(OR per additional dose ≥ 1.05) = 76% → Good evidence for at least a modest effect

P(OR per additional dose ≥ 1.10) = 35% → More limited evidence for a larger effect

# Density object
dens <- density(posterior_dose_linear_skeptical)
dens_df <- data.frame(x = dens$x, y = dens$y)
ymax <- max(dens_df$y)

# 89% credible interval
ci_lower <- quantile(posterior_dose_linear_skeptical, 0.055)
ci_upper <- quantile(posterior_dose_linear_skeptical, 0.945)

# Tail probabilities
p_gt_0   <- mean(posterior_dose_linear_skeptical > 0)
p_or_105 <- mean(posterior_dose_linear_skeptical >= beta_threshold_105)
p_or_110 <- mean(posterior_dose_linear_skeptical >= beta_threshold_110)

# Colors
col_gt_0 <- "#1B9E77"
col_005  <- "#1F78B4"
col_010  <- "#D95F02"
col_line <- "#08306B"

plot_linear_skeptical <-
  ggplot() +

  # Density outline
  geom_line(
    data = dens_df,
    aes(x = x, y = y),
    linewidth = 0.9,
    color = col_line
  ) +

  # Piecewise shading
  geom_area(
    data = subset(
      dens_df,
      x > 0 & x < beta_threshold_105
    ),
    aes(x = x, y = y),
    fill = col_gt_0,
    alpha = 0.35
  ) +
  geom_area(
    data = subset(
      dens_df,
      x >= beta_threshold_105 & x < beta_threshold_110
    ),
    aes(x = x, y = y),
    fill = col_005,
    alpha = 0.45
  ) +
  geom_area(
    data = subset(
      dens_df,
      x >= beta_threshold_110
    ),
    aes(x = x, y = y),
    fill = col_010,
    alpha = 0.50
  ) +

  # Reference lines
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "red",
    linewidth = 1
  ) +
  geom_vline(
    xintercept = c(ci_lower, ci_upper),
    linetype = "dotted",
    color = "grey40",
    linewidth = 0.8
  ) +

  # Labels
  annotate(
    "text",
    x = 0,
    y = ymax * 0.95,
    label = "Null",
    hjust = -0.2,
    color = "red"
  ) +
  annotate(
    "text",
    x = mean(posterior_dose_linear_skeptical),
    y = ymax * 0.40,
    label = paste0(
      "89% CrI: [",
      round(ci_lower, 2),
      ", ",
      round(ci_upper, 2),
      "]"
    ),
    size = 3.4
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax * 0.95,
    label = paste0(
      "P(β > 0 | ↑ dose → ↑ odds of response) = ",
      round(p_gt_0 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_gt_0
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax * 0.85,
    label = paste0(
      "P(OR per additional dose ≥ ",
      sprintf("%.2f", or_threshold_105),
      ") = ",
      round(p_or_105 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_005
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax * 0.75,
    label = paste0(
      "P(OR per additional dose ≥ ",
      sprintf("%.2f", or_threshold_110),
      ") = ",
      round(p_or_110 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_010
  ) +

  labs(
    title = "Posterior Distribution: Linear Dose Effect",
    subtitle = "Skeptical Prior",
    x = "β_dose (log-odds per dose)",
    y = "Density"
  ) +
  theme_minimal() +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.28))) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

plot_linear_skeptical

save_files(
  save_object = plot_linear_skeptical,
  filename = "Density Plot Linear Model Skeptical Priors",
  subD = "files/Modeling/Density Plot Linear Model Skeptical Priors"
)

save_files(
  save_object = plot_linear_skeptical,
  filename = "Density Plot Linear Model Skeptical Priors",
  subD = "files/Modeling/Density Plot Linear Model Skeptical Priors",
  extension = ".png",
  width = 7.6,
  height = 7.5,
  dpi = 300,
  units = "in"
)

---

4 Fit model with weakly informative prior

m_linear_real_weak <- ulam(
  alist(
    y ~ bernoulli(p),

    logit(p) <- alpha +
      beta_dose * dose +
      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_dose ~ 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_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_linear,
  chains = 4,
  cores = 4,
  iter = 2000,
  cmdstan = TRUE
)

#save_files(
#  save_object = m_linear_real_weak,
#  filename = "Model Linear Weak Prior",
#  subD = "files/Modeling/Model Linear Weak Prior"
#)

m_linear_real_weak@cstanfit$save_output_files(
  dir = here::here("files", "Modeling", "cmdstan_output", "model_linear_weak")
)

saveRDS(
  m_linear_real_weak,
  here::here(
    "files", "Modeling", "Model Linear Weak Prior",
    paste0("Model Linear Weak Prior-", format(Sys.time(), "%Y-%m-%d-%H-%M-%S"), ".rds")
  )
)

saveRDS(
  m_linear_real_weak,
  here::here(
    "files",
    "Modeling",
    "Model Linear Weak Prior",
    paste0("Model Linear Weak Prior-", format(Sys.time(), "%Y-%m-%d-%H-%M-%S"), ".rds")
  )
)

4.1 Posterior predtive check of model with weak prior

post_linear_weak <- extract.samples(m_linear_real_weak)


saveRDS(
  post_linear_weak,
  here::here(
    "files", "Modeling", "Posterior Dose Linear Weak Prior",
    paste0("Posterior Linear Weak Prior-", format(Sys.time(), "%Y-%m-%d-%H-%M-%S"), ".rds")
  )
)
p_linear_weak <- link(m_linear_real_weak)
sim_response_linear_weak <- sim(m_linear_real_weak)

posterior_alpha_linear_weak <- post_linear_weak$alpha

save_files(
  save_object = posterior_alpha_linear_weak,
  filename = "Posterior Alpha Linear Weak Prior",
  subD = "files/Modeling/Posterior Alpha Linear Weak Prior"
)

mean(df_analysis$response)

[1] 0.6455026

mean(sim_response_linear_weak)

[1] 0.797

hist(
  colMeans(sim_response_linear_weak),
  main = "Posterior Predictive Response Rate",
  xlab = "Simulated ORR"
)

abline(v = mean(df_analysis$response), col = "red", lwd = 2)

4.2 Posterior Predictive Rank (T-rank / PIT) Check

pred_rank_linear_weak <- sapply(
  seq_len(ncol(sim_response_linear_weak)),
  function(i) sum(sim_response_linear_weak[, i] < df_analysis$response[i])
)

hist(
  pred_rank_linear_weak,
  main = "Posterior Predictive Rank Distribution (Linear Model, Weak Prior)",
  sub = "Binary outcome: non-responders at rank 0; responders at higher ranks",
  xlab = "Rank of observed outcome within simulated outcomes",
  ylab = "Number of patients",
  col = "grey80",
  border = "white"
)

4.3 Posterior parameter interpretation

precis(m_linear_real_weak, depth = 2)

                          mean         sd        5.5%       94.5%      rhat
alpha              0.753497764 0.26776463  0.33595237  1.17666843 1.0004840
beta_dose          0.090362820 0.04384869  0.02337381  0.16117379 0.9997145
beta_age           0.002384067 0.01195762 -0.01723241  0.02092931 1.0010125
beta_imm          -0.910522447 0.18684516 -1.21282120 -0.61202249 1.0038241
beta_ecog         -0.312688990 0.26239021 -0.72365419  0.11234436 1.0004248
beta_conj         -2.300302179 0.43979953 -3.01930777 -1.60225513 1.0009581
beta_other        -0.156777576 0.26660060 -0.58589971  0.27290606 0.9996569
beta_unknown       0.040936234 0.29456179 -0.43204663  0.50422194 1.0014139
beta_stageIII     -0.006059324 0.23664102 -0.38835061  0.36729411 1.0020135
beta_stageIV      -0.318736535 0.23131880 -0.68954294  0.05393036 1.0007508
beta_stageUn       0.244667000 0.26751876 -0.18079053  0.66726161 1.0029004
beta_pembro_q3    -0.043551760 0.23238334 -0.41280773  0.33822749 1.0004781
beta_pembro_q6    -0.059511050 0.28876610 -0.51033559  0.40239627 1.0006898
beta_pembro_mixed  0.221496878 0.28090668 -0.22745664  0.67115780 1.0004570
beta_ipinivo       0.424702480 0.31010543 -0.06839071  0.91827093 1.0024773
                  ess_bulk
alpha             3872.565
beta_dose         5477.922
beta_age          8313.432
beta_imm          7496.241
beta_ecog         6770.133
beta_conj         8627.662
beta_other        7064.701
beta_unknown      7812.798
beta_stageIII     5267.926
beta_stageIV      5460.316
beta_stageUn      6470.538
beta_pembro_q3    7265.779
beta_pembro_q6    7005.360
beta_pembro_mixed 7784.908
beta_ipinivo      7732.356

plot(precis(m_linear_real_weak))

post_linear_weak <- extract.samples(m_linear_real_weak)

4.4 Intercept to baseline probability

baseline_prob_linear_weak <- plogis(post_linear_weak$alpha)

mean(baseline_prob_linear_weak)

[1] 0.6772046

quantile(baseline_prob_linear_weak, c(0.055, 0.945))

     5.5%     94.5% 
0.5832070 0.7643482 

4.5 Posterior predicted dose-response curve

doses <- 0:14

new_data_linear <- data.frame(
  dose = doses,
  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_linear_weak_new <- link(m_linear_real_weak, data = new_data_linear)

mean_p_linear_weak <- apply(p_linear_weak_new, 2, mean)
ci_linear_weak <- apply(p_linear_weak_new, 2, PI, 0.89)

plot_df_linear_weak <- tibble(
  dose_number = doses + 1,
  mean_p = mean_p_linear_weak,
  lower = ci_linear_weak[1, ],
  upper = ci_linear_weak[2, ]
)

save_files(
  save_object = plot_df_linear_weak,
  filename = "Plot DF Predicted Response Linear Weak Prior",
  subD = "files/Modeling/Plot DF Predicted Response Linear Weak Prior"
)

plot_predicted_response_linear_weak <-
  ggplot(plot_df_linear_weak, aes(x = dose_number, y = mean_p)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    fill = "steelblue",
    alpha = 0.22
  ) +
  geom_line(
    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_linear_weak$dose_number),
      max(plot_df_linear_weak$dose_number),
      by = 1
    ),
    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\nLinear dose 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_linear_weak

save_files(
  save_object = plot_predicted_response_linear_weak,
  filename = "Predicted Response Probability Linear Model Weak Prior",
  subD = "files/Modeling/Predicted Response Probability Linear Model Weak Prior"
)

save_files(
  save_object = plot_predicted_response_linear_weak,
  filename = "Predicted Response Probability Linear Model Weak Prior",
  subD = "files/Modeling/Predicted Response Probability Linear Model Weak Prior",
  extension = ".png",
  width = 8,
  height = 6,
  dpi = 300,
  units = "in"
)

4.6 Posterior dose-effect probabilities

post_linear_weak <- extract.samples(m_linear_real_weak)

posterior_dose_linear_weak <- post_linear_weak$beta_dose

save_files(
  save_object = posterior_dose_linear_weak,
  filename = "Posterior Dose Linear Weak Prior",
  subD = "files/Modeling/Posterior Dose Linear Weak Prior"
)

or_threshold_modest <- 1.05
or_threshold_larger <- 1.10

beta_threshold_modest <- log(or_threshold_modest)
beta_threshold_larger <- log(or_threshold_larger)

p_gt_0_linear_weak <- mean(posterior_dose_linear_weak > 0)

p_or_modest_linear_weak <- mean(
  posterior_dose_linear_weak >= beta_threshold_modest
)

p_or_larger_linear_weak <- mean(
  posterior_dose_linear_weak >= beta_threshold_larger
)

P(β > 0) = 99% → Strong evidence for a positive dose effect

P(OR per additional dose ≥ 1.05) = 82% → Good evidence for at least a modest effect

P(OR per additional dose ≥ 1.1) = 44% → More limited evidence for a larger effect

4.7 Posterior distribution plot

dens_linear_weak <- density(posterior_dose_linear_weak)
dens_df_linear_weak <- data.frame(x = dens_linear_weak$x, y = dens_linear_weak$y)
ymax_linear_weak <- max(dens_df_linear_weak$y)

ci_lower_linear_weak <- quantile(posterior_dose_linear_weak, 0.055)
ci_upper_linear_weak <- quantile(posterior_dose_linear_weak, 0.945)

# OR thresholds expressed on beta/log-odds scale
or_threshold_105_linear_weak <- 1.05
or_threshold_110_linear_weak <- 1.10

beta_threshold_105_linear_weak <- log(or_threshold_105_linear_weak)
beta_threshold_110_linear_weak <- log(or_threshold_110_linear_weak)

p_gt_0_linear_weak <- mean(posterior_dose_linear_weak > 0)

p_or_105_linear_weak <- mean(
  posterior_dose_linear_weak >= beta_threshold_105_linear_weak
)

p_or_110_linear_weak <- mean(
  posterior_dose_linear_weak >= beta_threshold_110_linear_weak
)

col_gt_0 <- "#1B9E77"
col_005  <- "#1F78B4"
col_010  <- "#D95F02"
col_line <- "#08306B"

plot_linear_weak <-
  ggplot() +
  geom_line(
    data = dens_df_linear_weak,
    aes(x = x, y = y),
    linewidth = 0.9,
    color = col_line
  ) +
  geom_area(
    data = subset(
      dens_df_linear_weak,
      x > 0 & x < beta_threshold_105_linear_weak
    ),
    aes(x = x, y = y),
    fill = col_gt_0,
    alpha = 0.35
  ) +
  geom_area(
    data = subset(
      dens_df_linear_weak,
      x >= beta_threshold_105_linear_weak &
        x < beta_threshold_110_linear_weak
    ),
    aes(x = x, y = y),
    fill = col_005,
    alpha = 0.45
  ) +
  geom_area(
    data = subset(
      dens_df_linear_weak,
      x >= beta_threshold_110_linear_weak
    ),
    aes(x = x, y = y),
    fill = col_010,
    alpha = 0.50
  ) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
  geom_vline(
    xintercept = c(ci_lower_linear_weak, ci_upper_linear_weak),
    linetype = "dotted",
    color = "grey40",
    linewidth = 0.8
  ) +
  annotate(
    "text",
    x = 0,
    y = ymax_linear_weak * 0.95,
    label = "Null",
    hjust = -0.2,
    color = "red"
  ) +
  annotate(
    "text",
    x = mean(posterior_dose_linear_weak),
    y = ymax_linear_weak * 0.40,
    label = paste0(
      "89% CrI: [",
      round(ci_lower_linear_weak, 2),
      ", ",
      round(ci_upper_linear_weak, 2),
      "]"
    ),
    size = 3.4
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_linear_weak * 0.95,
    label = paste0(
      "P(β > 0 | ↑ dose → ↑ odds of response) = ",
      round(p_gt_0_linear_weak * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_gt_0
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_linear_weak * 0.85,
    label = paste0(
      "P(OR per additional dose ≥ ",
      or_threshold_105_linear_weak,
      ") = ",
      round(p_or_105_linear_weak * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_005
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_linear_weak * 0.75,
    label = paste0(
      "P(OR per additional dose ≥ ",
      or_threshold_110_linear_weak,
      ") = ",
      round(p_or_110_linear_weak * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_010
  ) +
  labs(
    title = "Posterior Distribution: Linear Dose Effect",
    subtitle = "Weakly informative prior",
    x = "β_dose (log-odds per dose)",
    y = "Density"
  ) +
  theme_minimal() +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.28))) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

plot_linear_weak

save_files(
  save_object = plot_linear_weak,
  filename = "Density Plot Linear Model Weak Prior",
  subD = "files/Modeling/Density Plot Linear Model Weak Prior"
)

save_files(
  save_object = plot_linear_weak,
  filename = "Density Plot Linear Model Weak Prior",
  subD = "files/Modeling/Density Plot Linear Model Weak Prior",
  extension = ".png",
  width = 7.6,
  height = 7.5,
  dpi = 300,
  units = "in"
)

---

4.8 Compare dose effect based on prior

delta_beta_linear <- posterior_dose_linear_weak - posterior_dose_linear_skeptical

save_files(
  save_object = delta_beta_linear,
  filename = "Posterior Difference Linear Dose Effect Weak vs Skeptical",
  subD = "files/Modeling/Posterior Difference Linear Dose Effect Weak vs Skeptical"
)

prior_sensitivity_linear <- tibble(
  prior = c("Skeptical", "Weakly informative"),
  mean_beta = c(
    mean(posterior_dose_linear_skeptical),
    mean(posterior_dose_linear_weak)
  ),
  median_beta = c(
    median(posterior_dose_linear_skeptical),
    median(posterior_dose_linear_weak)
  ),
  ci_lower_89 = c(
    quantile(posterior_dose_linear_skeptical, 0.055),
    quantile(posterior_dose_linear_weak, 0.055)
  ),
  ci_upper_89 = c(
    quantile(posterior_dose_linear_skeptical, 0.945),
    quantile(posterior_dose_linear_weak, 0.945)
  ),
  p_gt_0 = c(
    mean(posterior_dose_linear_skeptical > 0),
    mean(posterior_dose_linear_weak > 0)
  ),
  p_gt_005 = c(
    mean(posterior_dose_linear_skeptical > 0.05),
    mean(posterior_dose_linear_weak > 0.05)
  ),
  p_gt_010 = c(
    mean(posterior_dose_linear_skeptical > 0.10),
    mean(posterior_dose_linear_weak > 0.10)
  )
)

prior_sensitivity_linear

4.9 Summarize difference directly

tibble(
  delta_mean = mean(delta_beta_linear),
  delta_median = median(delta_beta_linear),
  delta_lower_89 = quantile(delta_beta_linear, 0.055),
  delta_upper_89 = quantile(delta_beta_linear, 0.945),
  p_delta_gt_0 = mean(delta_beta_linear > 0)
)

delta_df <- tibble(delta_beta = delta_beta_linear)

ggplot(delta_df, aes(x = delta_beta)) +
  geom_density(fill = "steelblue", alpha = 0.4) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(
    title = "Difference in Posterior Dose Effect Between Priors",
    subtitle = "Weak prior minus skeptical prior",
    x = "Δβ_dose",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5),
    plot.subtitle = element_text(face = "bold", hjust = 0.5)
  )

Prior sensitivity was evaluated by comparing posterior dose-effect estimates under skeptical and weakly informative priors. The posterior difference in the dose coefficient (Δβ = β_{weak} − β_{skeptical}) was centered near zero (median Δβ ≈ 0.01; 89% CrI -0.09 to 0.11; P(Δβ > 0) = 56%), indicating that the estimated dose–response relationship was largely robust to prior specification.

---

5 Fit model with dose intensity parameter

This section presents a sensitivity analysis that augments the primary linear objective response rate (ORR) model with a dose-intensity term derived from treatment timing prior to response assessment. The purpose is to test whether the apparent association between dose count and response is materially altered after accounting for how tightly treatment was delivered over time.

Conceptually, this analysis separates how much treatment was received from how intensively it was delivered. The primary model focuses on dose count, whereas the present model adds a timing-based exposure measure to evaluate whether the estimated dose–response relationship is robust to this alternative representation of treatment delivery.

We first construct the analysis dataset for the dose-intensity-adjusted model, then fit the model, assess its posterior predictive behavior, and compare the posterior dose-effect estimate with the corresponding estimate from the primary linear specification.

df_analysis_di <- df_analysis |>
  filter(
    !is.na(response),
    !is.na(dose_minus_1),
    !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)
  )

5.1 Build data

dat_linear_di_uncapped <- list(
  N = nrow(df_analysis_di),
  y = df_analysis_di$response,

  dose = df_analysis_di$dose_minus_1,
  age = df_analysis_di$age_centered,
  imm = as.integer(df_analysis_di$immunosuppressed) - 1,
  ecog23 = df_analysis_di$ecog23,

  loc_conj = df_analysis_di$loc_conj,
  loc_other = df_analysis_di$loc_other,
  loc_unknown = df_analysis_di$loc_unknown,

  stage_III = df_analysis_di$stage_III,
  stage_IV = df_analysis_di$stage_IV,
  stage_unstaged = df_analysis_di$stage_unstaged,

  agent_pembro_q3 = df_analysis_di$agent_pembro_q3,
  agent_pembro_q6 = df_analysis_di$agent_pembro_q6,
  agent_pembro_mixed = df_analysis_di$agent_pembro_mixed,
  agent_ipinivo = df_analysis_di$agent_ipinivo,

  di = df_analysis_di$dose_intensity_pre_response_uncapped_centered,

  alpha_mu = qlogis(0.40)
)

5.2 Fit Dose Intensity Model

m_linear_real_weak_di_uncapped <- ulam(
  alist(
    y ~ bernoulli(p),

    logit(p) <- alpha +
      beta_dose * dose +
      beta_di * di +
      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_dose ~ normal(0.10, 0.12),
    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_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_linear_di_uncapped,
  chains = 4,
  cores = 4,
  iter = 2000
)

save_files(
  save_object = m_linear_real_weak_di_uncapped,
  filename = "Model Linear Weak Prior Dose Intensity Uncapped",
  subD = "files/Modeling/Model Linear Weak Prior Dose Intensity Uncapped"
)

5.3 Diagnostics

We assessed model convergence and sampling efficiency using standard MCMC diagnostics, including the potential scale reduction factor (R̂) and effective sample size (ESS). Well-behaved chains should yield R̂ values close to 1.0 and sufficiently large effective sample sizes to support stable posterior inference.

precis(m_linear_real_weak_di_uncapped, depth = 2)

                          mean         sd        5.5%       94.5%      rhat
alpha              0.763032351 0.28352307  0.32016856  1.21924792 1.0002966
beta_dose          0.089554022 0.04474110  0.02148311  0.16500592 0.9998073
beta_di            0.104524570 0.11758891 -0.07921349  0.28967106 1.0005927
beta_age           0.002400799 0.01189692 -0.01662758  0.02138751 1.0005155
beta_imm          -0.914022530 0.17687299 -1.19972618 -0.63610995 1.0034986
beta_ecog         -0.305779744 0.25860013 -0.72014299  0.10734743 1.0011874
beta_conj         -2.301106854 0.44821167 -3.01874272 -1.58301066 1.0064290
beta_other        -0.162788829 0.27496751 -0.60586474  0.28072157 1.0015295
beta_unknown       0.047734893 0.28828874 -0.41021057  0.51827551 1.0016453
beta_stageIII     -0.005630685 0.24084595 -0.39525188  0.37962096 1.0003037
beta_stageIV      -0.316445768 0.23574925 -0.68898227  0.06174781 1.0007057
beta_stageUn       0.239431559 0.27381587 -0.19755287  0.67819268 0.9997691
beta_pembro_q3    -0.047605096 0.23918961 -0.43559284  0.34930442 1.0015029
beta_pembro_q6    -0.053495984 0.29385024 -0.53102159  0.41622040 1.0003742
beta_pembro_mixed  0.226382285 0.28550771 -0.23160291  0.67982791 1.0002912
beta_ipinivo       0.421339739 0.29183885 -0.03889074  0.87995708 1.0000288
                  ess_bulk
alpha             4462.688
beta_dose         5928.060
beta_di           6677.863
beta_age          7495.346
beta_imm          7600.677
beta_ecog         6504.077
beta_conj         6330.339
beta_other        7835.267
beta_unknown      5933.160
beta_stageIII     5457.921
beta_stageIV      5718.174
beta_stageUn      6916.605
beta_pembro_q3    6571.191
beta_pembro_q6    6574.964
beta_pembro_mixed 7248.950
beta_ipinivo      5871.340

Across all model parameters, R̂ values were approximately 1.00, indicating no evidence of chain non-convergence or between-chain discrepancies. Effective sample sizes were uniformly high (generally >5,000), supporting stable estimation of posterior summaries.

Together with the rank trace plots, these diagnostics indicate that the model sampled efficiently and that posterior estimates are numerically reliable from an MCMC standpoint.

From a substantive perspective, the posterior summaries show a modest positive association between dose count (β_dose) and response probability, while the dose-intensity term (β_di) has greater uncertainty and includes zero within its credible interval. This suggests that, after accounting for dose count, there is limited evidence that dose intensity independently contributes to response in this specification.

Several covariates show stronger associations with response, including immunosuppression (β_imm) and concurrent therapy (β_conj), both of which are associated with lower response probability.

5.3.1 MCMC mixing diagnostic: rank trace plots

Before interpreting posterior estimates, we examined chain behavior to confirm that the fitted model sampled efficiently and mixed well across chains.

trankplot(m_linear_real_weak_di_uncapped)

These rank trace plots assess MCMC mixing across chains for each model parameter. Under good mixing, the colored chains should overlap without visible drift, separation, or persistent structure across the plotting window. In this model, the chains appear well mixed across parameters, with no obvious evidence of poor convergence or chain-specific pathologies. This interpretation is consistent with the near-1 rhat values and high effective sample sizes reported in the posterior summary.

plot(precis(m_linear_real_weak_di_uncapped))

Step 4: posterior predictive check

post_linear_weak_di_uncapped <- extract.samples(m_linear_real_weak_di_uncapped)

p_linear_weak_di_uncapped <- link(m_linear_real_weak_di_uncapped)
sim_response_linear_weak_di_uncapped <- sim(m_linear_real_weak_di_uncapped)

mean(df_analysis$response)

[1] 0.6455026

mean(sim_response_linear_weak_di_uncapped)

[1] 0.7961587

hist(
  colMeans(sim_response_linear_weak_di_uncapped),
  main = "Posterior Predictive Response Rate",
  xlab = "Simulated ORR"
)

abline(v = mean(df_analysis$response), col = "red", lwd = 2)

5.4 Posterior Predictive Rank Check

pred_rank_linear_weak_di_uncapped <- sapply(
  seq_len(ncol(sim_response_linear_weak_di_uncapped)),
  function(i) {
    sum(sim_response_linear_weak_di_uncapped[, i] < df_analysis$response[i])
  }
)

hist(
  pred_rank_linear_weak_di_uncapped,
  main = "Posterior Predictive Rank Distribution\nLinear Model + Dose Intensity (Weak Prior)",
  sub = "Binary outcome: non-responders at rank 0; responders at higher ranks",
  xlab = "Rank of observed outcome within simulated outcomes",
  ylab = "Number of patients",
  col = "grey80",
  border = "white",
  cex.main = 0.9,
  cex.sub = 0.8,
  cex.lab = 0.9
)

5.5 extract posterior dose and dose-intensity effects

post_linear_weak_di_uncapped <- extract.samples(m_linear_real_weak_di_uncapped)

posterior_dose_linear_weak_di_uncapped <- post_linear_weak_di_uncapped$beta_dose
posterior_di_linear_weak_uncapped <- post_linear_weak_di_uncapped$beta_di

# OR thresholds expressed on beta/log-odds scale
or_threshold_105 <- 1.05
or_threshold_110 <- 1.10

beta_threshold_105 <- log(or_threshold_105)
beta_threshold_110 <- log(or_threshold_110)

# Dose term after adjustment for uncapped pre-response dose intensity
p_gt_0_linear_weak_di <- mean(
  posterior_dose_linear_weak_di_uncapped > 0
)

p_or_105_linear_weak_di <- mean(
  posterior_dose_linear_weak_di_uncapped >= beta_threshold_105
)

p_or_110_linear_weak_di <- mean(
  posterior_dose_linear_weak_di_uncapped >= beta_threshold_110
)

# Dose-intensity term
p_di_gt_0 <- mean(
  posterior_di_linear_weak_uncapped > 0
)

p_di_or_105 <- mean(
  posterior_di_linear_weak_uncapped >= beta_threshold_105
)

p_di_or_110 <- mean(
  posterior_di_linear_weak_uncapped >= beta_threshold_110
)

Dose term after adjustment for uncapped pre-response dose intensity

P(β_dose > 0) = 99%

P(OR per additional dose ≥ 1.05) = 81%

P(OR per additional dose ≥ 1.1) = 42%

Uncapped pre-response dose-intensity term

P(β_DI > 0) = 81.4%

P(OR for dose-intensity term ≥ 1.05) = 69%

P(OR for dose-intensity term ≥ 1.1) = 53% Step 7: compare beta_dose with and without dose intensity

This is really the key analysis.

delta_beta_dose_di_adjustment <-
  posterior_dose_linear_weak_di_uncapped - posterior_dose_linear_weak

save_files(
  save_object = delta_beta_dose_di_adjustment,
  filename = "Posterior Difference Beta Linear Dose With vs Without Uncapped PreResponse DI",
  subD = "files/Modeling/Posterior Difference Beta Linear Dose With vs Without Uncapped PreResponse DI"
)

dose_di_comparison <- tibble(
  model = c(
    "Weak prior linear model",
    "Weak prior linear model + uncapped pre-response dose intensity"
  ),
  mean_beta_dose = c(
    mean(posterior_dose_linear_weak),
    mean(posterior_dose_linear_weak_di_uncapped)
  ),
  median_beta_dose = c(
    median(posterior_dose_linear_weak),
    median(posterior_dose_linear_weak_di_uncapped)
  ),
  ci_lower_89 = c(
    quantile(posterior_dose_linear_weak, 0.055),
    quantile(posterior_dose_linear_weak_di_uncapped, 0.055)
  ),
  ci_upper_89 = c(
    quantile(posterior_dose_linear_weak, 0.945),
    quantile(posterior_dose_linear_weak_di_uncapped, 0.945)
  ),
  p_gt_0 = c(
    mean(posterior_dose_linear_weak > 0),
    mean(posterior_dose_linear_weak_di_uncapped > 0)
  ),
  p_gt_005 = c(
    mean(posterior_dose_linear_weak > 0.05),
    mean(posterior_dose_linear_weak_di_uncapped > 0.05)
  ),
  p_gt_010 = c(
    mean(posterior_dose_linear_weak > 0.10),
    mean(posterior_dose_linear_weak_di_uncapped > 0.10)
  )
)

dose_di_comparison

# A tibble: 2 × 8
  model  mean_beta_dose median_beta_dose ci_lower_89 ci_upper_89 p_gt_0 p_gt_005
  <chr>           <dbl>            <dbl>       <dbl>       <dbl>  <dbl>    <dbl>
1 Weak …         0.0904           0.0883      0.0234       0.161  0.986    0.817
2 Weak …         0.0896           0.0870      0.0215       0.165  0.987    0.808
# ℹ 1 more variable: p_gt_010 <dbl>

Direct difference summary:

tibble(
  delta_mean = mean(delta_beta_dose_di_adjustment),
  delta_median = median(delta_beta_dose_di_adjustment),
  delta_lower_89 = quantile(delta_beta_dose_di_adjustment, 0.055),
  delta_upper_89 = quantile(delta_beta_dose_di_adjustment, 0.945),
  p_delta_gt_0 = mean(delta_beta_dose_di_adjustment > 0)
)

# A tibble: 1 × 5
  delta_mean delta_median delta_lower_89 delta_upper_89 p_delta_gt_0
       <dbl>        <dbl>          <dbl>          <dbl>        <dbl>
1  -0.000809     -0.00123         -0.101         0.0987        0.492

Density of the difference:

delta_beta_df <- tibble(delta_beta = delta_beta_dose_di_adjustment)

ggplot(delta_beta_df, aes(x = delta_beta)) +
  geom_density(fill = "steelblue", alpha = 0.4) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  labs(
    title = "Difference in Posterior Dose Effect After Adding Dose Intensity",
    subtitle = "β_dose(with uncapped pre-response DI) − β_dose(without DI)",
    x = expression(Delta * beta[dose]),
    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 linear ORR model and comparing the posterior dose-effect estimate before and after adjustment. The posterior difference in the dose coefficient was centered at -0.001, with an 89% credible interval from -0.101 to 0.099. This comparison shows how much inference on the primary dose effect changes after accounting for treatment timing intensity prior to response.

5.6 Posterior predicted response curve from the model including dose intensity

For the dose-response plot, hold di = 0, which corresponds to the average centered uncapped pre-response dose intensity in the cohort.

doses <- 0:14

new_data_linear_di_uncapped <- data.frame(
  dose = doses,
  di = 0,
  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_linear_weak_di_new <- link(
  m_linear_real_weak_di_uncapped,
  data = new_data_linear_di_uncapped
)

mean_p_linear_weak_di <- apply(p_linear_weak_di_new, 2, mean)
ci_linear_weak_di <- apply(p_linear_weak_di_new, 2, PI, 0.89)

plot_df_linear_weak_di <- tibble(
  dose_number = doses + 1,
  mean_p = mean_p_linear_weak_di,
  lower = ci_linear_weak_di[1, ],
  upper = ci_linear_weak_di[2, ]
)

save_files(
  save_object = plot_df_linear_weak_di,
  filename = "Plot DF Predicted Response Linear Weak Prior With Uncapped PreResponse DI",
  subD = "files/Modeling/Plot DF Predicted Response Linear Weak Prior With Uncapped PreResponse DI"
)

plot_predicted_response_linear_weak_di_uncapped <-
  ggplot(plot_df_linear_weak_di, aes(x = dose_number, y = mean_p)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    fill = "steelblue",
    alpha = 0.22
  ) +
  geom_line(
    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_linear_weak_di$dose_number),
                 max(plot_df_linear_weak_di$dose_number),
                 by = 1),
    expand = expansion(mult = c(0.01, 0.02))
  ) +
  labs(
    title = "Predicted Response Probability by Number of Doses",
    subtitle = "Weak prior linear 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_linear_weak_di_uncapped

save_files(
  save_object = plot_predicted_response_linear_weak_di_uncapped,
  filename = "Predicted Response Probability Linear Model Weak Prior With Uncapped PreResponse DI",
  subD = "files/Modeling/Predicted Response Probability Linear Model Weak Prior With Uncapped PreResponse DI"
)

save_files(
  save_object = plot_predicted_response_linear_weak_di_uncapped,
  filename = "Predicted Response Probability Linear Model Weak Prior With Uncapped PreResponse DI",
  subD = "files/Modeling/Predicted Response Probability Linear Model Weak Prior With Uncapped PreResponse DI",
  extension = ".png",
  width = 8,
  height = 6,
  dpi = 300,
  units = "in"
)

Step 9: posterior distribution plot for beta_dose after dose-intensity adjustment

dens_linear_weak_di <- density(posterior_dose_linear_weak_di_uncapped)
dens_df_linear_weak_di <- data.frame(
  x = dens_linear_weak_di$x,
  y = dens_linear_weak_di$y
)

ymax_linear_weak_di <- max(dens_df_linear_weak_di$y)

ci_lower_linear_weak_di <- quantile(posterior_dose_linear_weak_di_uncapped, 0.055)
ci_upper_linear_weak_di <- quantile(posterior_dose_linear_weak_di_uncapped, 0.945)

# OR thresholds expressed on the beta/log-odds scale
or_threshold_105 <- 1.05
or_threshold_110 <- 1.10

beta_threshold_105 <- log(or_threshold_105)
beta_threshold_110 <- log(or_threshold_110)

# Tail probabilities
p_gt_0_linear_weak_di <- mean(posterior_dose_linear_weak_di_uncapped > 0)

p_or_105_linear_weak_di <- mean(
  posterior_dose_linear_weak_di_uncapped >= beta_threshold_105
)

p_or_110_linear_weak_di <- mean(
  posterior_dose_linear_weak_di_uncapped >= beta_threshold_110
)

col_gt_0 <- "#1B9E77"
col_005  <- "#1F78B4"
col_010  <- "#D95F02"
col_line <- "#08306B"

plot_linear_weak_di_uncapped <-
  ggplot() +
  geom_line(
    data = dens_df_linear_weak_di,
    aes(x = x, y = y),
    linewidth = 0.9,
    color = col_line
  ) +
  geom_area(
    data = subset(
      dens_df_linear_weak_di,
      x > 0 & x < beta_threshold_105
    ),
    aes(x = x, y = y),
    fill = col_gt_0,
    alpha = 0.35
  ) +
  geom_area(
    data = subset(
      dens_df_linear_weak_di,
      x >= beta_threshold_105 & x < beta_threshold_110
    ),
    aes(x = x, y = y),
    fill = col_005,
    alpha = 0.45
  ) +
  geom_area(
    data = subset(
      dens_df_linear_weak_di,
      x >= beta_threshold_110
    ),
    aes(x = x, y = y),
    fill = col_010,
    alpha = 0.50
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "red",
    linewidth = 1
  ) +
  geom_vline(
    xintercept = c(ci_lower_linear_weak_di, ci_upper_linear_weak_di),
    linetype = "dotted",
    color = "grey40",
    linewidth = 0.8
  ) +
  annotate(
    "text",
    x = 0,
    y = ymax_linear_weak_di * 0.95,
    label = "Null",
    hjust = -0.2,
    color = "red"
  ) +
  annotate(
    "text",
    x = mean(posterior_dose_linear_weak_di_uncapped),
    y = ymax_linear_weak_di * 0.40,
    label = paste0(
      "89% CrI: [",
      round(ci_lower_linear_weak_di, 2),
      ", ",
      round(ci_upper_linear_weak_di, 2),
      "]"
    ),
    size = 3.4
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_linear_weak_di * 0.95,
    label = paste0(
      "P(β > 0 | ↑ dose → ↑ odds of response) = ",
      round(p_gt_0_linear_weak_di * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_gt_0
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_linear_weak_di * 0.85,
    label = paste0(
      "P(OR per additional dose ≥ ",
      or_threshold_105,
      ") = ",
      round(p_or_105_linear_weak_di * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_005
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_linear_weak_di * 0.75,
    label = paste0(
      "P(OR per additional dose ≥ ",
      or_threshold_110,
      ") = ",
      round(p_or_110_linear_weak_di * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_010
  ) +
  labs(
    title = "Posterior Distribution: Linear Dose Effect",
    subtitle = "Weak prior model adjusted for uncapped pre-response dose intensity",
    x = "β_dose (log-odds per dose)",
    y = "Density"
  ) +
  theme_minimal() +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.28))) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )
plot_linear_weak_di_uncapped

save_files(
  save_object = plot_linear_weak_di_uncapped,
  filename = "Density Plot Linear Model Weak Prior With Uncapped PreResponse DI",
  subD = "files/Modeling/Density Plot Linear Model Weak Prior With Uncapped PreResponse DI"
)

save_files(
  save_object = plot_linear_weak_di_uncapped,
  filename = "Density Plot Linear Model Weak Prior With Uncapped PreResponse DI",
  subD = "files/Modeling/Density Plot Linear Model Weak Prior With Uncapped PreResponse DI",
  extension = ".png",
  width = 7.6,
  height = 7.5,
  dpi = 300,
  units = "in"
)

Step 10: optional companion plot for the dose-intensity coefficient itself

dens_di <- density(posterior_di_linear_weak_uncapped)
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_linear_weak_uncapped, 0.055)
ci_upper_di <- quantile(posterior_di_linear_weak_uncapped, 0.945)

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

Diagnostic

df_analysis |>
  summarise(
    corr_dose_di = cor(
      dose_minus_1,
      dose_intensity_pre_response_uncapped_centered,
      use = "complete.obs"
    )
  )

# A tibble: 1 × 1
  corr_dose_di
         <dbl>
1      -0.0704

ggplot(
  df_analysis,
  aes(x = dose_minus_1, y = dose_intensity_pre_response_uncapped_centered)
) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", se = FALSE, color = "red") +
  labs(
    title = "Relationship Between Dose Count and Pre-Response Uncapped Dose Intensity",
    x = "Dose minus 1",
    y = "Centered uncapped pre-response dose intensity"
  ) +
  theme_minimal()

post_linear_weak_di_uncapped <- extract.samples(m_linear_real_weak_di_uncapped)

posterior_dose_linear_weak_di_uncapped <- post_linear_weak_di_uncapped$beta_dose
posterior_di_linear_weak_uncapped <- post_linear_weak_di_uncapped$beta_di

mean(posterior_dose_linear_weak_di_uncapped > 0)

[1] 0.98675

mean(posterior_dose_linear_weak_di_uncapped > 0.05)

[1] 0.8075

mean(posterior_dose_linear_weak_di_uncapped > 0.10)

[1] 0.38625

mean(posterior_di_linear_weak_uncapped > 0)

[1] 0.814

mean(posterior_di_linear_weak_uncapped > 0.05)

[1] 0.6815

mean(posterior_di_linear_weak_uncapped > 0.10)

[1] 0.51425

delta_beta_dose_di_adjustment <-
  posterior_dose_linear_weak_di_uncapped - posterior_dose_linear_weak

mean(delta_beta_dose_di_adjustment)

[1] -0.0008087979

median(delta_beta_dose_di_adjustment)

[1] -0.001232673

quantile(delta_beta_dose_di_adjustment, c(0.055, 0.945))

      5.5%      94.5% 
-0.1009558  0.0986611 

mean(delta_beta_dose_di_adjustment > 0)

[1] 0.4925

df_analysis |>
  summarise(
    corr = cor(
      dose_minus_1,
      dose_intensity_pre_response_uncapped_centered,
      use = "complete.obs"
    )
  )

# A tibble: 1 × 1
     corr
    <dbl>
1 -0.0704

Interpretation: Adjustment for uncapped pre-response dose intensity did not materially change the estimated dose effect (Δβ ≈ 0; 89% CrI −0.10 to 0.10), indicating that treatment timing intensity did not confound the observed dose–response relationship

DAG Interpretation: patient factors ──► dose count ──► response │ │ │ └──► dose intensity └────────────────────► response

Joint Posterior

df_joint <- tibble(
  beta_dose = posterior_dose_linear_weak_di_uncapped,
  beta_di   = posterior_di_linear_weak_uncapped
)

library(ggplot2)

joint_posterior_plot_filled <-
  ggplot(df_joint, aes(x = beta_dose, 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 Dose Effect and Dose Intensity",
    subtitle = "Weak prior linear model with uncapped pre-response dose intensity",
    x = expression(beta[dose]),
    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)
  )

save_files(
  save_object = joint_posterior_plot_filled,
  filename = "Joint Posterior Dose and DI Filled",
  subD = "files/Modeling/Joint Posterior Dose and DI Filled"
)

save_files(
  save_object = joint_posterior_plot_filled,
  filename = "Joint Posterior Dose and DI Filled",
  subD = "files/Modeling/Joint Posterior Dose and DI Filled",
  extension = ".png",
  width = 8,
  height = 6,
  dpi = 300,
  units = "in"
)

ggplot(df_joint, aes(beta_dose, beta_di)) +

  stat_ellipse(level = 0.89, linewidth = 1.2) +

  geom_point(alpha = 0.1) +

  geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +

  labs(
    title = "Joint Posterior Distribution",
    subtitle = "Dose count vs pre-response dose intensity",
    x = expression(beta[dose]),
    y = expression(beta[DI])
  ) +

  theme_classic()

The joint posterior distribution of the dose effect and pre-response dose intensity showed minimal correlation (r ≈ -0.0703838), indicating that the estimated benefit associated with additional doses was not explained by differences in treatment timing or dose intensity.

6 ——————————————————————————

7 Sensitivity analysis: linear ORR model using pre-response dose count

8 ——————————————————————————

df_analysis_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"),
    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(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)
  )

df_analysis_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),
    n_total_gt_pre_response = sum(num_doses > dose_pre_response, na.rm = TRUE),
    n_total_lt_pre_response = sum(num_doses < dose_pre_response, na.rm = TRUE)
  )

# A tibble: 1 × 6
      n observed_orr min_pre_response_doses max_pre_response_doses
  <int>        <dbl>                  <int>                  <int>
1   189        0.646                      1                     20
# ℹ 2 more variables: n_total_gt_pre_response <int>,
#   n_total_lt_pre_response <int>

dat_linear_pre_response <- list(
  N = nrow(df_analysis_pre_response),
  y = df_analysis_pre_response$response,

  # Key sensitivity exposure:
  # pre-response dose count minus 1
  dose = df_analysis_pre_response$dose_pre_response_minus_1,

  age = df_analysis_pre_response$age_centered,
  imm = as.integer(df_analysis_pre_response$immunosuppressed) - 1,
  ecog23 = df_analysis_pre_response$ecog23,

  loc_conj = df_analysis_pre_response$loc_conj,
  loc_other = df_analysis_pre_response$loc_other,
  loc_unknown = df_analysis_pre_response$loc_unknown,

  stage_III = df_analysis_pre_response$stage_III,
  stage_IV = df_analysis_pre_response$stage_IV,
  stage_unstaged = df_analysis_pre_response$stage_unstaged,

  agent_pembro_q3 = df_analysis_pre_response$agent_pembro_q3,
  agent_pembro_q6 = df_analysis_pre_response$agent_pembro_q6,
  agent_pembro_mixed = df_analysis_pre_response$agent_pembro_mixed,
  agent_ipinivo = df_analysis_pre_response$agent_ipinivo,

  alpha_mu = qlogis(0.40)
)

m_linear_real_weak_pre_response <- ulam(
  alist(
    y ~ bernoulli(p),

    logit(p) <- alpha +
      beta_dose * dose +
      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_dose ~ 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_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_linear_pre_response,
  chains = 4,
  cores = 4,
  iter = 2000,
  cmdstan = TRUE
)

m_linear_real_weak_pre_response@cstanfit$save_output_files(
  dir = here::here(
    "files",
    "Modeling",
    "cmdstan_output",
    "model_linear_weak_pre_response"
  )
)

saveRDS(
  m_linear_real_weak_pre_response,
  here::here(
    "files",
    "Modeling",
    "Model Linear Weak Prior PreResponse Dose",
    paste0(
      "Model Linear Weak Prior PreResponse Dose-",
      format(Sys.time(), "%Y-%m-%d-%H-%M-%S"),
      ".rds"
    )
  )
) 

8.0.1 Diagnostics

precis(m_linear_real_weak_pre_response, depth = 2)

                          mean         sd        5.5%       94.5%      rhat
alpha              0.736517851 0.28409190  0.29031067  1.20149600 1.0017937
beta_dose          0.115448764 0.05142088  0.03587391  0.19846506 1.0006780
beta_age           0.001625089 0.01215133 -0.01773266  0.02125541 0.9997816
beta_imm          -0.910723593 0.17974701 -1.20270670 -0.62274844 1.0007817
beta_ecog         -0.305390814 0.26281611 -0.71836259  0.11726613 1.0018796
beta_conj         -2.274927247 0.42243317 -2.96509143 -1.60484889 1.0000842
beta_other        -0.148546633 0.27007659 -0.59591095  0.28648773 1.0002946
beta_unknown       0.045265131 0.29214236 -0.42082789  0.50949780 1.0030077
beta_stageIII      0.001036432 0.24040404 -0.37680068  0.39038181 0.9998875
beta_stageIV      -0.315720681 0.23382095 -0.68217216  0.07253556 1.0012877
beta_stageUn       0.239385714 0.26521538 -0.18234210  0.66511197 1.0001085
beta_pembro_q3    -0.043951546 0.23317515 -0.42265855  0.32711945 1.0016095
beta_pembro_q6    -0.059934003 0.28341177 -0.51586760  0.38943500 1.0013596
beta_pembro_mixed  0.218465428 0.27436564 -0.21876770  0.66072113 1.0003667
beta_ipinivo       0.418728888 0.28393111 -0.03749181  0.86222850 1.0004804
                  ess_bulk
alpha             3960.510
beta_dose         5809.132
beta_age          6678.511
beta_imm          6676.527
beta_ecog         6583.509
beta_conj         6455.864
beta_other        7011.504
beta_unknown      7045.331
beta_stageIII     4707.972
beta_stageIV      4643.728
beta_stageUn      6064.495
beta_pembro_q3    6663.377
beta_pembro_q6    7511.292
beta_pembro_mixed 6704.574
beta_ipinivo      6671.517

plot(precis(m_linear_real_weak_pre_response))

# ------------------------------------------------------------------------------
# Manual posterior predictive check from fitted probabilities
# ------------------------------------------------------------------------------

set.seed(123)

p_hat_pre_response <- link(m_linear_real_weak_pre_response)

stopifnot(ncol(p_hat_pre_response) == nrow(df_analysis_pre_response))

# Use 1,000 posterior draws for the PPC
n_ppc_draws <- 1000

draw_ids <- sample(
  seq_len(nrow(p_hat_pre_response)),
  size = n_ppc_draws,
  replace = FALSE
)

p_hat_ppc <- p_hat_pre_response[draw_ids, , drop = FALSE]

# Simulate Bernoulli outcomes from the fitted probabilities
yrep_pre_response <- 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_pre_response <- mean(df_analysis_pre_response$response)

ppc_pre_response_manual <- tibble(
  draw = seq_len(nrow(yrep_pre_response)),
  simulated_orr = rowMeans(yrep_pre_response)
)

ppc_pre_response_manual |>
  summarise(
    observed_orr = observed_orr_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_pre_response),
    mean_fitted_probability = mean(p_hat_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.641                0.640    0.571    0.714
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>

ggplot(ppc_pre_response_manual, aes(x = simulated_orr)) +
  geom_histogram(
    bins = 50,
    fill = "steelblue",
    alpha = 0.7,
    color = "white"
  ) +
  geom_vline(
    xintercept = observed_orr_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 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")
  )

post_linear_weak_pre_response <- extract.samples(m_linear_real_weak_pre_response)

posterior_dose_linear_weak_pre_response <-
  post_linear_weak_pre_response$beta_dose

save_files(
  save_object = posterior_dose_linear_weak_pre_response,
  filename = "Posterior Dose Linear Weak Prior PreResponse Dose",
  subD = "files/Modeling/Posterior Dose Linear 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 Dose Linear Weak Prior PreResponse Dose/Posterior Dose Linear Weak Prior PreResponse Dose-2026-05-08-04-48-17.rds 

p_gt_0_pre_response   <- mean(posterior_dose_linear_weak_pre_response > 0)
p_gt_005_pre_response <- mean(posterior_dose_linear_weak_pre_response > 0.05)
p_gt_010_pre_response <- mean(posterior_dose_linear_weak_pre_response > 0.10)

tibble(
  parameter = "beta_dose_pre_response",
  mean_beta = mean(posterior_dose_linear_weak_pre_response),
  median_beta = median(posterior_dose_linear_weak_pre_response),
  ci_lower_89 = quantile(posterior_dose_linear_weak_pre_response, 0.055),
  ci_upper_89 = quantile(posterior_dose_linear_weak_pre_response, 0.945),
  p_gt_0 = p_gt_0_pre_response,
  p_gt_005 = p_gt_005_pre_response,
  p_gt_010 = p_gt_010_pre_response,
  median_or_per_additional_pre_response_dose =
    exp(median(posterior_dose_linear_weak_pre_response))
)

# A tibble: 1 × 9
  parameter        mean_beta median_beta ci_lower_89 ci_upper_89 p_gt_0 p_gt_005
  <chr>                <dbl>       <dbl>       <dbl>       <dbl>  <dbl>    <dbl>
1 beta_dose_pre_r…     0.115       0.114      0.0359       0.198  0.989    0.903
# ℹ 2 more variables: p_gt_010 <dbl>,
#   median_or_per_additional_pre_response_dose <dbl>

dens_pre_response <- density(posterior_dose_linear_weak_pre_response)
dens_df_pre_response <- data.frame(
  x = dens_pre_response$x,
  y = dens_pre_response$y
)

ymax_pre_response <- max(dens_df_pre_response$y)

ci_lower_pre_response <- quantile(
  posterior_dose_linear_weak_pre_response,
  0.055
)

ci_upper_pre_response <- quantile(
  posterior_dose_linear_weak_pre_response,
  0.945
)

col_gt_0 <- "#1B9E77"
col_005  <- "#1F78B4"
col_010  <- "#D95F02"
col_line <- "#08306B"

plot_linear_weak_pre_response <-
  ggplot() +
  geom_line(
    data = dens_df_pre_response,
    aes(x = x, y = y),
    linewidth = 0.9,
    color = col_line
  ) +
  geom_area(
    data = subset(dens_df_pre_response, x > 0 & x < 0.05),
    aes(x = x, y = y),
    fill = col_gt_0,
    alpha = 0.35
  ) +
  geom_area(
    data = subset(dens_df_pre_response, x >= 0.05 & x < 0.10),
    aes(x = x, y = y),
    fill = col_005,
    alpha = 0.45
  ) +
  geom_area(
    data = subset(dens_df_pre_response, x >= 0.10),
    aes(x = x, y = y),
    fill = col_010,
    alpha = 0.50
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "red",
    linewidth = 1
  ) +
  geom_vline(
    xintercept = c(ci_lower_pre_response, ci_upper_pre_response),
    linetype = "dotted",
    color = "grey40",
    linewidth = 0.8
  ) +
  annotate(
    "text",
    x = 0,
    y = ymax_pre_response * 0.95,
    label = "Null",
    hjust = -0.2,
    color = "red"
  ) +
  annotate(
    "text",
    x = mean(posterior_dose_linear_weak_pre_response),
    y = ymax_pre_response * 0.40,
    label = paste0(
      "89% CrI: [",
      round(ci_lower_pre_response, 2),
      ", ",
      round(ci_upper_pre_response, 2),
      "]"
    ),
    size = 3.4
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_pre_response * 0.95,
    label = paste0(
      "P(β > 0 | ↑ pre-response dose → ↑ odds of response) = ",
      round(p_gt_0_pre_response * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_gt_0
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_pre_response * 0.85,
    label = paste0(
      "P(β > 0.05 | ≥5% ↑ odds per additional pre-response dose) = ",
      round(p_gt_005_pre_response * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_005
  ) +
  annotate(
    "text",
    x = 0.17,
    y = ymax_pre_response * 0.75,
    label = paste0(
      "P(β > 0.10 | ≥11% ↑ odds per additional pre-response dose) = ",
      round(p_gt_010_pre_response * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_010
  ) +
  labs(
    title = "Posterior Distribution: Pre-Response Linear Dose Effect",
    subtitle = "Weakly informative prior",
    x = "β_dose (log-odds per additional pre-response dose)",
    y = "Density"
  ) +
  theme_minimal() +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.28))) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

plot_linear_weak_pre_response

save_files(
  save_object = plot_linear_weak_pre_response,
  filename = "Density Plot Linear Model Weak Prior PreResponse Dose",
  subD = "files/Modeling/Density Plot Linear Model Weak Prior PreResponse Dose"
)

save_files(
  save_object = plot_linear_weak_pre_response,
  filename = "Density Plot Linear Model Weak Prior PreResponse Dose",
  subD = "files/Modeling/Density Plot Linear Model Weak Prior PreResponse Dose",
  extension = ".png",
  width = 7.6,
  height = 7.5,
  dpi = 300,
  units = "in"
)

max_pre_response_dose_for_plot <- min(
  15,
  max(df_analysis_pre_response$dose_pre_response, na.rm = TRUE)
)

doses_pre_response <- 0:(max_pre_response_dose_for_plot - 1)

new_data_linear_pre_response <- data.frame(
  dose = doses_pre_response,
  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_linear_weak_pre_response_new <- link(
  m_linear_real_weak_pre_response,
  data = new_data_linear_pre_response
)

mean_p_linear_weak_pre_response <- apply(
  p_linear_weak_pre_response_new,
  2,
  mean
)

ci_linear_weak_pre_response <- apply(
  p_linear_weak_pre_response_new,
  2,
  PI,
  0.89
)

plot_df_linear_weak_pre_response <- tibble(
  dose_number = doses_pre_response + 1,
  mean_p = mean_p_linear_weak_pre_response,
  lower = ci_linear_weak_pre_response[1, ],
  upper = ci_linear_weak_pre_response[2, ]
)

plot_predicted_response_linear_weak_pre_response <-
  ggplot(plot_df_linear_weak_pre_response, aes(x = dose_number, y = mean_p)) +
  geom_ribbon(
    aes(ymin = lower, ymax = upper),
    fill = "steelblue",
    alpha = 0.22
  ) +
  geom_line(
    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_linear_weak_pre_response$dose_number),
      max(plot_df_linear_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 linear 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_linear_weak_pre_response

save_files(
  save_object = plot_df_linear_weak_pre_response,
  filename = "Plot DF Predicted Response Linear Weak Prior PreResponse Dose",
  subD = "files/Modeling/Plot DF Predicted Response Linear Weak Prior PreResponse Dose"
)

save_files(
  save_object = plot_predicted_response_linear_weak_pre_response,
  filename = "Predicted Response Probability Linear Model Weak Prior PreResponse Dose",
  subD = "files/Modeling/Predicted Response Probability Linear Model Weak Prior PreResponse Dose"
)

save_files(
  save_object = plot_predicted_response_linear_weak_pre_response,
  filename = "Predicted Response Probability Linear Model Weak Prior PreResponse Dose",
  subD = "files/Modeling/Predicted Response Probability Linear Model Weak Prior PreResponse Dose",
  extension = ".png",
  width = 8,
  height = 6,
  dpi = 300,
  units = "in"
)

linear_total_vs_pre_response_comparison <- tibble(
  model = c(
    "Total cumulative dose count",
    "Pre-response dose count"
  ),
  mean_beta_dose = c(
    mean(posterior_dose_linear_weak),
    mean(posterior_dose_linear_weak_pre_response)
  ),
  median_beta_dose = c(
    median(posterior_dose_linear_weak),
    median(posterior_dose_linear_weak_pre_response)
  ),
  ci_lower_89 = c(
    quantile(posterior_dose_linear_weak, 0.055),
    quantile(posterior_dose_linear_weak_pre_response, 0.055)
  ),
  ci_upper_89 = c(
    quantile(posterior_dose_linear_weak, 0.945),
    quantile(posterior_dose_linear_weak_pre_response, 0.945)
  ),
  p_gt_0 = c(
    mean(posterior_dose_linear_weak > 0),
    mean(posterior_dose_linear_weak_pre_response > 0)
  ),
  p_gt_005 = c(
    mean(posterior_dose_linear_weak > 0.05),
    mean(posterior_dose_linear_weak_pre_response > 0.05)
  ),
  p_gt_010 = c(
    mean(posterior_dose_linear_weak > 0.10),
    mean(posterior_dose_linear_weak_pre_response > 0.10)
  )
)

linear_total_vs_pre_response_comparison

# A tibble: 2 × 8
  model  mean_beta_dose median_beta_dose ci_lower_89 ci_upper_89 p_gt_0 p_gt_005
  <chr>           <dbl>            <dbl>       <dbl>       <dbl>  <dbl>    <dbl>
1 Total…         0.0904           0.0883      0.0234       0.161  0.986    0.817
2 Pre-r…         0.115            0.114       0.0359       0.198  0.989    0.903
# ℹ 1 more variable: p_gt_010 <dbl>

df_analysis_pre_response |>
  summarise(
    n = n(),
    n_total_differs_from_pre_response = sum(num_doses != 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),
    correlation_total_pre_response = cor(num_doses, dose_pre_response, use = "complete.obs")
  )

# A tibble: 1 × 5
      n n_total_differs_from_pre_respo…¹ mean_total_doses mean_pre_response_do…²
  <int>                            <int>            <dbl>                  <dbl>
1   189                               31             4.62                   3.88
# ℹ abbreviated names: ¹​n_total_differs_from_pre_response,
#   ²​mean_pre_response_doses
# ℹ 1 more variable: correlation_total_pre_response <dbl>

df_analysis_pre_response |>
  filter(num_doses != dose_pre_response) |>
  count(response, num_doses, dose_pre_response)

# A tibble: 27 × 4
   response num_doses dose_pre_response     n
      <dbl>     <int>             <int> <int>
 1        0         3                 2     3
 2        0         4                 1     1
 3        0         4                 3     1
 4        0         5                 1     1
 5        0         5                 3     1
 6        0         5                 4     2
 7        0         6                 2     1
 8        0         7                 3     1
 9        0         7                 4     1
10        0        16                 6     1
# ℹ 17 more rows

9 Conclusion

Across multiple Bayesian ORR model specifications, the central question was whether greater treatment exposure was associated with a higher probability of response in the observed cohort. To address that question rigorously, we evaluated the dose–response relationship under alternative functional forms and then examined whether the inference changed after incorporating treatment-timing variables such as pre-response dose intensity.

This workflow was designed to separate three related issues: the direction of the dose–response association, the sensitivity of that association to modeling assumptions, and the extent to which treatment timing might account for the observed pattern. Presenting the full modeling pathway—including posterior predictive checks, prior sensitivity analyses, and dose-intensity sensitivity analyses—helps clarify what aspects of the inference are stable and which remain uncertain.

As a companion to the manuscript, this page is meant to make the ORR analysis more transparent, reproducible, and interpretable for readers who want to see how the model-based conclusions were constructed from the observed data.