lme () i lmer () dają sprzeczne wyniki

Pracowałem z niektórymi danymi, które mają pewne problemy z powtarzanymi pomiarami. W ten sposób zauważyłem bardzo różne zachowanie między danymi testowymi lme()i ich lmer()używanie i chcę wiedzieć, dlaczego.

Fałszywy zestaw danych, który utworzyłem, zawiera pomiary wzrostu i masy ciała dla 10 osób, wykonane dwukrotnie. Ustawiłem dane tak, aby między badanymi istniał pozytywny związek między wzrostem a wagą, ale negatywny związek między powtarzanymi pomiarami w obrębie każdej osoby.

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

Oto wykres danych z liniami łączącymi dwa pomiary od każdego z osobna.

Więc wpadłem dwa modele, jeden z lme()z nlmeopakowania i jeden z lmer()od lme4. W obu przypadkach przeprowadziłem regresję masy w stosunku do wzrostu z losowym efektem ID w celu kontroli powtarzanych pomiarów każdego osobnika.

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

Te dwa modele często (choć nie zawsze w zależności od nasion) generowały zupełnie inne wyniki. Widziałem, gdzie generują nieco różne oszacowania wariancji, obliczają różne stopnie swobody itp., Ale tutaj współczynniki są w przeciwnych kierunkach.

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

Aby zilustrować wizualnie, modeluj za pomocą lme()

I modeluj z lmer()

Dlaczego te modele tak bardzo się różnią?

tl; dr, jeśli zmienisz optymalizator na „nloptwrap”, myślę, że pozwoli to uniknąć tych problemów (prawdopodobnie).

Gratulacje, znalazłeś jeden z najprostszych przykładów wielu optymów w problemie z estymacją statystyczną! Parametrem, który lme4wykorzystuje wewnętrznie (jest to wygodne do zilustrowania) jest skalowane odchylenie standardowe efektów losowych, tj. Odchylenie standardowe std w grupie podzielone przez resztkowe odchylenie standardowe.

Wyodrębnij te wartości dla oryginału lmei lmerpasuje:

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

Zainstaluj ponownie za pomocą innego optymalizatora (prawdopodobnie będzie to ustawienie domyślne w następnej wersji lme4):

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

Mecze lme... zobaczmy, co się dzieje. Funkcja odchylenia (prawdopodobieństwo logarytmu -2 *), lub w tym przypadku analogiczna funkcja kryterium REML, dla LMM z pojedynczym efektem losowym wymaga tylko jednego argumentu, ponieważ parametry efektu stałego są profilowane ; można je obliczyć automatycznie dla danej wartości odchylenia standardowego 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()

I nadal prześladować dalej nad tym i prowadził napady dla losowych nasion od 1 do 1000, montażu lme, lmeroraz lmer+ nloptwrap dla każdego przypadku. Oto liczby z 1000, gdzie dana metoda otrzymuje odpowiedzi, które są co najmniej 0,001 jednostek dewiacji gorsze niż inne ...

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

Innymi słowy: (1) nie ma metody, która zawsze działałaby najlepiej; (2) lmerprzy domyślnym optymalizatorze jest najgorszy (nie udaje się w 1/3 czasu); (3) lmerz „nloptwrap” jest najlepszy (gorszy niż lme4% czasu, rzadko gorszy niż lmer).

Aby być nieco uspokajającym, myślę, że ta sytuacja może być najgorsza w przypadku małych, źle określonych przypadków (tj. Błąd resztkowy jest tu raczej jednolity niż normalny). Interesujące byłoby jednak zbadanie tego bardziej systematycznie ...