Co robi komenda anova () z obiektem modelu Lmer?

Mam nadzieję, że jest to pytanie, na które ktoś tutaj może mi odpowiedzieć na temat natury rozkładania sum kwadratów z dopasowanego modelu mieszanego lmer(z pakietu lme4 R).

Po pierwsze powinienem powiedzieć, że jestem świadomy kontrowersji związanych z zastosowaniem tego podejścia, aw praktyce bardziej prawdopodobne byłoby użycie LRT z bootstrapem do porównywania modeli (jak sugeruje Faraway, 2006). Jednak zastanawiam się, jak odtworzyć wyniki, więc dla własnego zdrowia psychicznego pomyślałem, że zapytam tutaj.

Zasadniczo zajmuję się używaniem modeli z efektami mieszanymi dopasowanymi do lme4pakietu. Wiem, że możesz użyćanova() polecenia, aby uzyskać podsumowanie sekwencyjnego testowania efektów stałych w modelu. O ile wiem, to właśnie Faraway (2006) określa jako podejście „oczekiwanych średnich kwadratów”. Chcę wiedzieć, w jaki sposób obliczane są sumy kwadratów?

Wiem, że mogłem wziąć oszacowane wartości z konkretnego modelu (używając coef()), założyć, że są one ustalone, a następnie wykonać testy z wykorzystaniem sum kwadratów reszt modelu z i bez czynników zainteresowania. To dobrze w przypadku modelu zawierającego pojedynczy czynnik wewnątrz podmiotu. Jednak przy wdrażaniu projektu podzielonego wykresu sumy wartości kwadratów, które otrzymuję, są równoważne wartości wytworzonej przez R przy użyciu aov()odpowiedniego Error()oznaczenia. Nie jest to jednak to samo, co suma kwadratów wyprodukowanych przezanova() polecenie na obiekcie modelowym, mimo że współczynniki F są takie same.

Oczywiście ma to sens, ponieważ nie ma takiej potrzeby Error() warstw w modelu mieszanym. Musi to jednak oznaczać, że sumy kwadratów są w jakiś sposób karane w modelu mieszanym, aby zapewnić odpowiednie współczynniki F. Jak to się osiąga? I w jaki sposób model w jakiś sposób koryguje sumę kwadratów między wykresami, ale nie koryguje sumy kwadratów wewnątrz wykresów. Najwyraźniej jest to coś, co jest konieczne w przypadku klasycznej analizy wariancji ANOVA, która została osiągnięta poprzez wyznaczenie różnych wartości błędów dla różnych efektów, więc w jaki sposób pozwala na to model mieszany?

Zasadniczo chcę być w stanie samodzielnie odtworzyć wyniki anova()polecenia zastosowanego do obiektu modelu Lmer, aby zweryfikować wyniki i moje zrozumienie, ale obecnie mogę to osiągnąć dla normalnego projektu wewnątrz podmiotu, ale nie dla podziału projekt fabuły i nie mogę się dowiedzieć, dlaczego tak jest.

Jako przykład:

library(faraway)
library(lme4)
data(irrigation)

anova(lmer(yield ~ irrigation + variety + (1|field), data = irrigation))

Analysis of Variance Table
           Df Sum Sq Mean Sq F value
irrigation  3 1.6605  0.5535  0.3882
variety     1 2.2500  2.2500  1.5782

summary(aov(yield ~ irrigation + variety + Error(field/irrigation), data = irrigation))

Error: field
           Df Sum Sq Mean Sq F value Pr(>F)
irrigation  3  40.19   13.40   0.388  0.769
Residuals   4 138.03   34.51               

Error: Within
          Df Sum Sq Mean Sq F value Pr(>F)
variety    1   2.25   2.250   1.578  0.249
Residuals  7   9.98   1.426               

Jak widać powyżej, wszystkie współczynniki F są zgodne. Sumy kwadratów dla różnorodności również się zgadzają. Jednak sumy kwadratów do nawadniania nie są zgodne, jednak wydaje się, że wydajność lmera jest skalowana. Co właściwie robi komenda anova ()?

Użyj źródła, Luke. W ten sposób możemy zajrzeć do funkcji ANOVA getAnywhere(anova.Mermod). Pierwsza część tej funkcji polega na porównaniu dwóch różnych modeli. Anova na ustalone efekty pojawia się w dużym elsebloku w drugiej połowie:

 dc <- getME(object, "devcomp")
        X <- getME(object, "X")
        asgn <- attr(X, "assign")
        stopifnot(length(asgn) == (p <- dc$dims[["p"]]))
            ss <- as.vector(object@pp$RX() %*% object@beta)^2
        names(ss) <- colnames(X)
        terms <- terms(object)
        nmeffects <- attr(terms, "term.labels")[unique(asgn)]
        if ("(Intercept)" %in% names(ss)) 
            nmeffects <- c("(Intercept)", nmeffects)
        ss <- unlist(lapply(split(ss, asgn), sum))
        stopifnot(length(ss) == length(nmeffects))
        df <- vapply(split(asgn, asgn), length, 1L)
        ms <- ss/df
        f <- ms/(sigma(object)^2)
        table <- data.frame(df, ss, ms, f)
        dimnames(table) <- list(nmeffects, c("Df", "Sum Sq", 
            "Mean Sq", "F value"))
        if ("(Intercept)" %in% nmeffects) 
            table <- table[-match("(Intercept)", nmeffects), 
                ]
        attr(table, "heading") <- "Analysis of Variance Table"
        class(table) <- c("anova", "data.frame")
        table

objectjest wyjściem Lmer. Zaczynamy obliczać sumę kwadratów w linii 5: ss <- as.vector ...Kod mnoży ustalone parametry (in beta) przez górną macierz trójkątną; następnie kwadraty dla każdego terminu. Oto górna trójkątna matryca dla przykładu nawadniania. Każdy rząd odpowiada jednemu z pięciu parametrów efektów stałych (punkt przecięcia, 3 stopnie swobody dla nawadniania, 1 df dla odmiany).

zapsmall(irrigation.lmer@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]  [,5]
[1,] 0.813 0.203  0.203  0.203 0.407
[2,] 0.000 0.352 -0.117 -0.117 0.000
[3,] 0.000 0.000  0.332 -0.166 0.000
[4,] 0.000 0.000  0.000  0.287 0.000
[5,] 0.000 0.000  0.000  0.000 2.000

Pierwszy rząd daje sumę kwadratów dla punktu przecięcia, a ostatni daje SS dla efektu różnorodności wewnątrz pola. Rzędy 2-4 dotyczą tylko 3 parametrów poziomów nawadniania, więc pomnożenie wstępne daje trzy części SS do nawadniania.

Kawałki te nie są same w sobie interesujące, ponieważ pochodzą z domyślnego kontrastu leczenia w R, ale w linii ss <- unlist(lapply(split ....Bates zbiera kawałki sum kwadratów zgodnie z liczbą poziomów i do jakich czynników się odnoszą. Tutaj dzieje się dużo księgowości. Otrzymujemy również stopnie swobody (które są 3 dla nawadniania). Następnie otrzymuje średnie kwadraty, które pojawiają się na wydruku anova. W końcu dzieli wszystkie jego środków przez kwadraty wewnątrz grup resztkowej wariancji sigma(object)^2.

lmeraovlmerRX$R_{00}$$\sqrt{\sigma^2/\sigma_f^2}$$\sigma_f^2$

Asymptotycznie szacunki ustalonych efektów mają rozkład:

$R_{00}\hat{\beta}$$\beta=0$$\sigma^2$$\sigma^2$$\sigma^2$$R_{00}$$\sigma^2$

Pamiętaj, że nie uzyskałbyś takich samych statystyk F, gdyby dane były niezrównoważone. Nie uzyskałbyś tych samych statystyk F, gdybyś użył ML zamiast REML.

aov$\sigma^2$$\sigma_f^2$$\sigma^2$$\sigma_f^2$

Co ciekawe, Bates i Pinheiro zalecają użycie ANOVA zamiast dopasowania dwóch modeli i wykonania testu współczynnika wiarygodności. Ten ostatni ma tendencję do zachowania konserwatywności.

$R_{00}$

zapsmall(fit2@pp$RX(), digits = 3)
      [,1]  [,2]   [,3]   [,4]   [,5]
[1,] 0.816 0.205  0.205  0.205  0.457
[2,] 0.000 0.354 -0.119 -0.119 -0.029
[3,] 0.000 0.000  0.334 -0.168 -0.040
[4,] 0.000 0.000  0.000  0.288 -0.071
[5,] 0.000 0.000  0.000  0.000  1.874

Jak widać, sumy kwadratów dla parametrów irygacji zawierają teraz także pewien varietyefekt.