Bayesian ORR Modeling — Log Dose Specification

1 Overview

This page presents the Bayesian objective response rate (ORR) analysis for the log-dose specification. Within the broader dose–response modeling framework, this formulation asks whether increasing treatment exposure is associated with improved response while allowing for diminishing returns as the number of doses increases.

All models are fit in R using rethinking::ulam(). Before being applied to the real cohort, the priors and model structure were developed in a separate generative workflow using simulated data. That preliminary work was used to examine prior behavior, test recovery of known effects, and evaluate whether a log-dose formulation could represent nonlinear exposure–response patterns in a stable and interpretable way. The analyses shown here therefore reflect the application of a pre-specified and previously stress-tested model to the observed clinical dataset.

The primary adjusted model estimates the association between log-transformed dose exposure and response probability while accounting for key clinical covariates, including age, immunosuppression, tumor location, stage, treatment agent, and ECOG performance status.

The rationale for this specification is straightforward: if treatment benefit accumulates rapidly early and then tapers, a simple linear effect of dose count may be too rigid. A log-dose model provides a more flexible alternative by allowing early increases in exposure to have larger modeled effects than later increases.

This page also includes prior-sensitivity analyses and dose-intensity sensitivity analyses based on pre-response treatment timing. These complementary analyses help assess whether the estimated log-dose association remains stable across reasonable modeling choices and alternative representations of treatment delivery.

Taken together, this document is intended to complement the main manuscript by making the log-dose ORR workflow transparent, reproducible, and easier to interpret.

NoteNote on numerical values

This document is intended to show the modeling workflow, diagnostics, and sensitivity analyses supporting the manuscript. Because Bayesian models were fit using MCMC, exact posterior summaries may vary slightly across model runs. In addition, some HTML modeling documents may have been rendered at a different point in the workflow than the final manuscript.

Accordingly, small differences in posterior probabilities, credible intervals, or figure annotations may occur between this document and the manuscript. The final manuscript should be treated as the source of record for reported numerical results.

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_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 = factor(
      agent_schedule,
      levels = c("cemiplimab", "pembro_q3", "pembro_q6", "pembro_mixed", "Ipi-nivo")
    ),

    # -----------------------------
    # Surgery
    # -----------------------------
    surgery = as.integer(has_surg)
  ) |>
  filter(!is.na(num_doses), num_doses >= 1)

df_analysis <- df_analysis |>
  select(
    record_id,
    response,
    num_doses,
    log_dose,
    dose_group,
    dose_2,
    dose_3,
    dose_4,
    dose_5plus,
    age_centered,
    sex,
    immunosuppressed,
    location_condensed,
    stage_condensed,
    agent,
    agent_schedule,
    ecog_condensed,
    expected_interval_days,
    likely_schedule,
    num_doses_recorded,
    num_doses_pre_response,
    dose_intensity_final,
    dose_intensity_pre_response,
    dose_intensity_pre_response_uncapped,
    mean_interval_days,
    interval_deviation_days,
    extra_doses_after_response,
    dose_intensity_source,
    surgery,
    response_cutoff_date
  )

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.1 Inspect Modeling Dataset

colnames(df_analysis)

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

vars_main <- c(
  "response",
  "num_doses",
  "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: 10
$ 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…
$ 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…

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

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

dat_log <- list(

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

  dose = df_analysis$log_dose,
  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 prior for dose-response relationship

We first fit the log-dose ORR model under a skeptical prior for the treatment-exposure effect. In this specification, the coefficient for log-dose is centered at 0, reflecting prior uncertainty about whether greater treatment exposure meaningfully improves response probability. This prior does not rule out benefit, but it shrinks the model away from large positive dose effects unless the observed data provide support.

The outcome is binary (response vs non-response), so the model is fit using Bayesian logistic regression. The primary exposure term is log_dose, which allows the effect of additional treatment to diminish as dose accumulates. The model also adjusts for prespecified clinical covariates, including age, immunosuppression, ECOG performance status, tumor location, stage, and treatment agent.

This skeptical-prior fit serves as a conservative reference model. We then compare it with the corresponding weakly informative prior specification to assess the robustness of the inferred dose-response relationship.

#| label: fit-log-dose-model-skeptical
#| code-fold: false
#| results: hide

# Bayesian logistic regression for objective response (y = 0/1)
# Primary exposure: log-transformed dose count
# Prior setting: skeptical prior on beta_dose, centered at 0

m_log_real_skeptical <- ulam(
  alist(
    
    # Likelihood: binary response outcome
    y ~ bernoulli(p),

    # Linear predictor on the log-odds scale
    logit(p) <- alpha +
      beta_dose * dose +              # log-dose exposure effect
      beta_age * age +                # centered age
      beta_imm * imm +                # immunosuppression
      beta_ecog * ecog23 +            # ECOG 2-3 vs 0-1
      beta_conj * loc_conj +          # conjunctival primary
      beta_other * loc_other +        # other primary location
      beta_unknown * loc_unknown +    # unknown primary
      beta_stageIII * stage_III +     # stage III vs stage I-II
      beta_stageIV * stage_IV +       # stage IV vs stage I-II
      beta_stageUn * stage_unstaged + # unstaged vs stage I-II
      beta_pembro_q3 * agent_pembro_q3 +
      beta_pembro_q6 * agent_pembro_q6 +
      beta_pembro_mixed * agent_pembro_mixed +
      beta_ipinivo * agent_ipinivo,

    # Intercept prior:
    # centered on alpha_mu = qlogis(0.40), corresponding to a prior
    # baseline response probability of about 40% for the reference patient
    alpha ~ normal(alpha_mu, 1.5),

    # Skeptical prior on dose effect:
    # centered at 0, allowing benefit or harm but shrinking away
    # from large effects unless supported by the data
    beta_dose ~ normal(0, 0.2),

    # Priors for remaining covariates
    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_log,
  chains = 4,
  cores = 4,
  iter = 2000
)

3.1 Posterior predictive check: overall response rate

We first assessed model fit at the aggregate level by comparing the observed objective response rate (ORR) with the distribution of ORRs generated from posterior predictive simulations. This check asks whether the fitted model can reproduce the overall frequency of response observed in the cohort.

As an initial model-checking criterion, the observed ORR should fall within the central range of the simulated ORR distribution. In this analysis, the observed response rate was close to the center of the posterior predictive distribution, indicating that the model captured the aggregate response frequency reasonably well. This check supports aggregate calibration, although it does not by itself establish that the model fully captures patient-level heterogeneity or the true functional form of the dose–response relationship.

manual_binary_ppc <- function(fit, observed_y, n_ppc_draws = 1000, seed = 123) {
  
  # Set a seed so that the posterior predictive simulations are reproducible.
  # This does not change the fitted model; it only makes the random simulated
  # outcomes below repeatable.
  set.seed(seed)

  # link() returns the posterior expected response probabilities from the model.
  #
  # For a binary logistic model, these are probabilities on the response scale.
  # The resulting object is a matrix:
  #   rows    = posterior draws
  #   columns = patients/observations
  #
  # So p_hat[d, i] is the model-implied probability of response for patient i
  # under posterior draw d.
  p_hat <- link(fit)

  # Check that the number of prediction columns matches the number of observed
  # patient outcomes. This helps catch accidental dataset mismatch.
  stopifnot(ncol(p_hat) == length(observed_y))

  # Select a subset of posterior draws to use for the posterior predictive check.
  #
  # If the model has more than n_ppc_draws posterior draws, we randomly sample
  # n_ppc_draws of them. If it has fewer, we use all available draws.
  draw_ids <- sample(
    seq_len(nrow(p_hat)),
    size = min(n_ppc_draws, nrow(p_hat)),
    replace = FALSE
  )

  # Keep only the selected posterior draws.
  #
  # drop = FALSE ensures the result remains a matrix even if only one draw were
  # selected.
  p_hat_ppc <- p_hat[draw_ids, , drop = FALSE]

  # Generate posterior predictive replicated outcomes.
  #
  # For each fitted probability in p_hat_ppc, simulate a binary response:
  #   1 = response
  #   0 = non-response
  #
  # This creates a matrix yrep with:
  #   rows    = simulated posterior predictive datasets
  #   columns = patients/observations
  #
  # So yrep[d, i] is the simulated outcome for patient i under posterior draw d.
  yrep <- matrix(
    rbinom(
      n = length(p_hat_ppc),
      size = 1,
      prob = as.vector(p_hat_ppc)
    ),
    nrow = nrow(p_hat_ppc),
    ncol = ncol(p_hat_ppc)
  )

  # Return the key posterior predictive objects and summaries.
  list(
    
    # Posterior fitted probabilities for all posterior draws and patients.
    # These are probabilities, not simulated outcomes.
    p_hat = p_hat,
    
    # Simulated binary outcomes generated from the fitted probabilities.
    # These are the posterior predictive replicated datasets.
    yrep = yrep,
    
    # Observed objective response rate in the real data.
    observed_orr = mean(observed_y),
    
    # Simulated objective response rate for each posterior predictive dataset.
    # Each row of yrep is one simulated dataset, so rowMeans(yrep) gives the ORR
    # in each simulated dataset.
    simulated_orr = rowMeans(yrep),
    
    # A compact summary comparing the observed ORR with the simulated ORR
    # distribution.
    summary = tibble(
      
      # Observed ORR in the real cohort.
      observed_orr = mean(observed_y),
      
      # Mean simulated ORR across posterior predictive datasets.
      mean_simulated_orr = mean(rowMeans(yrep)),
      
      # Median simulated ORR across posterior predictive datasets.
      median_simulated_orr = median(rowMeans(yrep)),
      
      # 89% posterior predictive interval for the simulated cohort ORR.
      lower_89 = quantile(rowMeans(yrep), 0.055),
      upper_89 = quantile(rowMeans(yrep), 0.945),
      
      # Posterior predictive tail probability:
      # proportion of simulated ORRs less than or equal to the observed ORR.
      #
      # If this is very close to 0 or 1, the observed ORR is in the tail of what
      # the model tends to generate. Values closer to the middle suggest better
      # aggregate calibration.
      p_sim_le_observed = mean(rowMeans(yrep) <= mean(observed_y)),
      
      # Average fitted response probability across all posterior draws and
      # patients. This is a model-based expected ORR, before simulating binary
      # outcomes.
      mean_fitted_probability = mean(p_hat)
    )
  )
}

ggplot(
  tibble(simulated_orr = ppc_log_skeptical$simulated_orr),
  aes(x = simulated_orr)
) +
  geom_histogram(
    bins = 50,
    fill = "steelblue",
    alpha = 0.7,
    color = "white"
  ) +
  geom_vline(
    xintercept = ppc_log_skeptical$observed_orr,
    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: Skeptical Prior Log-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", size = 14),
    plot.subtitle = element_text(hjust = 0.5, face = "bold", size = 12)
  )

Agreement at the aggregate level does not guarantee that the model captures the distribution of outcomes across individual patients. We therefore also examined a rank-based posterior predictive diagnostic.

3.2 Posterior predictive rank check

For each patient, we compared the observed outcome to the corresponding simulated outcomes generated from the fitted model and recorded the number of simulated outcomes that were lower than the observed value.

Because the outcome is binary, the rank distribution is not expected to be uniform: non-response observations necessarily have rank 0, whereas response observations tend to have higher ranks. In this setting, the rank histogram is best interpreted as a qualitative check on how observed outcomes are positioned relative to the model’s predictive distribution, rather than as a strict uniform-calibration target.

ggplot(rank_df_log_skeptical, aes(x = pred_rank_scaled)) +
  geom_histogram(
    binwidth = 0.025,
    boundary = 0,
    closed = "left",
    fill = "grey80",
    color = "white"
  ) +
  coord_cartesian(xlim = c(0, 1)) +
  scale_x_continuous(
    breaks = seq(0, 1, by = 0.25),
    labels = scales::percent_format(accuracy = 1),
    expand = expansion(mult = c(0.02, 0.02))
  ) +
  labs(
    title = "Posterior Predictive Rank Distribution",
    subtitle = "Skeptical-prior log-dose model; binary outcomes produce expected boundary mass",
    x = "Scaled rank of observed outcome within simulated outcomes",
    y = "Number of patients"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5)
  )

rank_df_log_skeptical |>
  count(observed_response)

# A tibble: 2 × 2
  observed_response     n
              <dbl> <int>
1                 0    67
2                 1   122

rank_df_log_skeptical |>
  summarise(
    n = n(),
    n_observed_nonresponse = sum(observed_response == 0, na.rm = TRUE),
    n_rank_zero = sum(pred_rank == 0, na.rm = TRUE),
    prop_rank_zero = mean(pred_rank == 0, na.rm = TRUE),
    min_rank_scaled = min(pred_rank_scaled, na.rm = TRUE),
    max_rank_scaled = max(pred_rank_scaled, na.rm = TRUE)
  )

# A tibble: 1 × 6
      n n_observed_nonresponse n_rank_zero prop_rank_zero min_rank_scaled
  <int>                  <int>       <int>          <dbl>           <dbl>
1   189                     67          67          0.354               0
# ℹ 1 more variable: max_rank_scaled <dbl>

rank_df_log_skeptical |>
  filter(observed_response == 0) |>
  summarise(
    n = n(),
    unique_ranks = paste(unique(pred_rank_scaled), collapse = ", ")
  )

# A tibble: 1 × 2
      n unique_ranks
  <int> <chr>       
1    67 0           

The rank-based posterior predictive check showed the expected boundary mass at zero for observed non-responders, with observed responders distributed primarily in the lower-to-middle rank range. This pattern is expected for a high-response binary outcome and suggests no obvious gross misfit, although rank diagnostics for binary outcomes should be interpreted qualitatively rather than as uniform PIT histograms.

The key point: this is reassuring as a qualitative posterior predictive check, but not a proof that the log-dose shape is correct. It supports the idea that the model generates patient-level outcomes that are broadly compatible with the observed data structure.

3.3 Posterior parameter interpretation

## Model diagnostics

precis(m_log_real_skeptical, depth = 2)

                           mean         sd        5.5%       94.5%      rhat
alpha              0.8349073919 0.30157061  0.35448551  1.31606326 1.0000726
beta_dose          0.1767457242 0.15075104 -0.06522615  0.41904150 1.0002614
beta_age           0.0009980339 0.01176256 -0.01784645  0.02002435 1.0014856
beta_imm          -0.9073080598 0.17638657 -1.19155177 -0.62687320 1.0006232
beta_ecog         -0.3141110678 0.25763058 -0.72864820  0.10254590 0.9995992
beta_conj         -2.3064078068 0.43809139 -3.02184484 -1.59751827 1.0019747
beta_other        -0.1341862150 0.27577399 -0.57064587  0.30278798 1.0033913
beta_unknown       0.0555423473 0.28568123 -0.41034888  0.51422467 1.0001446
beta_stageIII     -0.0157902964 0.23617048 -0.38822690  0.35784485 0.9999771
beta_stageIV      -0.3200314720 0.23614955 -0.69873962  0.04650585 1.0005052
beta_stageUn       0.2578080836 0.27182215 -0.17135256  0.69854882 1.0009373
beta_pembro_q3    -0.0547246016 0.22999082 -0.43073655  0.31678220 1.0015219
beta_pembro_q6    -0.0524759140 0.28995486 -0.52515552  0.40283630 1.0006141
beta_pembro_mixed  0.2572041343 0.28366342 -0.21537811  0.70990471 1.0012748
beta_ipinivo       0.4209911627 0.30167392 -0.05811905  0.90305766 1.0033660
                  ess_bulk
alpha             4067.048
beta_dose         5848.324
beta_age          7563.725
beta_imm          7373.261
beta_ecog         6462.437
beta_conj         8485.499
beta_other        7597.898
beta_unknown      8740.159
beta_stageIII     5507.411
beta_stageIV      5176.039
beta_stageUn      7766.676
beta_pembro_q3    6980.767
beta_pembro_q6    7807.803
beta_pembro_mixed 7811.164
beta_ipinivo      6569.571

plot(precis(m_log_real_skeptical))

3.4 Extract paramaters from the model

post <- extract.samples(m_log_real_skeptical)

3.4.1 Intercept

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

# Summarize
mean(baseline_prob)

[1] 0.6937108

quantile(baseline_prob, c(0.055, 0.945))

     5.5%     94.5% 
0.5877049 0.7885260 

3.4.2 Posterior predicted dose–response curve

dose_number <- 1:15

new_data <- data.frame(
  dose = log(dose_number),
  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 <- link(m_log_real_skeptical, data = new_data)


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

plot_df_log_skeptical <- tibble(
  dose_number = dose_number,
  mean_p = mean_p,
  lower = ci[1, ],
  upper = ci[2, ]
)

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

plot_predicted_response_log_skeptical <-
  ggplot(plot_df_log_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_log_skeptical$dose_number),
    max(plot_df_log_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\nLog-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_log_skeptical

post_log_skeptical <- extract.samples(m_log_real_skeptical)

posterior_dose_log_skeptical <- post_log_skeptical$beta_dose


# OR thresholds for dose doubling
or_doubling_115 <- 1.15
or_doubling_130 <- 1.30

beta_doubling_115 <- log(or_doubling_115) / log(2)
beta_doubling_130 <- log(or_doubling_130) / log(2)

p_gt_0 <- mean(posterior_dose_log_skeptical > 0)

p_or_115_doubling <- mean(
  posterior_dose_log_skeptical >= beta_doubling_115
)

p_or_130_doubling <- mean(
  posterior_dose_log_skeptical >= beta_doubling_130
)

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

P(OR per dose doubling ≥ 1.15) = 43% → Good evidence for at least a modest effect

P(OR per dose doubling ≥ 1.30) = 9% → More limited evidence for a larger effect

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

# 89% CrI
ci_lower <- quantile(posterior_dose_log_skeptical, 0.055)
ci_upper <- quantile(posterior_dose_log_skeptical, 0.945)

# --- Doubling thresholds ---
ln2 <- log(2)
beta_15 <- log(1.15) / ln2   # ~0.2016
beta_30 <- log(1.30) / ln2   # ~0.3785

# Tail probabilities
p_gt_0  <- mean(posterior_dose_log_skeptical > 0)
p_gt_15 <- mean(posterior_dose_log_skeptical >= beta_15)
p_gt_30 <- mean(posterior_dose_log_skeptical >= beta_30)

# Colors
col_gt_0  <- "#1B9E77"
col_15    <- "#1F78B4"
col_30    <- "#D95F02"
col_line  <- "#08306B"
alpha_fill <- 0.45

plot_log_skeptical <-
  ggplot() +

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

  # Non-overlapping shaded regions
  geom_area(
    data = subset(dens_df, x > 0 & x < beta_15),
    aes(x = x, y = y),
    fill = col_gt_0,
    alpha = alpha_fill
  ) +
  geom_area(
    data = subset(dens_df, x >= beta_15 & x < beta_30),
    aes(x = x, y = y),
    fill = col_15,
    alpha = alpha_fill
  ) +
  geom_area(
    data = subset(dens_df, x >= beta_30),
    aes(x = x, y = y),
    fill = col_30,
    alpha = alpha_fill
  ) +

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

  # Top labels
  annotate(
    "text",
    x = 0,
    y = ymax * 1.2,
    label = "Null",
    hjust = -0.2,
    vjust = 1,
    color = "red",
    size = 3
  ) +
  annotate(
    "text",
    x = beta_15,
    y = ymax * 1.2,
    label = "≥15% ↑ odds per doubling\nof dose (e.g., 1→2 or 2→4)",
    hjust = -0.05,
    vjust = 1,
    color = col_15,
    size = 3
  ) +
  annotate(
    "text",
    x = beta_30 + 0.05,
    y = ymax * 1.2,
    label = "≥30% ↑ odds per doubling\nof dose (e.g., 1→2 or 2→4)",
    hjust = 0,
    vjust = 1,
    color = col_30,
    size = 3
  ) +

  # CrI label
  annotate(
    "text",
    x = mean(posterior_dose_log_skeptical),
    y = ymax * 0.40,
    label = paste0(
      "89% CrI: [",
      round(ci_lower, 2),
      ", ",
      round(ci_upper, 2),
      "]"
    ),
    size = 3.5
  ) +

  # Posterior probability labels
  annotate(
    "text",
    x = beta_30 + 0.05,
    y = ymax * 0.95,
    label = paste0(
      "P(β > 0 | ↑ log(dose) → ↑ odds of response) = ",
      round(p_gt_0 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_gt_0
  ) +
  annotate(
    "text",
    x = beta_30 + 0.05,
    y = ymax * 0.85,
    label = paste0(
      "P(β > ",
      round(beta_15, 2),
      " | ≥15% ↑ odds per doubling) = ",
      round(p_gt_15 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_15
  ) +
  annotate(
    "text",
    x = beta_30 + 0.05,
    y = ymax * 0.75,
    label = paste0(
      "P(β > ",
      round(beta_30, 2),
      " | ≥30% ↑ odds per doubling) = ",
      round(p_gt_30 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_30
  ) +

  scale_y_continuous(expand = expansion(mult = c(0.02, 0.01))) +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.38))) +
  labs(
    title = "Posterior Distribution: Log(Dose) Effect",
    subtitle = "Skeptical Prior",
    x = "β_log_dose (log-odds per log-unit increase in dose)",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

plot_log_skeptical

---

4 Fit model with weakly informative prior

This section fits the primary log-dose Bayesian response model using weakly informative priors. In contrast to the skeptical-prior version, these priors allow somewhat more flexibility for clinically plausible treatment effects while still regularizing estimation and discouraging extreme coefficients unsupported by the data.

The model estimates the probability of response as a function of immunotherapy dose exposure together with key clinical covariates, including age, immunosuppression status, performance status, tumor location, stage, and treatment regimen. The dose coefficient is assigned a prior centered modestly above zero, reflecting the possibility that greater exposure may improve response probability, but with enough uncertainty to allow the data to dominate. Most other coefficients are centered near zero unless there is a strong clinical rationale to encode directional prior information.

This analysis serves as a useful complement to the skeptical-prior model: if posterior inferences remain similar across the two specifications, that strengthens confidence that conclusions are being driven by the observed data rather than by prior choice alone.

# Fit Bayesian logistic regression model using weakly informative priors.
# Outcome:
#   y = binary response indicator
#
# Link function:
#   logit(p), where p is the probability of response
#
# Main exposure:
#   dose = cumulative immunotherapy dose exposure
#
# Covariates:
#   age, immunosuppression, ECOG performance status,
#   tumor location, stage, and treatment regimen

m_log_real_weak <- ulam(
  alist(
    
    # Likelihood: observed responses follow a Bernoulli distribution
    y ~ bernoulli(p),

    # Linear predictor for response probability on the log-odds scale
    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,

    # Intercept prior:
    # alpha_mu is defined earlier from the expected baseline response rate
    # and anchors the model to a clinically plausible baseline probability
    alpha ~ normal(alpha_mu, 1.5),

    # Weakly informative prior for dose effect:
    # centered slightly above zero to reflect a modest expected positive
    # association between greater dose exposure and response, while still
    # allowing the data to move the estimate meaningfully in either direction
    beta_dose ~ normal(0.3, 0.2),

    # Age effect:
    # centered at no effect, with narrow scale because very large changes in
    # response odds per one-year increase in age are not expected
    beta_age ~ normal(0, 0.02),

    # Immunosuppression:
    # centered below zero to reflect the prior belief that immunosuppressed
    # patients may have lower probability of response
    beta_imm ~ normal(-0.8, 0.2),

    # ECOG 2–3 versus reference:
    # mildly negative prior reflecting possible worse outcomes with poorer
    # performance status, but with substantial uncertainty
    beta_ecog ~ normal(-0.3, 0.3),

    # Tumor location priors:
    # conjunctival tumors assigned a strongly negative prior based on clinical
    # expectation of lower responsiveness; other location categories remain
    # centered near zero unless a strong directional belief exists
    beta_conj ~ normal(-2.3, 0.5),
    beta_other ~ normal(0, 0.3),
    beta_unknown ~ normal(0, 0.3),

    # Stage priors:
    # weakly regularizing and centered at zero, allowing the data to estimate
    # whether more advanced stage shifts response probability
    beta_stageIII ~ normal(0, 0.3),
    beta_stageIV ~ normal(0, 0.3),
    beta_stageUn ~ normal(0, 0.3),

    # Regimen priors:
    # centered at zero for pembrolizumab schedule indicators,
    # with a mildly positive prior for ipilimumab+nivolumab reflecting a
    # plausible expectation of higher activity but with meaningful uncertainty
    beta_pembro_q3 ~ normal(0, 0.3),
    beta_pembro_q6 ~ normal(0, 0.3),
    beta_pembro_mixed ~ normal(0, 0.3),
    beta_ipinivo ~ normal(0.4, 0.3)
  ),

  # Analysis dataset prepared earlier for the log-dose model
  data = dat_log,

  # MCMC settings
  chains = 4,
  cores = 4,
  iter = 2000
)

4.1 Posterior predtive check of model with weak prior

post_log_weak <- extract.samples(m_log_real_weak)

# Manual posterior predictive check
# Uses link() to obtain fitted response probabilities, then simulates
# binary response outcomes with rbinom().
ppc_log_weak <- manual_binary_ppc(
  fit = m_log_real_weak,
  observed_y = df_analysis$response,
  n_ppc_draws = 1000,
  seed = 123
)

# Keep these objects with intuitive names for downstream code
p_log_weak <- ppc_log_weak$p_hat
yrep_log_weak <- ppc_log_weak$yrep

observed_orr_log_weak <- ppc_log_weak$observed_orr
simulated_orr_log_weak <- ppc_log_weak$simulated_orr

ppc_log_weak$summary

# A tibble: 1 × 7
  observed_orr mean_simulated_orr median_simulated_orr lower_89 upper_89
         <dbl>              <dbl>                <dbl>    <dbl>    <dbl>
1        0.646              0.642                0.640    0.571    0.714
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>

ppc_log_weak_orr <- tibble(
  simulated_orr = simulated_orr_log_weak
)

ggplot(ppc_log_weak_orr, aes(x = simulated_orr)) +
  geom_histogram(
    bins = 50,
    fill = "steelblue",
    alpha = 0.7,
    color = "white"
  ) +
  geom_vline(
    xintercept = observed_orr_log_weak,
    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: Weak Prior Log-Dose Model",
    subtitle = "Red dashed line = observed ORR",
    x = "Simulated cohort ORR",
    y = "Count"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

# Use manually generated posterior predictive outcomes
# from manual_binary_ppc()
observed_y <- df_analysis$response

stopifnot(ncol(yrep_log_weak) == length(observed_y))

pred_rank_log_weak <- sapply(
  seq_len(ncol(yrep_log_weak)),
  function(i) {
    sum(yrep_log_weak[, i] < observed_y[i])
  }
)

n_sims_log_weak <- nrow(yrep_log_weak)

rank_df_log_weak <- tibble(
  record_id = df_analysis$record_id,
  observed_response = observed_y,
  pred_rank = pred_rank_log_weak,
  pred_rank_scaled = pred_rank / n_sims_log_weak
)

summary_rank_log_weak <- rank_df_log_weak |>
  summarise(
    prop_rank_zero = mean(pred_rank == 0, na.rm = TRUE),
    prop_mid_rank = mean(
      pred_rank_scaled >= 0.40 &
        pred_rank_scaled <= 0.60,
      na.rm = TRUE
    ),
    mean_rank_scaled = mean(pred_rank_scaled, na.rm = TRUE),
    median_rank_scaled = median(pred_rank_scaled, na.rm = TRUE)
  )

summary_rank_log_weak

# A tibble: 1 × 4
  prop_rank_zero prop_mid_rank mean_rank_scaled median_rank_scaled
           <dbl>         <dbl>            <dbl>              <dbl>
1          0.354        0.0741            0.203              0.225

ggplot(rank_df_log_weak, aes(x = pred_rank_scaled)) +
  geom_histogram(
    binwidth = 0.025,
    boundary = 0,
    closed = "left",
    fill = "grey80",
    color = "white"
  ) +
  coord_cartesian(xlim = c(0, 1)) +
  scale_x_continuous(
    breaks = seq(0, 1, by = 0.25),
    labels = scales::percent_format(accuracy = 1),
    expand = expansion(mult = c(0.02, 0.02))
  ) +
  labs(
    title = "Posterior Predictive Rank Distribution",
    subtitle = "Weak-prior log-dose model; binary outcomes produce expected boundary mass",
    x = "Scaled rank of observed outcome within simulated outcomes",
    y = "Number of patients"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

4.2 Posterior parameter interpretation

precis(m_log_real_weak, depth = 2)

                          mean         sd        5.5%       94.5%      rhat
alpha              0.636802266 0.29691912  0.14942875  1.10898488 1.0002252
beta_dose          0.346537838 0.15341145  0.10088202  0.59648412 1.0002361
beta_age           0.002206076 0.01190334 -0.01645847  0.02098780 1.0005818
beta_imm          -0.907910231 0.18116275 -1.19319986 -0.62590679 1.0016692
beta_ecog         -0.308308657 0.25898841 -0.71321374  0.10729852 1.0022858
beta_conj         -2.303397859 0.43338196 -3.00087029 -1.61520774 1.0002021
beta_other        -0.155362374 0.27615227 -0.59028542  0.28085392 1.0000517
beta_unknown       0.050284344 0.28564496 -0.40817277  0.50831930 1.0004514
beta_stageIII     -0.004189566 0.23461609 -0.38001288  0.37113451 0.9999095
beta_stageIV      -0.312139223 0.24075038 -0.69239891  0.06926491 0.9999715
beta_stageUn       0.253999323 0.26517537 -0.17126832  0.66927171 1.0013901
beta_pembro_q3    -0.065954501 0.23844435 -0.44418394  0.31836303 0.9999876
beta_pembro_q6    -0.057150777 0.28591885 -0.51405266  0.39877272 1.0048860
beta_pembro_mixed  0.229496094 0.28417727 -0.21706719  0.68011776 1.0042449
beta_ipinivo       0.426048729 0.29550457 -0.03970642  0.91220761 1.0012683
                  ess_bulk
alpha             3949.596
beta_dose         6479.632
beta_age          8034.289
beta_imm          8737.813
beta_ecog         8049.540
beta_conj         7769.609
beta_other        8317.741
beta_unknown      7324.966
beta_stageIII     4769.746
beta_stageIV      6251.290
beta_stageUn      6804.573
beta_pembro_q3    7311.832
beta_pembro_q6    8670.298
beta_pembro_mixed 8690.068
beta_ipinivo      8223.970

plot(precis(m_log_real_weak))

post_log_weak <- extract.samples(m_log_real_weak)

4.3 Intercept to baseline probability

# In the log-dose model, log_dose = 0 corresponds to 1 dose.
# Therefore plogis(alpha) is the model-implied response probability
# at 1 dose for the reference covariate profile, not the cohort-average ORR.
baseline_prob_log_weak <- plogis(post_log_weak$alpha)

mean(baseline_prob_log_weak)

[1] 0.6510791

quantile(baseline_prob_log_weak, c(0.055, 0.945))

     5.5%     94.5% 
0.5372878 0.7519398 

4.4 Posterior predicted dose-response curve

dose_number <- 1:15

new_data_log <- data.frame(
  dose = log(dose_number),
  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_log_weak_new <- link(m_log_real_weak, data = new_data_log)

mean_p_log_weak <- apply(p_log_weak_new, 2, mean)
ci_log_weak <- apply(p_log_weak_new, 2, PI, 0.89)

plot_df_log_weak <- tibble(
  dose_number = dose_number,
  mean_p = mean_p_log_weak,
  lower = ci_log_weak[1, ],
  upper = ci_log_weak[2, ]
)

plot_predicted_response_log_weak <-
  ggplot(plot_df_log_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 = 1:15,
    expand = expansion(mult = c(0.01, 0.02))
  ) +
  labs(
    title = "Predicted Response Probability by Number of Doses",
    subtitle = "Posterior mean predictions from Bayesian logistic regression\nLog-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_log_weak

4.5 Posterior dose-effect probabilities

post_log_weak <- extract.samples(m_log_real_weak)

posterior_dose_log_weak <- post_log_weak$beta_dose


# OR thresholds for dose doubling
or_doubling_115 <- 1.15
or_doubling_130 <- 1.30

beta_doubling_115 <- log(or_doubling_115) / log(2)
beta_doubling_130 <- log(or_doubling_130) / log(2)

p_gt_0_log_weak <- mean(posterior_dose_log_weak > 0)

p_or_115_doubling_log_weak <- mean(
  posterior_dose_log_weak >= beta_doubling_115
)

p_or_130_doubling_log_weak <- mean(
  posterior_dose_log_weak >= beta_doubling_130
)

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

P(OR per dose doubling ≥ 1.15) = 83% → Good evidence for at least a modest effect

P(OR per dose doubling ≥ 1.30) = 41% → Evidence for a larger effect

4.6 Posterior distribution plot

dens <- density(posterior_dose_log_weak)
dens_df <- data.frame(x = dens$x, y = dens$y)
ymax <- max(dens_df$y)

ci_lower <- quantile(posterior_dose_log_weak, 0.055)
ci_upper <- quantile(posterior_dose_log_weak, 0.945)

ln2 <- log(2)
beta_15 <- log(1.15) / ln2
beta_30 <- log(1.30) / ln2

p_gt_0  <- mean(posterior_dose_log_weak > 0)
p_gt_15 <- mean(posterior_dose_log_weak >= beta_15)
p_gt_30 <- mean(posterior_dose_log_weak >= beta_30)

col_gt_0  <- "#1B9E77"
col_15    <- "#1F78B4"
col_30    <- "#D95F02"
col_line  <- "#08306B"
alpha_fill <- 0.45

plot_log_weak <-
  ggplot() +
  geom_line(
    data = dens_df,
    aes(x = x, y = y),
    linewidth = 0.9,
    color = col_line
  ) +
  geom_area(
    data = subset(dens_df, x > 0 & x < beta_15),
    aes(x = x, y = y),
    fill = col_gt_0,
    alpha = alpha_fill
  ) +
  geom_area(
    data = subset(dens_df, x >= beta_15 & x < beta_30),
    aes(x = x, y = y),
    fill = col_15,
    alpha = alpha_fill
  ) +
  geom_area(
    data = subset(dens_df, x >= beta_30),
    aes(x = x, y = y),
    fill = col_30,
    alpha = alpha_fill
  ) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "red", linewidth = 1) +
  geom_vline(xintercept = beta_15, linetype = "dashed", color = col_15, linewidth = 0.9) +
  geom_vline(xintercept = beta_30, linetype = "dashed", color = col_30, linewidth = 0.9) +
  geom_vline(
    xintercept = c(ci_lower, ci_upper),
    linetype = "dotted",
    color = "grey40",
    linewidth = 0.8
  ) +
  annotate(
    "text",
    x = 0,
    y = ymax * 1.2,
    label = "Null",
    hjust = -0.2,
    vjust = 1,
    color = "red",
    size = 3
  ) +
  annotate(
    "text",
    x = beta_15,
    y = ymax * 1.2,
    label = "≥15% ↑ odds per doubling\nof dose (e.g., 1→2 or 2→4)",
    hjust = -0.05,
    vjust = 1,
    color = col_15,
    size = 3
  ) +
  annotate(
    "text",
    x = beta_30 + 0.05,
    y = ymax * 1.2,
    label = "≥30% ↑ odds per doubling\nof dose (e.g., 1→2 or 2→4)",
    hjust = 0,
    vjust = 1,
    color = col_30,
    size = 3
  ) +
  annotate(
    "text",
    x = mean(posterior_dose_log_weak),
    y = ymax * 0.40,
    label = paste0(
      "89% CrI: [",
      round(ci_lower, 2),
      ", ",
      round(ci_upper, 2),
      "]"
    ),
    size = 3.5
  ) +
  annotate(
    "text",
    x = beta_30 + 0.05,
    y = ymax * 0.95,
    label = paste0(
      "P(β > 0 | ↑ log(dose) → ↑ odds of response) = ",
      round(p_gt_0 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_gt_0
  ) +
  annotate(
    "text",
    x = beta_30 + 0.05,
    y = ymax * 0.85,
    label = paste0(
      "P(β > ",
      round(beta_15, 2),
      " | ≥15% ↑ odds per doubling) = ",
      round(p_gt_15 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_15
  ) +
  annotate(
    "text",
    x = beta_30 + 0.05,
    y = ymax * 0.75,
    label = paste0(
      "P(β > ",
      round(beta_30, 2),
      " | ≥30% ↑ odds per doubling) = ",
      round(p_gt_30 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_30
  ) +
  scale_y_continuous(expand = expansion(mult = c(0.02, 0.01))) +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.38))) +
  labs(
    title = "Posterior Distribution: Log(Dose) Effect",
    subtitle = "Weakly Informative Prior",
    x = "β_log_dose (log-odds per log-unit increase in dose)",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

plot_log_weak

---

4.7 Compare dose effect based on prior

delta_beta_log <- posterior_dose_log_weak - posterior_dose_log_skeptical

# OR thresholds for dose doubling
or_doubling_115 <- 1.15
or_doubling_130 <- 1.30

beta_doubling_115 <- log(or_doubling_115) / log(2)
beta_doubling_130 <- log(or_doubling_130) / log(2)

prior_sensitivity_log <- tibble(
  prior = c("Skeptical", "Weakly informative"),
  mean_beta = c(
    mean(posterior_dose_log_skeptical),
    mean(posterior_dose_log_weak)
  ),
  median_beta = c(
    median(posterior_dose_log_skeptical),
    median(posterior_dose_log_weak)
  ),
  median_or_doubling = c(
    median(exp(log(2) * posterior_dose_log_skeptical)),
    median(exp(log(2) * posterior_dose_log_weak))
  ),
  ci_lower_89 = c(
    quantile(posterior_dose_log_skeptical, 0.055),
    quantile(posterior_dose_log_weak, 0.055)
  ),
  ci_upper_89 = c(
    quantile(posterior_dose_log_skeptical, 0.945),
    quantile(posterior_dose_log_weak, 0.945)
  ),
  p_gt_0 = c(
    mean(posterior_dose_log_skeptical > 0),
    mean(posterior_dose_log_weak > 0)
  ),
  p_or_115_doubling = c(
    mean(posterior_dose_log_skeptical >= beta_doubling_115),
    mean(posterior_dose_log_weak >= beta_doubling_115)
  ),
  p_or_130_doubling = c(
    mean(posterior_dose_log_skeptical >= beta_doubling_130),
    mean(posterior_dose_log_weak >= beta_doubling_130)
  )
)

prior_sensitivity_log

# A tibble: 2 × 9
  prior  mean_beta median_beta median_or_doubling ci_lower_89 ci_upper_89 p_gt_0
  <chr>      <dbl>       <dbl>              <dbl>       <dbl>       <dbl>  <dbl>
1 Skept…     0.177       0.175               1.13     -0.0652       0.419  0.881
2 Weakl…     0.347       0.344               1.27      0.101        0.596  0.992
# ℹ 2 more variables: p_or_115_doubling <dbl>, p_or_130_doubling <dbl>

4.8 Summarize difference directly

tibble(
  delta_mean = mean(delta_beta_log),
  delta_median = median(delta_beta_log),
  delta_lower_89 = quantile(delta_beta_log, 0.055),
  delta_upper_89 = quantile(delta_beta_log, 0.945),
  p_delta_gt_0 = mean(delta_beta_log > 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.170        0.167         -0.178          0.527        0.784

delta_df <- tibble(delta_beta = delta_beta_log)

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.17; 89% CrI -0.18 to 0.53; P(Δβ > 0) = 78.4%), indicating that the estimated dose–response relationship was largely robust to prior specification.

4.9 Posterior predictive checks and prior sensitivity

Under the weakly informative prior, the estimated log-dose effect increased substantially (β_dose ≈ 0.33 vs ≈ 0.16 under the skeptical prior), while the intercept decreased, indicating that the model shifted explanatory weight from the baseline response probability toward treatment exposure. This corresponds to a clinically meaningful increase in the odds of response with increasing dose under the log-dose formulation.

Despite this change in parameter estimates, the posterior predictive distribution of the overall response rate remained centered near the observed value, indicating that both prior specifications yield models that are well-calibrated at the aggregate level.

However, the rank-based posterior predictive check continued to show a concentration of low ranks, consistent with systematic overprediction of response at the individual level. Importantly, this pattern persisted despite the change in prior and the stronger estimated dose effect.

Taken together, these results suggest that the observed miscalibration is not primarily driven by prior specification, but rather reflects limitations of the current model structure. While the log-dose formulation improves the representation of the dose–response relationship, it does not fully capture the variability in response probabilities across patients.

---

5 Fit model with dose intensity parameter

5.1 Build the data list for the dose-intensity-adjusted log ORR model

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

df_analysis |>
  summarise(
    across(
      c(agent_pembro_q3, agent_pembro_q6, agent_pembro_mixed, agent_ipinivo),
      ~ sum(is.na(.))
    )
  )

# A tibble: 1 × 4
  agent_pembro_q3 agent_pembro_q6 agent_pembro_mixed agent_ipinivo
            <int>           <int>              <int>         <int>
1               0               0                  0             0

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

  dose = df_analysis_di$log_dose,
  di = df_analysis_di$dose_intensity_pre_response_uncapped_centered,
  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,

  alpha_mu = qlogis(0.40)
)

tibble(
  variable = names(dat_log_di_uncapped),
  class = map_chr(dat_log_di_uncapped, ~ paste(class(.x), collapse = ",")),
  n = map_int(dat_log_di_uncapped, length),
  n_na = map_int(dat_log_di_uncapped, ~ sum(is.na(.))),
  n_inf = map_int(dat_log_di_uncapped, ~ sum(is.infinite(.)))
)

# A tibble: 18 × 5
   variable           class       n  n_na n_inf
   <chr>              <chr>   <int> <int> <int>
 1 N                  integer     1     0     0
 2 y                  numeric   189     0     0
 3 dose               numeric   189     0     0
 4 di                 numeric   189     0     0
 5 age                numeric   189     0     0
 6 imm                numeric   189     0     0
 7 ecog23             integer   189     0     0
 8 loc_conj           integer   189     0     0
 9 loc_other          integer   189     0     0
10 loc_unknown        integer   189     0     0
11 stage_III          integer   189     0     0
12 stage_IV           integer   189     0     0
13 stage_unstaged     integer   189     0     0
14 agent_pembro_q3    integer   189     0     0
15 agent_pembro_q6    integer   189     0     0
16 agent_pembro_mixed integer   189     0     0
17 agent_ipinivo      integer   189     0     0
18 alpha_mu           numeric     1     0     0

df_analysis_di |>
  summarise(
    n = n(),
    min_di = min(dose_intensity_pre_response_uncapped_centered, na.rm = TRUE),
    max_di = max(dose_intensity_pre_response_uncapped_centered, na.rm = TRUE),
    sd_di = sd(dose_intensity_pre_response_uncapped_centered, na.rm = TRUE)
  )

# A tibble: 1 × 4
      n min_di max_di sd_di
  <int>  <dbl>  <dbl> <dbl>
1   189 -0.624  0.356 0.113

5.1.1 Fit the weak-prior log ORR model with uncapped pre-response dose intensity

m_log_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_log_di_uncapped,
  chains = 4,
  cores = 4,
  iter = 2000
)

5.2 Diagnostics

precis(m_log_real_weak_di_uncapped, depth = 2)

                          mean         sd         5.5%       94.5%      rhat
alpha              0.838059076 0.27082420  0.402768622  1.28368332 0.9999936
beta_dose          0.167539072 0.10350549  0.002929216  0.33205914 1.0006542
beta_di            0.105237093 0.12077715 -0.090104961  0.29950207 1.0016908
beta_age           0.001242610 0.01166573 -0.017472444  0.01993892 1.0008494
beta_imm          -0.906334551 0.18251229 -1.203073937 -0.61214857 0.9997860
beta_ecog         -0.312293444 0.26629476 -0.737751423  0.12330153 1.0000489
beta_conj         -2.306290309 0.43064342 -3.004301685 -1.61493272 1.0001574
beta_other        -0.140743752 0.27213156 -0.571209280  0.29456653 1.0022576
beta_unknown       0.063839669 0.27965203 -0.389318108  0.51367286 1.0005094
beta_stageIII     -0.007234724 0.23852370 -0.392325966  0.35857860 1.0003060
beta_stageIV      -0.315754439 0.23597342 -0.699173012  0.05262411 1.0004205
beta_stageUn       0.260731167 0.25987109 -0.156313381  0.67024975 1.0003241
beta_pembro_q3    -0.063042711 0.23764092 -0.440877116  0.31911562 1.0005863
beta_pembro_q6    -0.053713119 0.28425247 -0.510414009  0.39357533 1.0039006
beta_pembro_mixed  0.249716365 0.28355424 -0.219112486  0.70374974 1.0007687
beta_ipinivo       0.422778396 0.30807619 -0.072099353  0.91805947 1.0017888
                  ess_bulk
alpha             4759.047
beta_dose         6411.780
beta_di           7589.799
beta_age          7564.930
beta_imm          7295.208
beta_ecog         6524.894
beta_conj         8533.960
beta_other        7124.942
beta_unknown      6823.257
beta_stageIII     6316.059
beta_stageIV      6031.708
beta_stageUn      7769.727
beta_pembro_q3    6578.617
beta_pembro_q6    7088.099
beta_pembro_mixed 8239.106
beta_ipinivo      7384.715

trankplot(m_log_real_weak_di_uncapped)

plot(precis(m_log_real_weak_di_uncapped))

5.3 Posterior predictive check

post_log_weak_di_uncapped <- extract.samples(m_log_real_weak_di_uncapped)

# Manual posterior predictive check
# Uses link() to obtain fitted response probabilities, then simulates
# binary response outcomes with rbinom().
#
# Important: this model was fit on df_analysis_di, so the observed outcome
# vector must come from df_analysis_di rather than df_analysis.
ppc_log_weak_di_uncapped <- manual_binary_ppc(
  fit = m_log_real_weak_di_uncapped,
  observed_y = df_analysis_di$response,
  n_ppc_draws = 1000,
  seed = 123
)

# Keep these objects with intuitive names for downstream code
p_log_weak_di_uncapped <- ppc_log_weak_di_uncapped$p_hat
yrep_log_weak_di_uncapped <- ppc_log_weak_di_uncapped$yrep

observed_orr_log_weak_di_uncapped <- ppc_log_weak_di_uncapped$observed_orr
simulated_orr_log_weak_di_uncapped <- ppc_log_weak_di_uncapped$simulated_orr

ppc_log_weak_di_uncapped$summary

# A tibble: 1 × 7
  observed_orr mean_simulated_orr median_simulated_orr lower_89 upper_89
         <dbl>              <dbl>                <dbl>    <dbl>    <dbl>
1        0.646              0.643                0.640    0.571    0.714
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>

tibble(
  observed_orr = observed_orr_log_weak_di_uncapped,
  mean_simulated_orr = mean(simulated_orr_log_weak_di_uncapped),
  median_simulated_orr = median(simulated_orr_log_weak_di_uncapped),
  lower_89 = quantile(simulated_orr_log_weak_di_uncapped, 0.055),
  upper_89 = quantile(simulated_orr_log_weak_di_uncapped, 0.945)
)

# A tibble: 1 × 5
  observed_orr mean_simulated_orr median_simulated_orr lower_89 upper_89
         <dbl>              <dbl>                <dbl>    <dbl>    <dbl>
1        0.646              0.643                0.640    0.571    0.714

ppc_log_weak_di_orr <- tibble(
  simulated_orr = simulated_orr_log_weak_di_uncapped
)

ggplot(ppc_log_weak_di_orr, aes(x = simulated_orr)) +
  geom_histogram(
    bins = 50,
    fill = "steelblue",
    alpha = 0.7,
    color = "white"
  ) +
  geom_vline(
    xintercept = observed_orr_log_weak_di_uncapped,
    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: Log-Dose Model + Dose Intensity",
    subtitle = "Red dashed line = observed ORR",
    x = "Simulated cohort ORR",
    y = "Count"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

# Use manually generated posterior predictive outcomes
# from manual_binary_ppc().
#
# Important: this model was fit on df_analysis_di, so the observed outcome
# vector must also come from df_analysis_di rather than df_analysis.
observed_y_log_di <- df_analysis_di$response

stopifnot(ncol(yrep_log_weak_di_uncapped) == length(observed_y_log_di))

pred_rank_log_weak_di_uncapped <- sapply(
  seq_len(ncol(yrep_log_weak_di_uncapped)),
  function(i) {
    sum(yrep_log_weak_di_uncapped[, i] < observed_y_log_di[i])
  }
)

n_sims_log_weak_di_uncapped <- nrow(yrep_log_weak_di_uncapped)

rank_df_log_weak_di_uncapped <- tibble(
  record_id = df_analysis_di$record_id,
  observed_response = observed_y_log_di,
  pred_rank = pred_rank_log_weak_di_uncapped,
  pred_rank_scaled = pred_rank / n_sims_log_weak_di_uncapped
)

summary_rank_log_weak_di_uncapped <- rank_df_log_weak_di_uncapped |>
  summarise(
    prop_rank_zero = mean(pred_rank == 0, na.rm = TRUE),
    prop_mid_rank = mean(
      pred_rank_scaled >= 0.40 &
        pred_rank_scaled <= 0.60,
      na.rm = TRUE
    ),
    mean_rank_scaled = mean(pred_rank_scaled, na.rm = TRUE),
    median_rank_scaled = median(pred_rank_scaled, na.rm = TRUE)
  )

summary_rank_log_weak_di_uncapped

# A tibble: 1 × 4
  prop_rank_zero prop_mid_rank mean_rank_scaled median_rank_scaled
           <dbl>         <dbl>            <dbl>              <dbl>
1          0.354        0.0688            0.203              0.247

ggplot(rank_df_log_weak_di_uncapped, aes(x = pred_rank_scaled)) +
  geom_histogram(
    binwidth = 0.025,
    boundary = 0,
    closed = "left",
    fill = "grey80",
    color = "white"
  ) +
  coord_cartesian(xlim = c(0, 1)) +
  scale_x_continuous(
    breaks = seq(0, 1, by = 0.25),
    labels = scales::percent_format(accuracy = 1),
    expand = expansion(mult = c(0.02, 0.02))
  ) +
  labs(
    title = "Posterior Predictive Rank Distribution",
    subtitle = "Weak-prior log-dose model + dose intensity; \nbinary outcomes produce expected boundary mass",
    x = "Scaled rank of observed outcome within simulated outcomes",
    y = "Number of patients"
  ) +
  theme_minimal(base_size = 13) +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

5.4 Extract posterior dose and dose-intensity effects

post_log_weak_di_uncapped <- extract.samples(m_log_real_weak_di_uncapped)

posterior_dose_log_weak_di_uncapped <- post_log_weak_di_uncapped$beta_dose
posterior_di_log_weak_uncapped <- post_log_weak_di_uncapped$beta_di

# OR thresholds for dose doubling
or_doubling_115 <- 1.15
or_doubling_130 <- 1.30

beta_doubling_115 <- log(or_doubling_115) / log(2)
beta_doubling_130 <- log(or_doubling_130) / log(2)

p_gt_0_log_weak_di <- mean(
  posterior_dose_log_weak_di_uncapped > 0
)

p_or_115_doubling_log_weak_di <- mean(
  posterior_dose_log_weak_di_uncapped >= beta_doubling_115
)

p_or_130_doubling_log_weak_di <- mean(
  posterior_dose_log_weak_di_uncapped >= beta_doubling_130
)

p_di_gt_0 <- mean(
  posterior_di_log_weak_uncapped > 0
)

Dose term after adjustment for uncapped pre-response dose intensity

P(β_dose > 0) = 95%

P(OR per dose doubling ≥ 1.15) = 38%

P(OR per dose doubling ≥ 1.30) = 2%

Uncapped pre-response dose-intensity term

P(β_DI > 0) = 81%

5.5 Compare beta_dose with and without dose intensity

delta_beta_dose_di_adjustment <-
  posterior_dose_log_weak_di_uncapped - posterior_dose_log_weak



dose_di_comparison <- tibble(
  model = c(
    "Weak prior log model",
    "Weak prior log model + uncapped pre-response dose intensity"
  ),
  mean_beta_dose = c(
    mean(posterior_dose_log_weak),
    mean(posterior_dose_log_weak_di_uncapped)
  ),
  median_beta_dose = c(
    median(posterior_dose_log_weak),
    median(posterior_dose_log_weak_di_uncapped)
  ),
  ci_lower_89 = c(
    quantile(posterior_dose_log_weak, 0.055),
    quantile(posterior_dose_log_weak_di_uncapped, 0.055)
  ),
  ci_upper_89 = c(
    quantile(posterior_dose_log_weak, 0.945),
    quantile(posterior_dose_log_weak_di_uncapped, 0.945)
  ),
  p_gt_0 = c(
    mean(posterior_dose_log_weak > 0),
    mean(posterior_dose_log_weak_di_uncapped > 0)
  ),
  p_gt_005 = c(
    mean(posterior_dose_log_weak > 0.05),
    mean(posterior_dose_log_weak_di_uncapped > 0.05)
  ),
  p_gt_010 = c(
    mean(posterior_dose_log_weak > 0.10),
    mean(posterior_dose_log_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.347            0.344     0.101         0.596  0.992    0.977
2 Weak …          0.168            0.168     0.00293       0.332  0.947    0.871
# ℹ 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.179       -0.175         -0.477          0.110        0.160

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 log ORR model and comparing the posterior dose-effect estimate before and after adjustment. The posterior difference in the dose coefficient was centered at -0.175, with an 89% credible interval from -0.477 to 0.11. 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_log_di_uncapped <- data.frame(
  dose = log(dose_number),
  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_log_weak_di_new <- link(
  m_log_real_weak_di_uncapped,
  data = new_data_log_di_uncapped
)

mean_p_log_weak_di <- apply(p_log_weak_di_new, 2, mean)
ci_log_weak_di <- apply(p_log_weak_di_new, 2, PI, 0.89)

plot_df_log_weak_di <- tibble(
  dose_number = doses + 1,
  mean_p = mean_p_log_weak_di,
  lower = ci_log_weak_di[1, ],
  upper = ci_log_weak_di[2, ]
)



plot_predicted_response_log_weak_di_uncapped <-
  ggplot(plot_df_log_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_log_weak_di$dose_number),
                 max(plot_df_log_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 log 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_log_weak_di_uncapped

5.7 Posterior distribution plot for beta_dose after dose-intensity adjustment

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

# 89% CrI
ci_lower <- quantile(posterior_dose_log_weak_di_uncapped, 0.055)
ci_upper <- quantile(posterior_dose_log_weak_di_uncapped, 0.945)

# --- Doubling thresholds ---
ln2 <- log(2)
beta_15 <- log(1.15) / ln2   # ~0.2016
beta_30 <- log(1.30) / ln2   # ~0.3785

# Tail probabilities
p_gt_0  <- mean(posterior_dose_log_weak_di_uncapped > 0)
p_gt_15 <- mean(posterior_dose_log_weak_di_uncapped > beta_15)
p_gt_30 <- mean(posterior_dose_log_weak_di_uncapped > beta_30)

# Colors
col_gt_0  <- "#1B9E77"
col_15    <- "#1F78B4"
col_30    <- "#D95F02"
col_line  <- "#08306B"
alpha_fill <- 0.45

plot_log_weak_di_uncapped <-
  ggplot() +

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

  # Non-overlapping shaded regions
  geom_area(
    data = subset(dens_df, x > 0 & x < beta_15),
    aes(x = x, y = y),
    fill = col_gt_0,
    alpha = alpha_fill
  ) +
  geom_area(
    data = subset(dens_df, x >= beta_15 & x < beta_30),
    aes(x = x, y = y),
    fill = col_15,
    alpha = alpha_fill
  ) +
  geom_area(
    data = subset(dens_df, x >= beta_30),
    aes(x = x, y = y),
    fill = col_30,
    alpha = alpha_fill
  ) +

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

  # Top labels
  annotate(
    "text",
    x = 0,
    y = ymax * 1.2,
    label = "Null",
    hjust = -0.2,
    vjust = 1,
    color = "red",
    size = 3
  ) +
  annotate(
    "text",
    x = beta_15,
    y = ymax * 1.2,
    label = "≥15% ↑ odds per doubling\nof dose (e.g., 1→2 or 2→4)",
    hjust = -0.05,
    vjust = 1,
    color = col_15,
    size = 3
  ) +
  annotate(
    "text",
    x = 0.5,
    y = ymax * 1.2,
    label = "≥30% ↑ odds per doubling\nof dose (e.g., 1→2 or 2→4)",
    hjust = 0,
    vjust = 1,
    color = col_30,
    size = 3
  ) +

  # CrI label
  annotate(
    "text",
    x = mean(posterior_dose_log_weak_di_uncapped),
    y = ymax * 0.40,
    label = paste0(
      "89% CrI: [",
      round(ci_lower, 2),
      ", ",
      round(ci_upper, 2),
      "]"
    ),
    size = 3.5
  ) +

  # Posterior probability labels
  annotate(
    "text",
    x = 0.4,
    y = ymax * 0.95,
    label = paste0(
      "P(β > 0 | ↑ log(dose) → ↑ odds of response) = ",
      round(p_gt_0 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_gt_0
  ) +
  annotate(
    "text",
    x = 0.4,
    y = ymax * 0.85,
    label = paste0(
      "P(β > ",
      round(beta_15, 2),
      " | ≥15% ↑ odds per doubling) = ",
      round(p_gt_15 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_15
  ) +
  annotate(
    "text",
    x = 0.4,
    y = ymax * 0.75,
    label = paste0(
      "P(β > ",
      round(beta_30, 2),
      " | ≥30% ↑ odds per doubling) = ",
      round(p_gt_30 * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_30
  ) +

  scale_y_continuous(expand = expansion(mult = c(0.02, 0.01))) +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.38))) +
  labs(
    title = "Posterior Distribution: Log(Dose) Effect",
    subtitle = "Weak prior model adjusted for uncapped pre-response dose intensity",
    x = "β_log_dose (log-odds per log-unit increase in dose)",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

plot_log_weak_di_uncapped

5.8 Companion plot for the dose-intensity coefficient itself

dens_di <- density(posterior_di_log_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_log_weak_uncapped, 0.055)
ci_upper_di <- quantile(posterior_di_log_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_log_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(
      log_dose,
      dose_intensity_pre_response_uncapped_centered,
      use = "complete.obs"
    )
  )

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

ggplot(
  df_analysis,
  aes(x = log_dose, 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_log_weak_di_uncapped <- extract.samples(m_log_real_weak_di_uncapped)

posterior_dose_log_weak_di_uncapped <- post_log_weak_di_uncapped$beta_dose
posterior_di_log_weak_uncapped <- post_log_weak_di_uncapped$beta_di

mean(posterior_dose_log_weak_di_uncapped > 0)

[1] 0.947

mean(posterior_dose_log_weak_di_uncapped > 0.05)

[1] 0.87125

mean(posterior_dose_log_weak_di_uncapped > 0.10)

[1] 0.73875

mean(posterior_di_log_weak_uncapped > 0)

[1] 0.807

mean(posterior_di_log_weak_uncapped > 0.05)

[1] 0.6775

mean(posterior_di_log_weak_uncapped > 0.10)

[1] 0.51975

delta_beta_dose_di_adjustment <-
  posterior_dose_log_weak_di_uncapped - posterior_dose_log_weak

mean(delta_beta_dose_di_adjustment)

[1] -0.1789988

median(delta_beta_dose_di_adjustment)

[1] -0.1749719

quantile(delta_beta_dose_di_adjustment, c(0.055, 0.945))

      5.5%      94.5% 
-0.4774642  0.1096967 

mean(delta_beta_dose_di_adjustment > 0)

[1] 0.15975

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

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

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_log_weak_di_uncapped,
  beta_di   = posterior_di_log_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 log 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)
  )

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.108477), indicating that the estimated benefit associated with additional doses was not explained by differences in treatment timing or dose intensity.

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

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

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

df_analysis_log_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(log_dose_pre_response),
    !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_log_pre_response |>
  summarise(
    n = n(),
    observed_orr = mean(response),
    min_pre_response_doses = min(dose_pre_response, na.rm = TRUE),
    max_pre_response_doses = max(dose_pre_response, na.rm = TRUE),
    mean_total_doses = mean(num_doses, na.rm = TRUE),
    mean_pre_response_doses = mean(dose_pre_response, na.rm = TRUE),
    n_total_differs_from_pre_response = sum(num_doses != dose_pre_response, na.rm = TRUE)
  )

# A tibble: 1 × 7
      n observed_orr min_pre_response_doses max_pre_response_doses
  <int>        <dbl>                  <int>                  <int>
1   189        0.646                      1                     20
# ℹ 3 more variables: mean_total_doses <dbl>, mean_pre_response_doses <dbl>,
#   n_total_differs_from_pre_response <int>

dat_log_pre_response <- list(
  N = nrow(df_analysis_log_pre_response),
  y = df_analysis_log_pre_response$response,

  # Key sensitivity exposure:
  # log of doses administered before response assessment
  dose = df_analysis_log_pre_response$log_dose_pre_response,

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

  loc_conj = df_analysis_log_pre_response$loc_conj,
  loc_other = df_analysis_log_pre_response$loc_other,
  loc_unknown = df_analysis_log_pre_response$loc_unknown,

  stage_III = df_analysis_log_pre_response$stage_III,
  stage_IV = df_analysis_log_pre_response$stage_IV,
  stage_unstaged = df_analysis_log_pre_response$stage_unstaged,

  agent_pembro_q3 = df_analysis_log_pre_response$agent_pembro_q3,
  agent_pembro_q6 = df_analysis_log_pre_response$agent_pembro_q6,
  agent_pembro_mixed = df_analysis_log_pre_response$agent_pembro_mixed,
  agent_ipinivo = df_analysis_log_pre_response$agent_ipinivo,

  alpha_mu = qlogis(0.40)
)

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

    # Same weakly informative prior as the total-dose log model
    beta_dose ~ normal(0.3, 0.2),

    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_log_pre_response,
  chains = 4,
  cores = 4,
  iter = 2000
)

precis(m_log_real_weak_pre_response, depth = 2)

                          mean         sd        5.5%       94.5%     rhat
alpha              0.611830370 0.29985735  0.13899634  1.10064167 1.001853
beta_dose          0.400252797 0.15706665  0.14892751  0.65451343 1.001829
beta_age           0.001462052 0.01158757 -0.01672945  0.01977698 1.001388
beta_imm          -0.901502040 0.18261906 -1.18779732 -0.60991110 1.002674
beta_ecog         -0.306419843 0.26200048 -0.72545984  0.11600365 1.004593
beta_conj         -2.300635599 0.43803303 -3.00794549 -1.60459418 1.000891
beta_other        -0.138185352 0.26493570 -0.55905170  0.28775353 1.001462
beta_unknown       0.051749681 0.28412354 -0.40709336  0.49915216 1.000996
beta_stageIII     -0.006854824 0.23743662 -0.37748521  0.37904926 1.000696
beta_stageIV      -0.316814468 0.24198798 -0.70105855  0.07449127 1.000570
beta_stageUn       0.246067934 0.26384072 -0.18151258  0.66565021 1.000295
beta_pembro_q3    -0.054548882 0.23554700 -0.43890945  0.32551691 1.001714
beta_pembro_q6    -0.061687254 0.28629515 -0.52517766  0.39873296 1.002262
beta_pembro_mixed  0.227561387 0.28116697 -0.22363204  0.67110327 1.002266
beta_ipinivo       0.428019969 0.27958753 -0.01299929  0.88037120 1.004211
                  ess_bulk
alpha             3202.560
beta_dose         5551.984
beta_age          8151.086
beta_imm          8615.921
beta_ecog         7955.251
beta_conj         8668.444
beta_other        8136.451
beta_unknown      8515.615
beta_stageIII     4976.034
beta_stageIV      4366.091
beta_stageUn      6577.853
beta_pembro_q3    8102.045
beta_pembro_q6    8623.000
beta_pembro_mixed 7686.950
beta_ipinivo      8568.276

plot(precis(m_log_real_weak_pre_response))

post_log_weak_pre_response <- extract.samples(m_log_real_weak_pre_response)

posterior_dose_log_weak_pre_response <-
  post_log_weak_pre_response$beta_dose

save_files(
  save_object = posterior_dose_log_weak_pre_response,
  filename = "Posterior Dose Log Weak Prior PreResponse Dose",
  subD = "files/Modeling/Posterior Dose Log 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 Log Weak Prior PreResponse Dose/Posterior Dose Log Weak Prior PreResponse Dose-2026-05-09-05-38-53.rds 

ln2 <- log(2)

beta_15 <- log(1.15) / ln2
beta_30 <- log(1.30) / ln2

p_gt_0_log_pre_response  <- mean(posterior_dose_log_weak_pre_response > 0)
p_gt_15_log_pre_response <- mean(posterior_dose_log_weak_pre_response > beta_15)
p_gt_30_log_pre_response <- mean(posterior_dose_log_weak_pre_response > beta_30)

log_pre_response_summary <- tibble(
  parameter = "beta_log_dose_pre_response",
  mean_beta = mean(posterior_dose_log_weak_pre_response),
  median_beta = median(posterior_dose_log_weak_pre_response),
  ci_lower_89 = quantile(posterior_dose_log_weak_pre_response, 0.055),
  ci_upper_89 = quantile(posterior_dose_log_weak_pre_response, 0.945),
  p_gt_0 = p_gt_0_log_pre_response,
  p_gt_15pct_odds_per_doubling = p_gt_15_log_pre_response,
  p_gt_30pct_odds_per_doubling = p_gt_30_log_pre_response,
  median_or_per_doubling =
    exp(median(posterior_dose_log_weak_pre_response) * log(2))
)

log_pre_response_summary

# A tibble: 1 × 9
  parameter                 mean_beta median_beta ci_lower_89 ci_upper_89 p_gt_0
  <chr>                         <dbl>       <dbl>       <dbl>       <dbl>  <dbl>
1 beta_log_dose_pre_respon…     0.400       0.401       0.149       0.655  0.994
# ℹ 3 more variables: p_gt_15pct_odds_per_doubling <dbl>,
#   p_gt_30pct_odds_per_doubling <dbl>, median_or_per_doubling <dbl>

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

set.seed(123)

p_hat_log_pre_response <- link(m_log_real_weak_pre_response)

stopifnot(ncol(p_hat_log_pre_response) == nrow(df_analysis_log_pre_response))

n_ppc_draws <- 1000

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

p_hat_ppc_log_pre_response <-
  p_hat_log_pre_response[draw_ids, , drop = FALSE]

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

observed_orr_log_pre_response <-
  mean(df_analysis_log_pre_response$response)

ppc_log_pre_response_orr <- tibble(
  draw = seq_len(nrow(yrep_log_pre_response)),
  simulated_orr = rowMeans(yrep_log_pre_response)
)

ppc_log_pre_response_orr |>
  summarise(
    observed_orr = observed_orr_log_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_log_pre_response),
    mean_fitted_probability = mean(p_hat_log_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.566    0.714
# ℹ 2 more variables: p_sim_le_observed <dbl>, mean_fitted_probability <dbl>

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

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

dose_number_pre_response <- 1:max_pre_response_dose_for_plot

new_data_log_pre_response <- data.frame(
  dose = log(dose_number_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_log_weak_pre_response_new <- link(
  m_log_real_weak_pre_response,
  data = new_data_log_pre_response
)

mean_p_log_weak_pre_response <- apply(
  p_log_weak_pre_response_new,
  2,
  mean
)

ci_log_weak_pre_response <- apply(
  p_log_weak_pre_response_new,
  2,
  PI,
  0.89
)

plot_df_log_weak_pre_response <- tibble(
  dose_number = dose_number_pre_response,
  mean_p = mean_p_log_weak_pre_response,
  lower = ci_log_weak_pre_response[1, ],
  upper = ci_log_weak_pre_response[2, ]
)

plot_predicted_response_log_weak_pre_response <-
  ggplot(plot_df_log_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_log_weak_pre_response$dose_number),
      max(plot_df_log_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 log-dose 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_log_weak_pre_response

dens_log_pre_response <- density(posterior_dose_log_weak_pre_response)

dens_df_log_pre_response <- data.frame(
  x = dens_log_pre_response$x,
  y = dens_log_pre_response$y
)

ymax_log_pre_response <- max(dens_df_log_pre_response$y)

ci_lower_log_pre_response <- quantile(
  posterior_dose_log_weak_pre_response,
  0.055
)

ci_upper_log_pre_response <- quantile(
  posterior_dose_log_weak_pre_response,
  0.945
)

ln2 <- log(2)
beta_15 <- log(1.15) / ln2
beta_30 <- log(1.30) / ln2

p_gt_0_log_pre_response  <- mean(posterior_dose_log_weak_pre_response > 0)
p_gt_15_log_pre_response <- mean(posterior_dose_log_weak_pre_response > beta_15)
p_gt_30_log_pre_response <- mean(posterior_dose_log_weak_pre_response > beta_30)

col_gt_0  <- "#1B9E77"
col_15    <- "#1F78B4"
col_30    <- "#D95F02"
col_line  <- "#08306B"
alpha_fill <- 0.45

plot_log_weak_pre_response <-
  ggplot() +
  geom_line(
    data = dens_df_log_pre_response,
    aes(x = x, y = y),
    linewidth = 0.9,
    color = col_line
  ) +
  geom_area(
    data = subset(dens_df_log_pre_response, x > 0 & x < beta_15),
    aes(x = x, y = y),
    fill = col_gt_0,
    alpha = alpha_fill
  ) +
  geom_area(
    data = subset(dens_df_log_pre_response, x >= beta_15 & x < beta_30),
    aes(x = x, y = y),
    fill = col_15,
    alpha = alpha_fill
  ) +
  geom_area(
    data = subset(dens_df_log_pre_response, x >= beta_30),
    aes(x = x, y = y),
    fill = col_30,
    alpha = alpha_fill
  ) +
  geom_vline(
    xintercept = 0,
    linetype = "dashed",
    color = "red",
    linewidth = 1
  ) +
  geom_vline(
    xintercept = beta_15,
    linetype = "dashed",
    color = col_15,
    linewidth = 0.9
  ) +
  geom_vline(
    xintercept = beta_30,
    linetype = "dashed",
    color = col_30,
    linewidth = 0.9
  ) +
  geom_vline(
    xintercept = c(ci_lower_log_pre_response, ci_upper_log_pre_response),
    linetype = "dotted",
    color = "grey40",
    linewidth = 0.8
  ) +
  annotate(
    "text",
    x = mean(posterior_dose_log_weak_pre_response),
    y = ymax_log_pre_response * 0.40,
    label = paste0(
      "89% CrI: [",
      round(ci_lower_log_pre_response, 2),
      ", ",
      round(ci_upper_log_pre_response, 2),
      "]"
    ),
    size = 3.5
  ) +
  annotate(
    "text",
    x = beta_30 + 0.15,
    y = ymax_log_pre_response * 0.95,
    label = paste0(
      "P(β > 0 | ↑ log(pre-response dose) → ↑ odds of response) = ",
      round(p_gt_0_log_pre_response * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_gt_0
  ) +
  annotate(
    "text",
    x = beta_30 + 0.15,
    y = ymax_log_pre_response * 0.85,
    label = paste0(
      "P(β > ",
      round(beta_15, 2),
      " | ≥15% ↑ odds per doubling) = ",
      round(p_gt_15_log_pre_response * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_15
  ) +
  annotate(
    "text",
    x = beta_30 + 0.15,
    y = ymax_log_pre_response * 0.75,
    label = paste0(
      "P(β > ",
      round(beta_30, 2),
      " | ≥30% ↑ odds per doubling) = ",
      round(p_gt_30_log_pre_response * 100, 0),
      "%"
    ),
    hjust = 0,
    size = 3.4,
    color = col_30
  ) +
  scale_y_continuous(expand = expansion(mult = c(0.02, 0.01))) +
  scale_x_continuous(expand = expansion(mult = c(0.02, 0.38))) +
  labs(
    title = "Posterior Distribution: Pre-Response Log(Dose) Effect",
    subtitle = "Weakly informative prior",
    x = "β_log_dose (log-odds per log-unit increase in pre-response dose)",
    y = "Density"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, face = "bold"),
    plot.subtitle = element_text(hjust = 0.5, face = "bold")
  )

plot_log_weak_pre_response

log_total_vs_pre_response_comparison <- tibble(
  model = c(
    "Total cumulative log-dose",
    "Pre-response log-dose"
  ),
  mean_beta_dose = c(
    mean(posterior_dose_log_weak),
    mean(posterior_dose_log_weak_pre_response)
  ),
  median_beta_dose = c(
    median(posterior_dose_log_weak),
    median(posterior_dose_log_weak_pre_response)
  ),
  ci_lower_89 = c(
    quantile(posterior_dose_log_weak, 0.055),
    quantile(posterior_dose_log_weak_pre_response, 0.055)
  ),
  ci_upper_89 = c(
    quantile(posterior_dose_log_weak, 0.945),
    quantile(posterior_dose_log_weak_pre_response, 0.945)
  ),
  p_gt_0 = c(
    mean(posterior_dose_log_weak > 0),
    mean(posterior_dose_log_weak_pre_response > 0)
  ),
  p_gt_15pct_odds_per_doubling = c(
    mean(posterior_dose_log_weak > beta_15),
    mean(posterior_dose_log_weak_pre_response > beta_15)
  ),
  p_gt_30pct_odds_per_doubling = c(
    mean(posterior_dose_log_weak > beta_30),
    mean(posterior_dose_log_weak_pre_response > beta_30)
  ),
  median_or_per_doubling = c(
    exp(median(posterior_dose_log_weak) * log(2)),
    exp(median(posterior_dose_log_weak_pre_response) * log(2))
  )
)

log_total_vs_pre_response_comparison

# A tibble: 2 × 9
  model           mean_beta_dose median_beta_dose ci_lower_89 ci_upper_89 p_gt_0
  <chr>                    <dbl>            <dbl>       <dbl>       <dbl>  <dbl>
1 Total cumulati…          0.347            0.344       0.101       0.596  0.992
2 Pre-response l…          0.400            0.401       0.149       0.655  0.994
# ℹ 3 more variables: p_gt_15pct_odds_per_doubling <dbl>,
#   p_gt_30pct_odds_per_doubling <dbl>, median_or_per_doubling <dbl>

9 Conclusion

This document evaluated the relationship between treatment exposure and objective response using a Bayesian ORR model with a log-dose specification. This formulation was chosen to test a clinically plausible nonlinear pattern in which early increases in dose exposure may matter more than later ones.

By fitting the model to the observed cohort, examining posterior predictive performance, comparing skeptical and weakly informative priors, and incorporating dose-intensity sensitivity analyses, this workflow separates three related questions: whether higher exposure is associated with response, whether that association is sensitive to assumptions about functional form and prior specification, and whether it is materially altered after accounting for treatment timing before response assessment.

The log-dose framework is especially useful because it moves beyond a simple constant-per-dose assumption and asks whether the exposure–response relationship appears to taper as treatment accumulates. The sensitivity analyses further clarify whether the modeled dose effect reflects dose count itself or is substantially explained by variation in treatment intensity.

As a companion to the manuscript, this page is intended to provide a transparent record of how the log-dose ORR analysis was constructed, checked, and interpreted in the observed data.