I've been playing around with both nlme::lme and lme4::lmer.
I fitted a simple random intercepts model using lme() and lmer(). As you can see below, I got completely different results from lmer() and lme(). Even signs of coefficients are different! Am I doing something wrong? I also fitted an empty model with the two packages. In this case, the results were practically the same (results not shown). Would you educate me to understand this issue? Unless I made a mistake, I think there is something wrong with the lme4 package.
multi<multi <- structure(list(x = c(4.9, 4.84, 4.91, 5, 4.95, 3.94, 3.88, 3.95,
4.04, 3.99, 2.97, 2.92, 2.99, 3.08, 3.03, 2.01, 1.96, 2.03, 2.12,
2.07, 1.05, 1, 1.07, 1.16, 1.11), y = c(3.2, 3.21, 3.256, 3.25,
3.256, 3.386, 3.396, 3.442, 3.436, 3.442, 3.572, 3.582, 3.628,
3.622, 3.628, 3.758, 3.768, 3.814, 3.808, 3.814, 3.944, 3.954,
4, 3.994, 4), pid = 1:25, gid = c(1L, 1L, 1L, 1L, 1L, 2L, 2L,
2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L,
5L, 5L)), class = "data.frame", row.names = c(NA, -25L))
#lme
> lme(y~x, random=~1|gid,data=multi,method="REML")
Linear mixed-effects model fit by REML
Data: multi
Log-restricted-likelihood: 41.76745
Fixed: y ~ x
(Intercept) x
4.1846756 -0.1928357
#lmer
lmer(y~x+(1|(gid)), data=multi, REML=T)
Linear mixed model fit by REML ['lmerMod']
Formula: y ~ x + (1 | (gid))
Data: multi
REML criterion at convergence: -78.4862
Random effects:
Groups Name Std.Dev.
(gid) (Intercept) 0.70325
Residual 0.02031
Number of obs: 25, groups: (gid), 5
Fixed Effects:
(Intercept) x
2.8152 0.2638
