Package {mlumr}

Title:	Multilevel Unanchored Meta-Regression for Indirect Treatment Comparisons
Version:	0.1.0
Description:	Bayesian multilevel unanchored meta-regression (ML-UMR) for indirect treatment comparisons using individual patient data (IPD) and aggregate data (AgD). Implements shared prognostic factor assumption (SPFA) and relaxed SPFA models for binary, continuous, and count outcomes via 'Stan'. Also provides simulated treatment comparison (STC) via parametric G-computation and naive unadjusted benchmarks. ML-UMR is an adaptation of the ML-NMR methodology (Phillippo et al. 2020, <doi:10.1111/rssa.12579>) implemented in the 'multinma' package (GPL-3) to the unanchored two-trial case; the public API deliberately mirrors multinma's so users can move between ML-NMR and ML-UMR with the same workflow.
License:	GPL-3
URL:	https://github.com/choxos/mlumr, https://choxos.github.io/mlumr/
BugReports:	https://github.com/choxos/mlumr/issues
Encoding:	UTF-8
Biarch:	true
Depends:	R (≥ 4.1.0)
Imports:	methods, stats, Rcpp (≥ 0.12.0), RcppParallel (≥ 5.0.1), rstan (≥ 2.26.0), rstantools (≥ 2.3.0), randtoolbox, copula
LinkingTo:	BH (≥ 1.66.0), Rcpp (≥ 0.12.0), RcppEigen (≥ 0.3.3.3.0), RcppParallel (≥ 5.0.1), rstan (≥ 2.26.0), StanHeaders (≥ 2.26.0)
Suggests:	testthat (≥ 3.0.0), knitr, rmarkdown, posterior, bayesplot, loo, Matrix, cmdstanr, withr
Additional_repositories:	https://stan-dev.r-universe.dev
SystemRequirements:	GNU make
Config/testthat/edition:	3
VignetteBuilder:	knitr
Config/roxygen2/version:	8.0.0
NeedsCompilation:	yes
Packaged:	2026-05-15 19:14:56 UTC; choxos
Author:	Ahmad Sofi-Mahmudi [aut, cre], Conor Chandler [aut]
Maintainer:	Ahmad Sofi-Mahmudi <a.sofimahmudi@gmail.com>
Repository:	CRAN
Date/Publication:	2026-05-20 08:30:02 UTC

mlumr: Multilevel Unanchored Meta-Regression for Indirect Treatment Comparisons

Description

Bayesian multilevel unanchored meta-regression (ML-UMR) for population- adjusted indirect comparisons between a single-arm individual patient data (IPD) study on the index treatment and aggregate data (AgD) on a comparator. Two model variants are provided: a shared-prognostic-factor (SPFA) model and a relaxed-SPFA model that allows treatment-specific covariate coefficients. Supports binary (binomial), continuous (normal), and count (Poisson) outcomes, each with appropriate link functions. Frequentist outcome regression (Simulated Treatment Comparison) and a naive unadjusted benchmark are included for side-by-side comparison.

Typical workflow

1. set_ipd() — declare the individual patient data and its covariates.

2. set_agd() — declare the aggregate-data arm.

3. combine_data() — join the two into an mlumr_data object.

4. add_integration() — build the Gaussian-copula QMC integration points used to marginalize over the AgD covariate distribution.

5. mlumr() — fit the Bayesian ML-UMR model (SPFA or relaxed).

6. prior_summary(), check_integration() — sanity-check the model before inferring from it.

7. predict(), marginal_effects(), conditional_effects(), conditional_predict() — extract population-level and profile-specific quantities.

8. calculate_dic(), calculate_loo(), calculate_waic(), compare_models() — compare SPFA vs relaxed or competing specifications.

9. prior_sensitivity() — verify the posterior is data-driven rather than prior-driven.

Priors

The prior constructors prior_normal(), prior_student_t(), prior_cauchy(), and prior_exponential() all plug into mlumr() via prior_intercept, prior_beta, and (normal family only) prior_sigma. Defaults follow the Stan community's prior-choice recommendations (Vehtari et al., 2025).

Alternative methods

stc() runs frequentist parametric G-computation (Simulated Treatment Comparison). naive() computes an unadjusted crude-rate benchmark. Both consume the same mlumr_data object as mlumr() so their estimands can be compared directly.

Author(s)

Maintainer: Ahmad Sofi-Mahmudi a.sofimahmudi@gmail.com (ORCID)

Authors:

• Ahmad Sofi-Mahmudi a.sofimahmudi@gmail.com (ORCID)

• Conor Chandler (ORCID)

See Also

Useful links:

• https://github.com/choxos/mlumr

• Report bugs at https://github.com/choxos/mlumr/issues

Adjust correlations to Gaussian-copula scale

Description

Adjust correlations to Gaussian-copula scale

Usage

.adjust_integration_cor(cor, cor_adjust, dtypes)

Build AgD covariate metadata

Description

Build AgD covariate metadata

Usage

.agd_cov_info(cov_names, cov_sds, cov_types)

Normalize AgD covariate specification

Description

Normalize AgD covariate specification

Usage

.agd_covariate_spec(cov_means, cov_sds = NULL, cov_types = NULL)

Summarize AgD outcomes

Description

Summarize AgD outcomes

Usage

.agd_outcome_summary(
  data,
  family,
  outcome_n,
  outcome_r,
  outcome_mean,
  outcome_se,
  outcome_E
)

AgD source columns required for setup

Description

AgD source columns required for setup

Usage

.agd_required_columns(
  treatment,
  cov_means,
  cov_sds,
  outcome_n = NULL,
  outcome_r = NULL,
  outcome_mean = NULL,
  outcome_se = NULL,
  outcome_E = NULL,
  study = NULL
)

Attach generated integration points to an mlumr_data object

Description

Attach generated integration points to an mlumr_data object

Usage

.attach_integration_points(
  data,
  X_int_array,
  n_int,
  cor,
  copula_cor,
  cor_adjust
)

Bound probabilities to the open unit interval

Description

Bound probabilities to the open unit interval

Usage

.bound_unit_interval(p)

Bound a Wald interval to a valid numerical range

Description

Bound a Wald interval to a valid numerical range

Usage

.bounded_wald_interval(center, se, z, lower = -Inf, upper = Inf)

Validate that required columns are present

Description

Validate that required columns are present

Usage

.check_required_columns(data, required_cols)

Check whether CmdStan is configured for cmdstanr

Description

Check whether CmdStan is configured for cmdstanr

Usage

.cmdstan_available()

Resolve a writable cmdstanr cache root

Description

Resolve a writable cmdstanr cache root

Usage

.cmdstanr_cache_root()

Cache directory for cmdstanr-compiled model executables

Description

Cache directory for cmdstanr-compiled model executables

Usage

.cmdstanr_compile_dir(model_name, stan_file)

Use user-supplied comparison names when available

Description

Use user-supplied comparison names when available

Usage

.comparison_names(models, fallback)

Compute mean linear predictors from parameter draws

Description

Internal helper for predict.mlumr_fit with type="link". Computes mean(eta_i) directly from parameter draws instead of applying the link to marginalized response-scale draws, avoiding Jensen's inequality bias.

Usage

.compute_mean_lp(object, pred_cols)

Arguments

Value

Data frame with the same columns as pred_cols, on the link scale

Conditional effect choices by family

Description

Conditional effect choices by family

Usage

.conditional_effect_choices(family)

Compute conditional linear predictors for one profile

Description

Compute conditional linear predictors for one profile

Usage

.conditional_eta(params, x)

Extract posterior parameter draws for conditional summaries

Description

Extract posterior parameter draws for conditional summaries

Usage

.conditional_parameters(object, covariates)

Build covariate profiles for conditional summaries

Description

Build covariate profiles for conditional summaries

Usage

.conditional_profiles(object, newdata = NULL)

Read a non-negative scalar diagnostic count

Description

Read a non-negative scalar diagnostic count

Usage

.diagnostic_count(x)

Format a diagnostic setting for warning messages

Description

Format a diagnostic setting for warning messages

Usage

.diagnostic_value(x)

Drop rows with missing setup inputs

Description

Drop rows with missing setup inputs

Usage

.drop_missing_rows(data, required_cols)

Ensure adjusted integration correlation is positive definite

Description

Ensure adjusted integration correlation is positive definite

Usage

.ensure_positive_definite_cor(copula_cor, n_cov)

Return finite numeric values from a summary column

Description

Return finite numeric values from a summary column

Usage

.finite_numeric_values(x)

Generate correlated uniform quasi-Monte Carlo points

Description

Generate correlated uniform quasi-Monte Carlo points

Usage

.generate_copula_uniforms(n_int, n_cov, copula_cor)

Pairwise-correlation diagnostics for integration points

Description

Pairwise-correlation diagnostics for integration points

Usage

.int_cor_stats(X_orig, X_double, cov_names, n_agd, cor_target = NULL)

Compute summary statistics for integration points

Description

Compute summary statistics for integration points

Usage

.int_stats(X_int, cov_names, n_agd)

Summarize standardized IPD outcomes

Description

Summarize standardized IPD outcomes

Usage

.ipd_outcome_summary(ipd_data, family)

Required source columns for IPD setup

Description

Required source columns for IPD setup

Usage

.ipd_required_columns(
  treatment,
  outcome,
  covariates,
  exposure = NULL,
  study = NULL
)

Check whether values are valid non-negative whole-number counts

Description

Check whether values are valid non-negative whole-number counts

Usage

.is_whole_number_count(x, allow_missing = FALSE)

Extract integer indexes from Stan vector column names

Description

Extract integer indexes from Stan vector column names

Usage

.log_lik_column_index(cols, source)

Report that the current engine was not changed

Description

Report that the current engine was not changed

Usage

.message_engine_unchanged()

Build the Stan data list for mlumr()

Description

Build the Stan data list for mlumr()

Usage

.mlumr_build_stan_data(
  data,
  family,
  link_info,
  prior_intercept,
  prior_beta,
  prior_sigma
)

Dispatch mlumr() sampling to the selected backend

Description

Dispatch mlumr() sampling to the selected backend

Usage

.mlumr_fit_backend(
  engine,
  model_name,
  stan_data,
  chains,
  iter,
  warmup,
  seed,
  adapt_delta,
  max_treedepth,
  refresh,
  ...
)

Log mlumr() fit metadata

Description

Log mlumr() fit metadata

Usage

.mlumr_log_fit_start(model_name, family, link, stan_data, engine, verbose)

Human-readable model label for an mlumr fit or DIC object

Description

Human-readable model label for an mlumr fit or DIC object

Usage

.mlumr_model_label(object)

Store user priors plus the resolved Stan-scale beta prior

Description

Store user priors plus the resolved Stan-scale beta prior

Usage

.mlumr_prior_metadata(
  data,
  family,
  prior_intercept,
  prior_beta,
  prior_sigma,
  beta_fields,
  sd_x
)

Naive comparison for binomial outcomes

Description

Naive comparison for binomial outcomes

Usage

.naive_binomial(data, ipd, agd, link, conf_level, z)

Naive comparison for normal outcomes

Description

Naive comparison for normal outcomes

Usage

.naive_normal(data, ipd, agd, conf_level, z)

Naive comparison for Poisson outcomes

Description

Naive comparison for Poisson outcomes

Usage

.naive_poisson(data, ipd, agd, conf_level, z)

Truncate tiny negative variance estimates caused by numerical noise

Description

Truncate tiny negative variance estimates caused by numerical noise

Usage

.nonnegative_variance(x, name = "variance", tol = 1e-10)

Ordered pointwise log-likelihood columns for one data source

Description

Ordered pointwise log-likelihood columns for one data source

Usage

.ordered_log_lik_columns(draws, source)

Compute stable relative effective sample sizes from log-likelihoods

Description

Compute stable relative effective sample sizes from log-likelihoods

Usage

.relative_eff_from_log_lik(log_lik, chain_id)

Require posterior draw columns

Description

Require posterior draw columns

Usage

.require_draw_columns(draws, columns, context)

Rescale a prior_beta's scale, preserving family / mean / df.

Description

Rescale a prior_beta's scale, preserving family / mean / df.

Usage

.rescale_prior_beta(prior, new_scale)

Resolve raw and copula correlation matrices for integration

Description

Resolve raw and copula correlation matrices for integration

Usage

.resolve_integration_cor(data, cov_names, n_cov, cor, cor_adjust, ds, verbose)

Resolve and validate mlumr() backend engine

Description

Resolve and validate mlumr() backend engine

Usage

.resolve_mlumr_engine(engine)

Format a homogeneous resolved beta prior

Description

Format a homogeneous resolved beta prior

Usage

.resolved_prior_broadcast_label(br, digits)

Label a resolved Stan prior family code

Description

Label a resolved Stan prior family code

Usage

.resolved_prior_family_label(dist)

Ratio with the same positive-denominator guard used by Stan generated quantities

Description

Ratio with the same positive-denominator guard used by Stan generated quantities

Usage

.safe_positive_ratio(num, denom)

Square root of a variance estimate with numerical guarding

Description

Square root of a variance estimate with numerical guarding

Usage

.sqrt_variance(x, name = "variance", tol = 1e-10)

Standardize AgD to mlumr's internal column contract

Description

Standardize AgD to mlumr's internal column contract

Usage

.standardize_agd_data(
  data,
  treatment,
  family,
  outcome_n,
  outcome_r,
  outcome_mean,
  outcome_se,
  outcome_E,
  cov_means,
  cov_sds,
  cov_names,
  study = NULL
)

Standardize IPD to mlumr's internal column contract

Description

Standardize IPD to mlumr's internal column contract

Usage

.standardize_ipd_data(
  data,
  treatment,
  outcome,
  covariates,
  family,
  exposure = NULL,
  study = NULL
)

Build comparator-population covariates from AgD means

Description

Build comparator-population covariates from AgD means

Usage

.stc_agd_mean_newdata(agd, cov_names)

Binomial STC estimator

Description

Binomial STC estimator

Usage

.stc_binomial(
  data,
  fit,
  ipd,
  agd,
  newdata,
  link_resolved,
  conf_level,
  z,
  beta_hat,
  V,
  n_int
)

Comparator-population delta-method terms for binomial STC

Description

Comparator-population delta-method terms for binomial STC

Usage

.stc_binomial_comparator_delta(
  fit,
  newdata,
  weights,
  beta_hat,
  V,
  link_resolved,
  p_A,
  p_B,
  n_B
)

Index-population delta-method terms for binomial STC

Description

Index-population delta-method terms for binomial STC

Usage

.stc_binomial_index_delta(
  fit,
  ipd,
  beta_hat,
  V,
  link_resolved,
  p_A_comp,
  p_B_comp,
  comp_delta
)

Build comparator-population covariates for STC prediction

Description

Build comparator-population covariates for STC prediction

Usage

.stc_comparator_data(data, cov_names, family)

Build the STC GLM formula without assuming syntactic covariate names

Description

Build the STC GLM formula without assuming syntactic covariate names

Usage

.stc_formula(cov_names, family)

Validate fitted STC GLM parameters before delta-method calculations

Description

Validate fitted STC GLM parameters before delta-method calculations

Usage

.stc_glm_parameters(fit)

Build a model matrix aligned with fitted GLM coefficients

Description

Build a model matrix aligned with fitted GLM coefficients

Usage

.stc_model_matrix(fit, newdata)

Normal-outcome STC estimator

Description

Normal-outcome STC estimator

Usage

.stc_normal(
  data,
  fit,
  ipd,
  agd,
  newdata,
  link_resolved,
  conf_level,
  z,
  beta_hat,
  V,
  n_int
)

Poisson-outcome STC estimator

Description

Poisson-outcome STC estimator

Usage

.stc_poisson(
  data,
  fit,
  ipd,
  agd,
  newdata,
  link_resolved,
  conf_level,
  z,
  beta_hat,
  V,
  n_int
)

Strip AgD covariate suffixes

Description

Strip AgD covariate suffixes

Usage

.strip_agd_cov_suffix(cov_means)

Summarize a draws matrix column-wise into a tidy data frame.

Description

Applies .summarize_draw_vector() across columns of draws and renames the quantile columns to qNN form (e.g., q2.5, q50, q97.5).

Usage

.summarize_draw_matrix(draws, probs)

Arguments

Value

Data frame with columns mean, sd, and one qNN column per element of probs.

Summarize a single draw vector into mean / sd / quantiles.

Description

Internal helper. Centralizes the (mean, sd, quantile) triplet used by predict.mlumr_fit(), marginal_effects(), conditional_effects() and conditional_predict() so that any later change to the canonical posterior summary only needs to happen in one place.

Usage

.summarize_draw_vector(x, probs)

Arguments

Value

Named numeric vector: ⁠c(mean, sd, <named quantiles>)⁠.

Summarize a sensitivity refit.

Description

Summarize a sensitivity refit.

Usage

.summarize_sensitivity(fit, scale, probs)

Supported Stan engines

Description

Supported Stan engines

Usage

.supported_engines()

Transform uniform integration points to covariate scales

Description

Transform uniform integration points to covariate scales

Usage

.transform_integration_points(u_cor, data, ds, cov_names, n_int, n_cov)

Validate binary AgD covariate summaries

Description

Validate binary AgD covariate summaries

Usage

.validate_agd_binary_covariates(data, cov_means, cov_sds, cov_types)

Validate binomial AgD outcomes

Description

Validate binomial AgD outcomes

Usage

.validate_agd_binomial_outcomes(r_vals, n_vals)

Validate AgD covariate type labels

Description

Validate AgD covariate type labels

Usage

.validate_agd_cov_types(cov_types)

Validate AgD covariate names after suffix stripping

Description

Validate AgD covariate names after suffix stripping

Usage

.validate_agd_covariate_names(cov_means)

Validate AgD covariate summaries

Description

Validate AgD covariate summaries

Usage

.validate_agd_covariates(data, cov_means, cov_sds)

Validate normal AgD outcomes

Description

Validate normal AgD outcomes

Usage

.validate_agd_normal_outcomes(y_vals, se_vals)

Validate AgD outcome column arguments

Description

Validate AgD outcome column arguments

Usage

.validate_agd_outcome_args(
  family,
  outcome_n,
  outcome_r,
  outcome_mean,
  outcome_se,
  outcome_E
)

Validate AgD outcome summaries

Description

Validate AgD outcome summaries

Usage

.validate_agd_outcomes(
  data,
  family,
  outcome_n,
  outcome_r,
  outcome_mean,
  outcome_se,
  outcome_E
)

Validate Poisson AgD outcomes

Description

Validate Poisson AgD outcomes

Usage

.validate_agd_poisson_outcomes(r_vals, E_vals)

Validate optional AgD sample sizes

Description

Validate optional AgD sample sizes

Usage

.validate_agd_sample_size(n_vals)

Validate that complete-case filtering left rows to analyze

Description

Validate that complete-case filtering left rows to analyze

Usage

.validate_complete_rows_remain(data, label)

Validate an IPD-derived correlation matrix before copula adjustment

Description

Validate an IPD-derived correlation matrix before copula adjustment

Usage

.validate_computed_integration_cor(cor)

Validate a confidence level

Description

Validate a confidence level

Usage

.validate_conf_level(conf_level)

Validate a diagnostics/model-comparison choice

Description

Validate a diagnostics/model-comparison choice

Usage

.validate_diagnostic_choice(x, choices, name)

Validate a scalar marginal-effect choice

Description

Validate a scalar marginal-effect choice

Usage

.validate_effect_choice(effect)

Validate a Stan engine name

Description

Validate a Stan engine name

Usage

.validate_engine_name(engine)

Validate top-level integration arguments

Description

Validate top-level integration arguments

Usage

.validate_integration_args(data, n_int, cor_adjust, verbose, ds)

Validate a user-supplied integration correlation matrix

Description

Validate a user-supplied integration correlation matrix

Usage

.validate_integration_cor(cor, n_cov, cov_names = NULL)

Validate integration distributions against the combined data

Description

Validate integration distributions against the combined data

Usage

.validate_integration_distributions(ds, data)

Validate IPD covariates

Description

Validate IPD covariates

Usage

.validate_ipd_covariates(data, covariates)

Validate finite IPD outcome and exposure values

Description

Validate finite IPD outcome and exposure values

Usage

.validate_ipd_finite_columns(data, outcome, exposure, family)

Validate IPD outcome and exposure columns

Description

Validate IPD outcome and exposure columns

Usage

.validate_ipd_outcome(data, outcome, family, exposure = NULL)

Validate a scalar link/family string

Description

Validate a scalar link/family string

Usage

.validate_link_string(x, name)

Validate a pointwise log-likelihood matrix

Description

Validate a pointwise log-likelihood matrix

Usage

.validate_log_lik_matrix(log_lik)

Validate an mlumr_data object

Description

Validate an mlumr_data object

Usage

.validate_mlumr_data_object(data)

Validate an mlumr fit object

Description

Validate an mlumr fit object

Usage

.validate_mlumr_fit_object(object)

Validate an integer-like mlumr() argument

Description

Validate an integer-like mlumr() argument

Usage

.validate_mlumr_integer(x, name, lower)

Validate mlumr() sampler controls before backend dispatch

Description

Validate mlumr() sampler controls before backend dispatch

Usage

.validate_mlumr_sampling_args(
  chains,
  iter,
  warmup,
  seed,
  adapt_delta,
  max_treedepth,
  refresh
)

Validate that model comparison has at least two candidates

Description

Validate that model comparison has at least two candidates

Usage

.validate_model_count(models)

Validate that input data contain at least one row

Description

Validate that input data contain at least one row

Usage

.validate_non_empty_data(data, label)

Validate a numeric vector

Description

Validate a numeric vector

Usage

.validate_numeric_vector(x, name)

Validate positive finite numeric input

Description

Validate positive finite numeric input

Usage

.validate_positive_numeric(x, name)

Validate an exact scalar prediction argument

Description

Validate an exact scalar prediction argument

Usage

.validate_predict_choice(x, choices, name)

Validate prior_summary digits

Description

Validate prior_summary digits

Usage

.validate_prior_summary_digits(digits)

Validate quantile probabilities

Description

Validate quantile probabilities

Usage

.validate_probs(probs)

Validate covariate argument shape before column lookup

Description

Validate covariate argument shape before column lookup

Usage

.validate_required_covariates(covariates, arg_name)

Reject user names that would overwrite standardized internal columns

Description

Reject user names that would overwrite standardized internal columns

Usage

.validate_reserved_internal_names(user_names, reserved, label)

Validate resolved beta-prior metadata stored on a fit

Description

Validate resolved beta-prior metadata stored on a fit

Usage

.validate_resolved_beta_prior(br)

Validate a data source contains only one treatment

Description

Validate a data source contains only one treatment

Usage

.validate_single_treatment(data, treatment, label)

Validate a scalar summary flag

Description

Validate a scalar summary flag

Usage

.validate_summary_flag(summary)

Warn when IPD covariates have no empirical variation

Description

Warn when IPD covariates have no empirical variation

Usage

.warn_constant_ipd_covariates(data, covariates)

Warn when integration resolution is low for the covariate dimension

Description

Warn when integration resolution is low for the covariate dimension

Usage

.warn_integration_size(n_int, n_cov)

Convert a confidence level to a two-sided normal critical value

Description

Convert a confidence level to a two-sided normal critical value

Usage

.z_from_conf_level(conf_level)

Add numerical integration points

Description

Generate quasi-Monte Carlo integration points using Sobol sequences and a Gaussian copula to account for correlations between covariates in the AgD.

Usage

add_integration(
  data,
  n_int = 64,
  cor = NULL,
  cor_adjust = NULL,
  verbose = TRUE,
  ...
)

Arguments

Value

An mlumr_data object with integration points added

Examples

## Not run: 
dat <- add_integration(
  dat,
  n_int = 64,
  x1 = distr(qnorm, mean = x1_mean, sd = x1_sd),
  x2 = distr(qbern, prob = x2_mean)
)

## End(Not run)

Delta-method variance for a transformed binomial proportion

Description

Delta-method variance for a transformed binomial proportion

Usage

binomial_link_variance(p, n, link = c("logit", "probit", "cloglog"))

Bound probabilities away from 0 and 1

Description

Bound probabilities away from 0 and 1

Usage

bound_probability(p, n, min_count = 0.5)

Calculate DIC for model comparison

Description

Computes the Deviance Information Criterion using the variance-based effective-parameters formula from Gelman et al. (2004): pD = 0.5 * Var(D). This is more stable than the plug-in alternative for multimodal posteriors.

Usage

calculate_dic(object)

Arguments

Details

DIC is retained for backward compatibility and rough comparison. For principled Bayesian model comparison, prefer calculate_loo() or calculate_waic() (Vehtari, Gelman, Gabry 2017).

Value

A list of class mlumr_dic with components DIC, pD, D_bar

Examples

## Not run: 
dic_spfa <- calculate_dic(fit_spfa)
dic_relaxed <- calculate_dic(fit_relaxed)

## End(Not run)

Calculate LOO-CV for an mlumr_fit

Description

Computes approximate leave-one-out cross-validation (PSIS-LOO, Vehtari, Gelman, Gabry 2017) using the pointwise log-likelihoods stored by the Stan models. Returns a loo object from the loo package.

Usage

calculate_loo(object, ...)

Arguments

Details

Pareto-k diagnostics: values > 0.7 indicate observations for which the PSIS approximation is unreliable; the printed output flags these. Typical remedies are running more iterations, using moment_match = TRUE, or (for highly influential AgD rows) refitting without the offending observation to check sensitivity.

Value

An object of class psis_loo (see loo::loo()).

Note

AgD rows are treated as independent observations. Each AgD row contributes one column to the pointwise log_lik matrix. If two or more AgD rows come from the same study (e.g. subgroup summaries within a single trial) the PSIS-LOO approximation does not account for the within-study clustering; effective sample sizes are inflated and Pareto-k warnings are understated. For clustered AgD, corroborate with prior_sensitivity() or refit omitting suspect rows to check the influence on the posterior.

Examples

## Not run: 
loo_spfa <- calculate_loo(fit_spfa)
print(loo_spfa)

## End(Not run)

Calculate WAIC for an mlumr_fit

Description

Watanabe-Akaike Information Criterion (Watanabe 2010) based on the pointwise log-likelihoods. WAIC is asymptotically equivalent to LOO-CV; prefer calculate_loo() when Pareto-k is well-behaved.

Usage

calculate_waic(object, ...)

Arguments

Value

An object of class waic (see loo::waic()).

Note

As with calculate_loo(), each AgD row is treated as an independent observation. WAIC will be optimistic for AgD rows that share a study (see the note on calculate_loo()).

Examples

## Not run: 
waic_spfa <- calculate_waic(fit_spfa)

## End(Not run)

Check MCMC diagnostics and warn if issues found

Description

Check MCMC diagnostics and warn if issues found

Usage

check_diagnostics(fit)

Arguments

Check integration point adequacy

Description

Compare integration results at the current n_int against a doubled resolution to assess numerical accuracy. Large discrepancies indicate that n_int should be increased.

Usage

check_integration(
  data,
  ...,
  cor = NULL,
  cor_adjust = NULL,
  check_joint = TRUE,
  verbose = TRUE
)

Arguments

Value

A list with components marginals (the original data frame returned by previous versions) and, if check_joint = TRUE, correlations – a data frame of pairwise covariate correlations at the current and doubled n_int for each AgD row. Printed with a pass/warn verdict.

Validate and resolve link function for a given family

Description

Checks that link is valid for family and returns the resolved link name plus an integer code for Stan. Accepts canonical family names ("binomial", "normal", "poisson") and data-type aliases ("binary", "count", "rate", "continuous").

Usage

check_link(family, link = NULL)

Arguments

Details

The V1 likelihood/link matrix is:

*"count" refers to count/total (binomial denominator) data, not Poisson event counts. For Poisson rate or count outcomes, use "poisson" or "rate".

Value

List with components:

family

Canonical family name (e.g. "binomial")

link

Resolved link name (e.g. "probit")

code

Integer code for Stan data block

Combine IPD and AgD for unanchored comparison

Description

Combine IPD and AgD for unanchored comparison

Usage

combine_data(ipd, agd)

Arguments

Value

An object of class mlumr_data

Examples

## Not run: 
dat <- combine_data(ipd, agd)

## End(Not run)

Compare fitted ML-UMR models

Description

Compare two or more mlumr_fit objects by DIC (default), LOO, or WAIC. For LOO/WAIC, loo::loo_compare() is used under the hood; the output is the standard loo_compare table. For DIC the return is a data frame ordered by DIC.

Usage

compare_models(..., criterion = c("dic", "loo", "waic"))

Arguments

Details

DIC is the default for backward compatibility and because it has no additional package dependencies. For principled Bayesian model comparison, LOO (Vehtari, Gelman, Gabry 2017) is preferred and requires the optional loo package.

Value

For "loo" / "waic": a compare.loo table from loo::loo_compare(). For "dic": a data frame (invisibly) with columns Model, DIC, pD, Delta_DIC.

Conditional treatment effects

Description

Compute conditional (individual-level) treatment effects at specific covariate values from a fitted ML-UMR model. Unlike marginal_effects(), which averages over a population's covariate distribution, conditional effects evaluate the treatment effect at a particular covariate profile.

Usage

conditional_effects(
  object,
  newdata = NULL,
  effect = "all",
  summary = TRUE,
  probs = c(0.025, 0.5, 0.975)
)

Arguments

Details

For SPFA models, the conditional link-scale treatment effect is constant across all covariate values because the shared beta cancels in the treatment contrast on the fitted link scale. However, risk difference (RD) and risk ratio (RR) still vary with covariates because they depend on absolute probability levels. For relaxed models, all conditional effects vary with covariate values because the index and comparator treatments have different regression coefficients.

The conditional link-scale effect is computed directly as eta_index - eta_comparator. This avoids numerical distortion from transforming extreme response-scale probabilities back through the link function.

Conditional vs marginal on non-identity links. Conditional effects are evaluated at a single covariate profile, so there is no averaging over a population and no Jensen's-inequality gap between the conditional and marginal response. Compare with marginal_effects() and predict.mlumr_fit(), which average over either the IPD individuals (index population) or the AgD integration points (comparator population) and therefore return E[g^{-1}(eta)], not g^{-1}(E[eta]).

Value

A data frame. If summary = TRUE, contains columns profile, effect, mean, sd, and quantile columns. If summary = FALSE, returns a single combined data frame of full posterior draws with a profile column indicating which covariate profile each draw belongs to.

See Also

marginal_effects() for population-averaged treatment effects; conditional_predict() for absolute predictions at specific profiles; predict.mlumr_fit() for population-level predictions.

Examples

## Not run: 
# Conditional effects at IPD covariate means (default)
conditional_effects(fit)

# At specific covariate values
conditional_effects(fit, newdata = data.frame(age = 60, sex = 1))

# Multiple profiles
profiles <- data.frame(age = c(50, 60, 70), sex = c(0, 0, 1))
conditional_effects(fit, newdata = profiles)

## End(Not run)

Conditional predictions

Description

Generate absolute predictions at specific covariate values for both treatments.

Usage

conditional_predict(
  object,
  newdata = NULL,
  type = c("response", "link"),
  summary = TRUE,
  probs = c(0.025, 0.5, 0.975)
)

Arguments

Value

A data frame with predictions for each treatment at each profile

See Also

conditional_effects() for covariate-conditional treatment effects; predict.mlumr_fit() for population-level predictions.

Examples

## Not run: 
conditional_predict(fit)
conditional_predict(fit, newdata = data.frame(age = 60, sex = 1))

## End(Not run)

Convert Pearson correlations to Gaussian copula correlations

Description

For continuous-continuous pairs, Pearson rho equals the Gaussian copula parameter only under normality of both margins. For non-normal continuous covariates, this no-adjustment assumption introduces approximation error. For binary and mixed pairs:

• Binary-binary: rho_copula = sin(pi * rho_P / 2)

• Continuous-binary: rho_copula = sqrt(pi/2) * rho_P

Usage

cor_adjust_pearson(X, types)

Arguments

Value

Adjusted correlation matrix for Gaussian copula

Convert Spearman correlations to Gaussian copula correlations

Description

Applies theoretical relationships between Spearman's rho and the Gaussian copula parameter (Kurowicka & Cooke, 2006; Lebrun & Dutfoy, 2009):

• Continuous-continuous: rho_copula = 2 * sin(pi * rho_S / 6) (exact)

• Binary-binary: rho_copula = sin(pi * rho_S / 2) (heuristic; the exact relationship depends on marginal prevalences, not accounted for here)

• Continuous-binary: rho_copula = sqrt(2) * sin(pi * rho_S / (2*sqrt(3))) (heuristic)

Usage

cor_adjust_spearman(X, types)

Arguments

Value

Adjusted correlation matrix for Gaussian copula

Bernoulli PMF

Description

Bernoulli PMF

Usage

dbern(x, prob, log = FALSE)

Arguments

Value

Numeric vector

Default priors used by mlumr()

Description

These accessors return the current default priors used by mlumr(), tagged with ⁠$default = TRUE⁠ and the package version. prior_summary() prints the version so cross-release reproducibility is diagnosable: if a later release changes a default, fits produced with an older version will still carry the correct ⁠$version⁠ tag.

Usage

default_prior_intercept()

default_prior_beta()

default_prior_sigma()

Value

A prior list (see prior_normal()).

Examples

default_prior_intercept()
default_prior_beta()
default_prior_sigma()

Specify a marginal distribution

Description

Used to specify marginal distributions for covariates when adding integration points. Wraps an inverse CDF (quantile) function with its parameters.

Usage

distr(qfun, ...)

Arguments

Value

A list with class "mlumr_distr" containing the distribution specification

Examples

# Normal distribution
distr(qnorm, mean = 0, sd = 1)

# Bernoulli distribution with probability 0.3
distr(qbern, prob = 0.3)

Evaluate a mlumr_distr object with data context

Description

Evaluate a mlumr_distr object with data context

Usage

eval_distr(d, p, data = list())

Arguments

Value

Numeric vector of quantiles

Evaluate a single mlumr_distr argument expression

Description

Evaluate a single mlumr_distr argument expression

Usage

eval_distr_arg(expr, data)

Arguments

Value

Evaluated value

Extract the full pointwise log-likelihood matrix from an mlumr_fit

Description

Combines the IPD and AgD per-observation log-likelihood draws into a single matrix with one row per posterior draw and one column per observation. The result is suitable for direct use with loo::loo() / loo::waic().

Usage

extract_log_lik(object)

Arguments

Value

A numeric matrix of dimension ⁠n_draws x (n_ipd + n_agd_rows)⁠. IPD columns come first, then AgD rows.

Family metadata registry

Description

Internal single-source-of-truth for family-specific Stan model names, AgD weighting, prediction-variable prefixes, and supported links and effect measures. Every R-side call site that hard-coded a per-family branch now looks up the relevant field here. This file is deliberately pure-R and makes no Stan calls; keeping it a data registry makes rstantools regeneration of R/stanmodels.R independent.

Details

Fields:

stan_prefix

Prefix for the Stan model name (the full name is ⁠<stan_prefix>_{spfa,relaxed}⁠).

predict_prefix

Column prefix for generated-quantity variables in predict.mlumr_fit() (e.g. "p", "y", "rate").

link_default

The default link when the user passes link = NULL.

links

Vector of supported links. Should match the branches in check_link().

effect_measures

Supported values of the effect argument in marginal_effects() (excluding "all").

marginal_effect_vars

Generated-quantity column names for each effect measure, per population. Expanded in marginal_effects().

comp_weight_field

Name of the Stan-data field used to weight the comparator-population marginal predictions (n_agd for binomial, E_agd for poisson, NULL for normal = equal weights).

Fit a Stan model using cmdstanr

Description

Fit a Stan model using cmdstanr

Usage

fit_cmdstanr(
  model_name,
  stan_data,
  chains,
  iter,
  warmup,
  seed,
  adapt_delta,
  max_treedepth,
  refresh,
  ...
)

Fit a Stan model using rstan

Description

Fit a Stan model using rstan

Usage

fit_rstan(
  model_name,
  stan_data,
  chains,
  iter,
  warmup,
  seed,
  adapt_delta,
  max_treedepth,
  refresh,
  ...
)

Get distribution type (continuous, discrete, or binary)

Description

Get distribution type (continuous, discrete, or binary)

Usage

get_distribution_type(..., data = list())

Arguments

Value

Named character vector

Get the current Stan engine (internal)

Description

Get the current Stan engine (internal)

Usage

get_engine()

Lookup helper for the family registry

Description

Fails fast with an informative message when a caller requests an unregistered family, rather than returning NULL and letting a downstream $ chain produce a cryptic error.

Usage

get_family_config(family)

Arguments

Value

The corresponding entry from family_config.

Inverse link function (linear predictor -> response scale)

Description

Inverse link function (linear predictor -> response scale)

Usage

inverse_link(x, link = c("identity", "log", "logit", "probit", "cloglog"))

Arguments

Value

Numeric vector on response scale

Derivative of inverse link with respect to linear predictor

Description

Derivative of inverse link with respect to linear predictor

Usage

inverse_link_derivative(
  eta,
  p = inverse_link(eta, link),
  link = c("logit", "probit", "cloglog")
)

Is this object a single prior (vs a list of priors)?

Description

Is this object a single prior (vs a list of priors)?

Usage

is_single_prior(x)

Derivative of a binomial link with respect to probability

Description

Derivative of a binomial link with respect to probability

Usage

link_derivative_response(p, link = c("logit", "probit", "cloglog"))

Link function (response scale -> linear predictor)

Description

Link function (response scale -> linear predictor)

Usage

link_fun(x, link = c("identity", "log", "logit", "probit", "cloglog"))

Arguments

Value

Numeric vector on linear predictor scale

Marginal treatment effects

Description

Extract marginal treatment effects from a fitted ML-UMR model. For binomial: log odds ratio, risk difference, risk ratio. For normal: mean difference. For poisson: rate ratio.

Usage

marginal_effects(
  object,
  population = c("both", "index", "comparator"),
  effect = "all",
  summary = TRUE,
  probs = c(0.025, 0.5, 0.975)
)

Arguments

Value

A data frame

See Also

predict.mlumr_fit() for absolute predictions; conditional_effects() for covariate-conditional effects at specific profiles; prior_sensitivity() to check how strongly the marginal effect depends on prior_beta.

Examples

## Not run: 
# All effect measures for both populations
marginal_effects(fit)

# Only the log odds ratio in the index population
marginal_effects(fit, population = "index", effect = "lor")

# Full posterior draws rather than summary statistics
marginal_effects(fit, summary = FALSE)

## End(Not run)

Fit ML-UMR Model

Description

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

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 model — beta_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.

Value

An object of class mlumr_fit

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

## Not run: 
# Binary SPFA model
fit_spfa <- mlumr(dat, model = "spfa")

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

## End(Not run)

Numerical guards in the Stan models

Description

mlumr's Stan models clip a few quantities to avoid undefined arithmetic in ⁠generated quantities⁠. This page documents the exact behavior so users can assess when the guard matters for their analysis.

Guards currently in place

safe_logit(p) (binomial models)

Returns logit(clamp(p, 1e-10, 1 - 1e-10)). The logit and log-odds-ratio posterior summaries use this rather than the raw logit(p) so that finite samples with p at the boundaries do not propagate Inf/NaN into the posterior.

safe_divide(num, denom) (binomial and Poisson models)

Returns num / max(denom, 1e-10). Used for risk ratios and rate ratios to avoid division-by-zero when a comparator posterior draw has a near-zero mean outcome.

Stan ⁠<lower=0>⁠ constraints on sigma (normal family)

Truncate the prior to the positive half-line — the reason why prior_sigma = prior_normal(0, 2.5) is interpreted as a half-normal and why prior_exponential() is also valid.

When the guards introduce bias

The guards are one-sided: they replace invalid values with extreme but finite ones rather than with a symmetric distribution. In practice:

• For binomial outcomes with event rates in the range ⁠[1e-4, 1 - 1e-4]⁠ the clipped output is indistinguishable from the unclipped value at double precision.

• For event rates outside that range — or Poisson models with expected rates below 1e-4 events per unit exposure — the summary of the posterior for LOR, RR, and related ratios may be biased toward the clip boundary. The point estimate for the treatment effect on the link scale (delta_index on log-odds, log-rate, or mean-difference scale) is not affected since it does not go through the guards.

• Direction of safe_divide bias. When the comparator rate (the denominator) is smaller than 1e-10, the ratio is clipped to num / 1e-10. Because rates and probabilities in mlumr are non-negative, the clipped value is always larger (more positive) than the mathematical limit — so risk ratios and rate ratios in near-zero-comparator regimes are upward-biased, never downward. Expect the clipped posterior summary to overstate, not understate, the relative treatment effect in those draws.

• For very-rare-event data, consider reframing the model in terms of rates: an event count with exposure via set_agd() family = "poisson" is numerically better behaved than Bernoulli / Binomial at extreme probabilities.

How to diagnose guard activation

Posterior draws of ⁠p_*_index⁠ / ⁠p_*_comparator⁠ or ⁠rate_*_*⁠ are returned by the Stan models and are available in fit$draws. If any posterior draw sits below 5e-10 or above 1 - 5e-10 for binomial probabilities (respectively below 5e-10 for Poisson rates), the guards are active for that draw. A quick check:

  any(fit$draws[, "p_comparator_comparator"] < 5e-10)

Get or Set the Stan Engine

Description

Control which Stan backend mlumr uses for model fitting. The default engine is "rstan" (compiled C++ models via rstantools). Users who prefer cmdstanr can switch engines after installation.

Usage

mlumr_engine(engine = NULL)

Arguments

Details

When switching to "cmdstanr", this function checks whether the cmdstanr package and CmdStan toolchain are installed. If either is missing, it offers to install them interactively.

Engine names must be matched exactly. Partial strings such as "c" are not accepted.

The engine preference is stored as options(mlumr.stan_engine = ...) and persists for the current R session. To set a permanent default, add to your .Rprofile:

options(mlumr.stan_engine = "cmdstanr")

Value

The current engine (character), returned invisibly when setting.

Examples

# Check current engine
mlumr_engine()

# Switch to cmdstanr (interactive)
## Not run: 
mlumr_engine("cmdstanr")

## End(Not run)

Emit package progress messages when enabled

Description

Emit package progress messages when enabled

Usage

mlumr_message(..., verbose = TRUE)

Naive unadjusted indirect comparison

Description

Compute an unadjusted (naive) indirect treatment comparison by comparing crude outcomes from the IPD and AgD without any covariate adjustment. Returns treatment effect on the link scale plus absolute predictions (event probabilities, risk difference, risk ratio) for both populations.

Usage

naive(data, link = NULL, conf_level = 0.95)

Arguments

Details

For binomial outcomes, event-probability intervals use Wald standard errors and are bounded to ⁠[0, 1]⁠. For Poisson outcomes, the log-rate contrast uses a 0.5 continuity correction when an observed event count is zero.

Value

An object of class mlumr_naive

Examples

## Not run: 
result <- naive(dat)
print(result)

## End(Not run)

Bernoulli CDF

Description

Bernoulli CDF

Usage

pbern(q, prob, lower.tail = TRUE, log.p = FALSE)

Arguments

Value

Numeric vector

Predictions from ML-UMR model

Description

Generate population-average absolute-outcome predictions in the index and comparator populations.

Usage

## S3 method for class 'mlumr_fit'
predict(
  object,
  population = c("both", "index", "comparator"),
  type = c("response", "link"),
  summary = TRUE,
  probs = c(0.025, 0.5, 0.975),
  ...
)

Arguments

Details

Marginalization on non-identity links. For type = "response" the reported values are E[g^{-1}(eta)] — the posterior expectation of the inverse-link-transformed linear predictor — not g^{-1}(E[eta]). The two differ whenever g is non-linear (logit, probit, cloglog, log) by Jensen's inequality. In the index population the expectation is taken over IPD individuals; in the comparator population it is taken over the integration points constructed by add_integration() from the AgD moments. This is the correct population-average prediction for an individual randomly drawn from that population, and it matches what the Stan ⁠generated quantities⁠ block computes. For the link scale (type = "link") the reported value is E[eta], a linear functional, and the two interpretations coincide.

Value

A data frame with predictions. When type = "link", values are mean linear predictors computed directly from parameter draws (avoiding Jensen's inequality bias).

See Also

marginal_effects() for treatment-effect summaries; conditional_predict() and conditional_effects() for predictions at specific covariate profiles.

Specify a Cauchy prior

Description

Cauchy is Student-t with df = 1; this constructor is a convenience wrapper around prior_student_t(). It has very heavy tails and should be used with care — modern recommendations generally prefer ⁠prior_student_t(df in 3:7, ...)⁠ over Cauchy for regression coefficients to keep sampling well-behaved (see Piironen & Vehtari on the horseshoe; Ghosh et al. 2015).

Usage

prior_cauchy(mean = 0, sd = 2.5, autoscale = FALSE)

Arguments

Value

A list with components distribution = "student_t", df = 1, mean, sd, autoscale.

Examples

prior_cauchy(mean = 0, sd = 2.5)

Specify an exponential prior

Description

Exponential prior on a positive scalar. Currently supported for prior_sigma (normal-family residual SD) only; rejected for unconstrained intercepts and regression coefficients. prior_exponential(rate = 1) has mean 1 and is a reasonable weakly-informative choice when the outcome is standardized.

Usage

prior_exponential(rate = 1)

Arguments

Value

A list with components distribution = "exponential", rate, and placeholders (mean = 0, sd = 1 / rate) so the same Stan-field translation works.

Examples

prior_exponential(rate = 1)

Specify a normal prior

Description

Construct a normal prior for passing to mlumr() via prior_intercept, prior_beta, or prior_sigma. The resulting list is consumed by the Stan models.

Usage

prior_normal(mean = 0, sd = 10, autoscale = FALSE)

Arguments

Value

A list with components distribution, mean, sd, df, autoscale.

Choosing a scale

The Stan community's prior-choice wiki (Vehtari et al., 2025) describes five broad categories, from least to most informative:

1. Flat prior — not recommended.

2. Super-vague proper prior, e.g., normal(0, 1e6) — not recommended.

3. Weakly informative, very weak, e.g., normal(0, 10).

4. Generic weakly informative, e.g., normal(0, 1).

5. Specific informative, e.g., normal(0.4, 0.2).

Those scales assume parameters are on roughly unit scale. In ML-UMR models the natural scales are:

Treatment intercepts

On the linear-predictor (link) scale. For a binary outcome with logit link, the intercept is a baseline log odds; normal(0, 10) spans ±20 log-odds at 95 percent and is "very weak". It is the default because the data usually constrain the intercept strongly. Tightening to normal(0, 5) is reasonable when the expected event rate is far from the extremes.

Regression coefficients (beta)

On the link scale, per unit of covariate. For logistic regression with predictors on unit scale, Gelman et al. (2008) and the Stan wiki recommend student_t(df, 0, 2.5) with df in 3:7, or — as a practical approximation — normal(0, 2.5). That is the default used by mlumr(). Use normal(0, 1) if you expect small effects (e.g., standardized predictors in a normal-outcome model). If predictors are on very different scales, set autoscale = TRUE so the scale is divided by each covariate's SD.

Residual SD (sigma, normal family only)

prior_sigma is interpreted as a half-normal via the Stan ⁠<lower=0>⁠ constraint. The default normal(0, 2.5) (i.e., half-normal(0, 2.5)) is weakly informative for residual SDs on the scale of the outcome. Scale to the outcome if it is far from unit scale, or use prior_exponential().

Prior sensitivity is especially important for the relaxed model, where beta_comparator is identified only by the AgD likelihood. Run prior_sensitivity() to quantify how much conclusions move under alternative scales; see vignette("mlumr-models").

References

Gelman, A., Jakulin, A., Pittau, M. G., & Su, Y.-S. (2008). A weakly informative default prior distribution for logistic and other regression models. Annals of Applied Statistics, 2(4), 1360–1383.

Vehtari, A. et al. Prior Choice Recommendations (Stan wiki): https://github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations.

Examples

# Default weakly-very-weak intercept prior
prior_normal(mean = 0, sd = 10)

# Gelman 2008 default for logistic-regression coefficients
prior_normal(mean = 0, sd = 2.5)

# Autoscaled coefficient prior (dividing 2.5 by each covariate's SD)
prior_normal(mean = 0, sd = 2.5, autoscale = TRUE)

Prior sensitivity analysis for an ML-UMR fit

Description

Refit an mlumr() model across a grid of prior_beta scales (keeping the family, mean, and df fixed) and summarize how the posterior for the marginal treatment effects (delta_index, delta_comparator) moves. This is the workflow recommended by Vehtari et al.'s prior-choice wiki for judging how much of the posterior is driven by the data versus the prior.

Usage

prior_sensitivity(
  fit,
  prior_beta_scales = c(0.5, 1, 2.5, 5, 10),
  probs = c(0.025, 0.5, 0.975),
  verbose = TRUE,
  ...
)

Arguments

Details

Only the scale of the prior_beta family is varied; its distribution (normal / student_t) and mean are preserved so comparisons are apples to apples. prior_intercept and prior_sigma are carried through unchanged from the original fit. Each value in prior_beta_scales is used as the absolute scale for every coefficient at that refit — if the original fit used per-coefficient priors, all coefficients are set to the same scale (the sweep is deliberately homogeneous so the grid reflects a single level of prior informativeness per refit, not a rescaling of existing relative differences). If the original prior_beta used an exponential family, it is swapped for a prior_normal(0, scale) at each grid point since exponential has no scale parameter to vary.

Value

A data frame (tibble-style) with one row per (scale, population, quantile) combination and columns scale, parameter, mean, sd, and the requested quantiles. Side effect: prints a summary table at the end.

See Also

prior_summary() for a one-shot description of the priors on a fit; marginal_effects() for the posterior summary quantities this sweep tracks.

Examples

## Not run: 
sens <- prior_sensitivity(fit_spfa, prior_beta_scales = c(1, 2.5, 5))

## End(Not run)

Specify a Student-t prior

Description

Heavier-tailed alternative to prior_normal(). For logistic-regression coefficients with unit-scale predictors, the Stan community prior-choice recommendations suggest Student-t with df between 3 and 7 as a robust weakly informative prior (Gelman et al., 2008).

Usage

prior_student_t(df = 5, mean = 0, sd = 2.5, autoscale = FALSE)

Arguments

Value

A list with components distribution = "student_t", df, mean, sd, autoscale.

Examples

# Gelman et al. 2008 recommendation for logistic-regression coefficients
prior_student_t(df = 5, mean = 0, sd = 2.5)

Summary of priors used by a fitted ML-UMR model

Description

Print a human-readable summary of every prior that was used to fit an mlumr() model, including the effective per-coefficient scales after autoscaling. Mirrors the spirit of rstanarm::prior_summary().

Usage

prior_summary(object, ...)

## Default S3 method:
prior_summary(object, ...)

## S3 method for class 'mlumr_fit'
prior_summary(object, digits = 3, ...)

Arguments

Value

Invisibly returns a list describing the priors; the side effect is printing a formatted summary.

See Also

prior_sensitivity() to quantify how much the posterior moves under alternative prior_beta scales; prior_normal(), prior_student_t(), prior_cauchy(), prior_exponential() for the prior constructors themselves.

Examples

## Not run: 
fit <- mlumr(dat)
prior_summary(fit)

## End(Not run)

Bernoulli quantile function

Description

Bernoulli quantile function

Usage

qbern(p, prob, lower.tail = TRUE, log.p = FALSE)

Arguments

Value

Integer vector of 0s and 1s

Set up aggregate data (AgD)

Description

Prepare AgD from the comparator treatment for an unanchored indirect comparison.

Usage

set_agd(
  data,
  treatment,
  family = c("binomial", "normal", "poisson"),
  outcome_n = NULL,
  outcome_r = NULL,
  outcome_mean = NULL,
  outcome_se = NULL,
  outcome_E = NULL,
  cov_means,
  cov_sds = NULL,
  cov_types = NULL,
  study = NULL
)

Arguments

Details

Scale assumptions for 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 outcome_mean and outcome_se must be on the arithmetic (original, untransformed) scale. If a publication reports only a log-scale mean / SD or a geometric mean, back-transform before calling set_agd() and propagate uncertainty via the delta method; passing log-scale or geometric summaries silently misspecifies the likelihood and biases the posterior.

Scale assumptions for family = "poisson". outcome_r is the total count in each AgD row and outcome_E is the total person-time (or other exposure). The Stan likelihood uses log(E_agd) as an offset, so rates are modeled on the log scale regardless of how outcome_r is tabulated.

Scale assumptions for family = "binomial". outcome_r / outcome_n are counts of events and trials. The log-odds (or probit / cloglog under alternative links) are formed from outcome_r / outcome_n, so no scale conversion is required.

Value

An object of class mlumr_agd

Examples

## Not run: 
# Binary outcome
agd <- set_agd(
  data = trial_b,
  treatment = "trt",
  outcome_n = "n_total",
  outcome_r = "n_events",
  cov_means = c("age_mean", "sex_prop"),
  cov_sds = c("age_sd", NA),
  cov_types = c("continuous", "binary")
)

# Continuous outcome
agd <- set_agd(
  data = trial_b,
  treatment = "trt",
  family = "normal",
  outcome_mean = "mean_score",
  outcome_se = "se_score",
  outcome_n = "n_total",
  cov_means = c("age_mean", "sex_prop")
)

# Count outcome
agd <- set_agd(
  data = trial_b,
  treatment = "trt",
  family = "poisson",
  outcome_r = "n_events",
  outcome_E = "person_years",
  cov_means = c("age_mean", "sex_prop")
)

## End(Not run)

Set up individual patient data (IPD)

Description

Prepare IPD from the index treatment for an unanchored indirect comparison.

Usage

set_ipd(
  data,
  treatment,
  outcome,
  covariates,
  family = c("binomial", "normal", "poisson"),
  exposure = NULL,
  study = NULL
)

Arguments

Value

An object of class mlumr_ipd

Examples

## Not run: 
# Binary outcome
ipd <- set_ipd(
  data = trial_a,
  treatment = "trt",
  outcome = "response",
  covariates = c("age", "sex")
)

# Continuous outcome
ipd <- set_ipd(
  data = trial_a,
  treatment = "trt",
  outcome = "score",
  covariates = c("age", "sex"),
  family = "normal"
)

# Count outcome with exposure
ipd <- set_ipd(
  data = trial_a,
  treatment = "trt",
  outcome = "events",
  covariates = c("age", "sex"),
  family = "poisson",
  exposure = "person_years"
)

## End(Not run)

Translate a scalar prior spec to the Stan data fields

Description

Stan scalar-prior fields: ⁠prior_*_mean⁠, ⁠prior_*_sd⁠, ⁠prior_*_dist⁠ (0 = normal, 1 = student_t, 2 = exponential for sigma only), ⁠prior_*_df⁠ (used only when dist == 1; positive placeholder otherwise).

Usage

stan_prior_fields(prior)

Arguments

Value

A list with mean, sd, dist, df.

Translate a prior_beta spec to vector-valued Stan data fields

Description

The Stan models declare prior_beta_mean and prior_beta_sd as vector[n_cov] so the same data contract supports:

Usage

stan_prior_fields_beta(prior, n_cov, sd_x = NULL, covariate_names = NULL)

Arguments

Details

• (i) a single prior broadcast to all covariates (scalar user input),

• (ii) per-coefficient priors (a list of prior lists, length n_cov),

• (iii) autoscaling: each covariate's scale is divided by sd(x_j) so the prior is weakly-informative on the standardized scale.

For a list of per-coefficient priors, all elements must use the same distribution family and df (Stan branches on a single dist code).

Value

A list with numeric vectors mean and sd (length n_cov) and scalars dist, df, and a logical vector autoscale recording which elements were autoscaled (for prior_summary()).

Simulated treatment comparison via G-computation

Description

Perform an unanchored simulated treatment comparison (STC) using parametric G-computation. Fits a regression on IPD, predicts counterfactual outcomes in both the index and comparator populations, and computes marginal treatment effects with delta-method standard errors.

Usage

stc(data, link = NULL, conf_level = 0.95)

Arguments

Details

For binomial outcomes, returns the treatment effect on the link scale plus event probabilities, risk difference, and log risk ratio with SEs and CIs for both populations. Event-probability intervals use Wald standard errors and are bounded to ⁠[0, 1]⁠. For Poisson outcomes, the comparator log rate uses a 0.5 continuity correction when the observed event count is zero.

The STC procedure is:

1. Fit a GLM on IPD (binomial/gaussian/poisson as appropriate).

2. Predict on comparator-population covariates (from integration points or AgD covariate means).

3. Marginalize predictions over the comparator population.

4. Predict on index-population covariates (IPD).

5. Compute treatment effects and SEs via the delta method.

STC is a parametric G-computation benchmark. It relies on the IPD outcome model being correctly specified and transportable to the comparator population. It does not model posterior uncertainty in population covariate distributions or relax treatment-specific covariate effects. When clinically meaningful effect modification is plausible, prefer mlumr(..., model = "relaxed") as the primary analysis and use STC as a sensitivity or benchmarking analysis.

For binomial outcomes, the comparator-population treatment contrast is transported to the index population assuming the treatment contrast is constant on the fitted link scale (i.e., no effect modification on that scale). Under this assumption, the index-population comparator probability is computed as inv_link(link(p_A_index) - (link(p_A_comp) - link(p_B))), and its uncertainty is propagated through the delta method. If effect modification is expected, fit a Bayesian relaxed model with mlumr(..., model = "relaxed") and use predict(..., population = "index") instead, which does not require this assumption.

Value

An object of class mlumr_stc

Examples

## Not run: 
result <- stc(dat)
print(result)

## End(Not run)

Expand integration points into a long-format data frame

Description

Expand integration points into a long-format data frame

Usage

unnest_integration(data)

Arguments

Value

A data frame with columns for each covariate plus .int_id and .agd_row

Validate prior specification

Description

Validate prior specification

Usage

validate_prior(prior, param_name = "parameter")

Arguments