---
title: "PD-(L)1 immunotherapy MBMA in metastatic NSCLC (Franzese 2026)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{PD-(L)1 immunotherapy MBMA in metastatic NSCLC (Franzese 2026)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4.5)
library(nlmixr2lib)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
set.seed(20260724)
```

# Model and source

```{r model-load, include = FALSE}
mod_fn <- readModelDb("Franzese_2026_pdl1_nsclc_mbma")
mod_meta <- mod_fn()
mod <- rxode2::rxode(mod_fn)
```

* Citation: `r mod_meta$reference`
* Description: `r mod_meta$description`
* Article: [CPT Pharmacometrics Syst Pharmacol. 2026](https://doi.org/10.1002/psp4.70196)
* Supplement: Appendix S1, Tables S1-S4, Figures S1-S7

This is a model-based meta-analysis (MBMA) of Objective Response Rate
(ORR), Overall Survival (OS), and Progression-Free Survival (PFS) for
programmed cell death protein 1 (PD-1) and programmed cell death
ligand 1 (PD-L1) inhibitors in metastatic non-small cell lung cancer
(mNSCLC), assembled from 114 studies (46 unique treatments) in the
Certara CODEX PD1ST database. The three sub-models are fit
sequentially: first the ORR mixed-effects logistic regression, then
the OS and PFS mixed-effects semi-parametric proportional hazards
models with predicted ORR as an input covariate. The packaged model
in `inst/modeldb/therapeuticArea/oncology/Franzese_2026_pdl1_nsclc_mbma.R`
integrates all three sub-models in a single `rxUi` object so that a
simulation of a treatment arm returns per-arm ORR + per-arm S_OS(t) +
per-arm S_PFS(t) survival curves in one call.

# Population

The analysis pools 114 mNSCLC studies (Phase 1, n = 38; Phase 2,
n = 38; Phase 3, n = 33; undefined, n = 5), from which:

* **ORR**: 197 arms, 284 strata arms across 114 studies.
* **OS**: 147 arms, 187 strata arms across 87 studies.
* **PFS**: 154 arms, 215 strata arms across 88 studies.

Median (range) follow-up: 13 (1.5-54) months for ORR; 21 (6-63) months
maturity for PFS; 29 (6-75) months maturity for OS (source paper
Section 3.1 and Table S3).

Six standing covariates are consumed by the model:

| Covariate | Type | Description |
|-----------|------|-------------|
| `TRT` | integer 1-46 | Specific treatment intercept for ORR + treatment category (1-5) for OS/PFS. Coding in `mod_meta$covariateData$TRT$notes`. |
| `PDL1_TUM` | continuous 0-100 | Per-arm average tumor PD-L1 expression (Franzese 2026 Equation S1). |
| `LINE_1L` | binary 0/1 | 1 = first-line arm, 0 = second-line-or-later arm. |
| `TUMTP_SQUAM_PCT` | continuous 0-100 | Per-arm percent squamous histology (100 - percent non-squamous). |
| `PS_ECOG_0_PCT` | continuous 0-100 | Per-arm percent participants with ECOG PS score of 0. |
| `RACE_ASIAN_PCT` | continuous 0-100 | Per-arm percent Asian race. |

The metadata is available programmatically:

```{r population, message = FALSE}
str(mod_meta$population)
```

**MBMA scope disclaimer.** All simulations are per study-strata-arm
mean quantities: ORR is the arm's predicted objective-response
fraction; S_OS(t) and S_PFS(t) are the arm's typical survival curves
under the paper's proportional-hazards framework. The random effects
`eta_study_*` are between-study-strata (not between-subject), so a
subject-level Monte Carlo simulation is NOT the correct use of this
model. Instead, treat each "id" in the rxode2 event dataset as one
study strata arm.

# Source trace

Every value in `ini()` is drawn from the on-disk paper + supplement.
The table below lists the source for every equation and every
parameter estimate.

| Component | Source location |
|-----------|-----------------|
| ORR mixed-effects logistic regression form | Franzese 2026 Table 1 ORR row |
| 46 treatment intercepts | Franzese 2026 Table S4 |
| ORR treatment-line effect (-0.697) | Table 1 ORR row |
| ORR PD-1 monotherapy PD-L1 quadratic (1.736) and linear (0.274) effects | Table 1 ORR row; Table S4 |
| ORR PD-L1 monotherapy PD-L1 linear effect (1.642) | Table 1 ORR row; Table S4 |
| ORR any-PD-(L)1-combination PD-L1 linear effect (1.074) | Table 1 ORR row; Table S4 |
| ORR chemotherapy squamous effect (0.282) | Table 1 ORR row; Table S4 |
| ORR strata random effect variance (0.346^2) | Table 1 ORR row |
| OS mixed-effects semi-parametric PH form | Table 1 OS row; Appendix S1 Section 2.3 |
| OS treatment-category ORR-slope effects (5 estimates) | Table 1 OS row; Table 2 |
| OS chemotherapy squamous shift (0.213) | Table 2 |
| OS non-chemotherapy ECOG-PS-0 shift (-0.400) | Table 2 |
| OS Asian-race interaction on ORR slope (-0.595) | Table 2 |
| OS strata random effects (2x2 diagonal, 0.011^2, 0.673^2) | Table 1 OS row |
| OS monthly baseline log(P_m) months 1-60 individually + 61-75 in 3-month groups | Table S2 OS column |
| PFS mixed-effects semi-parametric PH form + PD-(L)1 shift | Table 1 PFS row |
| PFS PD-(L)1 monotherapy hazard-intercept shift (0.410) | Table 1 PFS row; Table 2 |
| PFS treatment-category ORR-slope effects (5 estimates) | Table 1 PFS row; Table 2 |
| PFS chemotherapy squamous shift (0.186) | Table 2 |
| PFS ECOG-PS-0 shift (-0.293) | Table 2 |
| PFS strata random effects (2x2 correlated, 0.199^2, 0.251^2, r=-0.642) | Table 1 PFS row |
| PFS monthly baseline log(P_m) months 1-40 individually + 41-58 in 3-month groups | Table S2 PFS column |
| PFS chemotherapy monthly shift months 2-29 individually + 30-40 grouped (0.805) + 41-56 grouped (1.234) | Table S2 chemotherapy shift column |

# Virtual cohort

We simulate two four-arm cohorts based on Franzese 2026 Figure 2:

* **First-line high-PD-L1** matched population: pembrolizumab
  monotherapy (TRT = 5) versus atezolizumab monotherapy (TRT = 35),
  both at first line, PDL1_TUM = 70% (high), matched squamous /
  ECOG / race fractions.
* **First-line low/negative-PD-L1** matched population: pembrolizumab
  + chemotherapy (TRT = 12) versus atezolizumab + chemotherapy
  (TRT = 36), both at first line, PDL1_TUM = 15% (low), matched
  covariates.

The reference-population covariate values approximate the typical
mNSCLC trial cohort in the analysis dataset (paper Section 3.1 does
not tabulate the pooled covariate medians, so we adopt the following
values consistent with the model's centred reference behaviour):

```{r cohort-refs}
cohort_refs <- tibble::tribble(
  ~cov,               ~value, ~scale,
  "PDL1_TUM (high)",   70L,   "% average PD-L1 expression",
  "PDL1_TUM (low)",    15L,   "% average PD-L1 expression",
  "TUMTP_SQUAM_PCT",   25L,   "% squamous histology",
  "PS_ECOG_0_PCT",     35L,   "% ECOG PS = 0",
  "RACE_ASIAN_PCT",    20L,   "% Asian race"
)
knitr::kable(cohort_refs)
```

We use 200 study-strata "subjects" per arm (each rxode2 "id"
corresponds to one hypothetical study strata arm) so the paper's
between-study random effects (`eta_study_*`) generate an ensemble of
per-arm survival curves.

```{r nsub}
N_PER_ARM <- 200L
```

# Simulation: high-PD-L1 pembrolizumab vs atezolizumab monotherapy

```{r sim-high}
build_cohort <- function(trt, n, pdl1, squam = 25, ecog0 = 35, asian = 20, line_1L = 1L) {
  times <- seq(0, 60, by = 1)
  n_t <- length(times)
  ids <- seq_len(n)
  do.call(rbind, lapply(ids, function(i) {
    data.frame(
      id = rep(i, n_t),
      time = times,
      evid = rep(0L, n_t),
      amt  = rep(0, n_t),
      cmt  = rep("cumhaz_os", n_t),
      TRT = rep(trt, n_t),
      PDL1_TUM = rep(pdl1, n_t),
      LINE_1L = rep(line_1L, n_t),
      TUMTP_SQUAM_PCT = rep(squam, n_t),
      PS_ECOG_0_PCT = rep(ecog0, n_t),
      RACE_ASIAN_PCT = rep(asian, n_t),
      stringsAsFactors = FALSE
    )
  }))
}

# Two-arm high-PD-L1 cohort: TRT = 5 (pembrolizumab), TRT = 35 (atezolizumab).
high_cohort <- rbind(
  build_cohort(trt = 5,  n = N_PER_ARM, pdl1 = 70),
  build_cohort(trt = 35, n = N_PER_ARM, pdl1 = 70)
)
high_cohort$id_ext <- with(high_cohort, ifelse(TRT == 5, id, id + N_PER_ARM))
high_cohort$id <- high_cohort$id_ext

sim_high <- rxode2::rxSolve(mod, high_cohort, returnType = "data.frame") %>%
  as_tibble() %>%
  mutate(arm = ifelse(TRT == 5, "Pembrolizumab", "Atezolizumab"))
```

## ORR distribution (per-arm predicted objective-response fraction)

```{r orr-high, fig.height = 3.5}
sim_high %>%
  filter(time == 0) %>%
  distinct(id, arm, Cc) %>%
  ggplot(aes(x = Cc, fill = arm)) +
  geom_histogram(alpha = 0.6, position = "identity", bins = 40) +
  scale_fill_manual(values = c("Pembrolizumab" = "#2c7fb8", "Atezolizumab" = "#e6550d")) +
  labs(
    x = "Predicted per-arm ORR (proportion)",
    y = "Count of simulated arms",
    fill = "Arm",
    title = "Franzese 2026 Figure 2 reproduction: high-PD-L1 first-line monotherapy",
    subtitle = paste0(
      "PDL1_TUM = 70%; 1L; ", N_PER_ARM, " simulated strata per arm; per-arm eta_study_orr ~ N(0, 0.346^2)"
    )
  ) +
  theme_minimal(base_size = 11)
```

## Median and 95% CI of predicted ORR / OS(24 mo) / PFS(24 mo)

```{r summary-high}
horizons <- c(6, 12, 24, 36, 48)

surv_summary <- sim_high %>%
  filter(time %in% horizons) %>%
  group_by(arm, time) %>%
  summarise(
    S_OS_median  = median(S_os),
    S_OS_q025    = quantile(S_os,  0.025),
    S_OS_q975    = quantile(S_os,  0.975),
    S_PFS_median = median(S_pfs),
    S_PFS_q025   = quantile(S_pfs, 0.025),
    S_PFS_q975   = quantile(S_pfs, 0.975),
    .groups = "drop"
  )
knitr::kable(surv_summary, digits = 3,
             caption = "Franzese 2026 high-PD-L1 first-line monotherapy: simulated S_OS(t) and S_PFS(t) medians and 95% CIs across 200 simulated strata per arm.")
```

## Median survival curves

```{r curves-high}
curve_dat <- sim_high %>%
  group_by(arm, time) %>%
  summarise(
    S_OS_median  = median(S_os),
    S_OS_q025    = quantile(S_os,  0.025),
    S_OS_q975    = quantile(S_os,  0.975),
    S_PFS_median = median(S_pfs),
    S_PFS_q025   = quantile(S_pfs, 0.025),
    S_PFS_q975   = quantile(S_pfs, 0.975),
    .groups = "drop"
  )

p_os <- ggplot(curve_dat, aes(x = time, y = S_OS_median, colour = arm, fill = arm)) +
  geom_ribbon(aes(ymin = S_OS_q025, ymax = S_OS_q975), alpha = 0.15, colour = NA) +
  geom_line(size = 0.9) +
  scale_colour_manual(values = c("Pembrolizumab" = "#2c7fb8", "Atezolizumab" = "#e6550d")) +
  scale_fill_manual  (values = c("Pembrolizumab" = "#2c7fb8", "Atezolizumab" = "#e6550d")) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(x = "Time (months)", y = "S_OS(t)", colour = "Arm", fill = "Arm",
       title = "Overall survival, high-PD-L1 first-line monotherapy",
       subtitle = "Median (line) and 2.5-97.5 percentile ribbon across 200 simulated strata per arm") +
  theme_minimal(base_size = 11)
print(p_os)

p_pfs <- ggplot(curve_dat, aes(x = time, y = S_PFS_median, colour = arm, fill = arm)) +
  geom_ribbon(aes(ymin = S_PFS_q025, ymax = S_PFS_q975), alpha = 0.15, colour = NA) +
  geom_line(size = 0.9) +
  scale_colour_manual(values = c("Pembrolizumab" = "#2c7fb8", "Atezolizumab" = "#e6550d")) +
  scale_fill_manual  (values = c("Pembrolizumab" = "#2c7fb8", "Atezolizumab" = "#e6550d")) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(x = "Time (months)", y = "S_PFS(t)", colour = "Arm", fill = "Arm",
       title = "Progression-free survival, high-PD-L1 first-line monotherapy",
       subtitle = "Median (line) and 2.5-97.5 percentile ribbon across 200 simulated strata per arm") +
  theme_minimal(base_size = 11)
print(p_pfs)
```

# Reproducing Franzese 2026 Table 1 / Table 2 head-to-head HR

Franzese 2026 Figure 2 reports pembrolizumab-vs-atezolizumab summary
estimates for a hypothetical first-line high-PD-L1 population:

| Endpoint | Reported estimate (95% CI) | Simulation |
|----------|---------------------------|-----------|
| ORR OR (pembro vs atezo) | 1.39 (0.94-2.05) | see below |
| OS HR (pembro vs atezo) | 0.87 (0.69-1.02) | see below |
| PFS HR (pembro vs atezo) | 0.83 (0.67-1.03) | see below |

Compute the ORR odds ratio from typical predictions (zeroRe):

```{r or-hr-typical}
typ_dat <- rbind(
  build_cohort(trt = 5,  n = 1, pdl1 = 70),
  build_cohort(trt = 35, n = 1, pdl1 = 70)
)
typ_dat$id <- c(1, 2)

typ <- rxode2::rxSolve(rxode2::zeroRe(mod), typ_dat, returnType = "data.frame") %>%
  as_tibble()

typ_final <- typ %>%
  filter(time == 0) %>%
  distinct(id, TRT, Cc) %>%
  mutate(arm = ifelse(TRT == 5, "pembrolizumab", "atezolizumab"))

typ_final
odds <- setNames(typ_final$Cc / (1 - typ_final$Cc), typ_final$arm)
odds
or_typical <- odds[["pembrolizumab"]] / odds[["atezolizumab"]]
message(sprintf("Typical-value ORR odds ratio pembro/atezo = %.3f (Franzese 2026 Figure 2: 1.39)", or_typical))
```

The typical-value ORR-OR estimate above is the point estimate; the
95% CI reported by Franzese 2026 (0.94-2.05) captures uncertainty in
the fixed-effect parameters (Table S4 intercept RSEs) which is not
propagated in the packaged model (all parameters are `fixed()`).
Reproducing the paper's 95% CI would require simulating from a
joint sampling distribution over the treatment intercepts and the
PD-L1 covariate coefficients (paper Section 2.4 uses 10,000 trial
simulations); the packaged model is sufficient for the typical-value
point estimate.

Compute the OS and PFS HR from typical predictions at t = 24 months
via `HR = log(S_2(t)) / log(S_1(t))`:

```{r hr-typical}
typ_horizon <- typ %>%
  filter(time == 24) %>%
  distinct(id, TRT, S_os, S_pfs) %>%
  mutate(arm = ifelse(TRT == 5, "pembrolizumab", "atezolizumab"))

typ_horizon
hr_os  <- log(typ_horizon$S_os [typ_horizon$arm == "pembrolizumab"]) /
          log(typ_horizon$S_os [typ_horizon$arm == "atezolizumab"])
hr_pfs <- log(typ_horizon$S_pfs[typ_horizon$arm == "pembrolizumab"]) /
          log(typ_horizon$S_pfs[typ_horizon$arm == "atezolizumab"])
message(sprintf("Typical-value OS HR pembro/atezo at 24 mo = %.3f (paper Figure 2: 0.87)", hr_os))
message(sprintf("Typical-value PFS HR pembro/atezo at 24 mo = %.3f (paper Figure 2: 0.83)", hr_pfs))
```

The typical-value HR calculations above use a single-time-point ratio;
because the paper's proportional-hazards model imposes a constant
log-HR across the reference baseline hazard, the HR is
time-independent at typical values (the ratio would come out the same
at any horizon).

# Simulation: low-PD-L1 pembrolizumab + chemo vs atezolizumab + chemo

The second published head-to-head compares chemotherapy combinations
in a low-PD-L1 first-line population:

```{r sim-lowchemo, fig.height = 3.5}
low_cohort <- rbind(
  build_cohort(trt = 12, n = N_PER_ARM, pdl1 = 15),
  build_cohort(trt = 36, n = N_PER_ARM, pdl1 = 15)
)
low_cohort$id_ext <- with(low_cohort, ifelse(TRT == 12, id, id + N_PER_ARM))
low_cohort$id <- low_cohort$id_ext

sim_low <- rxode2::rxSolve(mod, low_cohort, returnType = "data.frame") %>%
  as_tibble() %>%
  mutate(arm = ifelse(TRT == 12, "Pembro + chemo", "Atezo + chemo"))

curve_low <- sim_low %>%
  group_by(arm, time) %>%
  summarise(
    S_OS_median  = median(S_os),
    S_PFS_median = median(S_pfs),
    .groups = "drop"
  )

p_low <- curve_low %>%
  pivot_longer(cols = c(S_OS_median, S_PFS_median), names_to = "endpoint",
               values_to = "S") %>%
  mutate(endpoint = ifelse(endpoint == "S_OS_median", "OS", "PFS")) %>%
  ggplot(aes(x = time, y = S, colour = arm, linetype = endpoint)) +
  geom_line(size = 0.9) +
  scale_colour_manual(values = c("Pembro + chemo" = "#2c7fb8", "Atezo + chemo" = "#e6550d")) +
  coord_cartesian(ylim = c(0, 1)) +
  labs(x = "Time (months)", y = "S(t)", colour = "Arm", linetype = "Endpoint",
       title = "Low-PD-L1 first-line +chemotherapy: pembrolizumab vs atezolizumab",
       subtitle = "Median across 200 simulated strata per arm; median S_OS(t) and S_PFS(t)") +
  theme_minimal(base_size = 11)
print(p_low)
```

Franzese 2026 Figure 2 reports for this comparison:

| Endpoint | Reported estimate (95% CI) | Simulation |
|----------|---------------------------|-----------|
| ORR OR (pembro+chemo vs atezo+chemo) | 1.08 (0.81-1.44) | see below |
| OS HR (pembro+chemo vs atezo+chemo) | 0.97 (0.86-1.08) | see below |
| PFS HR (pembro+chemo vs atezo+chemo) | 0.98 (0.89-1.07) | see below |

```{r or-hr-lowchemo}
typ_dat_low <- rbind(
  build_cohort(trt = 12, n = 1, pdl1 = 15),
  build_cohort(trt = 36, n = 1, pdl1 = 15)
)
typ_dat_low$id <- c(1, 2)

typ_low <- rxode2::rxSolve(rxode2::zeroRe(mod), typ_dat_low, returnType = "data.frame") %>%
  as_tibble()

typ_low_orr <- typ_low %>% filter(time == 0) %>% distinct(id, TRT, Cc)
odds_low <- setNames(typ_low_orr$Cc / (1 - typ_low_orr$Cc),
                     ifelse(typ_low_orr$TRT == 12, "pembro_chemo", "atezo_chemo"))
odds_low
or_low <- odds_low[["pembro_chemo"]] / odds_low[["atezo_chemo"]]
message(sprintf("Typical-value low-PD-L1 +chemo ORR OR = %.3f (paper Figure 2: 1.08)", or_low))

typ_low_24 <- typ_low %>% filter(time == 24) %>% distinct(id, TRT, S_os, S_pfs)
hr_os_low  <- log(typ_low_24$S_os [typ_low_24$TRT == 12]) /
              log(typ_low_24$S_os [typ_low_24$TRT == 36])
hr_pfs_low <- log(typ_low_24$S_pfs[typ_low_24$TRT == 12]) /
              log(typ_low_24$S_pfs[typ_low_24$TRT == 36])
message(sprintf("Typical-value low-PD-L1 +chemo OS HR at 24 mo = %.3f (paper Figure 2: 0.97)", hr_os_low))
message(sprintf("Typical-value low-PD-L1 +chemo PFS HR at 24 mo = %.3f (paper Figure 2: 0.98)", hr_pfs_low))
```

# Verifying the monthly baseline hazard against Table S2

The reference OS survival curve at typical values with all covariates
= 0 and all etas = 0 is `S_0(t) = product from m = 1 to floor(t) of
(1 - exp(P_m))`. We can verify the model reproduces the tabulated
Table S2 baseline hazards by solving with zeroed covariates:

```{r baseline-verify}
# The paper's baseline hazard represents "0% ORR" and a reference arm; setting
# ORR to 0 requires zeroing the treatment intercepts too. Because ORR enters
# multiplicatively in the log(HR) equations, the simplest reproduction is to
# select the paper's chemotherapy-only baseline (TRT = 1) at LINE_1L = 1 and
# all covariates zero -- the OS log(HR) then reduces to
# (eta - 1.651) * orr_pred (with orr_pred = 0 at typical value if we could
# force logit_orr = -inf). In practice we use the typical value at the
# chemotherapy reference and interpret the resulting S_0(t) as the arm's
# baseline curve modulo the small residual HR from the non-zero orr_pred.
base_dat <- data.frame(
  id = 1L,
  time = seq(0, 24, by = 1),
  evid = 0L,
  amt = 0,
  cmt = "cumhaz_os",
  TRT = 1L, PDL1_TUM = 0, LINE_1L = 0L,
  TUMTP_SQUAM_PCT = 0, PS_ECOG_0_PCT = 0, RACE_ASIAN_PCT = 0
)
base_sim <- rxode2::rxSolve(rxode2::zeroRe(mod), base_dat, returnType = "data.frame") %>%
  as_tibble()

# Paper reference: S_0(t = 1) = 1 - exp(-2.743) = 0.9356; S_0(t = 2) = S_0(t = 1) * (1 - exp(-2.502)) = 0.9356 * 0.9181 = 0.8590
tab_s2 <- tibble::tribble(
  ~month, ~logP_os,
  1,  -2.743,
  2,  -2.502,
  3,  -2.511,
  4,  -2.449,
  5,  -2.379,
  6,  -2.402,
  12, -2.426,
  24, -2.814
)
tab_s2$P_m <- exp(tab_s2$logP_os)
tab_s2$expected_S0_monthly <- 1 - tab_s2$P_m

knitr::kable(tab_s2, digits = 4,
             caption = "Franzese 2026 Table S2 (OS): monthly conditional death probability at each month; expected S_0 contribution 1 - P_m.")
```

# Assumptions and deviations

* **All parameters wrapped in `fixed()`.** The model is a downstream
  consumer of the published fit, not a re-estimation of it. Users
  reproducing the paper's 10,000-trial simulations (Section 2.4) that
  propagate parameter-estimation uncertainty need to inject the
  Table 1 / Table 2 %RSE covariance externally; the packaged model
  reproduces the point-estimate predictions.
* **`fixed()` on both diagonal and off-diagonal random effects.** The
  three random-effect blocks (ORR, OS 2x2 diagonal, PFS 2x2 correlated)
  encode the reported variances / covariance verbatim.
* **`eta_study_*` naming.** These are between-study-strata random
  effects (MBMA convention), not between-subject. They are declared
  via the model's `paper_specific_etas` metadata field so
  `checkModelConventions()` does not flag them as missing a 1-to-1
  fixed-effect pairing (the study-strata random effect naturally sits
  on the treatment-category ORR-slope terms, not on a single
  structural theta).
* **Line-of-therapy sign flip.** The paper encodes the effect as
  `-0.697 * (1 if >=2L)`; the canonical `LINE_1L` (1 = first-line,
  0 = >=2L) reverses the reference category, so the packaged model
  encodes it as `+0.697 * LINE_1L`. Effect direction is preserved
  (first-line arms have +0.697 log-odds relative to the >=2L
  reference; equivalently, >=2L arms have -0.697 relative to the 1L
  reference).
* **Monthly discrete-time baseline hazard.** The paper's semi-parametric
  proportional-hazards model uses monthly discrete conditional event
  probabilities as the reference baseline hazard. The packaged model
  implements this via a piecewise-constant hazard function h(t)
  reading from 65 nested `ifelse` selectors (OS: months 1-60
  individual + 5 grouped 3-month tails 61-63, 64-66, 67-69, 70-72,
  73-75; PFS: months 1-40 individual + 6 grouped 3-month tails
  41-43, 44-46, 47-49, 50-52, 53-55, 56-58; PFS chemotherapy shift:
  months 2-29 individual + 2 pooled 0.805 for months 30-40 and 1.234
  for months 41-56). The integration `d/dt(cumhaz_*) = h(t) * exp(log_HR)`
  reproduces the paper's `S(t) = S_0(t)^exp(log_HR)` proportional-hazards
  formulation.
* **Simulation scope: study-strata mean only.** Each rxode2 "id" is
  one hypothetical study strata arm, not one individual patient. The
  model is NOT suitable for individual-subject Monte Carlo simulation.
  Uses supporting the packaged model: reproducing published
  head-to-head trial simulations (Figures 1-5); benchmarking emerging
  Phase 1/2 ORR data against the historical PD-(L)1 monotherapy /
  combination landscape; predicting OS / PFS from a stated per-arm
  ORR.
* **Beyond the paper's time support.** Model beyond OS month 75 or
  PFS month 58 holds the last grouped estimate constant per the
  paper's estimator pooling. Extrapolation is nonetheless outside
  the supported time-window; use with caution.
* **No PKNCA validation.** This model has no PK compartment (input
  is a treatment integer + arm-level covariates, not a dose), so
  the standard PKNCA cross-check does not apply. Validation focuses
  on reproducing the published head-to-head simulation summary
  statistics from Figure 2.
* **Nominal residual SD.** `addSd = 0.01` on `Cc`; the paper's actual
  residual structure is the normal-approximation binomial weighting
  (Equations 3 and 5) which does not produce a single scalar SD.
  Users doing forward simulation should either set etas to zero
  (deterministic per-arm prediction) or rely on the `eta_study_orr`
  random effect to generate an ensemble of arm-level ORR values.

# Errata

* The paper's Table 1 typesets the OS and PFS random-effect matrices
  with the correlation coefficient (r equals -0.642) shown as
  "r =- 0.642". We interpret this as (r = -0.642) and encode the
  covariance as (-0.642 times 0.199 times 0.251 equals -0.032063).
* Franzese 2026 Table S4 lists 46 unique treatments. The packaged
  model encodes all 46 verbatim; a treatment absent from the paper's
  dataset (e.g., a newly-approved anti-PD-(L)1 combination) is not
  supported by any of the 1-46 TRT integers and would require adding a
  new intercept extending the register.
