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 


  [1]: https://i.sstatic.net/l4WPG.png