big differences between brms and rstanarm

[Bruno Nicenboim]

I'm seeing some strange differences in the population level effects between brms and rstanarm. I'm attaching the dataset.

I'm using quite weak priors, so I wouldn't expect them to change my posterior so much. I think that the posteriors of rstanarm make more sense, so I might be doing something wrong with brms.

This is the output of brms:

df_F05 <- read.table("Data/Fuchs_2005_complete.tsv", sep="\t", header=T)
df_F05 %<>% mutate(cvoicing=if_else(voicing=="voiced",.5,-.5))

head(df_F05)

#   publication subject experiment    dur voicing stimuli pair cvoicing
# 1  Fuchs_2005      cg        ema 189.87  voiced  gedaad   aa      0.5
# 2  Fuchs_2005      cg        ema 178.56  voiced  gedaad   aa      0.5
# 3  Fuchs_2005      cg        ema 184.19  voiced  gedaad   aa      0.5
# 4  Fuchs_2005      cg        ema 177.50  voiced  gedaad   aa      0.5
# 5  Fuchs_2005      cg        ema 200.31  voiced  gedaad   aa      0.5
# 6  Fuchs_2005      cg        ema 184.56  voiced  gedaad   aa      0.5


quantile(df_F05$dur, c(.275,.5,.975))
 #  27.5%     50%   97.5%
 # 91.000 113.895 193.000

# Linear mixed model:
fit_F05 <- brm(dur ~  cvoicing + (cvoicing|pair)  +(cvoicing|subject), data = df_F05,
    prior= c(prior(normal(0,300), class = Intercept),
             prior(normal(0,100), class = b, coef = cvoicing),
             prior(normal(0,100), class = sd),
             prior(normal(0,100), class = sigma),
             prior(lkj(2), class = cor)
             ),iter=iter,
                control=list(adapt_delta=.9999, max_treedepth=17))


fit_F05
#  Family: gaussian
#   Links: mu = identity; sigma = identity
# Formula: dur ~ cvoicing + (cvoicing | pair) + (cvoicing | subject)
#    Data: df_F05 (Number of observations: 270)
# Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#          total post-warmup samples = 4000
#     ICs: LOO = NA; WAIC = NA; R2 = NA
 
# Group-Level Effects:
# ~pair (Number of levels: 3)
#                         Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
# sd(Intercept)              57.74     35.21    18.68   153.36       1779 1.00
# sd(cvoicing)               10.60     17.17     0.22    59.60       1003 1.01
# cor(Intercept,cvoicing)     0.05      0.45    -0.78     0.84       4000 1.00

# ~subject (Number of levels: 3)
#                         Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
# sd(Intercept)              49.39     33.75    14.23   140.22       1808 1.00
# sd(cvoicing)               19.95     22.81     1.95    88.47       1119 1.00
# cor(Intercept,cvoicing)    -0.02      0.45    -0.82     0.78       3324 1.00

# Population-Level Effects:
#           Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
# Intercept   115.71     49.18    19.63   221.39       1178 1.01
# cvoicing     31.83     20.78   -16.44    69.75       1294 1.00

# Family Specific Parameters:
#       Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
# sigma    14.54      0.63    13.38    15.83       4000 1.00

# Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
# is a crude measure of effective sample size, and Rhat is the potential
# scale reduction factor on split chains (at convergence, Rhat = 1).


I think that the population level-effects are too wide. Given the quantiles of the variable dur (above), the posterior of the Intercept seems unreasonable. And I also think that the posterior of cvoicing is too wide.


If I try rstanarm these posteriors make more sense (to me). I know the models aren't exactly the same (because of the centered intercept), but I wouldn't expect such a big difference:

fit_F05s <- stan_glmer(dur ~  cvoicing  + (cvoicing|pair)+(cvoicing|subject), data = df_F05, prior = normal(0,100), prior_intercept = normal(0,300),  prior_aux = normal(0,100),prior_covariance=decov(2), adapt_delta=.9999, family = gaussian(),iter=iter)
summary(fit_F05s, probs=c(.0275,.975))
# Model Info:

#  function:     stan_glmer
#  family:       gaussian [identity]
#  formula:      dur ~ cvoicing + (cvoicing | pair) + (cvoicing | subject)
#  algorithm:    sampling
#  priors:       see help('prior_summary')
#  sample:       4000 (posterior sample size)
#  observations: 270
#  groups:       pair (3), subject (3)

# Estimates:
#                                          mean    sd      2.75%   97.5%
# (Intercept)                              118.1    22.7    73.4   164.7
# cvoicing                                  33.6    11.8    10.1    58.3
# b[(Intercept) pair:aa]                    33.7    17.8    -1.0    70.5
# b[cvoicing pair:aa]                        0.9     8.0   -14.5    17.8
# b[(Intercept) pair:ie]                   -25.2    17.8   -60.3    10.6
# b[cvoicing pair:ie]                       -2.9     7.9   -19.4    13.3
# b[(Intercept) pair:uh]                   -10.9    17.8   -45.8    25.9
# b[cvoicing pair:uh]                        2.0     8.0   -13.0    19.6
# b[(Intercept) subject:cg]                 18.7    14.5   -10.4    46.9
# b[cvoicing subject:cg]                     2.8     9.3   -15.2    21.4
# b[(Intercept) subject:jd]                  4.8    14.5   -23.9    33.0
# b[cvoicing subject:jd]                    -6.6     9.4   -26.1    10.6
# b[(Intercept) subject:sf]                -24.1    14.5   -53.4     4.4
# b[cvoicing subject:sf]                     3.4     9.3   -15.1    22.1
# sigma                                     14.5     0.6    13.4    15.8
# Sigma[pair:(Intercept),(Intercept)]      896.1   836.2   202.0  3073.1
# Sigma[pair:cvoicing,(Intercept)]          17.4   187.6  -318.2   403.7
# Sigma[pair:cvoicing,cvoicing]            173.8   334.3     1.7  1207.4
# Sigma[subject:(Intercept),(Intercept)]   607.9   643.2   119.6  2276.2
# Sigma[subject:cvoicing,(Intercept)]      -13.7   219.7  -407.5   382.1
# Sigma[subject:cvoicing,cvoicing]         253.0   475.3     9.6  1345.8
# mean_PPD                                 117.1     1.2   114.7   119.5
# log-posterior                          -1157.4     4.3 -1166.8 -1150.2

# Diagnostics:
#                                        mcse Rhat n_eff
# (Intercept)                             0.6  1.0 1401
# cvoicing                                0.3  1.0 1970
# b[(Intercept) pair:aa]                  0.4  1.0 1577
# b[cvoicing pair:aa]                     0.2  1.0 2338
# b[(Intercept) pair:ie]                  0.5  1.0 1559
# b[cvoicing pair:ie]                     0.2  1.0 2135
# b[(Intercept) pair:uh]                  0.5  1.0 1556
# b[cvoicing pair:uh]                     0.2  1.0 2217
# b[(Intercept) subject:cg]               0.3  1.0 1870
# b[cvoicing subject:cg]                  0.2  1.0 2389
# b[(Intercept) subject:jd]               0.3  1.0 1883
# b[cvoicing subject:jd]                  0.2  1.0 2248
# b[(Intercept) subject:sf]               0.3  1.0 1877
# b[cvoicing subject:sf]                  0.2  1.0 2313
# sigma                                   0.0  1.0 4000
# Sigma[pair:(Intercept),(Intercept)]    17.7  1.0 2243
# Sigma[pair:cvoicing,(Intercept)]        3.6  1.0 2646
# Sigma[pair:cvoicing,cvoicing]           8.0  1.0 1764
# Sigma[subject:(Intercept),(Intercept)] 12.7  1.0 2558
# Sigma[subject:cvoicing,(Intercept)]     4.2  1.0 2764
# Sigma[subject:cvoicing,cvoicing]       10.9  1.0 1899
# mean_PPD                                0.0  1.0 3915
# log-posterior                           0.1  1.0 1208

# For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

As a sanity check, I've also tried a similar model with lme4:

library(lme4)
summary(lmer(dur ~  cvoicing + (cvoicing||pair)  +(cvoicing||subject), data = df_F05))
# Linear mixed model fit by REML ['lmerMod']
# Formula: dur ~ cvoicing + ((1 | pair) + (0 + cvoicing | pair)) + ((1 |  
#     subject) + (0 + cvoicing | subject))
#    Data: df_F05

# REML criterion at convergence: 2228.7

# Scaled residuals:
#     Min      1Q  Median      3Q     Max
# -3.2698 -0.6221  0.0035  0.5158  3.0776

# Random effects:
#  Groups    Name        Variance Std.Dev.
#  pair      (Intercept) 950.056  30.823  
#  pair.1    cvoicing      1.806   1.344  
#  subject   (Intercept) 479.101  21.888  
#  subject.1 cvoicing     33.869   5.820  
#  Residual              209.340  14.469  
# Number of obs: 270, groups:  pair, 3; subject, 3

# Fixed effects:
#             Estimate Std. Error t value
# (Intercept)  117.054     21.844   5.359
# cvoicing      33.509      3.872   8.654

# Correlation of Fixed Effects:
#          (Intr)
# cvoicing 0.000

[Bruno Nicenboim]

By the way, this dataset is quite small, with similar models but more data, I  only see minor differences between rstanarm and brms.

[Paul Buerkner]

I wouldn't call those "big" differences ;-)

They are driven by rstanarm's default priors, which can end up being rather informative for small data sets.

See also https://groups.google.com/forum/#!topic/brms-users/gj5-NvAUlEU

--
You received this message because you are subscribed to the Google Groups "brms-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to brms-users+unsubscribe@googlegroups.com.
To post to this group, send email to brms-...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/brms-users/c77810c3-04cc-4373-b711-7431de8a86c0%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

[Paul Buerkner]

I wouldn't call those "big" differences ;-)

They are driven by rstanarm's default priors, which can end up being rather informative for small data sets.

See also https://groups.google.com/forum/#!topic/brms-users/gj5-NvAUlEU