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 nlme
opakowania 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ą?
źródło
Odpowiedzi:
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
lme4
wykorzystuje 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
lme
ilmer
pasuje:Zainstaluj ponownie za pomocą innego optymalizatora (prawdopodobnie będzie to ustawienie domyślne w następnej wersji
lme4
):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.I nadal prześladować dalej nad tym i prowadził napady dla losowych nasion od 1 do 1000, montażu
lme
,lmer
orazlmer
+ 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 ...Innymi słowy: (1) nie ma metody, która zawsze działałaby najlepiej; (2)
lmer
przy domyślnym optymalizatorze jest najgorszy (nie udaje się w 1/3 czasu); (3)lmer
z „nloptwrap” jest najlepszy (gorszy niżlme
4% 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 ...
źródło