Hierarchical random intercept model - Stan vs brms

Modeling
,
1

I am fitting a hierarchical random intercept model in Stan and brms but getting different results. Perhaps I’ve incorrectly specified the brms model, but I can’t see the problem.

Data: Y_{ij} \mid \mu, a_i, \sigma \stackrel{indep.}{\sim} N(\mu, \sigma), i = 1,2,\ldots,20; j = 1,2.

Priors:
\mu \sim N(0, 30)
\sigma \sim t_{4}(0, 1)
a_1, \ldots, a_{20} \stackrel{i.i.d}{\sim} N(0, \sigma_a)
\sigma_a \sim t_{4}(0, 1)

Here is reproducible code for the Stan version of the model:

library(rstan)
library(brms)
# data
y <- c(108, 98, 91, 94, 93, 96, 104, 99, 99, 97, 95, 98, 93, 97, 99,
96, 90, 100, 92, 95, 101, 89, 97, 97, 97, 100, 96, 95, 106, 100, 100,
98, 90, 99, 88, 98, 92, 92, 100, 101)
subject <- rep(1:20, each = 2)

stan_code <- "
data {
  vector[40] y;  // response values
}
parameters {
  real mu;      // common mean
  real<lower=0> sigma_a; // sd of between-subject random effects
  real<lower=0> sigma;  // sd of within-subject random effects
  vector[20] a;  // random effect for mean
}
model {
  // prior distributions
  mu ~ normal(0, 30);
  sigma_a ~ student_t(4, 0, 1);
  sigma ~ student_t(4,0, 1);
  // distribution of random effects for mean
  for(i in 1:20) a[i] ~ normal(0, sigma_a);
  // data distribution
  for(i in 1:20) {
    y[(i - 1) * 2 + 1] ~ normal(mu + a[i], sigma);
    y[(i - 1) * 2 + 2] ~ normal(mu + a[i], sigma);
  }
}
"
stan_mod <- stan(model_code = stan_code, data = list(y = y),
                 chains = 4, iter = 10000, refresh = 0, seed = 23)
stan_sum <- summary(stan_mod, prob = c(0.025, 0.975))
round(stan_sum$summary[1:3,c(1, 3:5)], 3)

The posterior results are

mean sd 2.5% 97.5%
\mu 96.693 0.723 95.254 98.105
\sigma_a 1.174 0.730 0.181 2.851
\sigma 4.111 0.530 3.176 5.276

Here is my attempt at fitting the same model using brms:

bdf <- data.frame(y = y, subject = as.factor(subject))
bmod <- brm(y ~ 1 + (1 | subject), data = bdat, family = gaussian(),
            prior = c(set_prior("normal(0, sqrt(1000))", class = "Intercept"),
                      set_prior("student_t(4, 0, 1)", class = "sd"),
                      set_prior("student_t(4, 0, 1)", class = "sd",
                                group = "subject")),
            chains = 4, iter = 10000, refresh = 0, seed = 23)
b_sum <- summary(bmod, prob = 0.95)
round(rbind(b_sum$fixed[1:4], b_sum$random$subject[1:4], b_sum$spec_pars[1:4]), 3)

This produces the following posterior results, which are similar, but not close enough to be the same.

Estimate Est.Error l-95% CI u-95% CI
Intercept 96.695 0.736 95.227 98.142
sd(Intercept) 0.928 0.703 0.035 2.586
sigma 4.344 0.541 3.395 5.512

I looked at the Stan code used by brms, and there is one line that seems to be causing a problem:

stancode(bmod)

The problematic line is

lprior += student_t_lpdf(sigma | 3, 0, 4.4) - 1 * student_t_lccdf(0 | 3, 0, 4.4);

It seems that brms is using t_3(0, 4.4) prior for \sigma instead of t_4(0, 3), as is described in the Stan model and in the brms code.

Any insight into how to produce the Stan results using the brms package is appreciated.