Reproducing simple rstanarm meta-analysis in brms

Interfaces brms
1

Hi All,

I’m attempting to carry out a simple meta-analysis via brms. I’ve started with a very simple model, but the results I get depending on whether I use rstanarm or brms have been different. I wonder what’s the source of this difference, and how to set up the model in brms to reproduce the result from rstanarm. Here’s a simple example:

library(tidyverse)
library(rstanarm)
library(brms)


test_data <- tibble(yi = c(0.3, 0.8, -0.05, -0.2), dataset_name = c("A", "B", "C", "D"))

rstanarm_test <- stan_lmer(
  yi ~ (1 | dataset_name),
  data = test_data,
  adapt_delta = .999,
  cores = 4
)

## prior summary includes an exponential prior and covariance
rstanarm_test %>% summary
rstanarm_test %>% prior_summary

## Extracting samples...
rstanarm_samples <- rstanarm_test %>% as.matrix %>% as_tibble %>% select(1:5) %>% rename_all(function (x) x %>% str_remove(".*dataset_name:") %>% str_remove("\\]")) %>% rename(Summary = `(Intercept)`) %>% mutate_at(2:5, function (x) {.$Summary + x})

## And this is the summary of the sample:
rstanarm_samples %>% gather(key = "Collection") %>% group_by(Collection) %>% summarize(es_mean = mean(value), es_se = sd(value))


## Now trying to get a similar model with brms:
brms_test <- brm(
  yi ~ (1 | dataset_name),
  data = test_data,
  prior = c(
    set_prior("normal(0, 4)", class = "Intercept"), set_prior("exponential(0.44)", class = "sigma"), set_prior("normal(1, 1)", class = "sd") ## normal(0, 4.44) and exponential(0.44) since those was the "scaled" prior summary outcomes from rstanarm
  ),
  control = list(adapt_delta = 0.999),
  cores = 4
)
brms_test %>% summary
brms_test %>% prior_summary

## Extracting samples...
brms_samples <- brms_test %>% as.matrix %>% as_tibble %>% select(4:8) %>% rename_all(function (x) x %>% str_remove("r_dataset_name\\[") %>% str_remove(",Intercept\\]")) %>% rename(Summary = Intercept) %>% mutate_at(2:5, function (x) {.$Summary + x})

## Summary of the sample:
brms_samples %>% gather(key = "Collection") %>% group_by(Collection) %>% summarize(es_mean = mean(value), es_se = sd(value))

The samples drawn from the two models are significantly different in terms of the posterior estimate distributions of effect sizes. Is there a way to reproduce the rstanarm result in brms?

  • Operating System: Ubuntu 16.04
  • brms Version: 2.11.1
  • rstanarm Version: 2.18.2

Thanks!

2

Default priors are different in brms and rstanarm. Thus, for very small data sets, these default priors play a relevant role in the posterior and hence results may be different.