Why do Stan model and rstanarm/brms model with same priors have different posterior?

Modeling
1

The rstanarm/brms model have same posterior, but are different with Stan code.

Thank you very much for your help and suggestions

library(tidyverse)
library(cmdstanr)
check_cmdstan_toolchain(fix = TRUE, quiet = TRUE)

penguins <- palmerpenguins::penguins %>% tidyr::drop_na()

1. stan code

stan_program <- write_stan_file("
data {
  int N;
  vector[N] y;
  vector[N] x;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  y ~ normal(alpha + beta * x, sigma);

  alpha  ~ normal(0, 5);
  beta   ~ normal(0, 5);
  sigma  ~ exponential(1);
}

")


stan_data <- list(
  N = nrow(penguins),
  x = penguins$flipper_length_mm,
  y = penguins$body_mass_g
)

mod <- cmdstan_model(stan_file = stan_program)
fit <- mod$sample(data = stan_data)


fit$summary(variables = c("alpha", "beta"), "mean", "median", "sd")
# A tibble: 2 × 4
  variable  mean median    sd
  <chr>    <dbl>  <dbl> <dbl>
1 alpha    -1.52  -1.56 4.99 
2 beta     21.1   21.1  0.108

2. brms code

library(brms)

mod_brms <- brm(
  data = penguins,
  body_mass_g ~ 1 + flipper_length_mm,
  family = gaussian(link = "identity"),
  prior = c(
    prior(normal(0, 5), class = Intercept),
    prior(normal(0, 5), class = b),
    prior(exponential(1), class = sigma)
  ),
  backend = "cmdstanr"
)

mod_brms %>% broom.mixed::tidy(effects = "fixed")
# A tibble: 2 × 7
  effect component term              estimate std.error conf.low conf.high
  <chr>  <chr>     <chr>                <dbl>     <dbl>    <dbl>     <dbl>
1 fixed  cond      (Intercept)        -3588.     817.   -5159.     -1967. 
2 fixed  cond      flipper_length_mm     17.9      4.07     9.86      25.7

3. cmdstanr

library(rstanarm)

mod_rstanarm <- stan_glm(
  data = penguins,
  body_mass_g ~ 1 + flipper_length_mm,
  family = gaussian(link = "identity"),
  prior_intercept = normal(0, 5),
  prior = normal(0, 5),
  prior_aux = exponential(1)
)

mod_rstanarm %>% broom.mixed::tidy()
# A tibble: 2 × 3
  term              estimate std.error
  <chr>                <dbl>     <dbl>
1 (Intercept)        -3583.     796.  
2 flipper_length_mm     17.9      3.96

The full code is available here
three_models.Rmd (1.8 KB)

1 Like
2

This is due to the fact that brms places a prior not on the intercept itself (i.e. the parameter alpha in the stan program), but rather on the intercept after centering the columns of the design matrix (i.e. on the value when all predictors are held at their means). To avoid this behavior, you can use the following syntax if desired:

library(brms)

mod_brms <- brm(
  data = penguins,
  body_mass_g ~ 0 + Intercept + flipper_length_mm,
  family = gaussian(link = "identity"),
  prior = c(
    prior(normal(0, 5), class = b),
    prior(exponential(1), class = sigma)
  ),
  backend = "cmdstanr"
)
2 Likes
3

I see. Thank you for your answer