Binary outcomes: an unanchored PASI 75
comparison
library(mlumr)
library(ggplot2)
options(mc.cores = parallel::detectCores())
This vignette is a complete worked example of an
unanchored indirect comparison for a
binary endpoint, using plaque-psoriasis trial data
bundled with mlumr (psoriasis_ipd /
psoriasis_agd, copied from the GPL-3 multinma
package (Phillippo
2024)). It is the unanchored, two-study analogue of
multinma’s plaque_psoriasis ML-NMR example (Phillippo et al.
2020). For the shared data-preparation machinery see
vignette("data-preparation"); for sampler/prior/diagnostic
detail see vignette("fitting-and-diagnostics"); for how to
choose between methods see
vignette("choosing-a-method").
About these data. psoriasis_ipd /
psoriasis_agd are copied verbatim from
multinma; they were not re-simulated for mlumr. multinma’s
patient records are themselves simulated to resemble
the published UNCOVER-2 (Griffiths et al. 2015) and FIXTURE
(Langley et al.
2014) trials. mlumr bundles one arm from each
trial: ixekizumab Q4W from UNCOVER-2 and secukinumab 300 mg
from FIXTURE. Both trials also randomized patients to placebo
and etanercept, and those arms are deliberately
omitted here. Keeping them would connect the two trials through
a common comparator and make an anchored comparison possible,
which is not the problem this vignette is about. The unanchored framing
is a construction for illustration.
Hypothetical single-arm trials. UNCOVER-2 and
FIXTURE are randomized trials. Dropping their common reference arms
(placebo and etanercept) turns them into hypothetical single-arm trials
and leaves fully disconnected evidence; we do this only to illustrate
ML-UMR. In practice, randomized trials should never be analyzed this
way: ML-NMR (for example with multinma) or another
appropriate method should be used for them. ML-UMR should only be used
for fully unanchored, single-arm comparisons.
The clinical question
We want the relative efficacy of two interleukin inhibitors on
PASI 75 response (≥ 75% improvement in the Psoriasis
Area and Severity Index) at week 12, when:
- Index (IPD). We hold individual patient data for
ixekizumab Q4W (
IXE_Q4W) from the
UNCOVER-2 trial (Griffiths et al. 2015).
- Comparator (AgD). Only published aggregate data are
available for secukinumab 300 mg (
SEC_300)
from the FIXTURE trial (Langley et al. 2014).
The two trials share no common arm, so the comparison is
unanchored and relies on the shared-effect-modifier /
conditional-constancy assumptions discussed in
vignette("introduction"). ML-UMR adjusts for cross-trial
differences in prognostic factors by standardizing the IPD model to the
comparator population (and vice versa).
The ML-UMR model
On the index IPD, ML-UMR fits an individual-level logistic regression
for the probability of PASI 75 response, \[\operatorname{logit}\Pr(y_{ik}=1\mid
x_i)=\mu_k+x_i^\top\beta_k,\] where \(\mu_k\) is a treatment-specific intercept
and \(\beta_k\) the vector of
prognostic-factor coefficients. The aggregate comparator does not
contribute patient rows; instead its likelihood is the individual model
integrated over the comparator covariate distribution
\(f_{\text{AgD}}(x)\) implied by the
published moments, \[\bar p_k=\Pr(Y=1\mid
\text{AgD
population})=\int\operatorname{logit}^{-1}(\mu_k+x^\top\beta_k)\,f_{\text{AgD}}(x)\,dx,\]
which mlumr approximates with quasi-Monte Carlo integration points (see
Setting up the ML-UMR data). This is the same
population-adjustment device as ML-NMR (Phillippo et al. 2020), specialized
to the unanchored two-study case.
Two variants differ only in how the prognostic effects are shared
across treatments:
- SPFA (shared prognostic factors) assumes \(\beta_{\text{index}}=\beta_{\text{comparator}}=\beta\),
the conditional-constancy / shared-effect-modifier assumption.
- Relaxed frees them, \(\beta_{\text{index}}\neq\beta_{\text{comparator}}\),
allowing effect modification at the cost of identifiability (the
comparator coefficients are informed only by the single AgD likelihood
term).
Population-standardized treatment effects are formed by averaging the
\(\operatorname{logit}^{-1}\)
predictions over a target population and contrasting the treatments,
yielding effects in both the index and comparator
populations.
SPFA is the unanchored analogue of the shared-effect-modifier
assumption used in anchored comparisons, but it is a
stronger assumption, so transportability is
correspondingly harder to justify (Chandler and Ishak 2025). A
useful consequence, reported from simulation (Chandler and Ishak 2025): in
the scenarios that study examined, the
comparator-population effect was recovered with little
bias even when SPFA was violated, because the prognostic coefficients
are anchored by the index IPD, whereas the
index-population effect was biased when SPFA was
wrongly imposed under strong effect modification. Read that as evidence
from the data-generating mechanisms, effect sizes and sample sizes that
were simulated, not as a guarantee that the comparator-population
estimand tolerates any violation. (Restricted to the comparator
population, the SPFA estimand has the same standardization target as
unanchored STC. The estimators still differ: ML-UMR fits a joint
likelihood in which the aggregate outcome informs the comparator
intercept, whereas STC regresses on the IPD alone.)
# IPD + AgD bundled with mlumr (copied from multinma's plaque_psoriasis, GPL-3).
data("psoriasis_ipd")
data("psoriasis_agd")
covs <- c("age", "bsa", "weight", "prevsys") # adjustment set for this example
# --- Index IPD: UNCOVER-2, ixekizumab Q4W -----------------------------------
ipd <- psoriasis_ipd
ipd$bsa <- ipd$bsa / 100 # body-surface area: % -> proportion
ipd <- ipd[ipd$study == "UNCOVER-2" & ipd$treatment == "IXE_Q4W", ]
ipd <- ipd[stats::complete.cases(ipd[, c("pasi75", covs)]), ]
# --- Comparator AgD: FIXTURE, secukinumab 300 mg ----------------------------
agd <- psoriasis_agd
agd$bsa_mean <- agd$bsa_mean / 100
agd$bsa_sd <- agd$bsa_sd / 100
agd <- agd[agd$study == "FIXTURE" & agd$treatment == "SEC_300", ]
A glimpse of the individual patient data (one row per patient) and
the single aggregate comparator row:
knitr::kable(head(ipd[, c("study", "treatment", "pasi75", covs)]),
caption = "Index IPD (UNCOVER-2, ixekizumab Q4W), first rows")
Index IPD (UNCOVER-2, ixekizumab Q4W), first rows
| UNCOVER-2 |
IXE_Q4W |
1 |
49 |
0.25 |
88.7 |
1 |
| UNCOVER-2 |
IXE_Q4W |
1 |
33 |
0.29 |
96.3 |
1 |
| UNCOVER-2 |
IXE_Q4W |
1 |
25 |
0.40 |
93.6 |
1 |
| UNCOVER-2 |
IXE_Q4W |
0 |
60 |
0.20 |
112.0 |
1 |
| UNCOVER-2 |
IXE_Q4W |
1 |
49 |
0.20 |
119.3 |
0 |
| UNCOVER-2 |
IXE_Q4W |
1 |
58 |
0.16 |
92.0 |
1 |
knitr::kable(agd[, c("study", "treatment", "pasi75_r", "pasi75_n",
"age_mean", "bsa_mean", "weight_mean", "prevsys_prop")],
caption = "Comparator AgD (FIXTURE, secukinumab 300 mg)")
Comparator AgD (FIXTURE, secukinumab 300 mg)
| FIXTURE |
SEC_300 |
249 |
323 |
44.5 |
0.343 |
83 |
0.63 |
Covariate balance
Population adjustment is motivated by the fact that the two trial
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)", "Body-surface area (prop.)",
"Weight (kg)", "Previous systemic (prop.)"),
Index_IXE = c(mean(ipd$age), mean(ipd$bsa), mean(ipd$weight), mean(ipd$prevsys)),
Comparator_SEC = c(agd$age_mean, agd$bsa_mean, agd$weight_mean, agd$prevsys_prop)
)
knitr::kable(balance, caption = "Prognostic-factor balance: index vs comparator population")
Prognostic-factor balance: index vs comparator
population
| Age (years) |
45.7275362 |
44.500 |
| Body-surface area (prop.) |
0.2761159 |
0.343 |
| Weight (kg) |
94.0605797 |
83.000 |
| Previous systemic (prop.) |
0.6782609 |
0.630 |
Setting up the ML-UMR data
set_ipd() takes the patient-level data;
set_agd() takes the comparator event count
(outcome_r), sample size (outcome_n), and
covariate summaries. Covariate column suffixes (_mean,
_sd) are stripped to match the IPD names.
ipd_obj <- set_ipd(ipd, treatment = "treatment", outcome = "pasi75", covariates = covs)
agd_obj <- set_agd(agd, treatment = "treatment",
outcome_n = "pasi75_n", outcome_r = "pasi75_r",
cov_means = c("age_mean", "bsa_mean", "weight_mean", "prevsys_prop"),
cov_sds = c("age_sd", "bsa_sd", "weight_sd", NA),
cov_types = c("continuous", "continuous", "continuous", "binary"))
dat <- combine_data(ipd_obj, agd_obj)
dat
#> Unanchored Comparison Data (Binary)
#> ====================================
#>
#> Index treatment (IPD): IXE_Q4W
#> N = 345
#> Events = 267 (77.4%)
#>
#> Comparator treatment (AgD): SEC_300
#> N = 323
#> Events = 249 (77.1%)
#>
#> Covariates ( 4 ): age, bsa, weight, prevsys
#> Integration points: not yet added (use add_integration())
The comparator covariate distribution is represented by quasi-Monte
Carlo integration points drawn from each covariate’s marginal,
correlated via a Gaussian copula calibrated from the IPD. The marginals
are the ones multinma’s own plaque-psoriasis example uses: gammas for
the right-skewed age and weight, a logit-normal for body-surface area (a
proportion, so it must stay inside (0, 1)), and a Bernoulli
for previous systemic therapy. We use a modest n_int = 64
here so the vignette builds quickly; a real analysis would use several
hundred or more (see vignette("data-preparation")).
dat <- add_integration(
dat, n_int = 64,
# Marginals follow multinma's own plaque-psoriasis example: gamma for the
# right-skewed continuous covariates, logit-normal for body surface area
# (a proportion), Bernoulli for the binary one. mlumr's qgamma() and
# qlogitnorm() accept a mean and sd directly, so the published baseline
# table can be used as printed.
age = distr(qgamma, mean = age_mean, sd = age_sd),
bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
prevsys = distr(qbern, prob = prevsys_mean)
)
check_integration() confirms the quasi-Monte Carlo
points reproduce the requested covariate moments, run it whenever
covariates are correlated or n_int is small:
check_integration(
dat,
age = distr(qgamma, mean = age_mean, sd = age_sd),
bsa = distr(qlogitnorm, mean = bsa_mean, sd = bsa_sd),
weight = distr(qgamma, mean = weight_mean, sd = weight_sd),
prevsys = distr(qbern, prob = prevsys_mean)
)
#> Integration check: n_int = 64 vs 128
#> Resolution heuristic, max relative difference: 0.0348
#> Caution: 1-5% marginal relative difference. Consider increasing n_int.
#> Declared-target fidelity, max relative difference: 0.0712
#> Warning: grid moments differ from declared AgD moments by >5%.
#> Joint: max |cor(current) - cor(doubled)|: 0.1025
#> Warning: pairwise correlations differ by > 0.05 between resolutions.
#> Target (spearman): max |cor(doubled) - cor_target|: 0.0390
Frequentist benchmarks
The unadjusted (naive) comparison ignores the covariate imbalance;
STC adjusts for it by parametric G-computation. Both are fast and serve
as reference points.
res_naive <- naive(dat)
res_stc <- stc(dat)
res_naive
#> Naive Unadjusted Indirect Comparison
#> =====================================
#>
#> Treatments: IXE_Q4W vs SEC_300
#>
#> Population basis: index-study outcome versus comparator-population outcome; no common standardized target.
#>
#> Event rates:
#> Index (IPD): 0.774 (267/345)
#> Comparator (AgD): 0.771 (249/323)
#>
#> Log Odds Ratio: 0.0172 (SE: 0.1846)
#> 95% CI: [-0.3448, 0.3791]
#>
#> All effect measures (95% CI):
#> Log odds ratio 0.0172 (SE 0.1846) [-0.3448, 0.3791]
#> Odds ratio 1.0173 [0.7084, 1.4609]
#> Risk difference 0.0030 (SE 0.0325) [-0.0606, 0.0666]
#> Risk ratio 1.0039 [0.9245, 1.0901]
res_stc
#> Simulated Treatment Comparison (G-computation)
#> ===============================================
#>
#> Treatments: IXE_Q4W vs SEC_300
#>
#> 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.
#>
#> Marginalized P(Y=1|index trt, comp pop): 0.8106
#> Observed P(Y=1|comp trt, comp pop): 0.7709
#>
#> Log Odds Ratio: 0.2403 (SE: 0.2029)
#> 95% CI: [-0.1575, 0.6380]
#>
#> All effect measures (95% CI):
#> Log odds ratio 0.2403 (SE 0.2029) [-0.1575, 0.6380]
#> Odds ratio 1.2716 [0.8543, 1.8928]
#> Risk difference 0.0397 (SE 0.0332) [-0.0255, 0.1048]
#> Risk ratio 1.0515 [0.9683, 1.1418]
#>
#> Outcome model coefficients:
#> (Intercept) age bsa weight prevsys
#> 4.3879 -0.0404 0.0676 -0.0139 0.1003
Fitting the ML-UMR models
The binomial family supports "logit" (default),
"probit", and "cloglog" links. We fit both the
shared-prognostic-factor model (SPFA) and the
relaxed model (treatment-specific prognostic effects,
allowing effect modification). Because age is on a ~50-unit scale,
autoscale = TRUE keeps the coefficient prior’s contribution
on a common SD scale. The numerical scale is still an
application-specific choice and should be checked predictively.
fit_spfa <- mlumr(dat, model = "spfa", link = "logit",
prior_beta = prior_normal(0, 2.5, autoscale = TRUE),
chains = 4, iter = 2000, warmup = 1000,
seed = 2026, refresh = 0)
#> Running MCMC with 4 parallel chains...
#>
#> Chain 3 finished in 3.0 seconds.
#> Chain 1 finished in 3.2 seconds.
#> Chain 2 finished in 3.3 seconds.
#> Chain 4 finished in 3.2 seconds.
#>
#> All 4 chains finished successfully.
#> Mean chain execution time: 3.2 seconds.
#> Total execution time: 3.6 seconds.
fit_relaxed <- mlumr(dat, model = "relaxed", link = "logit",
prior_beta = prior_normal(0, 2.5, autoscale = TRUE),
chains = 4, iter = 2000, warmup = 1000,
seed = 2026, refresh = 0)
#> Running MCMC with 4 parallel chains...
#>
#> Chain 3 finished in 12.6 seconds.
#> Chain 2 finished in 14.2 seconds.
#> Chain 1 finished in 15.7 seconds.
#> Chain 4 finished in 16.0 seconds.
#>
#> All 4 chains finished successfully.
#> Mean chain execution time: 14.6 seconds.
#> Total execution time: 16.3 seconds.
summary(fit_spfa)
#> ML-UMR Model Summary
#> ====================
#>
#> Model: SPFA
#> Family: Binary
#> Link: logit
#> Engine: cmdstanr
#> Treatments: IXE_Q4W (IPD) vs SEC_300 (AgD)
#>
#> MCMC Diagnostics:
#> Divergent transitions: 0
#> Max treedepth hits: 0
#> Max Rhat: 1.004
#> Min ESS: 1546
#>
#> Intercepts (logit scale):
#> variable mean sd 2.5% 97.5% Rhat
#> mu_index 1.447602 0.1526277 1.1576075 1.742907 1.000540
#> mu_comparator 1.180924 0.1513114 0.8840591 1.473595 1.000812
#>
#> Regression Coefficients:
#> variable mean sd 2.5% 97.5% Rhat
#> beta[age] -0.04150236 0.01130282 -0.06426055 -0.019883407 1.001242
#> beta[bsa] 0.11811085 0.75278434 -1.33547013 1.634856495 1.001090
#> beta[weight] -0.01422774 0.00634066 -0.02624840 -0.001559571 1.003903
#> beta[prevsys] 0.09495758 0.30194582 -0.50583182 0.687454176 1.000363
#>
#> Marginal Treatment Effects:
#> Log Odds Ratios:
#> variable mean sd 2.5% 97.5%
#> lor_index 0.2463067 0.2043813 -0.1612142 0.6518974
#> lor_comparator 0.2472410 0.2052714 -0.1614905 0.6557578
#> Risk Differences:
#> variable mean sd 2.5% 97.5%
#> rd_index 0.04626266 0.03863495 -0.02944331 0.1231594
#> rd_comparator 0.04043621 0.03363255 -0.02668746 0.1069070
#> Risk Ratios:
#> variable mean sd 2.5% 97.5%
#> rr_index 1.065625 0.05615489 0.9623604 1.181424
#> rr_comparator 1.053508 0.04515453 0.9667412 1.146098
Priors
These fits use \(\mu_k\sim\mathrm{N}(0,10^2)\) for the
intercepts and \(\beta\sim\mathrm{N}(0,2.5^2)\) for the
coefficients, the latter autoscaled by each covariate’s
standard deviation. Autoscaling puts contributions on a common SD scale;
it does not establish that 2.5 is appropriate for a particular endpoint.
prior_summary() shows exactly what was passed to Stan,
including the autoscaled scales:
prior_summary(fit_spfa)
#> Priors for ML-UMR Fit
#> =====================
#>
#> Intercepts (mu_index, mu_comparator):
#> normal(0, 10)
#> (package default, mlumr 0.2.0)
#>
#> Regression coefficients (beta):
#> Family: normal
#> coefficient mean scale autoscaled sd_x
#> age 0 0.196 TRUE 12.784
#> bsa 0 14.055 TRUE 0.178
#> weight 0 0.120 TRUE 20.794
#> prevsys 0 5.344 TRUE 0.468
#> (scale = user_scale / sd_x for autoscaled rows)
The data are substantially more informative than the priors, the
posterior intercepts are much tighter than the \(\mathrm{N}(0,10)\) prior:
plot_prior_posterior(fit_spfa, pars = c("mu_index", "mu_comparator"))
Handling quasi-complete separation
With small or sparse binary data, rare events, or a covariate that
almost perfectly predicts the outcome, the maximum-likelihood logistic
fit can diverge (quasi-complete separation), and a flat or very
wide prior lets the posterior chase it to implausibly large
coefficients. A proper, substantively calibrated coefficient prior
regularizes separated directions. Gelman et al. (2008) propose a
Cauchy prior after a particular predictor scaling; a moderate Student-t
is often easier to sample, but neither family makes 2.5 universally
appropriate. No separate model is needed: the same
prior_beta argument accepts a Student-t (or Cauchy)
prior:
# Student-t coefficients (df 3-7 is a robust default), autoscaled as above:
fit_sep <- mlumr(dat, model = "spfa",
prior_beta = prior_student_t(df = 5, 0, 2.5, autoscale = TRUE))
# prior_cauchy(0, 2.5) is the df = 1 special case (heaviest tails).
Use prior-predictive checks to set plausible scales and
prior_sensitivity() to report how the result changes across
the scales examined.
Convergence
mlumr() warns automatically about divergences, treedepth
saturation, high \(\hat R\), and low
effective sample size. They are clean here:
data.frame(
n_divergent = fit_spfa$diagnostics$n_divergent,
max_treedepth = fit_spfa$diagnostics$n_max_treedepth,
max_Rhat = round(max(fit_spfa$summary$Rhat, na.rm = TRUE), 3),
min_ESS = round(min(fit_spfa$summary$n_eff, na.rm = TRUE))
)
#> n_divergent max_treedepth max_Rhat min_ESS
#> 1 0 0 1.004 1546
Treatment effects
marginal_effects() reports three population-standardized
binary effect measures in each population: the log odds
ratio (LOR), the risk difference
(RD), and the risk ratio
(RR). Because ML-UMR targets both populations, you
get a population-specific effect for each.
me_spfa <- marginal_effects(fit_spfa, effect = "all")
knitr::kable(me_spfa, caption = "SPFA population-standardized effects (IXE_Q4W vs SEC_300)")
SPFA population-standardized effects (IXE_Q4W vs
SEC_300)
| lor_index |
LOR |
Index |
0.2463067 |
0.2043813 |
-0.1612142 |
0.2495905 |
0.6518974 |
| lor_comparator |
LOR |
Comparator |
0.2472410 |
0.2052714 |
-0.1614905 |
0.2503917 |
0.6557578 |
| rd_index |
RD |
Index |
0.0462627 |
0.0386350 |
-0.0294433 |
0.0465180 |
0.1231594 |
| rd_comparator |
RD |
Comparator |
0.0404362 |
0.0336325 |
-0.0266875 |
0.0409191 |
0.1069070 |
| rr_index |
RR |
Index |
1.0656252 |
0.0561549 |
0.9623604 |
1.0637364 |
1.1814243 |
| rr_comparator |
RR |
Comparator |
1.0535080 |
0.0451545 |
0.9667412 |
1.0530110 |
1.1460981 |
A forest plot makes the value of adjustment visible, on the
log odds ratio. ML-UMR is shown in
both target populations. The index
population is normally the decision-relevant one for HTA, since
cost-effectiveness models are built for the population a reimbursement
decision is about, which is usually the index trial’s (Chandler and
Ishak 2026). STC is standardized to the comparator
population; the naive contrast is not standardized to either population.
Showing the comparator ML-UMR row provides the like-for-like STC
comparison.
me_spfa_b <- marginal_effects(fit_spfa, effect = "lor", population = "both")
me_rel_b <- marginal_effects(fit_relaxed, effect = "lor", population = "both")
# One row per (model, population); `pop()` pulls the requested one.
pop <- function(d, which) d[d$population == which, ]
forest_df <- data.frame(
label = c("Naive (unstandardized)", "STC (comparator)",
"ML-UMR SPFA (index)", "ML-UMR SPFA (comparator)",
"ML-UMR relaxed (index)", "ML-UMR relaxed (comparator)"),
est = c(res_naive$link_effect, res_stc$link_effect,
pop(me_spfa_b, "Index")$mean, pop(me_spfa_b, "Comparator")$mean,
pop(me_rel_b, "Index")$mean, pop(me_rel_b, "Comparator")$mean),
lo = c(res_naive$ci_lower, res_stc$ci_lower,
pop(me_spfa_b, "Index")$q2.5, pop(me_spfa_b, "Comparator")$q2.5,
pop(me_rel_b, "Index")$q2.5, pop(me_rel_b, "Comparator")$q2.5),
hi = c(res_naive$ci_upper, res_stc$ci_upper,
pop(me_spfa_b, "Index")$q97.5, pop(me_spfa_b, "Comparator")$q97.5,
pop(me_rel_b, "Index")$q97.5, pop(me_rel_b, "Comparator")$q97.5)
)
mlumr_forest(forest_df, ref_line = 0,
x = "Log odds ratio",
title = "PASI 75: ixekizumab vs secukinumab",
subtitle = "Unadjusted vs population-adjusted, in both target populations")
The full set of standardized effects in both populations, as a
table:
knitr::kable(marginal_effects(fit_relaxed, effect = "all", population = "both"),
caption = "Relaxed-model standardized effects, both populations")
Relaxed-model standardized effects, both populations
| lor_index |
LOR |
Index |
0.0028398 |
0.5849539 |
-1.1027831 |
0.0207378 |
1.0585269 |
| lor_comparator |
LOR |
Comparator |
0.2444256 |
0.1994718 |
-0.1371652 |
0.2421185 |
0.6322585 |
| rd_index |
RD |
Index |
0.0144135 |
0.1020933 |
-0.1437769 |
0.0035884 |
0.2314025 |
| rd_comparator |
RD |
Comparator |
0.0400623 |
0.0326564 |
-0.0234464 |
0.0398917 |
0.1027967 |
| rr_index |
RR |
Index |
1.0388531 |
0.1558008 |
0.8393618 |
1.0046632 |
1.4242271 |
| rr_comparator |
RR |
Comparator |
1.0529897 |
0.0437571 |
0.9703815 |
1.0518831 |
1.1403727 |
The population-standardized log odds ratios, in both populations,
plot straight from marginal_effects() with its
plot() method (point estimate and 95% credible interval,
faceted by measure):
plot(marginal_effects(fit_spfa, effect = "lor"))
Absolute predictions
Absolute predicted response probabilities for each treatment in each
population, as a table and as a plot. plot() on a
predict() result draws the point-intervals directly, with
the two populations distinguished by color:
knitr::kable(predict(fit_spfa, population = "both", type = "response"),
caption = "Standardized PASI 75 response probabilities")
Standardized PASI 75 response probabilities
| IXE_Q4W |
Index |
0.7742054 |
0.0215267 |
0.7311367 |
0.7744590 |
0.8156391 |
| SEC_300 |
Index |
0.7279427 |
0.0317215 |
0.6626771 |
0.7288011 |
0.7869427 |
| IXE_Q4W |
Comparator |
0.8110121 |
0.0230334 |
0.7637451 |
0.8119185 |
0.8531067 |
| SEC_300 |
Comparator |
0.7705759 |
0.0235685 |
0.7221162 |
0.7712480 |
0.8138316 |
plot(predict(fit_spfa, population = "both", type = "response"))
Conditional effects
Conditional effects evaluate the contrast at specific covariate
profiles rather than averaging over a population, useful for checking
whether the effect is plausibly constant across the prognostic range
(the SPFA assumption) or varies with the covariates:
# Weight stays in kilograms in this vignette (IPD mean about 94), unlike
# `vignette("choosing-a-method")`, which divides it by 10. Profiles must be on
# the scale the model was fitted on, or the contrast is evaluated outside the
# covariate support.
profiles <- data.frame(age = c(40, 55, 70), bsa = c(0.15, 0.30, 0.45),
weight = c(75, 90, 110), prevsys = c(0, 1, 1))
knitr::kable(conditional_effects(fit_spfa, newdata = profiles),
caption = "Conditional log odds ratios at three covariate profiles")
Conditional log odds ratios at three covariate
profiles
| 1 |
LINK_EFFECT |
0.2666779 |
0.2222286 |
-0.1697892 |
0.2697886 |
0.7096585 |
| 1 |
RD |
0.0360436 |
0.0308723 |
-0.0236670 |
0.0355675 |
0.1003421 |
| 1 |
RR |
1.0454193 |
0.0401692 |
0.9714840 |
1.0431503 |
1.1330295 |
| 2 |
LINK_EFFECT |
0.2666779 |
0.2222286 |
-0.1697892 |
0.2697886 |
0.7096585 |
| 2 |
RD |
0.0545626 |
0.0458718 |
-0.0336336 |
0.0541125 |
0.1460648 |
| 2 |
RR |
1.0831630 |
0.0727151 |
0.9548934 |
1.0782510 |
1.2388811 |
| 3 |
LINK_EFFECT |
0.2666779 |
0.2222286 |
-0.1697892 |
0.2697886 |
0.7096585 |
| 3 |
RD |
0.0646590 |
0.0537742 |
-0.0413750 |
0.0653554 |
0.1729787 |
| 3 |
RR |
1.1526318 |
0.1400779 |
0.9261354 |
1.1364555 |
1.4762950 |
Model comparison
Compare SPFA against the relaxed model with leave-one-out
cross-validation and DIC. A markedly better relaxed fit would
suggest effect modification, but the comparator coefficients
are weakly identified from a single AgD row, so treat it as suggestive
and check prior_sensitivity() (see
vignette("fitting-and-diagnostics")).
compare_models(SPFA = fit_spfa, Relaxed = fit_relaxed, criterion = "loo")
#>
#> Model Comparison (LOO)
#> ======================
#>
#> model elpd_diff se_diff p_worse diag_diff diag_elpd
#> Relaxed 0.0 0.0 NA 1 k_psis > 0.7
#> SPFA -0.2 0.1 0.97 |elpd_diff| < 4 1 k_psis > 0.7
#>
#> 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.
compare_models(SPFA = fit_spfa, Relaxed = fit_relaxed, criterion = "dic")
#>
#> Model Comparison (DIC)
#> ======================
#>
#> Model DIC pD Delta_DIC
#> Relaxed 367.87 5.71 0.00
#> SPFA 368.60 6.22 0.73
#>
#> Lower DIC = better fit. Delta_DIC > 5 is a rough heuristic for
#> meaningful difference, not a formally calibrated threshold.
#> DIC should not be the sole basis for model selection.
Interpretation
- The naive log odds ratio mixes the treatment
contrast with between-study differences, including prognostic-factor
imbalance.
- In this example STC and ML-UMR
SPFA both standardize over the named covariates and broadly
agree. ML-UMR reports posterior uncertainty and effects in both
populations, which matters when the decision population is not the
comparator’s.
- The relaxed model relaxes the shared-effect
assumption but pays for it in precision here, because the
comparator-specific coefficients lean on one aggregate row. Its
comparator-population estimate stays close to SPFA, but its
index-population estimate falls back to about the naive value, with an
interval more than two log odds ratios wide (see the both-populations
table above). The index-population estimand extrapolates comparator
coefficients that a single aggregate row cannot identify, so read that
divergence as a statement about identification, not about effect
modification.
Data provenance. psoriasis_ipd /
psoriasis_agd are bundled with mlumr, copied verbatim from
the GPL-3 multinma package (Phillippo 2024)
(plaque_psoriasis_ipd / plaque_psoriasis_agd),
which provides simulated IPD resembling the UNCOVER-2 (Griffiths et al.
2015) and FIXTURE (Langley et al. 2014) trials. This
unanchored two-study framing is for illustration; the trials were not
designed for a head-to-head indirect comparison.
References
Chandler, C., and K. J. Ishak. 2025.
Anchors Away: Navigating
Unanchored Indirect Comparisons with Multilevel Unanchored
Meta-Regression. ISPOR Europe, Glasgow, UK; abstract MSR28.
https://www.valueinhealthjournal.com/article/S1098-3015(25)05944-3/abstract.
Chandler, Conor, and K. Jack Ishak. 2026.
Reframing
Population-Adjusted Indirect Comparisons as a Transportability Problem:
An Estimand-Based Perspective and Implications for Health Technology
Assessment. Preprint.
https://doi.org/10.48550/arXiv.2602.17041.
Gelman, Andrew, Aleks Jakulin, Maria Grazia Pittau, and Yu-Sung Su.
2008.
“A Weakly Informative Default Prior Distribution for
Logistic and Other Regression Models.” The Annals of Applied
Statistics 2 (4): 1360–83.
https://doi.org/10.1214/08-AOAS191.
Griffiths, C. E. M., K. Reich, M. Lebwohl, et al. 2015.
“Comparison of Ixekizumab with Etanercept or Placebo in
Moderate-to-Severe Psoriasis (UNCOVER-2 and
UNCOVER-3): Results from Two Phase 3 Randomised
Trials.” The Lancet 386 (9993): 541–51.
https://doi.org/10.1016/s0140-6736(15)60125-8.
Langley, R. G., B. E. Elewski, M. Lebwohl, et al. 2014.
“Secukinumab in Plaque Psoriasis — Results of Two Phase 3
Trials.” New England Journal of Medicine 371 (4):
326–38.
https://doi.org/10.1056/nejmoa1314258.
Phillippo, D. M., S. Dias, A. E. Ades, et al. 2020.
“Multilevel
Network Meta-Regression for Population-Adjusted Treatment
Comparisons.” Journal of the Royal Statistical Society:
Series A (Statistics in Society) 183 (3): 1189–210.
https://doi.org/10.1111/rssa.12579.
Phillippo, David M. 2024.
multinma:
Bayesian Network Meta-Analysis of Individual and Aggregate
Data.
https://doi.org/10.5281/zenodo.3904454.
