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

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

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\beta + Zu + e \def\R{\mathbb{R}}$$ $X\in\R^{n\times p}$$Z\in\R^{n\times q}$$\beta \in \R^p$$u \sim \mathcal N(0, \tau I_q)$$e \sim \mathcal N(0, \sigma^2 I_n)$. 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 $n-p$$C \in \R^{(n-p)\times n}$$CX = 0$$CC' = I_{n-p}$$C'$$X$ z ϵ ∼ N ( 0 , σ 2 I n - p ) , a prawdopodobieństwo dla τ , σ 2 dla C Y jest ograniczonym prawdopodobieństwem. $$CY = CZ u + \epsilon$$ $\epsilon \sim \mathcal N(0, \sigma^2 I_{n-p})$$\tau, \sigma^2$$CY$

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 Σ $\Sigma_g$$\Sigma_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 Σ$LR=-2 (log(\mathcal{L}_s(\hat{\Sigma}_s)) - log(\mathcal{L}_g(\hat{\Sigma}_g) )$$\hat{\Sigma}_s$$\hat{\Sigma}_g$ są estymatory maksymalnego prawdopodobieństwa.

Statystyka jest asymptotycznie (!) Χ $LR$ $\chi^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 ? )$\hat{\Sigma}_s$$\hat{\Sigma}_g$$n-1$

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$$\chi^2$$\hat{\Sigma}_s$$\hat{\Sigma}_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 .$\hat{\Sigma}_s$$\hat{\Sigma}_g$$LR$$\chi^2$$\Sigma_s$$\Sigma_g$$LR$$\chi^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.

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.