Survival (time-to-event) outcomes

The unanchored survival problem

In an unanchored comparison of two single-arm studies we have individual patient data (IPD) for the index treatment and only published, aggregate evidence for the comparator. For survival endpoints that evidence is almost always a Kaplan-Meier (KM) curve plus a baseline-characteristics table.

mlumr handles this exactly as for its other families: the comparator arm is represented by reconstructed pseudo-IPD, event/censoring times digitized from the published KM curve (e.g. with the Guyot algorithm), and the model integrates the individual-level survival likelihood over the comparator’s covariate distribution, transporting the index-derived survival model to the comparator (or any) population. The unanchored setting removes the need for a common arm but still relies on conditional constancy / the shared-prognostic-factor assumption (SPFA).

The ML-UMR survival model

For a proportional-hazards distribution ML-UMR models the individual hazard \[h(t\mid x)=h_0(t)\exp(\mu_k+x^\top\beta_k),\] with a parametric or M-spline baseline \(h_0(t)\); accelerated-failure-time distributions instead scale time by \(\exp(\mu_k+x^\top\beta_k)\). The index IPD enters through the usual survival likelihood; the comparator’s reconstructed pseudo-IPD contributes a likelihood integrated over its covariate distribution \(f_{\text{AgD}}(x)\), \[\ell_{\text{AgD}}=\sum_j\log\!\int h(t_j\mid x)^{\delta_j}\,S(t_j\mid x)\,f_{\text{AgD}}(x)\,dx,\] approximated by the quasi-Monte Carlo integration points. Population-standardized quantities follow from the marginal survival \(\bar S(t)=\mathbb{E}_x[S(t\mid x)]\) and marginal hazard \(\bar h(t)=\mathbb{E}_x[h(t\mid x)S(t\mid x)]/\mathbb{E}_x[S(t\mid x)]\). Because hazard ratios are non-collapsible, the marginal log HR is time-varying even under conditional PH, so we report it as a curve (predict(type = "loghr")) and lean on the collapsible RMST estimands for scalar summaries.

Survival distributions

The distribution argument of mlumr() selects the survival model, nine parametric forms (proportional hazards, PH; or accelerated failure time, AFT) plus two flexible baselines:

marginal_effects() reports these on the natural scale, null 1, as for the Poisson rate ratio. For AFT distributions a population time ratio requires both a shared baseline shape and a shared coefficient vector; under the default aux_by = ".study", or under the relaxed model, marginal_effects() reports EXP_DELTA_ETA instead and an explicit request for "tr" is an error. predict(type = "loghr") gives the time-varying marginal log hazard ratio curve (null 0). The "gengamma" form is the generalized gamma restricted to positive Lawless shape Q (Q = 1 / sqrt(aux2) > 0); it nests the Weibull, gamma, and (as a limit) log-normal, but does not cover negative-Q shapes (use an "mspline" baseline for a hazard shape outside that family).

The flexible baselines estimate the baseline-hazard shape from the data using a normalized M-spline (the piecewise exponential is the degree-0 case) with a random-walk smoothing prior (Phillippo et al. 2026). Knots are placed by make_knots() (control with n_knots).

Censoring

Outcomes are given to set_ipd() / set_agd_surv() either as a survival::Surv() object or as time / status / entry_time columns. All standard schemes are supported (internally mapped to status 0 right, 1 event, 2 left, 3 interval; delayed entry handled as left truncation):

library(survival)
Surv(time, status)                      # right censoring (status 1 = event)
Surv(entry_time, time, status)          # delayed entry / left truncation
Surv(lower, upper, type = "interval2")  # interval censoring

The data: progression-free survival in NDMM

• Index (IPD). Patient-level progression-free survival for lenalidomide (Len) maintenance from the McCarthy2012 study.
• Comparator (AgD). A reconstructed KM curve plus covariate moments for thalidomide (Thal) from the Morgan2012 study.

We adjust for age, ISS stage III, complete/very-good-partial response, and sex.

data("ndmm_ipd")        # bundled with mlumr (copied from multinma, GPL-3)
data("ndmm_agd")        # reconstructed pseudo-IPD (KM)
data("ndmm_agd_covs")   # covariate moments for the AgD arm

covs <- c("age", "iss_stage3", "response_cr_vgpr", "male")

# Index IPD: McCarthy2012, lenalidomide. Age stays in years, as in multinma's
# own NDMM example, so the coefficient reads as a per-year log hazard ratio.
ipd_one <- ndmm_ipd[ndmm_ipd$study == "McCarthy2012" & ndmm_ipd$treatment == "Len", ]

# Comparator AgD: Morgan2012, thalidomide (pseudo-IPD + covariate moments).
agd_one <- ndmm_agd[ndmm_agd$study == "Morgan2012" & ndmm_agd$treatment == "Thal", ]
cv <- ndmm_agd_covs[ndmm_agd_covs$study == "Morgan2012" & ndmm_agd_covs$treatment == "Thal", ]
agd_one$age_mean         <- cv$age_mean
agd_one$age_sd           <- cv$age_sd
agd_one$iss_stage3_prop       <- cv$iss_stage3_prop        # proportions (binary covariates)
agd_one$response_cr_vgpr_prop <- cv$response_cr_vgpr_prop
agd_one$male_prop             <- cv$male_prop

Covariate balance

Population adjustment is motivated by the fact that the two single-arm populations differ on prognostic factors. Comparing the index sample means against the published comparator means shows the imbalance we must adjust for:

balance <- data.frame(
  Covariate = c("Age (years)", "ISS stage III (prop.)",
                "CR/VGPR response (prop.)", "Male (prop.)"),
  Index_Len = c(mean(ipd_one$age), mean(ipd_one$iss_stage3),
                mean(ipd_one$response_cr_vgpr), mean(ipd_one$male)),
  Comparator_Thal = c(agd_one$age_mean[1], agd_one$iss_stage3_prop[1],
                      agd_one$response_cr_vgpr_prop[1], agd_one$male_prop[1])
)
knitr::kable(balance, caption = "Prognostic-factor balance: McCarthy2012 Len vs Morgan2012 Thal")

The observed Kaplan-Meier curves (index IPD and the reconstructed comparator pseudo-IPD) show the crude, unadjusted survival difference:

km_obs <- rbind(
  data.frame(eventtime = ipd_one$eventtime, status = ipd_one$status,
             arm = "Lenalidomide"),
  data.frame(eventtime = agd_one$eventtime, status = agd_one$status,
             arm = "Thalidomide")
)
# The displayed observed curves use ggsurvfit: step lines with censoring marks
# and a numbers-at-risk table, the standard survival-reporting layout. (The
# predicted-vs-observed overlay further down draws its KM with mlumr's geom_km().)
arm_cols <- c("Lenalidomide" = "#00BA38",   # green (index)
              "Thalidomide" = "#F8766D")    # red (comparator)
ggsurvfit::survfit2(survival::Surv(eventtime, status) ~ arm, data = km_obs) |>
  ggsurvfit::ggsurvfit(linewidth = 0.7) +
  ggsurvfit::add_censor_mark(shape = 3, size = 1.6, alpha = 0.7) +
  ggsurvfit::add_risktable(risktable_stats = "n.risk", size = 3.2) +
  scale_color_manual(values = arm_cols) +
  scale_fill_manual(values = arm_cols) +
  scale_x_continuous(breaks = seq(0, 60, 10)) +
  scale_y_continuous(
    breaks = seq(0, 1, 0.1),                       # gridline at every 0.1
    labels = function(b) ifelse(round(b * 10) %% 2 == 0,  # label 0, 0.2, .., 1
                                sprintf("%.1f", b), ""),
    limits = c(0, 1), expand = expansion(mult = c(0, 0.02))) +
  labs(x = "Progression-free survival (months)", y = "Survival probability",
       title = "Observed Kaplan-Meier curves") +
  theme_minimal(base_size = 11) + theme(legend.position = "bottom")

plot of chunk km

Setting up the ML-UMR data

The index IPD uses family = "survival" with time/status columns; the comparator uses set_agd_surv(), which takes the reconstructed pseudo-IPD plus the covariate moments from the baseline table.

ipd <- set_ipd(ipd_one, treatment = "treatment", covariates = covs,
               family = "survival", time = "eventtime", status = "status")

agd <- set_agd_surv(agd_one, treatment = "treatment",
                    time = "eventtime", status = "status",
                    cov_means = c("age_mean", "iss_stage3_prop", "response_cr_vgpr_prop", "male_prop"),
                    cov_sds   = c("age_sd", NA, NA, NA),
                    cov_types = c("continuous", "binary", "binary", "binary"))

dat <- suppressWarnings(combine_data(ipd, agd))
dat <- suppressWarnings(add_integration(
  dat, n_int = 64,
  # Age is right-skewed, so a gamma marginal, matching multinma's NDMM example.
  age              = distr(qgamma, mean = age_mean, sd = age_sd),
  iss_stage3       = distr(qbern, prob = iss_stage3_mean),
  response_cr_vgpr = distr(qbern, prob = response_cr_vgpr_mean),
  male             = distr(qbern, prob = male_mean)))
dat
#> Unanchored Comparison Data (Time-to-event)
#> ====================================
#> 
#> Index treatment (IPD): Len 
#>   N = 231 
#>   Events = 133 (57.6%), censored = 98
#> 
#> Comparator treatment (AgD): Thal 
#>   Reconstructed pseudo-IPD: 408 row(s)
#>   Events = 263 (64.5%), censored = 145
#> 
#> Covariates ( 4 ): age, iss_stage3, response_cr_vgpr, male 
#> Integration points: 64 (QMC with Gaussian copula)

Frequentist benchmarks

The naive comparison is an unadjusted Cox model (a log hazard ratio), ignoring the covariate imbalance. Survival STC fits a parametric model to the index IPD and transports it to the comparator population by G-computation, returning a restricted-mean-survival-time (RMST) difference via flexsurv. Both are fast frequentist reference points for the Bayesian fit below (the two estimands, log HR and RMST difference, are on different scales, so each is compared like for like later).

RMST is an integral to a restriction time, so two RMST differences computed to different horizons are different estimands and must not be placed side by side. The defaults do not agree by construction: a stratified flexible baseline defaults to the follow-up both studies observed, because each study’s baseline is extrapolated as a constant hazard past its own last event, whereas a parametric fit and stc() default to the pooled maximum. We therefore fix one horizon here and pass it to every estimator below.

# The follow-up both arms in THIS analysis actually observed, taken from the two
# subsets that were combined above rather than from the full bundled datasets.
# Every RMST estimate in this vignette is integrated to this time, so the forest
# plot compares like with like.
tau_rmst <- min(max(ipd_one$eventtime), max(agd_one$eventtime))
tau_rmst
#> [1] 64.8

res_naive <- naive(dat)                            # unadjusted Cox log HR
res_stc   <- stc(dat, distribution = "weibull", seed = 2026,
                 rmst_horizon = tau_rmst)          # G-computation RMST diff
res_naive
#> Naive Unadjusted Indirect Comparison
#> =====================================
#> 
#> Treatments: Len vs Thal 
#> 
#> Population basis: index-study outcome versus comparator-population outcome; no common standardized target.
#> 
#> Median survival:
#>   Index (IPD):      54.526 (133 events / 231)
#>   Comparator (AgD): 25.499 (263 events / 408)
#> 
#> Log Hazard Ratio (Cox): -0.5541 (SE: 0.1087)
#> 95% CI: [-0.7672, -0.3410]
#> 
#> All effect measures (95% CI):
#>   Log hazard ratio (Cox)        -0.5541 (SE 0.1087) [-0.7672, -0.3410]
#>   Hazard ratio (Cox)             0.5746 [0.4643, 0.7111]
res_stc
#> Simulated Treatment Comparison (G-computation)
#> ===============================================
#> 
#> Treatments: Len vs Thal 
#> 
#> Estimand population: comparator
#> Treating this as the index-population effect requires a separate effect-equality assumption; this calculation does not transport to the index population.
#> 
#> Method note: package-specific parametric survival extension.
#> Distribution: weibull | RMST horizon: 64.800
#> Marginalized RMST (index trt, comp pop): 36.8432
#> Fitted RMST (comp trt, comp pop):        33.6594
#> 
#> RMST Difference: 3.1837 (bootstrap SE: 2.6586, 200 reps)
#> 95% CI: [-2.0271, 8.3946]
#> 
#> All effect measures (95% CI):
#>   RMST difference                3.1837 (SE 2.6586) [-2.0271, 8.3946]
#>   Log cumulative-hazard ratio (at horizon)   0.0461 (SE 0.1222) [-0.1935, 0.2857]
#>   Cumulative-hazard ratio (at horizon)   1.0472 [0.8241, 1.3306]

Fitting the ML-UMR models

Following multinma’s NDMM survival example, the primary model uses a cubic M-spline baseline hazard (distribution = "mspline", mspline_degree = 3 by default) with seven internal knots placed at event-time quantiles, adjusting for all four prognostic factors. We fit the shared-effect (SPFA) model and, mirroring multinma’s ~ (covariates) * .trt interaction specification, a relaxed model that lets the prognostic effects differ by treatment (effect modification). A Weibull proportional-hazards model is fit alongside as a parametric comparison. The settings (64 integration points, 4 chains of 2000 iterations, cubic M-spline with 7 internal knots) match multinma’s NDMM example. prior_aux (shape/scale) and prior_smooth (M-spline smoothing SD) supplement the usual priors; see vignette("fitting-and-diagnostics").

The baseline hazard is estimated separately for each study (aux_by = ".study", the default), as in multinma, where .study is always part of the stratification. Two single-arm trials rarely share a hazard shape, and forcing one on them would impose proportional hazards across studies, which nothing in an unanchored design supports; a shared shape is available as aux_by = "none" when the two Kaplan-Meier curves plainly agree. (As in multinma, aux_by = NULL also means a baseline per study, so "none" is a mlumr-specific spelling for a model multinma cannot fit.) Identification survives because each study gets its own knots over its own observed support (as in multinma’s default type = "quantile") and its spline coefficients are a simplex, which pins that study’s cumulative hazard to 1 at a time the study actually observed. Both halves matter: a single pooled basis spanning the longer study would leave the shorter one with basis functions it never observes, and the likelihood would then be exactly flat along a direction that trades the spline scale against the study intercept.

One consequence matters for reading everything below. The marginal hazard ratio is time-varying in any case, because it weights the covariate distribution by each arm’s own survival and the two risk sets diverge; with two baseline shapes it additionally picks up the ratio of the baselines. So the scalar hazard ratio is only its value at one instant, and marginal_effects() records which one in an at_time column: the first prediction time here, since the prediction grid excludes zero. Changing pred_times moves that number, which is why it should never be read as “the” hazard ratio. predict(type = "loghr") gives the whole curve, and the RMST difference is the estimand that does not depend on when you look.

# Primary: cubic M-spline baseline (the multinma NDMM survival model)
fit_mspline <- mlumr(dat, model = "spfa", distribution = "mspline", n_knots = 7,
                     rmst_horizon = tau_rmst,
                     chains = 4, iter = 2000, warmup = 1000, seed = 2026, refresh = 0)
#> Running MCMC with 4 parallel chains...
#> 
#> Chain 3 finished in 1026.5 seconds.
#> Chain 4 finished in 1106.9 seconds.
#> Chain 1 finished in 1214.1 seconds.
#> Chain 2 finished in 1256.3 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 1151.0 seconds.
#> Total execution time: 1256.7 seconds.

# Relaxed M-spline: treatment-specific prognostic effects (effect modification),
# the unanchored analogue of multinma's ~ (covariates) * .trt. beta_comparator is
# identified only by the reconstructed comparator likelihood, so we regularize it
# with a tighter autoscaled prior (see the caveat at the end of this vignette).
fit_mspline_relaxed <- mlumr(dat, model = "relaxed", distribution = "mspline",
                             n_knots = 7, rmst_horizon = tau_rmst,
                             prior_beta_comparator = prior_normal(0, 1, autoscale = TRUE),
                             chains = 4, iter = 2000, warmup = 1000, seed = 2026,
                             refresh = 0)
#> Running MCMC with 4 parallel chains...
#> 
#> Chain 4 finished in 1472.6 seconds.
#> Chain 3 finished in 1492.9 seconds.
#> Chain 1 finished in 1525.4 seconds.
#> Chain 2 finished in 1547.0 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 1509.5 seconds.
#> Total execution time: 1547.7 seconds.

# Parametric comparison: Weibull proportional hazards
fit_weibull <- mlumr(dat, model = "spfa", distribution = "weibull",
                     rmst_horizon = tau_rmst,
                     chains = 4, iter = 2000, warmup = 1000, seed = 2026, refresh = 0)
#> Running MCMC with 4 parallel chains...
#> 
#> Chain 2 finished in 1346.0 seconds.
#> Chain 1 finished in 1411.4 seconds.
#> Chain 3 finished in 1413.4 seconds.
#> Chain 4 finished in 1478.4 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 1412.3 seconds.
#> Total execution time: 1479.3 seconds.

summary(fit_mspline)
#> ML-UMR Model Summary
#> ====================
#> 
#> Model: SPFA 
#> Family: Time-to-event [mspline] 
#> Link: log 
#> Engine: cmdstanr 
#> Treatments: Len (IPD) vs Thal (AgD)
#> 
#> MCMC Diagnostics:
#>   Divergent transitions: 0 
#>   Max treedepth hits: 0 
#>   Max Rhat: 1.003  
#>   Min ESS: 1295  
#> 
#> Intercepts (log scale):
#>       variable      mean         sd        2.5%     97.5%     Rhat
#>       mu_index 0.3278324 0.11767435  0.09338423 0.5571708 1.000259
#>  mu_comparator 0.1072612 0.08812043 -0.07031650 0.2766931 1.001350
#> 
#> Regression Coefficients:
#>                variable       mean         sd        2.5%     97.5%      Rhat
#>               beta[age] 0.06993144 0.01627319  0.03846601 0.1015201 1.0006454
#>        beta[iss_stage3] 0.40244365 0.18953509  0.01672759 0.7567338 0.9998211
#>  beta[response_cr_vgpr] 0.06266713 0.18381897 -0.30799778 0.4183324 1.0005757
#>              beta[male] 0.07541778 0.18129370 -0.27756031 0.4240964 1.0024914
#> 
#> Baseline hazard: stratified by study (aux_by = ".study")
#> 
#> Shape / smoothing parameters:
#>         variable      mean        sd      2.5%    97.5%     Rhat
#>  sigma_smooth[1] 0.9861364 0.4947026 0.2604557 2.157619 1.000875
#>  sigma_smooth[2] 1.2443250 0.4832417 0.4561599 2.348745 1.000971
#> 
#> Marginal Treatment Effects:
#>   Log Hazard Ratios at t = 1.4:
#>          variable       mean        sd      2.5%     97.5%
#>       delta_index -0.4853513 0.4247787 -1.389658 0.2756818
#>  delta_comparator -0.4825707 0.4218878 -1.383165 0.2693954
#>   RMST Differences (to t = 64.8):
#>              variable     mean       sd      2.5%    97.5%
#>       rmst_diff_index 3.014535 2.433137 -1.665791 7.953528
#>  rmst_diff_comparator 3.846938 2.635502 -1.282804 9.223502

Priors

Beyond the usual intercept/coefficient priors, survival models add prior_aux (shape/scale parameters) and, for flexible baselines, prior_smooth (the M-spline random-walk smoothing SD). prior_summary() shows exactly what was passed to Stan:

prior_summary(fit_mspline)
#> Priors for ML-UMR Fit
#> =====================
#> 
#> Intercepts (mu_index, mu_comparator):
#>   normal(0, 10)
#>   (package default, mlumr 0.2.0)
#> 
#> Regression coefficients (beta):
#>   normal(0, 2.5) applied to all 4 covariate(s)
#>   (package default, mlumr 0.2.0)
#> 
#> Survival baseline smoothing (random-walk SD, half-normal via <lower=0>):
#>   normal(0, 1)
#>   (package default, mlumr 0.2.0)

The posterior intercepts (log-baseline scale) versus their \(\mathrm{N}(0,10)\) prior:

plot_prior_posterior(fit_mspline, pars = c("mu_index", "mu_comparator"))

plot of chunk prior-post

data.frame(
  n_divergent  = fit_mspline$diagnostics$n_divergent,
  max_treedepth = fit_mspline$diagnostics$n_max_treedepth,
  max_Rhat     = round(max(fit_mspline$summary$Rhat, na.rm = TRUE), 3),
  min_ESS      = round(min(fit_mspline$summary$n_eff, na.rm = TRUE))
)
#>   n_divergent max_treedepth max_Rhat min_ESS
#> 1           0             0    1.003    1295

Treatment effects

marginal_effects() reports, in each population, the marginal hazard ratio (HR, natural scale, null 1), the RMST difference (RMSTD, in months of follow-up, null 0), and the RMST ratio (RMSTR, natural scale, null 1). For accelerated failure time distributions the scalar is a time ratio (TR) only when the two studies share one baseline shape and one coefficient vector; otherwise it is the exponentiated linear-predictor contrast (EXP_DELTA_ETA). The hazard ratio is exponentiated from the log hazard ratio, while predict(type = "loghr") remains on the log scale. Each RMST row carries the restriction horizon it was integrated to, since RMST at a different horizon is a different estimand:

knitr::kable(marginal_effects(fit_mspline, effect = "all"),
   caption = "M-spline population-standardized effects (Len vs Thal)")

Because hazard ratios are non-collapsible, the marginal hazard ratio is time-varying even under conditional PH. Plot predict(type = "loghr") directly; its plot() method draws the credible band, facets by population, and adds the null line at 0:

plot(predict(fit_mspline, type = "loghr"))

plot of chunk loghr

A forest plot of the collapsible RMST difference. The Bayesian estimate is shown in the index population, which is the decision-relevant target for HTA: cost- effectiveness models are built for the population a reimbursement decision is about, and that is normally the index trial’s, not the comparator’s (Conor Chandler and Ishak 2026). The STC benchmark is a comparator-population estimand by construction, so ML-UMR is shown in both populations: the comparator row is the like-for-like comparison with STC, and the index row is the decision-relevant estimand. (The naive Cox estimate is a log hazard ratio, a different scale, so it is reported under Frequentist benchmarks above rather than on this axis.)

# Every ML-UMR fit, in both target populations, on the collapsible RMST scale.
# Each row must be integrated to the SAME restriction time or the plot compares
# different estimands, so check that before drawing rather than assuming it.
rd <- function(fit) marginal_effects(fit, effect = "rmstd", population = "both")
pop <- function(d, which) d[d$population == which, ]
rd_ms  <- rd(fit_mspline)
rd_rel <- rd(fit_mspline_relaxed)
rd_wb  <- rd(fit_weibull)
stopifnot(all.equal(unique(c(rd_ms$horizon, rd_rel$horizon, rd_wb$horizon,
                             res_stc$horizon)), tau_rmst))
row <- function(d, which) pop(d, which)[, c("mean", "q2.5", "q97.5")]
bayes <- rbind(row(rd_ms, "Index"),  row(rd_ms, "Comparator"),
               row(rd_rel, "Index"), row(rd_rel, "Comparator"),
               row(rd_wb, "Index"),  row(rd_wb, "Comparator"))
forest_df <- data.frame(
  label = c("STC (flexsurv, comparator pop.)",
            "ML-UMR M-spline SPFA (index pop.)",
            "ML-UMR M-spline SPFA (comparator pop.)",
            "ML-UMR M-spline relaxed (index pop.)",
            "ML-UMR M-spline relaxed (comparator pop.)",
            "ML-UMR Weibull SPFA (index pop.)",
            "ML-UMR Weibull SPFA (comparator pop.)"),
  est = c(res_stc$rmst_diff, bayes$mean),
  lo  = c(res_stc$ci_lower,  bayes$q2.5),
  hi  = c(res_stc$ci_upper,  bayes$q97.5)
)
mlumr_forest(forest_df, ref_line = 0,
             x = "RMST difference (months)",
             title = "Restricted mean survival time difference",
             subtitle = sprintf(paste0("Every ML-UMR fit and the STC benchmark,",
                                       " integrated to tau = %.4g months"),
                                tau_rmst))

plot of chunk forest

Absolute predictions

Absolute predicted survival summaries (RMST and median) for each treatment, then the standardized survival and hazard curves. Median survival is obtained by linear interpolation on the fitted pred_times grid. If the true median falls before the first grid point there is still an exact known point to interpolate from, S(0) = 1, so the interpolation runs between (0, 1) and the first fitted time rather than returning that time as if it were the median. The estimate is still coarser there than a denser grid would give, so pass a pred_times with more resolution near zero when very early medians are expected.

knitr::kable(predict(fit_mspline, type = "rmst"),    caption = "Restricted mean survival time")

knitr::kable(predict(fit_mspline, type = "median"),  caption = "Median survival")

The model-based standardized survival curves. predict() returns both populations and its plot() method facets by population, so the index population (the decision-relevant target, see above) is shown alongside the comparator. predict() anchors the model curves at (time = 0, survival = 1); one manual scale keeps the index treatment green and the comparator red:

surv_cols <- c(Len = "#00BA38", Thal = "#F8766D")   # green index, red comparator
plot(predict(fit_mspline, type = "survival")) +
  scale_color_manual(values = surv_cols, aesthetics = c("colour", "fill"))

plot of chunk surv-curves

Each arm’s own observed Kaplan-Meier is a check on the fitted curve for that arm, and it can only be read in the population the arm was observed in. geom_km() (mlumr’s analogue of multinma’s) draws them; the index arm’s KM against the index-population curves, and the comparator arm’s against the comparator- population curves:

plot(predict(fit_mspline, population = "index", type = "survival")) +
  geom_km(dat, dat$index_treatment, marks = FALSE) +
  scale_color_manual(values = surv_cols, aesthetics = c("colour", "fill")) +
  ggplot2::labs(subtitle = "Index population, with the observed Len KM")

plot of chunk surv-curves-km

plot(predict(fit_mspline, population = "comparator", type = "survival")) +
  geom_km(dat, dat$comparator_treatment, marks = FALSE) +
  scale_color_manual(values = surv_cols, aesthetics = c("colour", "fill")) +
  ggplot2::labs(subtitle = "Comparator population, with the observed Thal KM")

plot of chunk surv-curves-km

The standardized marginal hazard, again in both populations. It is time-varying, and the two treatments’ marginal hazards are non-proportional after integrating over covariates even though the conditional model is PH:

plot(predict(fit_mspline, type = "hazard")) +
  scale_color_manual(values = surv_cols, aesthetics = c("colour", "fill"))

plot of chunk hazard-curve

Effect modification (relaxed M-spline model)

multinma’s NDMM example specifies covariate-by-treatment interactions (~ (covariates) * .trt), i.e. effect modification. The unanchored analogue is the relaxed model fitted above (fit_mspline_relaxed), which gives each treatment its own prognostic coefficients:

knitr::kable(marginal_effects(fit_mspline_relaxed, effect = "rmstd"),
   caption = "Relaxed M-spline RMST difference (effect modification, collapsible summary)")

Under the relaxed model the conditional hazard ratio depends on the covariates, so there is no single conditional HR. The scalar HR is then only the marginal hazard ratio at one instant, and the at_time column says which: the \(t \rightarrow 0\) limit for an unstratified fit, and otherwise the first prediction time, which is positive and set by pred_times. Prefer predict(type = "loghr") for the time-varying contrast, or the collapsible RMST estimands (effect = "rmstd" / "rmstr"), whose horizon is reported beside them. Caveat (unanchored two-study setting). Unlike multinma’s connected network, here the comparator coefficients are identified only through the single reconstructed comparator likelihood, so the index-population relaxed estimand extrapolates them over the IPD covariate distribution and is correspondingly wide. We regularize beta_comparator with a tighter autoscaled prior (set in the fit above), report both populations so the extra width of the index-population estimand is visible rather than hidden, and check prior_sensitivity(); see the Relaxed-model identification notes in NEWS and ?mlumr.

Model comparison

The choice of survival baseline (parametric versus flexible) is itself a model comparison. We compare the Weibull and M-spline fits by leave-one-out cross-validation; a lower looic favors the more flexible baseline when the hazard shape departs from the Weibull form:

compare_models(Weibull = fit_weibull, MSpline = fit_mspline, criterion = "loo")
#> 
#> Model Comparison (LOO)
#> ======================
#> 
#>    model elpd_diff se_diff p_worse diag_diff diag_elpd
#>  MSpline       0.0     0.0      NA                    
#>  Weibull      -5.0     3.9    0.90                    
#> 
#> elpd_diff is the difference in expected log pointwise predictive
#> density vs the best model, and se_diff is its standard error: the
#> uncertainty about that difference, not evidence for it. Read the two
#> together. A difference small relative to se_diff is not distinguished
#> from zero by this comparison, whatever se_diff itself is.
#> Treat any ratio as a heuristic, not a decision rule, and check the
#> PSIS diagnostics and whether the difference matters for the
#> prediction you care about.

The Weibull fit’s own standardized effects, in both populations, so the baseline-shape choice can be judged on the estimand and not only on looic:

knitr::kable(marginal_effects(fit_weibull, effect = "all"),
   caption = "Weibull population-standardized effects (Len vs Thal), both populations")

plot(predict(fit_weibull, type = "survival")) +
  scale_color_manual(values = surv_cols, aesthetics = c("colour", "fill"))

plot of chunk weibull-surv

Log hazard ratios against the naive estimate

The naive Cox estimate under Frequentist benchmarks is a log hazard ratio, so it is set against ML-UMR on that scale here, with ML-UMR in both target populations.

# ML-UMR reports a log HR in each population; pull both from the fit summary.
lhr_row <- function(fit, var) fit$summary[fit$summary$variable == var, ]

# `fit$summary` stores coefficients under Stan's positional names, `beta[1]`
# and so on; only `print()` relabels them as `beta[age]`. Look them up by
# position so the prose below cannot silently cite an empty result.
beta_mean <- function(fit, cov) {
  k <- match(cov, fit$data$covariates)
  # Fail loudly. An unmatched name gives `beta[NA]`, which matches nothing and
  # renders as blank prose with no error, which is the failure this helper
  # exists to avoid. A relaxed fit stores `beta_index[k]` and
  # `beta_comparator[k]`, so `beta[k]` would be empty there too.
  stopifnot(length(k) == 1L, !is.na(k), fit$model == "spfa")
  out <- fit$summary$mean[fit$summary$variable == sprintf("beta[%d]", k)]
  stopifnot(length(out) == 1L)
  out
}

lhr_i     <- lhr_row(fit_mspline, "delta_index")
lhr_c     <- lhr_row(fit_mspline, "delta_comparator")
lhr_rel_i <- lhr_row(fit_mspline_relaxed, "delta_index")
lhr_rel_c <- lhr_row(fit_mspline_relaxed, "delta_comparator")

knitr::kable(data.frame(
  Method = c("Naive Cox (unadjusted)",
             "ML-UMR M-spline (SPFA, prognostic only)",
             "ML-UMR M-spline (SPFA, prognostic only)",
             "ML-UMR M-spline (relaxed, effect modification)",
             "ML-UMR M-spline (relaxed, effect modification)"),
  # The naive Cox row is standardized to no population; it contrasts the two
  # crude arms. `vignette("choosing-a-method")` uses the same label.
  Population = c("unstandardized", "index", "comparator", "index", "comparator"),
  `log HR` = c(res_naive$estimate,
               lhr_i$mean, lhr_c$mean, lhr_rel_i$mean, lhr_rel_c$mean),
  `2.5%`   = c(res_naive$ci_lower,
               lhr_i[["2.5%"]], lhr_c[["2.5%"]],
               lhr_rel_i[["2.5%"]], lhr_rel_c[["2.5%"]]),
  `97.5%`  = c(res_naive$ci_upper,
               lhr_i[["97.5%"]], lhr_c[["97.5%"]],
               lhr_rel_i[["97.5%"]], lhr_rel_c[["97.5%"]]),
  check.names = FALSE), digits = 3,
  caption = "Naive and ML-UMR log hazard ratios, in both target populations (Len vs Thal)")

The SPFA and relaxed fits separate in an informative way. In the comparator population they agree almost exactly (-0.483 and -0.490). In the index population SPFA is close to its own comparator value (-0.485), but the relaxed model moves to -1.047, roughly twice as far from the null, with an interval to match.

This divergence should not be read as evidence of effect modification; it is attributable to weak identification. beta_comparator is informed only by the single reconstructed comparator curve, which constrains beta'X but not the direction of beta. The comparator-population estimand does not require that direction to be resolved, whereas the index-population estimand extrapolates it across the IPD covariate distribution. marginal_effects() reports the marginal posterior-to-prior variance change for each comparator coefficient for this reason. That quantity is descriptive: a small or negative value can signal sensitivity worth investigating, but it is not a fraction of information learned from the data and does not establish the source of the index-population result.

Note where SPFA and relaxed agree and where they do not. In the comparator population, freeing the comparator coefficients costs precision without moving the estimate, which says the shared-effect assumption is not straining the data it can actually see. In the index population they part company. Read that as a statement about what two single-arm trials can identify rather than as a finding about the treatments: the honest summary here is the comparator-population contrast and the collapsible RMST difference, with the index-population relaxed value reported alongside its width rather than in place of them.

A note on estimands, since two kinds of log hazard ratio sit in one table. The naive figure is an unadjusted marginal Cox estimate, and ML-UMR reports a population-standardized marginal effect. The two ML-UMR SPFA values above are close (delta_index -0.485, delta_comparator -0.483), which says the standardization population barely matters here, but they are not identical and neither is a conditional contrast. With a baseline hazard per study the marginal hazard ratio also varies with time, so these scalars are its value at the first prediction time. Compare the methods on the collapsible RMST difference, or on the predict(type = "loghr") curve, when the distinction matters.

ML-UMR does not substitute for randomization: where a common comparator exists it should be used. Where none exists, both the SPFA and relaxed fits should be reported so that the contribution of the modeling assumptions is visible.

Interpretation

• The naive Cox log hazard ratio ignores a real prognostic-factor imbalance between McCarthy2012 (Len) and Morgan2012 (Thal): mean age 57.9 against 65.6 years, and 27% against 32% ISS stage 3. Adjusting for it moves the SPFA estimate toward the null, from -0.554 to -0.485 in the index population.
• The prognostic coefficients are the part to check first, because they are what the adjustment is built from. Here they have plausible signs and sizes: age 0.07 per year, ISS stage III 0.40, and the response and sex coefficients small, with intervals spanning zero.
• ML-UMR adjusts for the imbalance (SPFA -0.485 in the index population and -0.483 in the comparator, against the naive -0.554). Its intervals are much wider than the naive one, which reflects sampling error alone and carries none of the uncertainty of adjusting. Read the SPFA and relaxed fits together: they agree in the comparator population (-0.490) but diverge in the index population (-1.047), which measures how much of the index-population answer rests on comparator coefficients the aggregate curve cannot identify. ML-UMR adds what the frequentist benchmarks cannot: full posterior uncertainty, and effects in both populations.
• Because the two studies have their own baseline hazards (aux_by = ".study", the default) the marginal hazard ratio genuinely varies with time, and because the hazard ratio is also non-collapsible the marginal and conditional contrasts differ under effect modification. Prefer the time-varying log HR curve (predict(type = "loghr")) and the collapsible RMST estimands for scalar summaries; the scalar marginal HR is only the contrast at the first prediction time.
• The relaxed model frees the shared-effect assumption but, with the comparator coefficients identified from one reconstructed likelihood, pays for it in precision. Reserve ML-UMR for genuinely unanchored evidence; when a common comparator exists, prefer an anchored network meta-analysis.

Data provenance. ndmm_ipd / ndmm_agd / ndmm_agd_covs are bundled with mlumr, copied verbatim from the GPL-3 multinma package (Phillippo 2024), which provides simulated IPD and reconstructed Kaplan-Meier pseudo-IPD for the aggregate arm (Leahy and Walsh 2019). This unanchored two-study framing is for illustration.