lme () dan lmer () memberikan hasil yang bertentangan

Saya telah bekerja dengan beberapa data yang memiliki beberapa masalah dengan pengukuran berulang. Dengan melakukan itu saya memperhatikan perilaku yang sangat berbeda antara lme()dan lmer()menggunakan data pengujian saya dan ingin tahu mengapa.

Kumpulan data palsu yang saya buat memiliki pengukuran tinggi dan berat badan untuk 10 subjek, masing-masing diambil dua kali. Saya mengatur data sehingga antara subjek akan ada hubungan positif antara tinggi dan berat badan, tetapi hubungan negatif antara tindakan berulang dalam setiap individu.

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

Berikut adalah sebidang data, dengan garis yang menghubungkan dua pengukuran dari masing-masing individu.

Jadi saya menjalankan dua model, satu dengan lme()dari nlmepaket dan satu dengan lmer()dari lme4. Dalam kedua kasus saya menjalankan regresi berat terhadap tinggi badan dengan efek acak ID untuk mengontrol pengukuran berulang setiap individu.

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

Kedua model ini seringkali (walaupun tidak selalu tergantung pada benih) menghasilkan hasil yang sangat berbeda. Saya telah melihat di mana mereka menghasilkan estimasi varians yang sedikit berbeda, menghitung derajat kebebasan yang berbeda, dll., Tetapi di sini koefisien berada di arah yang berlawanan.

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

Untuk menggambarkan secara visual, modelkan dengan lme()

Dan model dengan lmer()

Mengapa model-model ini sangat berbeda?

tl; dr jika Anda mengubah pengoptimal menjadi "nloptwrap" Saya pikir itu akan menghindari masalah ini (mungkin).

Selamat, Anda telah menemukan salah satu contoh paling sederhana dari beberapa optima dalam masalah estimasi statistik! Parameter yang lme4menggunakan internal (sehingga mudah untuk ilustrasi) adalah standar deviasi skala efek acak, yaitu std dev antar-kelompok dibagi dengan sisa std dev.

Ekstrak nilai-nilai ini untuk yang asli lmedan lmercocok:

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

Pasang kembali dengan pengoptimal lain (ini mungkin akan menjadi default di rilis berikutnya lme4):

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

Cocok lme... mari kita lihat apa yang terjadi. Fungsi penyimpangan (-2 * log likelihood), atau dalam hal ini fungsi kriteria REML analog, untuk LMM dengan efek acak tunggal hanya memerlukan satu argumen, karena parameter efek tetap diprofilkan ; mereka dapat dihitung secara otomatis untuk nilai deviasi standar RE yang diberikan.

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()

Aku terus terobsesi lebih lanjut atas ini dan berlari cocok untuk bibit random 1 sampai 1000, pas lme, lmerdan lmer+ nloptwrap untuk setiap kasus. Berikut adalah angka-angka dari 1000 di mana metode yang diberikan mendapat jawaban yang setidaknya 0,001 unit penyimpangan lebih buruk daripada yang lain ...

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

Dengan kata lain, (1) tidak ada metode yang selalu bekerja paling baik; (2) lmerdengan pengoptimal default yang terburuk (gagal sekitar 1/3 waktu); (3) lmerdengan "nloptwrap" adalah yang terbaik (lebih buruk dari lme4% dari waktu, jarang lebih buruk daripada lmer).

Untuk sedikit meyakinkan, saya pikir situasi ini kemungkinan terburuk untuk kasus-kasus kecil yang tidak ditentukan (mis. Kesalahan residual di sini seragam daripada Normal). Akan menarik untuk mengeksplorasi ini lebih sistematis meskipun ...