Что делает команда anova () с объектом модели lmer?

Надеюсь, что это вопрос, который кто-то здесь может ответить для меня о природе разложения сумм квадратов из модели смешанных эффектов lmer(из пакета lme4 R).

Прежде всего я должен сказать, что мне известно о противоречиях с использованием этого подхода, и на практике я бы с большей вероятностью использовал загрузочный LRT для сравнения моделей (как это было предложено Faraway, 2006). Тем не менее, я озадачен тем, как воспроизвести результаты, и поэтому для собственного здравомыслия я подумал, что я спрошу здесь.

По сути, я справляюсь с использованием моделей со смешанными эффектами, подходящих под lme4пакет. Я знаю, что вы можете использовать anova()команду, чтобы получить сводную информацию о последовательном тестировании фиксированных эффектов в модели. Насколько я знаю, именно это Faraway (2006) называет подходом «ожидаемых средних квадратов». Я хочу знать, как рассчитываются суммы квадратов?

Я знаю, что могу взять оценочные значения из конкретной модели (используя coef()), предположить, что они являются фиксированными, а затем провести тесты, используя суммы квадратов невязок модели с интересующими факторами и без них. Это хорошо для модели, содержащей один внутрисубъектный фактор. Тем не менее, при реализации сплит-графика значение суммы квадратов, которое я получаю, эквивалентно значению, полученному с помощью R aov()с соответствующим Error()обозначением. Однако это не то же самое, что суммы квадратов, полученных anova()командой на модельном объекте, несмотря на то, что F-отношения одинаковы.

Конечно, это имеет полный смысл, так как Error()в смешанной модели нет необходимости в стратах. Однако это должно означать, что суммы квадратов как-то штрафуются в смешанной модели, чтобы обеспечить соответствующие F-отношения. Как это достигается? И как модель каким-то образом исправляет сумму квадратов между участками, но не исправляет сумму квадратов в пределах сюжета. Очевидно, что это то, что необходимо для классического сплит-графика ANOVA, который был достигнут путем назначения разных значений ошибок для разных эффектов, так как модель смешанного эффекта позволяет это сделать?

По сути, я хочу иметь возможность реплицировать результаты из anova()команды, примененной к объекту модели lmer, самостоятельно, чтобы проверить результаты и мое понимание, однако в настоящее время я могу добиться этого для обычного внутрисубъектного дизайна, но не для разделения. дизайн сюжета, и я не могу понять, почему это так.

В качестве примера:

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               

Как видно выше, все F-соотношения согласны. Суммы квадратов для разнообразия также согласуются. Тем не менее, суммы квадратов для орошения не согласуются, однако, похоже, что объем производства лмера масштабируется. Так что же на самом деле делает команда anova ()?

Используйте Источник, Люк. Мы можем заглянуть внутрь функции ANOVA, выполнив getAnywhere(anova.Mermod). Первая часть этой функции предназначена для сравнения двух разных моделей. Анова на фиксированные эффекты приходит во elseвтором блоке большим блоком:

 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

objectэто выход лмера. Мы начинаем вычислять сумму квадратов в строке 5: ss <- as.vector ...код умножает фиксированные параметры (in beta) на верхнюю треугольную матрицу; затем возводит в квадрат каждый член. Вот та верхняя треугольная матрица для примера орошения. Каждый ряд соответствует одному из пяти фиксированных параметров эффектов (перехват, 3 степени свободы для орошения, 1 df для разнообразия).

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

В первой строке дается сумма квадратов для перехвата, а в последней - SS для эффекта разнообразия внутри полей. Строки 2-4 включают только 3 параметра для уровней орошения, поэтому предварительное умножение дает вам три части СС для орошения.

Эти фрагменты сами по себе не интересны, поскольку они основаны на контрасте обработки по умолчанию в R, но в строке ss <- unlist(lapply(split ....Бейтс собирает биты сумм квадратов в зависимости от количества уровней и факторов, к которым они относятся. Здесь много бухгалтерского учета. Мы также получаем степени свободы (3 для орошения). Затем он получает средние квадраты, которые отображаются на распечатке anova. Наконец, он делит все свои средние квадраты на остаточную дисперсию внутри группы sigma(object)^2.

Итак, что происходит? Философия lmerне имеет ничего общего с методом подхода моментов, который используется aov. Идея lmerсостоит в том, чтобы максимизировать предельную вероятность, полученную путем интегрирования невидимых случайных эффектов. В этом случае случайный уровень рождаемости каждого поля. Глава 2 Пиньейру и Бейтса описывает уродство этого процесса. RXМатрица используется для получения суммы квадратов их матрица из уравнения 2.17, стр 70 текста. Эта матрица получена из набора QR-разложений на (среди прочего) матрицу дизайна случайных эффектов и , где - дисперсия эффекта поля , Это тот недостающий фактор$R_{00}$$\sqrt{\sigma^2/\sigma_f^2}$$\sigma_f^2$ Вы спрашивали о, но это не входит в решение прозрачным или простым способом.

Асимптотически оценки фиксированных эффектов имеют распределение:

Это означает, что компоненты распределены независимо. Если , квадраты этих членов (разделенные на ) следуют распределению хи-квадрат. Разделение на оценку (еще один хи-квадрат при делении на ) дает F-статистику. Мы не делим среднеквадратичную ошибку для анализа между группами, потому что это не имеет ничего общего с тем, что здесь происходит. Нам нужно сравнить суммы квадратов, полученные через , сравнив их с оценкой .$R_{00}\hat{\beta}$$\beta=0$$\sigma^2$$\sigma^2$$\sigma^2$$R_{00}$$\sigma^2$

Обратите внимание, что вы не получили бы ту же F статистику, если бы данные были несбалансированными. Также вы не получили бы такую ​​же статистику F, если бы использовали ML вместо REML.

Идея aovзаключается в том, что ожидаемые средние квадраты для орошения являются функцией , и эффектов орошения. Ожидаемый средний квадрат для остатков поля является функцией и . Когда ирригационные эффекты равны 0, обе эти величины оценивают одно и то же, а их отношение соответствует распределению F.$\sigma^2$$\sigma_f^2$$\sigma^2$$\sigma_f^2$

Интересно, что Бейтс и Пинейро рекомендуют использовать ANOVA для подгонки двух моделей и проведения теста отношения правдоподобия. Последний имеет тенденцию быть антиконсервативным.

Конечно, если данные не сбалансированы, уже не ясно, какую именно гипотезу проверяет ANOVA. Я удалил первое наблюдение из данных орошения и восстановил. Вот снова матрица :$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

Как вы можете видеть, суммы квадратов для параметров орошения теперь также содержат некоторые varietyэффекты.