delta <- -4 - (-3)
delta_se <- sqrt(0.6^2 + 0.2^2)
ci_lo <- delta - delta_se * qnorm(0.975)
ci_up <- delta + delta_se * qnorm(0.975)
round(data.frame(delta, ci_lo, ci_up), 1) delta ci_lo ci_up
1 -1 -2.2 0.2From Bucher to Bayes
A Brief Introduction to Bayesian Model-Based Network Meta-Analysis for Indirect Treatment Comparisons using R
Introduction
It’s common for meta-analyses in the clinical trial literature to lump together studies with different doses, follow-up times, or populations. One common workaround is to stratify the meta-analysis into subgroups. But often we can do better than that and instead model the dependencies or discrepancies explicitly. One potential payoff is increased statistical precision.
In this post, I’ll give a brief introduction to Bayesian model-based network meta-analysis (MBNMA) to model studies on different dose levels of the same drug and show how it can be used to inform an indirect treatment comparison between competitive drugs that have not been studied in the same trial. We simulate data to keep the example agnostic of therapy area.
To follow along, I assume some basic familiarity with Bayesian statistics and R (R Core Team 2024a). I’ll use the rethinking package (McElreath 2020), an interface to Stan (Carpenter et al. 2017), because its syntax mirrors the formal model very closely and makes the structure of MBNMA easy to see.
The evidence network
Suppose we have three trials:
- Trial A1: Drug A at dose 1 mg vs. placebo. Estimated treatment difference: \(-4\) (SE \(\pm 0.6\)). Smaller trial, \(N \approx 100\).
- Trial A2: Drug A at dose 2 mg vs. placebo. Estimated treatment difference: \(-4.5\) (SE \(\pm 0.3\)). Larger trial, \(N \approx 400\).
- Trial B1: Drug B at dose 1 mg vs. placebo. Estimated treatment difference: \(-3\) (SE \(\pm 0.2\)).
The outcome is some continuous endpoint where lower (more negative) is better. There is no head-to-head trial of Drug A vs. Drug B but let’s say they are expected to compete in the same market on their respective 1 mg dose.
The evidence network looks like this:
Drug A, 2 mg
|
placebo
/ \
Drug A, 1 mg Drug B, 1 mg
\\------//
(indirect comparison)Solid edges are the direct randomised comparisons against placebo. The dashed edge between Drug A 1 mg and Drug B 1 mg is the indirect comparison we’re after.
The Bucher method
The “Bucher method” (Bucher et al. 1997) takes the dose 1 mg estimate of Drug A and the 1 mg estimate of Drug B and subtracts them, and the confidence interval is constructed by adding the variances of the two trial estimates.
The point estimate is \(-1\) and the 95% CI roughly \(-2.2\) to \(+0.2\). The estimate favors Drug A, but the interval just crosses zero – it’s not “statistically significant.” This is where we can exploit the dose-response structure of the evidence network sketched above.
The Drug A programme studied Drug A at both 1 mg and 2 mg, providing two data points on the same dose-response curve. The 2 mg estimate is substantially more precise than the 1 mg estimate, because the 2 mg trial had four times the sample size compared to the 1 mg trial. Both doses are estimating effects that arise from the same underlying pharmacology. A model that exploits this structure can allow the precision of the 2 mg estimate to flow back along the curve and improve our estimate of 1 mg.
Taken together, this setup – multiple trials sharing a common comparator – is the structure of a network meta-analysis: by treating the trials as a connected network rather than a set of isolated pairwise comparisons, we can estimate effects between treatments that have not been directly compared in the same trial (Lu and Ades 2004; Lumley 2002). The “model-based” bit in MBNMA refers to the additional structure we layer on top (Mawdsley et al. 2016) – here, a dose–response curve linking the two Drug A trials.
Let’s start by collecting the data in a list.
dlist <- list(D = c(1.0, 2.0),
Y_obs = c(-4, -4.5),
Y_se = c(0.6, 0.3),
N = 2)A Bayesian model-based network meta-analysis
A common pharmacometric dose–response model is the Emax model (see e.g. Danielle Navarro’s and Kristian Brock’s smooth introductions to the Emax model), which says the response at dose \(D\) follows
\[\begin{aligned} \mu(D) = \frac{E_{\max} \times D}{ED_{50} + D} \end{aligned}\]where \(E_{\max}\) is the asymptotic effect at a very large dose and \(ED_{50}\) is the dose at which half of \(E_{\max}\) is achieved. We assume that zero dose gives zero placebo-adjusted effect (i.e., the intercept \(E_0\) is fixed at 0).
Combining the Emax structure with a meta-analytic likelihood gives the model:
\[\begin{aligned} Y_{\text{obs}, i} &\sim \text{Normal}(Y_{\text{true}, i},\ Y_{\text{se}, i}) \\ Y_{\text{true}, i} &\sim \text{Normal}(\mu_i,\ \sigma) \\ \mu_i &= \frac{E_{\max} \times D_i}{ED_{50} + D_i} \\ E_{\max} &\sim \text{Normal}(-5.5,\ 0.75) \\ \log ED_{50} &\sim \text{Normal}(-1.5,\ 1) \\ \sigma &\sim \text{Exponential}(4) \end{aligned}\]The first two lines are a standard random-effects meta-analysis: For the effect of the two dose levels of Drug A, the observed point estimates \(Y_{\text{obs}, i}\) are normally distributed around unobserved true study means \(Y_{\text{true}, i}\), with each trial’s reported standard error \(Y_{\text{se}, i}\) reflecting sampling variability. The true means themselves come from a population of true means with mean \(\mu\) and between-trial heterogeneity \(\sigma\). This notation comes from McElreath (2020, ch. 15), who frames it as a measurement error model: the observed estimates are measurements of the true study-level effects but they are measured with error. In the third line, \(\mu\) is given its own structural shape, namely the Emax model.
This gives us three parameters – the \(E_{\max}\) and \(ED_{50}\) as well as the between-trial heterogeneity \(\sigma\) – to estimate from two data points. But this is no problem for Bayes. The minimum sample size for a Bayesian analysis is 0 (we’ll just get the prior back). But it does mean that the priors will exert some influence.
Setting priors
We have only two observations on the Drug A dose–response curve, so the priors will matter. In practice, we’d need some clinical guidance and/or previous studies to qualify these choices, but since this is a simulated example, we’ll just go with the flow.
Recall that \(E_{\max}\) is the asymptotic effect at very large doses. It should be at least as large in magnitude as the largest finite-dose effect observed. The largest observed estimate in our data is \(-4.5\), so a prior mean of \(-5.5\) allows the asymptote to sit modestly beyond observed effects. The SD of 0.75 acknowledges that the true ceiling could be somewhat larger or smaller, but probably not much smaller than what we’ve already observed.
\(ED_{50}\), the dose that achieves half of \(E_{\max}\), has to be positive (we cannot have negative doses), so we model it on the log scale. A prior mean of \(-1.5\) on the log scale puts most mass near \(\exp(-1.5) \approx 0.22\) mg, which is roughly a fifth of our lowest studied dose – a reasonable starting guess for a drug whose two studied doses already produce “statistically significant” effects against placebo. The SD of 1 keeps this prior fairly weak.
Finally, \(\sigma\). An exponential keeps \(\sigma\) non-negative and a rate of 4 gives a prior mean of 0.25 on the effect-measure scale, encoding our expectation that two trials of the same drug at well-characterised doses should not disagree much once dose is accounted for.
While we can quibble about these exact prior choices, let’s simulate from them to find out if they give reasonable prior predictions. In a real analysis we would want to do a sensitivity analysis to the priors, especially the prior on \(\sigma\) given the very limited data on between-trial variability.
Prior predictive check
Before fitting, it’s worth seeing what these priors actually imply for the dose–response curve. Below I draw 30 Emax curves with parameters sampled from the priors:
Show the code
library(ggplot2)
set.seed(2)
n_draws <- 30
D_grid <- seq(0, 3, length.out = 200)
draws <- lapply(1:n_draws, function(i) {
emax <- rnorm(1, -5.5, 0.75)
log_ed50 <- rnorm(1, -1.5, 1)
mu <- (emax * D_grid) / (exp(log_ed50) + D_grid)
data.frame(D = D_grid, y = mu, draw = i)
})
prior_df <- do.call(rbind, draws)
mbnma_prior <- ggplot(prior_df, aes(D, y, group = draw)) +
geom_line(alpha = 0.3) +
labs(x = "Dose (mg)", y = "Placebo-adjusted effect") +
scale_y_continuous(breaks = seq(0, -6, by = -1)) +
coord_cartesian(ylim = c(-6, 0)) +
theme_classic()
mbnma_prior
The curves cover a fairly wide range of plausible Emax shapes, but they’re mostly in a sensible ballpark for the effect-size scale.
Fitting the model
library(rethinking)
m_mbnma <- ulam(
alist(
Y_obs ~ normal(Y_true, Y_se),
vector[N]:Y_true ~ normal(mu, sigma),
mu <- (Emax * D) / (exp(logED50) + D),
Emax ~ normal(-5.5, 0.75),
logED50 ~ normal(-1.5, 1),
sigma ~ exponential(4)
),
data = dlist,
cores = 4, chains = 4, iter = 4000, seed = 1,
control = list(adapt_delta = 0.99, max_treedepth = 12),
constraints = list(Emax = "upper=0", sigma = "lower=0"))Note how the model syntax in alist() satisfyingly mirrors the formal model from above. Now, because we have only two data points, I help the sampler along by raising adapt_delta and max_treedepth and constraining Emax to be negative and sigma positive. The model gives a warning about some divergent transitions, which we’ll come back to.
divergent(m_mbnma)[1] 109First, a quick look at the fit:
precis(m_mbnma, depth = 2, prob = 0.95) mean sd 2.5% 97.5% rhat ess_bulk
Y_true[1] -4.0502476 0.3981446 -4.80726125 -3.2771823 1.002179 1985.2710
Y_true[2] -4.5328263 0.2603689 -5.04302150 -4.0334485 1.001766 2314.2516
Emax -5.2437998 0.5360207 -6.38889150 -4.3385252 1.000851 1627.1591
logED50 -1.4398116 0.7752004 -3.17093025 -0.1906518 1.004403 1604.3959
sigma 0.1988313 0.1894714 0.01331812 0.7166844 1.004692 809.3079Sampling looks fine: \(\hat R\) close to 1 and ESS (“effective sample size” for the sampler) in the high hundreds or thousands. We get a posterior estimate for each study mean (Y_true[1] for the 1 mg trial, Y_true[2] for the 2 mg trial) plus the Emax parameters.
Prior vs. posterior
A useful diagnostic is to compare priors against posteriors for each parameter. Below, the dashed curve is the prior, and the filled blue density is the posterior.
Show the code
post <- extract.samples(m_mbnma)
n <- length(post$Emax)
prior_Emax <- pmin(rnorm(n, -5.5, 0.75), 0) # respect upper=0 constraint
prior_ED50 <- exp(rnorm(n, -1.5, 1))
prior_sigma <- rexp(n, 4)
plot_prior_post <- function(prior, post_samples, title, xlim, add_legend = FALSE) {
d_prior <- density(prior, adjust = 1.5)
d_post <- density(post_samples, adjust = 1.5)
ylim <- c(0, max(d_prior$y, d_post$y) * 1.1)
plot(NA, xlim = xlim, ylim = ylim,
main = title, xlab = "", ylab = "Density",
bty = "l")
# Posterior: filled steelblue polygon, no border
polygon(c(d_post$x[1], d_post$x, d_post$x[length(d_post$x)]),
c(0, d_post$y, 0),
col = adjustcolor("steelblue", alpha.f = 0.4),
border = NA)
# Prior: dashed line
lines(d_prior, col = "black", lwd = 2, lty = 2)
# Free-form labels on first panel only
if (add_legend) {
idx_prior <- which.max(d_prior$y)
text(d_prior$x[idx_prior] * 1.1,
d_prior$y[idx_prior] * 1.02,
"Prior", col = "black", cex = 1.1)
idx_post <- which.max(d_post$y)
text(d_post$x[idx_post],
d_post$y[idx_post] * 1.075,
"Posterior", col = "steelblue", font = 2, cex = 1.1)
}
}
par(mfrow = c(1, 3), mar = c(4, 4, 3, 1))
plot_prior_post(prior_Emax, post$Emax, expression(bold(E[max])),
xlim = range(c(prior_Emax, post$Emax)), add_legend = TRUE)
plot_prior_post(prior_ED50, exp(post$logED50), expression(bold(ED[50])),
xlim = c(0, 1.5))
plot_prior_post(prior_sigma, post$sigma, expression(bold(sigma)),
xlim = c(0, 1))
The model has learned about \(E_{\max}\) and \(\log ED_{50}\) – the posteriors are tighter than the priors. It has learned less about \(\sigma\), which is unsurprising: with only two data points, between-study heterogeneity is hard to pin down.
The Emax curve and shrinkage
Let’s plot the posterior Emax curve with the observed data overlaid, alongside the prior predictive from above for comparison:
Show the code
library(patchwork)
dose_seq <- seq(0, 3, by = 0.025)
mu <- link(m_mbnma, data = list(D = dose_seq))
mu_median <- apply(mu, 2, median)
mu_PI95 <- apply(mu, 2, PI, prob = 0.95)
plot_data <- data.frame(
dose = dose_seq,
median = mu_median,
lower_95 = mu_PI95[1, ],
upper_95 = mu_PI95[2, ]
)
obs_data <- data.frame(
dose = c(1, 2),
Y_obs = c(-4, -4.5),
Y_sd = c(0.6, 0.3)
)
mbnma_post <- ggplot() +
# 95% CI ribbon
geom_ribbon(data = plot_data,
aes(x = dose, ymin = lower_95, ymax = upper_95),
fill = "steelblue", alpha = 0.4) +
# Posterior median line
geom_line(data = plot_data,
aes(x = dose, y = median),
color = "black", linewidth = 1) +
# Observed data with error bars
geom_errorbar(data = obs_data,
aes(x = dose,
ymin = Y_obs - 1.96 * Y_sd,
ymax = Y_obs + 1.96 * Y_sd),
width = 0.08, linewidth = 0.5) +
geom_point(data = obs_data,
aes(x = dose, y = Y_obs),
size = 2, shape = 19) +
labs(x = "Dose (mg)", y = "Placebo-adjusted effect") +
scale_y_continuous(breaks = seq(0, -6, by = -1)) +
scale_x_continuous(breaks = c(0, 1, 2)) +
coord_cartesian(ylim = c(-6, 0)) +
theme_classic()
mbnma_prior + ggtitle("Prior predictive") + mbnma_post + ggtitle("Posterior")
Two things are worth noting. First, we see that we did manage to increase precision on the target dose of 1 mg. The posterior interval at 1 mg is narrower than the empirical estimates, reflecting the information borrowed from the dose-response curve. Second, our inference on 2 mg remains largely unchanged – borrowing information from 1 mg does not help us much in getting a more precise estimate for 2 mg. This is because the 2 mg trial was a much larger trial than the 1 mg trial and so there’s less to learn from the latter.
An improved indirect treatment comparison
Now, let’s query the model about the posterior predictions for Drug A on 1 mg, so that we can plug this updated estimate into our indirect treatment comparison against Drug B 1 mg.
mu_1mg <- link(m_mbnma, data = list(D = 1))The object mu_1mg holds the posterior distribution of predictions for the mean placebo-adjusted effect of 1 mg. We can summarise the distribution using precis().
precis(data.frame(mu_1mg), prob = 0.95) mean sd 2.5% 97.5% histogram
mu_1mg -4.059585 0.3930425 -4.800853 -3.284729 ▁▁▁▂▇▅▁▁▁Since we now have a full posterior distribution at our disposal, we can do better than the Bucher above. Rather than subtracting point estimates and adding variances algebraically, we approximate the posterior of the Drug B 1 mg estimate as a normal distribution and subtract it sample-by-sample from our posterior draws of the updated Drug A 1 mg estimate. This propagates the full uncertainty from both sources of information and gives us a proper posterior interval for the difference. It also handles any potential skewness in mu_1mg.
delta_mbnma <- mu_1mg - rnorm(length(mu_1mg), mean = -3, sd = 0.2)
data.frame(
delta = round(mean(delta_mbnma), 1),
ci_lo = round(quantile(delta_mbnma, 0.025), 1),
ci_up = round(quantile(delta_mbnma, 0.975), 1)) delta ci_lo ci_up
2.5% -1.1 -1.9 -0.2The point estimate has barely moved – the Bucher gave \(-1\), the MBNMA gives roughly \(-1.1\) – but the intervals are narrower, such that the upper bound of the interval is below 0. We’re now more certain that Drug A is more effective than Drug B on average. How certain are we? With a full posterior on hand, we can easily compute the probability that the indirect treatment estimate favors Drug A.
round(mean(delta_mbnma < 0), 2)[1] 0.99Based on this model and data, the posterior probability that Drug A is more effective than Drug B is around 99%.
Non-centered parameterisation
Now, we have unfinished business with a few divergent transitions. Divergent transitions are basically signals that the sampler had difficulties exploring parts of the posterior. We are ready to tackle those now.
The model outlined above is what is known as a “centered” parameterisation: each study’s true mean \(Y_{\text{true}, i}\) is sampled directly from the population distribution:
\[\begin{aligned} Y_{\text{true}, i} \sim \text{Normal}(\mu, \sigma) \end{aligned}\]This means \(Y_{\text{true}, 1}\) and \(Y_{\text{true}, 2}\) are not independent – they both depend on \(\sigma\), so draws of one tend to move with draws of the other. A pairs plot confirms this: the two study-level parameters show a moderate correlation of around \(r = 0.5\) in the centered model. This correlated geometry is harder for the sampler to explore efficiently (Gelman et al. 2013, ch. 13 and 15). This is the source of our divergent transitions.
Show the code
pairs(m_mbnma)
The fix is what is sometimes referred to as “non-centered” parameterisation (McElreath 2020, ch. 14). Instead of sampling \(Y_{\text{true}, i}\) directly, we introduce a standardised offset \(z_i\):
\[\begin{aligned} z_i \sim \text{Normal}(0, 1) \end{aligned}\]and reconstruct \(Y_{\text{true}, i}\) using that offset:
\[\begin{aligned} Y_{\text{true}, i} = \mu_i + \sigma \times z_i \end{aligned}\]Now \(z_1\) and \(z_2\) are independent by construction – \(\sigma\) only enters later, when computing \(Y_{\text{true}, i}\). Why is this called “non-centered” parameterisation? One way to think of it is with reference to the normal distribution. In the centered parameterisation, the study-level parameters are placed directly at their natural scale – they are centered on their mean \(\mu\) with spread \(\sigma\). In the non-centered version, we instead sample dimensionless standard normal offsets \(z \sim \text{Normal}(0, 1)\) and reconstruct the original parameters as \(\mu_i + \sigma \times z_i\), effectively removing the center and scale from the sampling geometry. Anyway, the important thing is that it works. The full non-centered model becomes:
\[\begin{aligned} Y_{\text{obs}, i} &\sim \text{Normal}(Y_{\text{true}, i},\ Y_{\text{se}, i}) \\ Y_{\text{true}, i} &= \mu_i + \sigma \times z_i \\ \mu_i &= \frac{E_{\max} \times D_i}{ED_{50} + D_i} \\ z_i &\sim \text{Normal}(0, 1) \\ E_{\max} &\sim \text{Normal}(-5.5,\ 0.75) \\ \log ED_{50} &\sim \text{Normal}(-1.5,\ 1) \\ \sigma &\sim \text{Exponential}(4) \end{aligned}\]The key change from the centered model is the second line: \(Y_{\text{true}, i}\) goes from being sampled (\(\sim\)) to being deterministically defined (\(=\)), as well as the addition of the vector of z’s that are given a standardised normal prior. In code:
m_mbnma_nc <- ulam(
alist(
Y_obs ~ normal(Y_true, Y_se),
# non-centered: Y_true is now deterministic
vector[N]:Y_true <- mu + sigma * z[i],
mu <- (Emax * D) / (exp(logED50) + D),
vector[N]:z ~ normal(0, 1),
Emax ~ normal(-5.5, 0.75),
logED50 ~ normal(-1.5, 1),
sigma ~ exponential(4)
),
data = m_mbnma@data,
cores = 4, chains = 4, iter = 4000, seed = 1,
control = list(adapt_delta = 0.99, max_treedepth = 12),
constraints = list(Emax = "upper=0", sigma = "lower=0"))We’ll first ensure that parameter estimates are the same for the two models – looks good.
precis(m_mbnma) mean sd 5.5% 94.5% rhat ess_bulk
Emax -5.2437998 0.5360207 -6.18276045 -4.4672451 1.000851 1627.1591
logED50 -1.4398116 0.7752004 -2.82573045 -0.3523874 1.004403 1604.3959
sigma 0.1988313 0.1894714 0.01850621 0.5557643 1.004692 809.3079precis(m_mbnma_nc) mean sd 5.5% 94.5% rhat ess_bulk
Emax -5.2562210 0.5415828 -6.20665815 -4.4865684 1.001386 3945.348
logED50 -1.4051758 0.7634635 -2.77494120 -0.3594211 1.001368 3606.681
sigma 0.1983168 0.1892859 0.01289558 0.5406912 1.000882 3676.308Let’s then check the number of divergent transitions with the non-centered model.
divergent(m_mbnma_nc)[1] 0Zero divergences. And if we check the pairs plot for the non-centered model, this shows why: the correlation between \(z_1\) and \(z_2\) is now roughly zero.
Show the code
pairs(m_mbnma_nc)
The non-centered form of a multilevel model tends to work better when data are sparse, there are few groups, and \(\sigma\) is small – like in our example here. On the other hand, when the data are rich and informative, the centered form works just fine. The non-centered model is very much worth knowing about, and it can be an elegant solution to thorny sampling issues – McElreath (2020, ch. 14) goes into quite a bit of detail, including a clear illustration of the funnel geometry and the conditions under which each parameterisation is preferable. The non-centered parameterisation is the default in other popular R interfaces to Stan, such as rstanarm (Goodrich et al. 2025) and brms (Bürkner 2017).
Summary
Exploiting a dose-response structure rather than ignoring it is one way to get more out of an evidence network. The same principle extends naturally to other dimensions along which trials differ. One such example is follow-up time: if the same treatment has been evaluated at different follow-up times across a set of trials, a longitudinal dose-response-style model can borrow information across time points in the same way we borrowed it across doses (Strathe, Bøg, and Gorst-Rasmussen 2026). Another extension is “multilevel network meta-regression,” which combines individual patient data with aggregate trial data to adjust indirect comparisons for differences in patient-level effect modifiers across trials (Phillippo et al. 2020).
All in all, Bayes is a natural fit for these kinds of models, in that it allows uncertainty to propagate in a principled manner through evidence networks. The resulting posterior allows for direct, interpretable, and decision-relevant probability statements about the relative efficacy and safety of competing clinical interventions.
If you liked this post, you might also want to check out this brief introduction to “Bayesian dynamic borrowing”.
Software
library(grateful)
pkgs <- cite_packages(output = "table", out.dir = ".", cite.tidyverse = TRUE, dependencies = TRUE)
knitr::kable(pkgs)| Package | Version | Citation |
|---|---|---|
| base | 4.4.1 | R Core Team (2024b) |
| base64enc | 0.1.3 | Urbanek (2015) |
| bslib | 0.8.0 | Sievert, Cheng, and Aden-Buie (2024) |
| cachem | 1.1.0 | Chang (2024a) |
| colorspace | 2.1.1 | Zeileis, Hornik, and Murrell (2009); Stauffer et al. (2009); Zeileis et al. (2020) |
| digest | 0.6.37 | Eddelbuettel (2024) |
| evaluate | 0.24.0 | H. Wickham and Xie (2024) |
| fansi | 1.0.6 | Gaslam (2023) |
| farver | 2.1.2 | Pedersen, Nicolae, and François (2024) |
| fastmap | 1.2.0 | Chang (2024b) |
| fontawesome | 0.5.2 | Iannone (2023) |
| fs | 1.6.4 | Hester, Wickham, and Csárdi (2024) |
| glue | 1.8.0 | Hester and Bryan (2024) |
| gtable | 0.3.6 | H. Wickham and Pedersen (2024) |
| highr | 0.11 | Xie and Qiu (2024) |
| htmltools | 0.5.8.1 | Cheng, Sievert, et al. (2024) |
| isoband | 0.2.7 | H. Wickham, Wilke, and Pedersen (2022) |
| jquerylib | 0.1.4 | Sievert and Cheng (2021) |
| knitr | 1.48 | Xie (2014); Xie (2015); Xie (2024a) |
| labeling | 0.4.3 | Justin Talbot (2023) |
| lifecycle | 1.0.4 | Henry and Wickham (2023) |
| memoise | 2.0.1 | H. Wickham et al. (2021) |
| mime | 0.12 | Xie (2021) |
| munsell | 0.5.1 | C. Wickham (2024) |
| patchwork | 1.2.0 | Pedersen (2024) |
| pkgconfig | 2.0.3 | Csárdi (2019) |
| R6 | 2.6.1 | Chang (2025) |
| rappdirs | 0.3.3 | Ratnakumar, Mick, and Davis (2021) |
| RColorBrewer | 1.1.3 | Neuwirth (2022) |
| remotes | 2.5.0 | Csárdi et al. (2024) |
| renv | 1.0.7 | Ushey and Wickham (2024) |
| rethinking | 2.42 | McElreath (2024) |
| rmarkdown | 2.27 | Xie, Allaire, and Grolemund (2018); Xie, Dervieux, and Riederer (2020); Allaire et al. (2024) |
| sass | 0.4.9 | Cheng, Mastny, et al. (2024) |
| scales | 1.4.0 | H. Wickham, Pedersen, and Seidel (2025) |
| tidyverse | 2.0.0 | H. Wickham et al. (2019) |
| tinytex | 0.52 | Xie (2019); Xie (2024b) |
| utf8 | 1.2.6 | Perry (2025) |
| vctrs | 0.6.5 | H. Wickham, Henry, and Vaughan (2023) |
| viridisLite | 0.4.2 | Garnier et al. (2023) |
| withr | 3.0.2 | Hester et al. (2024) |
| xfun | 0.46 | Xie (2024c) |
| yaml | 2.3.10 | Garbett et al. (2024) |
References
Allaire, JJ, Yihui Xie, Christophe Dervieux, Jonathan McPherson, Javier Luraschi, Kevin Ushey, Aron Atkins, et al. 2024. rmarkdown: Dynamic Documents for r. https://github.com/rstudio/rmarkdown.
Bucher, Heiner C, Gordon H Guyatt, Lauren E Griffith, and Stephen D Walter. 1997. “The Results of Direct and Indirect Treatment Comparisons in Meta-Analysis of Randomized Controlled Trials.” Journal of Clinical Epidemiology 50 (6): 683–91.
Bürkner, Paul-Christian. 2017. “Brms: An r Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80: 1–28.
Carpenter, Bob, Andrew Gelman, Matthew D. Hoffman, Daniel Lee, Ben Goodrich, Michael Betancourt, Marcus Brubaker, Jiqiang Guo, Peter Li, and Allen Riddell. 2017. “Stan: A Probabilistic Programming Language.” Journal of Statistical Software 76: 1–32.
Chang, Winston. 2024a. cachem: Cache r Objects with Automatic Pruning. https://cachem.r-lib.org/.
———. 2024b. fastmap: Fast Data Structures. https://r-lib.github.io/fastmap/.
———. 2025. R6: Encapsulated Classes with Reference Semantics. https://CRAN.R-project.org/package=R6.
Cheng, Joe, Timothy Mastny, Richard Iannone, Barret Schloerke, and Carson Sievert. 2024. sass: Syntactically Awesome Style Sheets (“Sass”). https://rstudio.github.io/sass/.
Cheng, Joe, Carson Sievert, Barret Schloerke, Winston Chang, Yihui Xie, and Jeff Allen. 2024. htmltools: Tools for HTML. https://github.com/rstudio/htmltools.
Csárdi, Gábor. 2019. pkgconfig: Private Configuration for “R” Packages. https://CRAN.R-project.org/package=pkgconfig.
Csárdi, Gábor, Jim Hester, Hadley Wickham, Winston Chang, Martin Morgan, and Dan Tenenbaum. 2024. remotes: R Package Installation from Remote Repositories, Including “GitHub”. https://remotes.r-lib.org.
Eddelbuettel, Dirk. 2024. digest: Create Compact Hash Digests of r Objects. https://CRAN.R-project.org/package=digest.
Garbett, Shawn P, Jeremy Stephens, Kirill Simonov, Yihui Xie, Zhuoer Dong, Hadley Wickham, Jeffrey Horner, et al. 2024. yaml: Methods to Convert r Data to YAML and Back. https://github.com/vubiostat/r-yaml/.
Garnier, Simon, Ross, Noam, Rudis, Robert, Camargo, et al. 2023. viridis(Lite) - Colorblind-Friendly Color Maps for r. https://doi.org/10.5281/zenodo.4678327.
Gaslam, Brodie. 2023. fansi: ANSI Control Sequence Aware String Functions. https://github.com/brodieG/fansi.
Gelman, Andrew, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari, and Donald B. Rubin. 2013. Bayesian Data Analysis. 3rd ed. Chapman; Hall/CRC.
Goodrich, Ben, Jonah Gabry, Imad Ali, and Sam Brilleman. 2025. Rstanarm: Bayesian Applied Regression Modeling via Stan. https://mc-stan.org/rstanarm/.
Henry, Lionel, and Hadley Wickham. 2023. lifecycle: Manage the Life Cycle of Your Package Functions. https://CRAN.R-project.org/package=lifecycle.
Hester, Jim, and Jennifer Bryan. 2024. glue: Interpreted String Literals. https://CRAN.R-project.org/package=glue.
Hester, Jim, Lionel Henry, Kirill Müller, Kevin Ushey, Hadley Wickham, and Winston Chang. 2024. withr: Run Code “With” Temporarily Modified Global State. https://CRAN.R-project.org/package=withr.
Hester, Jim, Hadley Wickham, and Gábor Csárdi. 2024. fs: Cross-Platform File System Operations Based on “libuv”. https://fs.r-lib.org.
Iannone, Richard. 2023. fontawesome: Easily Work with “Font Awesome” Icons. https://github.com/rstudio/fontawesome.
Justin Talbot. 2023. labeling: Axis Labeling. https://CRAN.R-project.org/package=labeling.
Lu, Guobing, and AE15449338 Ades. 2004. “Combination of Direct and Indirect Evidence in Mixed Treatment Comparisons.” Statistics in Medicine 23 (20): 3105–24.
Lumley, Thomas. 2002. “Network Meta-Analysis for Indirect Treatment Comparisons.” Statistics in Medicine 21 (16): 2313–24.
Mawdsley, David, Meg Bennetts, Sofia Dias, Martin Boucher, and Nicky J Welton. 2016. “Model-Based Network Meta-Analysis: A Framework for Evidence Synthesis of Clinical Trial Data.” CPT: Pharmacometrics & Systems Pharmacology 5 (8): 393–401.
McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in r and Stan. Second. CRC Press.
———. 2024. rethinking: Statistical Rethinking Book Package. https://github.com/rmcelreath/rethinking.
Neuwirth, Erich. 2022. RColorBrewer: ColorBrewer Palettes. https://CRAN.R-project.org/package=RColorBrewer.
Pedersen, Thomas Lin. 2024. patchwork: The Composer of Plots. https://patchwork.data-imaginist.com.
Pedersen, Thomas Lin, Berendea Nicolae, and Romain François. 2024. farver: High Performance Colour Space Manipulation. https://CRAN.R-project.org/package=farver.
Perry, Patrick O. 2025. Utf8: Unicode Text Processing. https://CRAN.R-project.org/package=utf8.
Phillippo, David M, Sofia Dias, AE Ades, Mark Belger, Alan Brnabic, Alexander Schacht, Daniel Saure, Zbigniew Kadziola, and Nicky J Welton. 2020. “Multilevel Network Meta-Regression for Population-Adjusted Treatment Comparisons.” Journal of the Royal Statistical Society Series A: Statistics in Society 183 (3): 1189–1210.
R Core Team. 2024a. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
———. 2024b. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.
Ratnakumar, Sridhar, Trent Mick, and Trevor Davis. 2021. rappdirs: Application Directories: Determine Where to Save Data, Caches, and Logs. https://rappdirs.r-lib.org.
Sievert, Carson, and Joe Cheng. 2021. jquerylib: Obtain “jQuery” as an HTML Dependency Object.
Sievert, Carson, Joe Cheng, and Garrick Aden-Buie. 2024. bslib: Custom “Bootstrap” “Sass” Themes for “shiny” and “rmarkdown”. https://rstudio.github.io/bslib/.
Stauffer, Reto, Georg J. Mayr, Markus Dabernig, and Achim Zeileis. 2009. “Somewhere over the Rainbow: How to Make Effective Use of Colors in Meteorological Visualizations.” Bulletin of the American Meteorological Society 96 (2): 203–16. https://doi.org/10.1175/BAMS-D-13-00155.1.
Strathe, Anders, Martin Bøg, and Anders Gorst-Rasmussen. 2026. “Model-Based Network Meta-Analysis: Joint Estimation of Dose–Response and Time–Course Relationships.” Research Synthesis Methods, 1–20. https://doi.org/10.1017/rsm.2026.10098.
Urbanek, Simon. 2015. Base64enc: Tools for Base64 Encoding. http://www.rforge.net/base64enc.
Ushey, Kevin, and Hadley Wickham. 2024. renv: Project Environments. https://rstudio.github.io/renv/.
Wickham, Charlotte. 2024. munsell: Utilities for Using Munsell Colours. https://cran.r-project.org/package=munsell.
Wickham, Hadley, Mara Averick, Jennifer Bryan, Winston Chang, Lucy D’Agostino McGowan, Romain François, Garrett Grolemund, et al. 2019. “Welcome to the tidyverse.” Journal of Open Source Software 4 (43): 1686. https://doi.org/10.21105/joss.01686.
Wickham, Hadley, Lionel Henry, and Davis Vaughan. 2023. vctrs: Vector Helpers. https://CRAN.R-project.org/package=vctrs.
Wickham, Hadley, Jim Hester, Winston Chang, Kirill Müller, and Daniel Cook. 2021. memoise: “Memoisation” of Functions. https://memoise.r-lib.org.
Wickham, Hadley, and Thomas Lin Pedersen. 2024. gtable: Arrange “Grobs” in Tables. https://CRAN.R-project.org/package=gtable.
Wickham, Hadley, Thomas Lin Pedersen, and Dana Seidel. 2025. scales: Scale Functions for Visualization. https://CRAN.R-project.org/package=scales.
Wickham, Hadley, Claus O. Wilke, and Thomas Lin Pedersen. 2022. isoband: Generate Isolines and Isobands from Regularly Spaced Elevation Grids. https://CRAN.R-project.org/package=isoband.
Wickham, Hadley, and Yihui Xie. 2024. evaluate: Parsing and Evaluation Tools That Provide More Details Than the Default. https://github.com/r-lib/evaluate.
Xie, Yihui. 2014. “knitr: A Comprehensive Tool for Reproducible Research in R.” In Implementing Reproducible Computational Research, edited by Victoria Stodden, Friedrich Leisch, and Roger D. Peng. Chapman; Hall/CRC.
———. 2015. Dynamic Documents with R and Knitr. 2nd ed. Boca Raton, Florida: Chapman; Hall/CRC. https://yihui.org/knitr/.
———. 2019. “TinyTeX: A Lightweight, Cross-Platform, and Easy-to-Maintain LaTeX Distribution Based on TeX Live.” TUGboat 40 (1): 30–32. https://tug.org/TUGboat/Contents/contents40-1.html.
———. 2021. mime: Map Filenames to MIME Types. https://github.com/yihui/mime.
———. 2024a. knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.
———. 2024b. tinytex: Helper Functions to Install and Maintain TeX Live, and Compile LaTeX Documents. https://github.com/rstudio/tinytex.
———. 2024c. xfun: Supporting Functions for Packages Maintained by “Yihui Xie”. https://github.com/yihui/xfun.
Xie, Yihui, J. J. Allaire, and Garrett Grolemund. 2018. R Markdown: The Definitive Guide. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown.
Xie, Yihui, Christophe Dervieux, and Emily Riederer. 2020. R Markdown Cookbook. Boca Raton, Florida: Chapman; Hall/CRC. https://bookdown.org/yihui/rmarkdown-cookbook.
Xie, Yihui, and Yixuan Qiu. 2024. highr: Syntax Highlighting for r Source Code. https://github.com/yihui/highr.
Zeileis, Achim, Jason C. Fisher, Kurt Hornik, Ross Ihaka, Claire D. McWhite, Paul Murrell, Reto Stauffer, and Claus O. Wilke. 2020. “colorspace: A Toolbox for Manipulating and Assessing Colors and Palettes.” Journal of Statistical Software 96 (1): 1–49. https://doi.org/10.18637/jss.v096.i01.
Zeileis, Achim, Kurt Hornik, and Paul Murrell. 2009. “Escaping RGBland: Selecting Colors for Statistical Graphics.” Computational Statistics & Data Analysis 53 (9): 3259–70. https://doi.org/10.1016/j.csda.2008.11.033.
Session info
sessionInfo()R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Rocky Linux 8.10 (Green Obsidian)
Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblasp-r0.3.15.so; LAPACK version 3.9.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Europe/Copenhagen
tzcode source: system (glibc)
attached base packages:
[1] parallel stats graphics grDevices utils datasets methods
[8] base
other attached packages:
[1] grateful_0.2.4 patchwork_1.2.0 digest_0.6.37
[4] rethinking_2.42 posterior_1.6.0 cmdstanr_0.9.0.9000
[7] ggplot2_3.5.1
loaded via a namespace (and not attached):
[1] tensorA_0.36.2.1 generics_0.1.4 renv_1.0.7
[4] KernSmooth_2.23-24 shape_1.4.6.1 lattice_0.22-6
[7] magrittr_2.0.4 evaluate_0.24.0 grid_4.4.1
[10] RColorBrewer_1.1-3 mvtnorm_1.2-5 fastmap_1.2.0
[13] jsonlite_2.0.0 processx_3.8.4 backports_1.5.0
[16] ps_1.7.7 scales_1.4.0 codetools_0.2-20
[19] abind_1.4-5 cli_3.6.5 rlang_1.1.6
[22] remotes_2.5.0 withr_3.0.2 yaml_2.3.10
[25] tools_4.4.1 checkmate_2.3.2 coda_0.19-4.1
[28] dplyr_1.1.4 curl_5.2.1 vctrs_0.6.5
[31] R6_2.6.1 matrixStats_1.3.0 lifecycle_1.0.4
[34] htmlwidgets_1.6.4 MASS_7.3-61 pkgconfig_2.0.3
[37] pillar_1.11.1 gtable_0.3.6 loo_2.8.0
[40] glue_1.8.0 data.table_1.17.8 xfun_0.46
[43] tibble_3.3.0 tidyselect_1.2.1 rstudioapi_0.17.1
[46] knitr_1.48 farver_2.1.2 htmltools_0.5.8.1
[49] rmarkdown_2.27 labeling_0.4.3 compiler_4.4.1
[52] distributional_0.4.0