Hi, I want to fit a Binomial-Normal (meta-analysis) model in brms. As further described below, this model has been applied in {lme4} and {MetaStan} R packages. Although all three packages are supposed to fit the same model, MetaStan doesn’t generate an Intercept and I would like to replicate this parameterization in brms.
Background
There is an article that describes 6 different Binomial-Normal models. More specifically, I want to fit their “Model 4”, which is an extension of their “Model 2”. The latter is described as:
while “Model 4” is described as:
lme4 model
Fortunately, to facilitate our understanding, the authors share a code on how to fit “Model 4” with the R package {lme4}:
lme4::glmer(
cbind(responders,sampleSize-responders)~
1 + factor(study) + factor(treat) + (treat12-1|study),
data=data,
family=binomial(link="logit")
)
where data would be a data frame with this structure:
MetaStan model
Further, I found another article fitting the same “Model 4” above (as the authors explicitly say in the article), but with Stan through the MetaStan package. In this 2019 article, they share the following Stan code:
The authors implemented the non-centered parameterization, as described below:
" In the presence of sparse data such as in the meta-analysis of few studies involving rare events, Betancourt et al44 showed that centered parametrization of a hierarchical model (such as the BNHM) brings computational issues compared with a noncentered parametrization. Thus, we use the noncentered reparametrized version of the BNHM for our implementations. Specifically, applying both location and scale reparametrization, 1 becomes μi + θ xij + ui τxij where ui ~ N(0, 1) and xij=+0.5 (experimental) or xij=−0.5 (control)."
After the publication of this article in 2019, the authors recently put their MetaStan package in CRAN and published an article last week specifically on this package and its multiple new models. Now, their Stan model is a little bit different (since they also support different models other than “Model 4” above). Nevertheless, it should still fit “Model 4” as the binomial option in the code below, as explicitly referenced in their package vignette:
/*
Author: Burak Kuersad Guenhan
A standard Meta-analysis model
Dataset need is one-row-per-arm format
*/
data {
// num of observations
int<lower=1> Nobs;
// treatment arm indicator (0=control, 1=experimental)
vector[Nobs] t;
// Study number for each observation
int<lower=1> st[Nobs];
// link function (1=normal, 2=binary, 3=poisson)
int<lower=1,upper=3> link;
// normal data, link=identity=1
vector[Nobs] y;
vector[Nobs] y_se;
// binomial data, link=logit=2
int<lower=0> r[Nobs];
int<lower=1> n[Nobs];
// count data, link=log=3
int<lower=0> count[Nobs];
vector[Nobs] exposure;
// Priors
vector[2] mu_prior; // Prior parameters for mu
vector[2] theta_prior; // Prior parameters for theta
real tau_prior; // Prior parameter for tau
int tau_prior_dist; // Indicator for distribution of prior for tau
vector[2] beta_prior; // Prior parameters for beta
// Fixed-effects or Random-effects model: (0: fe, 1: re)
int<lower=0,upper=1> re;
// non-centered parametrization: (0: NO, 1: YES)
int<lower=0,upper=1> ncp;
// met-regression: (0: NO, 1: YES)
int<lower=0,upper=1> mreg;
// number of covariates
int<lower=0> ncov;
// trial-level covariate information
matrix[ncov,max(st)] cov;
}
parameters {
vector[Nobs] mu; // baseline risks (log odds)
real theta; // relative treatment effect (log odds ratio)
vector[Nobs] u[re]; // individual treatment effects
real<lower=0> tau[re]; // heterogeneity stdev.
vector[ncov] beta[mreg]; // beta coeffients in meta-regression
}
transformed parameters {
vector[Nobs] d;
vector[max(st)] temp;
vector[max(st)] beta_cov;
if(mreg) {
for(i in 1:Nobs)
temp[st[i]] = 0.5 * beta[1,1] * cov[1,st[i]];
if(ncov == 1) {
for(i in 1:Nobs)
beta_cov[st[i]] = temp[st[i]];
} else {
for(j in 1:ncov) {
for(i in 1:Nobs) {
beta_cov[st[i]] = temp[st[i]] + 0.5 * beta[1,j] * cov[j,st[i]];
temp[st[i]] = beta_cov[st[i]];
}
}
}
} else {
for(i in 1:Nobs) {
beta_cov[st[i]] = 0;
temp[st[i]] = 0;
}
};
if(re) {
if (ncp) {
for(i in 1:Nobs) {
if (t[i] == 0) {
d[i] = mu[st[i]] - 0.5 * theta - 0.5 * u[1,st[i]] * tau[1] - beta_cov[st[i]];
} else {
d[i] = mu[st[i]] + 0.5 * theta + 0.5 * u[1,st[i]] * tau[1] + beta_cov[st[i]];
}
}
} else {
for(i in 1:Nobs) {
if (t[i] == 0) {
d[i] = mu[st[i]] - 0.5 * u[1,st[i]] - beta_cov[st[i]];
} else {
d[i] = mu[st[i]] + 0.5 * u[1,st[i]] + beta_cov[st[i]];
}
}
}
} else {
for(i in 1:Nobs) {
if (t[i] == 0) {
d[i] = mu[st[i]] - 0.5 * theta - beta_cov[st[i]];
} else {
d[i] = mu[st[i]] + 0.5 * theta + beta_cov[st[i]];
}
}
}
}
model {
if(re) {
if (ncp) {
// latent variable (random effects)
u[1] ~ normal(0, 1);
} else {
// random effect distribution
u[1] ~ normal(d, tau[1]);
}
}
// prior distributions
mu ~ normal(mu_prior[1], mu_prior[2]);
theta ~ normal(theta_prior[1], theta_prior[2]);
if(re) {
// prior for tau
if(tau_prior_dist == 1) tau[1] ~ normal(0, tau_prior)T[0,];
if(tau_prior_dist == 2) tau[1] ~ uniform(0, tau_prior);
if(tau_prior_dist == 3) tau[1] ~ cauchy(0, tau_prior)T[0,];
}
if(mreg) beta[1] ~ normal(beta_prior[1], beta_prior[2]);
// likelihood
if(link == 1) y ~ normal(d, y_se);
if(link == 2) r ~ binomial_logit(n, d);
if(link == 3) count ~ poisson_log(exposure + d);
}
generated quantities {
vector[Nobs] log_lik; // pointwise log-likelihood contribution
real theta_pred[re]; // predicted log-odds ratio for the new study
for (s in 1:Nobs) {
if(link == 1) log_lik[s] = normal_lpdf(y[s] | d[s], y_se[s]);
if(link == 2) log_lik[s] = binomial_logit_lpmf(r[s] | n[s], d[s]);
if(link == 3) log_lik[s] = poisson_log_lpmf(count[s] | exposure[s] + d[s]);
}
// Prediction for new study
if(re) { theta_pred[1] = normal_rng(theta, tau[1]);}
}
brms model
As shown in the code below, I am able to replicate the model fitted with lme4 in brms. However, the model is different from the MetaStan model and I don’t know why.
Here is the code to load packages, data, and fit lme4/brms models:
# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")
### Install/Load packages
pacman::p_load(dplyr,
lme4,
brms,
MetaStan)
pacman::p_load_gh("stan/cmdstanr")
### Data ----
data("dat.Crins2014", package = "MetaStan")
# Remove NAs
data = subset(dat.Crins2014, is.finite(exp.PTLD.events))
# Adapt dataset to MetaStan's format
dat_MetaStan = MetaStan::create_MetaStan_dat(
dat = data,
armVars = c(responders = "r",
sampleSize = "n"))
dat_MetaStan
# Adapt dataset to lme4 and brms's format
data_jackson =
dat_MetaStan$data_long |>
dplyr::mutate(treat = rep(c(0, 1), 3),
treat12 = rep(c(-0.5, 0.5), 3))
data_jackson
### Fit models
## Model 4 in Jackson et al 2018 ----
# https://onlinelibrary.wiley.com/doi/10.1002/sim.7588
jackson_model4 = lme4::glmer(
cbind(responders,sampleSize-responders)~
1 + factor(study) + factor(treat) + (treat12-1|study),
data=data_jackson,
family=binomial(link="logit")
)
fixef(jackson_model4)
# factor(treat)1
# -1.7149217
## brms model ----
# Formula
mf =
responders | trials(sampleSize) ~
1 + factor(study) + factor(treat) + (treat12 - 1 | study)
brms::get_prior(mf, family = binomial, data = data_jackson)
brms_model = brms::brm(
formula = mf,
family = binomial,
prior =
prior(normal(0, 100), class = Intercept) +
prior(normal(0, 100), class = b, coef = "factortreat1") +
prior(normal(0, 0.5), class = sd),
data = data_jackson,
backend = "cmdstanr"
)
brms_model
fixef(brms_model)
# Estimate Q2.5 Q97.5
# factortreat1 -1.7715522 -2.6682835 -0.903934675
As expected (given the similarities between their model syntax), both models yield the same parameters, including an Intercept:
Discrepancies with MetaStan
Now, here is the code to fit MetaStan’s model:
## MetaStan ----
# 3 sources about the same model
# 01: "Note that this model corresponds to Model 4 in Jackson et al (2018)."
# https://cran.r-project.org/web/packages/MetaStan/vignettes/MetaStan_BNHM.html
# 02: "2.2. Pairwise meta-analysis"
# https://arxiv.org/pdf/2202.00502.pdf
# 03: Corresponding Stan code
# Non-centered parametrization (ncp = TRUE);
# Random-effects (re = TRUE),
# Binomial
# No meta-regression
# https://github.com/gunhanb/MetaStan/blob/master/inst/stan/SMA.stan
MetaStan_model = MetaStan::meta_stan(data = dat_MetaStan,
likelihood = "binomial",
mu_prior = c(0, 10),
theta_prior = c(0, 100),
tau_prior = 0.5,
tau_prior_dist = "half-normal")
MetaStan_model
# Treatment effect (theta) estimates
# mean 2.5% 50% 97.5%
# -1.78 -2.66 -1.76 -0.97
Although fitting the corresponding model with MetaStan yields very similar results (“Treatment effect (theta) estimates”), there is no intercept:
Summary
Can anyone give me some help on how to fit MetaStan’s model in brms?
Thanks!
Here is the full code which I had broken down in separate pieces above:
# Ensures the package "pacman" is installed
if (!require("pacman")) install.packages("pacman")
### Install/Load packages
pacman::p_load(dplyr,
lme4,
brms,
MetaStan)
pacman::p_load_gh("stan/cmdstanr")
### Data ----
data("dat.Crins2014", package = "MetaStan")
# Remove NAs
data = subset(dat.Crins2014, is.finite(exp.PTLD.events))
# Adapt dataset to MetaStan's format
dat_MetaStan = MetaStan::create_MetaStan_dat(
dat = data,
armVars = c(responders = "r",
sampleSize = "n"))
dat_MetaStan
# Adapt dataset to lme4 and brms's format
data_jackson =
dat_MetaStan$data_long |>
dplyr::mutate(treat = rep(c(0, 1), 3),
treat12 = rep(c(-0.5, 0.5), 3))
data_jackson
### Fit models
## Model 4 in Jackson et al 2018 ----
# https://onlinelibrary.wiley.com/doi/10.1002/sim.7588
jackson_model4 = lme4::glmer(
cbind(responders,sampleSize-responders)~
1 + factor(study) + factor(treat) + (treat12-1|study),
data=data_jackson,
family=binomial(link="logit")
)
fixef(jackson_model4)
# factor(treat)1
# -1.7149217
## brms model ----
# Formula
mf =
responders | trials(sampleSize) ~
1 + factor(study) + factor(treat) + (treat12 - 1 | study)
brms::get_prior(mf, family = binomial, data = data_jackson)
brms_model = brms::brm(
formula = mf,
family = binomial,
prior =
prior(normal(0, 100), class = Intercept) +
prior(normal(0, 100), class = b, coef = "factortreat1") +
prior(normal(0, 0.5), class = sd),
data = data_jackson,
backend = "cmdstanr"
)
brms_model
fixef(brms_model)
# Estimate Q2.5 Q97.5
# factortreat1 -1.7715522 -2.6682835 -0.903934675
## MetaStan ----
# 3 sources about the same model
# 01: "Note that this model corresponds to Model 4 in Jackson et al (2018)."
# https://cran.r-project.org/web/packages/MetaStan/vignettes/MetaStan_BNHM.html
# 02: "2.2. Pairwise meta-analysis"
# https://arxiv.org/pdf/2202.00502.pdf
# 03: Corresponding Stan code
# Non-centered parametrization (ncp = TRUE);
# Random-effects (re = TRUE),
# Binomial
# No meta-regression
# https://github.com/gunhanb/MetaStan/blob/master/inst/stan/SMA.stan
MetaStan_model = MetaStan::meta_stan(data = dat_MetaStan,
likelihood = "binomial",
mu_prior = c(0, 10),
theta_prior = c(0, 100),
tau_prior = 0.5,
tau_prior_dist = "half-normal")
MetaStan_model
# Treatment effect (theta) estimates
# mean 2.5% 50% 97.5%
# -1.78 -2.66 -1.76 -0.97
MetaStan_model$fit
## MetaStan vs. brms
# 01:
# MetaStan does not have an "Intercept". In fact, it models
# 1 baseline odds per study k ("mu[k]")
MetaStan_model$fit
# Here we have an intercept (which is the baseline odds of study 1 )
# Then each "factostudy..." corresponds the deviation from this intercept
brms_model
