Skip to contents

Custom models

Source: vignettes/custom-models.qmd

Introduction

bayesma covers a broad range of meta-analytic models, but applied synthesis sometimes requires departures from the built-in options. The package supports three levels of customisation with increasing flexibility and complexity.

Level 1: Modify priors or add constraints

The simplest customisation changes only the prior distributions, without touching the model structure.

Example: Use an informative prior on τ\tau from a published prior predictive analysis.

fit_informative <- bayesma(
  data,
  model_type = "random_effect",
  tau_prior  = half_normal(0, 0.25),
  mu_prior   = normal(0, 0.5)
)

This requires only knowledge of the *_prior arguments. See ?priors for all available constructors:

normal(mu, sigma)
half_normal(mu, sigma)
half_cauchy(mu, sigma)
half_student_t(nu, mu, sigma)
exponential(rate)
uniform(lower, upper)

Constraints on the pooled effect can be added via mu_constraint:

fit_positive <- bayesma(
  data,
  model_type    = "random_effect",
  mu_constraint = "positive"
)

This restricts μ>0\mu > 0 — appropriate when the direction of the effect is known a priori (e.g., a toxin that cannot be beneficial).

Level 2: Extend an existing model

The second level modifies the Stan data or adds blocks to the generated model without replacing the core likelihood. The workflow is:

  1. Extract the autogenerated Stan code with bayesma(data, return_stage = "code").
  2. Edit the Stan code.
  3. Pass the modified code back to bayesma() via custom_model.

Example: Add a study-level covariate to the random-effects model.

base_code <- bayesma(data, model_type = "random_effect", return_stage = "code")

modified_code <- base_code |>
  stringr::str_replace(
    "vector\\[N\\] y;",
    "vector[N] y;\n  vector[N] x_cov;"
  ) |>
  stringr::str_replace(
    "target \\+= normal_lpdf\\(y \\| mu \\+ u\\[study\\], se\\);",
    "target += normal_lpdf(y | mu + beta * x_cov + u[study], se);"
  )

fit_extended <- bayesma(
  data,
  custom_model = modified_code,
  custom_data  = list(x_cov = data$covariate)
)

The custom_data argument passes additional data blocks to Stan. The Stan code must declare matching entries in the data {} block.

Level 3: Fully custom Stan model

The most flexible option passes a complete Stan program written from scratch. bayesma provides the data preparation pipeline and posterior extraction infrastructure; the user supplies the model.

Required data block structure

The data block must include at minimum:

data {
  int<lower=1> N;          // number of observations
  vector[N] y;             // effect estimates
  vector<lower=0>[N] se;   // standard errors
  array[N] int<lower=1> study;  // study indicator (for RE models)
}

For one-stage binary models, replace y and se with:

  array[N] int<lower=0> events;
  array[N] int<lower=1> n;
  vector[N] treat;

Required generated quantities block

For compatibility with bayesma_extract() and bayesma_output(), the generated quantities {} block must export b_Intercept (the pooled effect):

generated quantities {
  real b_Intercept = mu;
}

Example: Three-level model

A three-level (study > sample > effect) model for dependent effects:

data {
  int<lower=1> N;
  int<lower=1> K;        // number of studies
  int<lower=1> J;        // number of samples
  vector[N] y;
  vector<lower=0>[N] se;
  array[N] int<lower=1> study;
  array[N] int<lower=1> sample;
}

parameters {
  real mu;
  real<lower=0> tau_study;
  real<lower=0> tau_sample;
  vector[K] z_study;
  vector[J] z_sample;
}

transformed parameters {
  vector[K] u_study  = tau_study  * z_study;
  vector[J] u_sample = tau_sample * z_sample;
}

model {
  target += normal_lpdf(mu | 0, 1);
  target += cauchy_lpdf(tau_study  | 0, 0.5);
  target += cauchy_lpdf(tau_sample | 0, 0.5);
  target += std_normal_lpdf(z_study);
  target += std_normal_lpdf(z_sample);

  target += normal_lpdf(y | mu + u_study[study] + u_sample[sample], se);
}

generated quantities {
  real b_Intercept = mu;
}

Pass this to bayesma():

fit_three_level <- bayesma(
  data,
  custom_model = three_level_stan,
  custom_data  = list(J = n_samples, sample = data$sample_id)
)

Validation and diagnostics

Custom models should be checked for:

  • Identifiability: pairs() plots of correlated parameters reveal non-identifiable structure.
  • Prior coverage: Prior predictive simulation should produce plausible ranges for all parameters.
  • Convergence: R̂<1.01\hat{R} < 1.01 and bulk ESS > 400 for all parameters.
  • Divergent transitions: Zero divergences at adapt_delta = 0.95. Use the non-centred parameterisation for random effects (as in the examples above) to reduce divergences.

Limitations

Custom models bypass all built-in argument validation. The user is responsible for ensuring the model is correctly specified, identified, and that the custom_data list contains all variables referenced in the Stan data {} block.