O que o comando anova () faz com um objeto de modelo lmer?

Espero que esta seja uma pergunta que alguém aqui possa responder sobre a natureza da decomposição de somas de quadrados de um modelo de efeitos mistos adequado lmer(do pacote lme4 R).

Primeiramente, devo dizer que estou ciente da controvérsia ao usar essa abordagem e, na prática, seria mais provável que eu usasse um LRT com inicialização para comparar modelos (como sugerido por Faraway, 2006). No entanto, estou intrigado com a forma de replicar os resultados e, por minha própria sanidade, pensei em perguntar aqui.

Basicamente, estou começando a entender o uso de modelos de efeitos mistos adequados ao lme4pacote. Eu sei que você pode usar o anova()comando para fornecer um resumo do teste seqüencial dos efeitos fixos no modelo. Até onde eu sei, isso é o que Faraway (2006) chama de abordagem dos quadrados médios esperados. O que eu quero saber é como são calculadas as somas dos quadrados?

Eu sei que eu poderia pegar os valores estimados de um modelo específico (usando coef()), supor que eles são fixos e, em seguida, fazer testes usando a soma dos quadrados dos resíduos do modelo com e sem os fatores de interesse. Isso é bom para um modelo que contém um único fator dentro do assunto. No entanto, ao implementar um design de plotagem dividida, o valor da soma dos quadrados que recebo é equivalente ao valor produzido por R usando aov()uma Error()designação apropriada . No entanto, isso não é o mesmo que a soma dos quadrados produzidos pelo anova()comando no objeto de modelo, apesar do fato de que as relações F são as mesmas.

É claro que isso faz todo sentido, pois não há necessidade de Error()estratos em um modelo misto. No entanto, isso deve significar que as somas de quadrados são penalizadas de alguma forma em um modelo misto, a fim de fornecer índices F apropriados. Como isso é alcançado? E como o modelo de alguma forma corrige a soma dos quadrados entre as parcelas, mas não corrige a soma dos quadrados entre as parcelas. Evidentemente, isso é algo que é necessário para uma ANOVA de parcelas subdivididas clássica que foi alcançada designando diferentes valores de erro para os diferentes efeitos, então como um modelo de efeito misto permite isso?

Basicamente, eu quero poder replicar os resultados do anova()comando aplicado a um objeto de modelo mais antigo para verificar os resultados e meu entendimento; no entanto, no momento, posso conseguir isso para um design normal dentro do assunto, mas não para a divisão de divisão. desenho de enredo e não consigo descobrir por que esse é o caso.

Como um exemplo:

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               

Como pode ser visto acima, todos os índices F concordam. As somas de quadrados por variedade também concordam. No entanto, as somas de quadrados para irrigação não concordam, no entanto, parece que a produção mais recente é escalada. Então, o que o comando anova () realmente faz?

Use a fonte, Luke. Podemos espiar dentro da função ANOVA fazendo getAnywhere(anova.Mermod). A primeira parte dessa função é comparar dois modelos diferentes. A anova nos efeitos fixos vem no grande elsebloco no segundo tempo:

 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é a saída final. Começamos a calcular a soma dos quadrados na linha 5: ss <- as.vector ...O código multiplica os parâmetros fixos (in beta) por uma matriz triangular superior; depois esquadrinha cada termo. Aqui está a matriz triangular superior para o exemplo de irrigação. Cada linha corresponde a um dos cinco parâmetros de efeitos fixos (interceptação, 3 graus de liberdade para irrigação, 1 df para variedade).

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

A primeira linha fornece a soma dos quadrados para a interceptação e a última fornece o SS para o efeito de variedade dentro dos campos. As linhas 2-4 envolvem apenas os 3 parâmetros para os níveis de irrigação; portanto, a pré-multiplicação fornece três partes do SS para irrigação.

Essas peças não são interessantes por si só, pois vêm do contraste padrão de tratamento em R, mas na linha ss <- unlist(lapply(split ....Bates recolhe bits de somas de quadrados de acordo com o número de níveis e a quais fatores eles se referem. Há muita contabilidade acontecendo aqui. Também obtemos os graus de liberdade (que são 3 para irrigação). Então, ele obtém os quadrados médios que aparecem na impressão de anova. Finalmente, ele divide todos os quadrados de suas médias pela variação residual dentro dos grupos sigma(object)^2,.

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

Assintoticamente, as estimativas dos efeitos fixos têm distribuição:

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

Observe que você não teria as mesmas estatísticas F se os dados estivessem desequilibrados. Você também não obteria as mesmas estatísticas F se tivesse usado ML em vez de REML.

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

Curiosamente, Bates e Pinheiro recomendam o uso da ANOVA para ajustar dois modelos e fazer um teste de razão de verossimilhança. Este último tende a ser anti-conservador.

$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

Como você pode ver, as somas de quadrados para os parâmetros de irrigação agora também contêm alguns dos varietyefeitos.