0

I'm trying the following model with the lme4 package:

library(nmle) # for the data
data("Machines") # the data
library(lme4)

# the model:
fit1 <- lmer(score ~ -1 + Machine + (1|Worker), data=Machines)
summary(fit1)

> summary(fit1)
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

I now try to fit the same model using rstan through the glmer2stan package:

library(glmer2stan)
Machines$Machine_idx <- as.numeric(Machines$Machine)
Machines$Worker_idx <- as.numeric(as.character(Machines$Worker))
fit3 <- lmer2stan(score ~ -1 + Machine_idx + (1|Worker_idx),  data=Machines)

this is the result

> stanmer(fit3)
glmer2stan model: score ~ -1 + Machine_idx + (1 | Worker_idx) [gaussian]

Level 1 estimates:
        Expectation StdDev 2.5% 97.5%
Machine_idx        7.04   0.55 5.95  8.08
sigma              3.26   0.35 2.66  4.02

Level 2 estimates:
(Std.dev. and correlations)

Group: Worker_idx (6 groups / imbalance: 0)
  (Intercept)  55.09  (SE 15.82)

DIC: 287   pDIC: 7.9   Deviance: 271.3

I don't think that's the same model. Is my glmer2stan specification wrong? I know that glmer2stan is not actively developed any more but it should handle this simple model, shouldn't it?

UPDATE:

thanks to the tip by Roland I changed the Machine factor levels to dummies and it now works:

Machines$Worker <- as.numeric(as.character(Machines$Worker))
m <- model.matrix(~ 0 + ., Machines)
m <- as.data.frame(m)

fit3 <- lmer2stan(score ~ -1 + (1|Worker) + MachineA + MachineB + MachineC, data=m, chains=2)
4

0 回答 0