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)