lme() 和 lmer() 給出相互矛盾的結果

我一直在處理一些重複測量存在問題的數據。在這樣做的過程中，我注意到我的測試數據之間lme()和lmer()使用我的測試數據之間的行為非常不同，並且想知道為什麼。

我創建的假數據集包含 10 名受試者的身高和體重測量值，每人測量兩次。我設置了數據，以便受試者之間的身高和體重之間存在正相關關係，但每個個體內部的重複測量之間存在負相關關係。

set.seed(21)
Height=1:10; Height=Height+runif(10,min=0,max=3) #First height measurement
Weight=1:10; Weight=Weight+runif(10,min=0,max=3) #First weight measurement

Height2=Height+runif(10,min=0,max=1) #second height measurement
Weight2=Weight-runif(10,min=0,max=1) #second weight measurement

Height=c(Height,Height2) #combine height and wight measurements
Weight=c(Weight,Weight2)

DF=data.frame(Height,Weight) #generate data frame
DF$ID=as.factor(rep(1:10,2)) #add subject ID
DF$Number=as.factor(c(rep(1,10),rep(2,10))) #differentiate between first and second measurement

這是數據圖，用線連接每個人的兩次測量。

所以我運行了兩個模型，一個帶有lme()fromnlme包，一個帶有lmer()from lme4。在這兩種情況下，我使用 ID 的隨機效應對體重與身高進行回歸，以控制每個人的重複測量。

library(nlme)
Mlme=lme(Height~Weight,random=~1|ID,data=DF)
library(lme4)
Mlmer=lmer(Height~Weight+(1|ID),data=DF)

這兩個模型經常（儘管並不總是取決於種子）產生完全不同的結果。我已經看到它們在哪裡產生略有不同的方差估計，計算不同的自由度等，但這裡的係數是相反的方向。

coef(Mlme)
# (Intercept) Weight
#1 1.57102183 0.7477639
#2 -0.08765784 0.7477639
#3 3.33128509 0.7477639
#4 1.09639883 0.7477639
#5 4.08969282 0.7477639
#6 4.48649982 0.7477639
#7 1.37824171 0.7477639
#8 2.54690995 0.7477639
#9 4.43051687 0.7477639
#10 4.04812243 0.7477639

coef(Mlmer)
# (Intercept) Weight
#1 4.689264 -0.516824
#2 5.427231 -0.516824
#3 6.943274 -0.516824
#4 7.832617 -0.516824
#5 10.656164 -0.516824
#6 12.256954 -0.516824
#7 11.963619 -0.516824
#8 13.304242 -0.516824
#9 17.637284 -0.516824
#10 18.883624 -0.516824

為了直觀地說明，使用lme()

和模型lmer()

為什麼這些模型差異如此之大？

tl；博士如果您將優化器更改為“nloptwrap”，我認為它將避免這些問題（可能）。

恭喜，您在統計估計問題中找到了多重最優的最簡單示例之一！內部使用的參數lme4（因此便於說明）是隨機效應的標度標準差，即組間標準差除以殘差標準差。

提取這些原始值lme和lmer擬合值：

(sd1 <- sqrt(getVarCov(Mlme)[[1]])/sigma(Mlme))
## 2.332469
(sd2 <- getME(Mlmer,"theta")) ## 14.48926

用另一個優化器改裝（這可能是下一個版本的默認值lme4）：

Mlmer2 <- update(Mlmer,
 control=lmerControl(optimizer="nloptwrap"))
sd3 <- getME(Mlmer2,"theta")   ## 2.33247

比賽lme……讓我們看看發生了什麼。對於具有單個隨機效應的 LMM，偏差函數（-2log 似然性），或者在這種情況下是類似的 REML 標準函數，只需要一個參數，因為固定效應參數已被分析出來*；對於給定的 RE 標準偏差值，它們可以自動計算。

ff <- as.function(Mlmer)
tvec <- seq(0,20,length=101)
Lvec <- sapply(tvec,ff)
png("CV38425.png")
par(bty="l",las=1)
plot(tvec,Lvec,type="l",
    ylab="REML criterion",
    xlab="scaled random effects standard deviation")
abline(v=1,lty=2)
points(sd1,ff(sd1),pch=16,col=1)
points(sd2,ff(sd2),pch=16,col=2)
points(sd3,ff(sd3),pch=1,col=4)
dev.off()

我繼續對此更加著迷，並為每個案例運行了從 1 到 1000 的隨機種子的擬合，擬合lme,lmer和lmer+nloptwrap。以下是 1000 個數字，其中給定方法得到的答案至少比另一種方法差 0.001 個偏差單位……

         lme.dev lmer.dev lmer2.dev
lme.dev         0       64        61
lmer.dev      369        0       326
lmer2.dev      43        3         0

換句話說，（1）沒有一種方法總是最有效的；(2)lmer使用默認優化器最差（大約 1/3 的時間失敗）；(3)lmer最好使用“nloptwrap”（比lme4% 的時間差，很少比 差lmer）。

稍微讓人放心的是，我認為這種情況對於小的、錯誤指定的情況可能是最糟糕的（即這裡的殘差是統一的而不是正常的）。不過，更系統地探索這一點會很有趣……

引用自：https://stats.stackexchange.com/questions/384528