here's my model in lme4

library(nlme)   # for the data
data("Machines")  # the data

fit2 <- lmer(score ~  -1 + Machine + (1|Worker), data=Machines)

Linear mixed model fit by REML ['lmerMod']
Formula: score ~ -1 + Machine + (1 | Worker)
   Data: Machines

REML criterion at convergence: 286.9

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-2.7249 -0.5233  0.1328  0.6513  1.7559 

Random effects:
 Groups   Name        Variance Std.Dev.
 Worker   (Intercept) 26.487   5.147   
 Residual              9.996   3.162   
Number of obs: 54, groups:  Worker, 6

Fixed effects:
         Estimate Std. Error t value
MachineA   52.356      2.229   23.48
MachineB   60.322      2.229   27.06
MachineC   66.272      2.229   29.73

Correlation of Fixed Effects:
         MachnA MachnB
MachineB 0.888        
MachineC 0.888  0.888 

when I model this in STAN using rstanarm, I get a totally different output:

fit2a <- stan_glmer(score ~  -1 + Machine + (1|Worker), data=Machines)

stan_glmer(formula = score ~ -1 + Machine + (1 | Worker), data = Machines)

         Median MAD_SD
MachineA -6.3   24.4  
MachineB  1.6   24.3  
MachineC  7.5   24.4  
sigma     3.6    0.4  

Error terms:
 Groups   Name        Std.Dev.
 Worker   (Intercept) 38.2    
 Residual              3.6    
Num. levels: Worker 6 

Sample avg. posterior predictive 
distribution of y (X = xbar):
         Median MAD_SD
mean_PPD 59.6    0.7 

my interpretation is this: I get the fixed effects by adding mean_PPD to the estimates, e.g. 59.6-6.3 = 53.3 ~ 52.3 as in lme.

Are the error terms in the stan output the random effects? And why does it use the median for the estimates or is this a rstanarm specific thing?

  • 2
    $\begingroup$ The print method in the rstanarm package reports posterior medians because they always exist (whereas the mean may not, although that is unlikely for the models currently in rstanarm) and perhaps are better estimated than is the mean in cases where the Markov chain has issues sampling from the tails of the posterior distribution. $\endgroup$ Jan 18 '16 at 3:29
  • 1
    $\begingroup$ The output under Error terms in rstanarm is comparable to the output under Random effects in lme4. But since rstanarm is largely Bayesian, the phrases "fixed effects" and "random effects" are avoided, although there are fixef and ranef extractor functions for compatibility with lme4. Either way, you can think of the (1|Worker) term as an error in the sense that the Machine predictor is unable to explain some variation in score. $\endgroup$ Jan 18 '16 at 3:35
  • 1
    $\begingroup$ To your main question, I think you have encountered a bug in rstanarm. See github.com/stan-dev/rstanarm/issues/49 . $\endgroup$ Jan 18 '16 at 3:36

The main question was the result of a bug in rstanarm that has since been fixed on GitHub.

However, in general, we do not recommend rstanarm models that exclude the intercept. A better alternative is to place a tight prior with mean zero on the intercept. In this case, including the intercept yields a better fit and similar results between lme4::lmer and rstanarm::stan_lmer.

  • $\begingroup$ "We allow priors, use them!" But more seriously, I hadn't thought of this before and it makes a lot of sense. There are probably several CV discussions about excluding intercepts and the issues involved, and priors provide a better (continuous, partially-data-driven) way to handle this kind of thing. $\endgroup$
    – Wayne
    Nov 22 '16 at 20:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.