Skip to contents

Fit a Bayesian multilevel unanchored meta-regression model using individual patient data (IPD) and aggregate data (AgD). Supports binary, continuous, and count outcomes.

Usage

mlumr(
  data,
  model = c("spfa", "relaxed"),
  link = NULL,
  prior_intercept = default_prior_intercept(),
  prior_beta = default_prior_beta(),
  prior_sigma = default_prior_sigma(),
  chains = 4,
  iter = 2000,
  warmup = 1000,
  seed = NULL,
  adapt_delta = 0.95,
  max_treedepth = 15,
  refresh = 200,
  engine = NULL,
  verbose = TRUE,
  ...
)

Arguments

data

An mlumr_data object with integration points (from add_integration())

model

Model type: "spfa" (shared prognostic factor assumption) or "relaxed" (treatment-specific coefficients). Default "spfa".

Link function. For binomial: "logit" (default), "probit", or "cloglog". For normal: "identity" (default) or "log". For poisson: "log" (default, only option). If NULL, uses the canonical default for the family.

prior_intercept

Prior for treatment intercepts. Default from default_prior_intercept() (prior_normal(0, 10), weakly informative on the linear-predictor scale). See prior_normal() for guidance.

prior_beta

Prior for regression coefficients. May be a single prior broadcast to all covariates, or a list of priors of length n_cov for per-coefficient specification. All per-coefficient priors must share the same family and (for Student-t) df. Default from default_prior_beta() (prior_normal(0, 2.5); Gelman et al. 2008; Stan community prior-choice wiki). Set autoscale = TRUE on the prior to divide the scale by each covariate's empirical SD — useful when predictors are on very different scales.

prior_sigma

Prior for residual SD (normal family only). Default from default_prior_sigma() (prior_normal(0, 2.5), half-normal via the Stan <lower=0> constraint). prior_exponential() is also supported for sigma.

chains

Number of MCMC chains (default 4)

iter

Total iterations per chain (default 2000)

warmup

Number of warmup iterations (default 1000)

seed

Random seed for reproducibility

adapt_delta

Target acceptance rate (default 0.95)

max_treedepth

Maximum tree depth for NUTS (default 15)

refresh

How often to print progress (0 = silent, default 200)

engine

Stan backend: "rstan" (default) or "cmdstanr". If NULL, uses the engine set by mlumr_engine(). See mlumr_engine() for setup.

verbose

Logical; if FALSE, suppresses mlumr progress messages. Stan sampler progress is still controlled by refresh.

...

Additional arguments passed to the Stan sampling function (rstan::sampling() or cmdstanr's $sample() method)

Value

An object of class mlumr_fit

Details

The model assumes that all AgD rows come from the same comparator treatment and that, conditional on covariates, there is no between-study heterogeneity. If AgD rows come from multiple studies with different designs or unmeasured confounders, this assumption may not hold. No random effects for study-level heterogeneity are included.

AgD scale assumptions (family = "normal"). The AgD likelihood is y_agd ~ normal(E[exp(eta)], se_agd) under link = "log" and y_agd ~ normal(E[eta], se_agd) under link = "identity". In both cases set_agd() expects outcome_mean and outcome_se on the arithmetic (original, untransformed) scale, not log-scale or geometric. Passing log-scale summaries silently misspecifies the likelihood. See set_agd() for details.

Comparator-population weighting is family-dependent. Integrated marginal predictions in the comparator population (*_comparator generated quantities) are weighted by:

  • binomial: n_agd[k] (AgD sample size), so larger AgD rows contribute more to the marginal mean.

  • normal: equal weights across AgD rows (/ n_agd_rows), reflecting the normal likelihood's treatment of each row as a single summary statistic.

  • poisson: E_agd[k] (AgD exposure), matching the rate-based likelihood.

Each weighting is natural for the corresponding likelihood; users comparing marginal effects across families should be aware they are not identically weighted.

Weakly-identified coefficients in the relaxed modelbeta_comparator is identified only through AgD, so the relaxed model needs informative priors (or many AgD rows) to estimate effect modification reliably. prior_sensitivity() is the recommended diagnostic.

See also

prior_sensitivity() for sensitivity of the posterior to prior_beta; set_agd() for AgD scale requirements; prior_summary() for introspection of the priors actually used.

Examples

if (FALSE) { # \dontrun{
# Binary SPFA model
fit_spfa <- mlumr(dat, model = "spfa")

# Relaxed SPFA (allows effect modification)
fit_relaxed <- mlumr(dat, model = "relaxed")
} # }