Dlaczego do wybierania zagnieżdżonych modeli var-covar należy używać REML (zamiast ML)?

16

Różne opisy wyboru modeli losowych efektów liniowych modeli mieszanych instruują użycie REML. Znam różnicę między REML i ML na pewnym poziomie, ale nie rozumiem, dlaczego REML powinien być używany, ponieważ ML jest stronniczy. Na przykład, czy błędem jest przeprowadzanie LRT na parametrze wariancji normalnego modelu dystrybucji przy użyciu ML (patrz kod poniżej)? Nie rozumiem, dlaczego ważniejsze jest bycie obiektywnym niż ML w wyborze modelu. Myślę, że ostateczną odpowiedzią musi być „ponieważ wybór modelu działa lepiej z REML niż z ML”, ale chciałbym wiedzieć coś więcej. Nie przeczytałem pochodnych LRT i AIC (nie jestem wystarczająco dobry, aby je dokładnie zrozumieć), ale jeśli REML jest wyraźnie użyty w pochodnych, po prostu wiedząc, że to będzie w rzeczywistości wystarczające (np.

n <- 100
a <- 10
b <- 1
alpha <- 5
beta <- 1
x <- runif(n,0,10)
y <- rnorm(n,a+b*x,alpha+beta*x)

loglik1 <- function(p,x,y){
   a <- p[1]
   b <- p[2]
   alpha <- p[3]
  -sum(dnorm(y,a+b*x,alpha,log=T))
}

loglik2 <- function(p,x,y){
   a <- p[1]
   b <- p[2]
   alpha <- p[3]
   beta <- p[4]
  -sum(dnorm(y,a+b*x,alpha+beta*x,log=T))
}

m1 <- optim(c(a,b,alpha),loglik1,x=x,y=y)$value
m2 <- optim(c(a,b,alpha,beta),loglik2,x=x,y=y)$value
D <- 2*(m1-m2)
1-pchisq(D,df=1) # p-value
spierać się
źródło
1
O REML i AIC powinieneś spojrzeć na to pytanie .
Elvis

Odpowiedzi:

13

Bardzo krótka odpowiedź: REML jest ML, więc test oparty na REML i tak jest poprawny. Ponieważ oszacowanie parametrów wariancji za pomocą REML jest lepsze, naturalne jest jego użycie.

Dlaczego REML to ML? Rozważ np. Model z X R n × p , Z R n × q , a β R p jest wektorem stałych efektów, u N ( 0 , τ I q ) jest wektorem efektów losowych, a e N ( 0 , σ 2 I n )

Y=Xβ+Zu+e
XRn×pZRn×qβRpuN(0,τIq)eN(0,σ2In). Ograniczone prawdopodobieństwo można uzyskać, biorąc pod uwagę kontrasty celu „usunięcia” ustalonych efektów. Dokładniej, niech C R ( n - p ) × n , tak aby C X = 0 i C C = I n - p (to znaczy kolumny C ' są ortonormalną podstawą przestrzeni wektorowej ortogonalnej do przestrzeń generowana przez kolumny X ); następnie C Y = C Z u +npCR(np)×nCX=0CC=InpCX z ϵ N ( 0 , σ 2 I n - p ) , a prawdopodobieństwo dla τ , σ 2 dla C Y jest ograniczonym prawdopodobieństwem.
CY=CZu+ϵ
ϵN(0,σ2Inp)τ,σ2CY
Elvis
źródło
Dobra odpowiedź (+1), czy mam rację twierdząc, że macierz zależy od modelu średniej? Więc możesz porównywać tylko szacunki REML dla tej samej macierzy C ? CC
Tak, zależy od X (przeredaguję odpowiedź za minutę, aby było jasne), więc twoje modele zagnieżdżone muszą mieć te same zmienne z ustalonymi efektami. CX
Elvis
REML nie jest ML! ML jest jednoznacznie określony dla danego modelu prawdopodobieństwa ale REML zależy od parametryzacji efektów stałych. Zobacz np. Ten komentarz Douga Batesa (a także wiele historycznych na temat modeli mieszanych R-SIG).
Livius
1
@Livius Myślę, że moja odpowiedź wystarczająco jasno określa, jak skonstruowane jest ograniczone prawdopodobieństwo. Jest to prawdopodobieństwo, po prostu nie jest to prawdopodobne, biorąc pod uwagę obserwowane w modelu zapisanym w pierwszym wyświetlanym równaniu, ale biorąc pod uwagę rzutowany wektor C Y w modelu zapisanym w drugim wyświetlanym równaniu. REML to ML uzyskane z tego prawdopodobieństwa. YCY
Elvis
2
Myślę, że jest to w pewnym sensie protesty DBatesa w tej sprawie: jest to inny model i jest to model, dla którego porównania są trudne, ponieważ model i parametryzacja są ze sobą powiązane. Więc nie jesteś obliczeniowych w ml dla oryginalnego modelu, ale w ml dla innego modelu wynikającego z określonego parametryzacji oryginalnego modelu. Dlatego modele wyposażone w REML z zagnieżdżonymi strukturami efektów stałych nie są już modelami zagnieżdżonymi (jak wspomniano powyżej). Ale modele z dopasowaniem ML są nadal zagnieżdżone, ponieważ maksymalizujesz prawdopodobieństwo dla określonego modelu.
Livius
9

Testy ilorazu wiarygodności to statystyczne testy hipotez oparte na stosunku dwóch prawdopodobieństw. Ich właściwości są powiązane z oszacowaniem maksymalnego prawdopodobieństwa (MLE). (patrz np. szacunek maksymalnej wiarygodności (MLE) w kategoriach laika) ).

W twoim przypadku (patrz pytanie) chcesz „wybrać” spośród dwóch zagnieżdżonych modeli var-covar, powiedzmy, że chcesz wybrać między modelem, w którym var-covar wynosi a modelem, w którym var-covar jest Σ sΣgΣs , gdzie druga (prosty model) jest szczególnym przypadkiem pierwszej (ogólna jeden).

Test oparty jest na wskaźnik prawdopodobieństwa . W przypadku, gdy Σ s i Σ gLR=2(log(Ls(Σ^s))log(Lg(Σ^g))Σ^sΣ^g są estymatory maksymalnego prawdopodobieństwa.

Statystyka jest asymptotycznie (!) Χ 2LR χ2 .

Estymatory maksymalnego prawdopodobieństwa są znane jako spójne, jednak w wielu przypadkach są one tendencyjne. Jest to przypadek dla estymatorów MLE dla i Ď g , może być pokazać, że są tendencyjne. Wynika to z tego, że są one obliczane przy użyciu średniej uzyskanej z danych, tak że rozrzut wokół tej „szacowanej średniej” jest mniejszy niż rozrzut wokół prawdziwej średniej (patrz np. Intuicyjne wyjaśnienie dzielenia przez n - 1 przy obliczaniu odchylenia standardowego ? )Σ^sΣ^gn1

Statystykę powyżej χ 2 w dużych próbkach, to tylko ze względu na fakt, że w dużych próbkach, Σ s i Σ g zbiegają się ich wartościami rzeczywistymi (MLE zgodnych). (Uwaga: w powyższym linku, w przypadku bardzo dużych próbek, dzielenie przez n lub przez (n-1) nie będzie miało znaczenia)LRχ2Σ^sΣ^g

W przypadku mniejszych próbek, MLE szacunki Σ s i Σ g będą tendencyjne i dlatego rozkład L R będzie odbiegać od × 2 , podczas gdy szacunki REML da bezstronne szacunki Σ s i Ď g , więc jeśli używasz w odniesieniu do wyboru modelu vAR KOWARIANCJA The REML szacuje następnie L R się na mniejsze próbki lepiej aproksymowane × 2 .Σ^sΣ^gLRχ2ΣsΣgLRχ2

Zauważ, że REML powinien być używany tylko do wybierania spośród zagnieżdżonych struktur var-covar modeli o tej samej średniej, w przypadku modeli o różnych średnich, REML nie jest odpowiedni, w przypadku modeli o różnych środkach należy używać ML.


źródło
Stwierdzenie „Statystyka LR jest asymptotycznie (!) Χ2” nie jest w tym przypadku prawdziwe. Jest tak, ponieważ jeśli jest zagnieżdżony w Ď g , a następnie Σ y znajduje się na granicy Ď g . W tym przypadku rozkład χ 2 nie obowiązuje. Na przykład zobacz tutajΣsΣgΣsΣgχ2
Cliff AB
@Cliff AB, oto wyjaśnienie poniżej tego oświadczenia i jest to powód, dla którego musisz użyć REML.
-4

Mam odpowiedź, która ma więcej wspólnego ze zdrowym rozsądkiem niż ze statystyką. Jeśli spojrzysz na PROC MIXED w SAS, oszacowania można dokonać za pomocą sześciu metod:

http://support.sas.com/documentation/cdl/en/statug/63033/HTML/default/viewer.htm#statug_mixed_sect008.htm

ale REML jest ustawieniem domyślnym. Dlaczego? Najwyraźniej praktyczne doświadczenie pokazało, że ma najlepszą wydajność (np. Najmniejszą szansę na problemy z konwergencją). Dlatego jeśli twój cel jest osiągalny dzięki REML, wówczas warto zastosować REML w przeciwieństwie do pozostałych pięciu metod.

James
źródło
2
Musi to być związane z „teorią dużej próby” i tendencyjnością szacunków MLE, patrz moja odpowiedź.
1
„To domyślne w SAS” nie jest akceptowalną odpowiedzią na pytanie „dlaczego” na tej stronie.
Paul,
Wartości p dla modeli mieszanych dostarczanych domyślnie przez SAS nie są z założenia dostępne w bibliotece lme4 dla R, ponieważ są niewiarygodne ( stat.ethz.ch/pipermail/r-help/2006-May/094765.html ). Tak więc „domyślny SAS” może być nawet błędny.
Tim