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 two pakcages. In this case, the results were practically same (results not shwon). Would you educate me to understand this issue? Unless I made a mistake, I think there is something wrong with the lme4 package. Thank you.

     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