¿Qué hace el comando anova () con un objeto modelo lmer?

Con suerte, esta es una pregunta que alguien aquí puede responder por mí sobre la naturaleza de la descomposición de sumas de cuadrados a partir de un modelo de efectos mixtos lmer(del paquete lme4 R).

En primer lugar, debo decir que soy consciente de la controversia con el uso de este enfoque, y en la práctica sería más probable que use un LRT de arranque para comparar modelos (como lo sugiere Faraway, 2006). Sin embargo, estoy desconcertado sobre cómo replicar los resultados, y por mi propia cordura, pensé en preguntar aquí.

Básicamente, me estoy acostumbrando al uso de modelos de efectos mixtos que se ajustan al lme4paquete. Sé que puede usar el anova()comando para dar un resumen de la prueba secuencial de los efectos fijos en el modelo. Hasta donde yo sé, esto es a lo que Faraway (2006) se refiere como el enfoque de "cuadrados medios esperados". Lo que quiero saber es cómo se calculan las sumas de cuadrados.

Sé que podría tomar los valores estimados de un modelo en particular (usando coef()), asumir que son fijos y luego hacer pruebas usando las sumas de cuadrados de los residuos del modelo con y sin los factores de interés. Esto está bien para un modelo que contiene un único factor dentro del sujeto. Sin embargo, al implementar un diseño de parcela dividida, el valor de las sumas de cuadrados que obtengo es equivalente al valor producido por R usando aov()una Error()designación apropiada . Sin embargo, esto no es lo mismo que las sumas de cuadrados producidas por el anova()comando en el objeto modelo, a pesar de que las proporciones F son las mismas.

Por supuesto, esto tiene mucho sentido ya que no hay necesidad de Error()estratos en un modelo mixto. Sin embargo, esto debe significar que las sumas de cuadrados se penalizan de alguna manera en un modelo mixto para proporcionar proporciones F apropiadas. ¿Cómo se logra esto? ¿Y cómo el modelo de alguna manera corrige la suma de cuadrados entre parcelas pero no corrige la suma de cuadrados dentro de la parcela? Evidentemente, esto es algo necesario para un ANOVA de parcela dividida clásico que se logró al designar diferentes valores de error para los diferentes efectos, entonces, ¿cómo lo permite un modelo de efectos mixtos?

Básicamente, quiero poder replicar los resultados del anova()comando aplicado a un objeto de modelo más pequeño para verificar los resultados y mi comprensión, sin embargo, actualmente puedo lograr esto para un diseño normal dentro del tema, pero no para la división. diseño de la trama y parece que no puedo descubrir por qué este es el caso.

Como ejemplo:

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 se puede ver sobre todo, las relaciones F están de acuerdo. Las sumas de cuadrados para la variedad también están de acuerdo. Sin embargo, las sumas de cuadrados para riego no están de acuerdo, sin embargo, parece que la salida de lmer está escalada. Entonces, ¿qué hace realmente el comando anova ()?

Usa la fuente, Luke. Podemos echar un vistazo dentro de la función ANOVA haciendo getAnywhere(anova.Mermod). La primera parte de esa función es para comparar dos modelos diferentes. La anova sobre los efectos fijos viene en el elsebloque grande en la segunda mitad:

 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

objectes la salida de lmer. Comenzamos a calcular la suma de cuadrados en la línea 5: ss <- as.vector ...el código multiplica los parámetros fijos (in beta) por una matriz triangular superior; luego cuadra cada término. Aquí está esa matriz triangular superior para el ejemplo de riego. Cada fila corresponde a uno de los cinco parámetros de efectos fijos (intercepción, 3 grados de libertad para riego, 1 df para variedad).

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

La primera fila le da la suma de cuadrados para la intersección y la última le da el SS para el efecto de variedad dentro de los campos. Las filas 2-4 solo involucran los 3 parámetros para los niveles de riego, por lo que la multiplicación previa le proporciona tres piezas de SS para el riego.

Estas piezas no son interesantes en sí mismas, ya que provienen del contraste de tratamiento predeterminado en R, pero en línea ss <- unlist(lapply(split ....Bates recoge trozos de sumas de cuadrados de acuerdo con el número de niveles y a qué factores se refieren. Hay mucha contabilidad aquí. También obtenemos los grados de libertad (que son 3 para riego). Luego, obtiene los cuadrados medios que aparecen en la impresión de anova. Por último, se divide a todos sus medios cuadrados por la varianza residual intra-grupos, sigma(object)^2.

Entonces, ¿qué está pasando? La filosofía de lmernada tiene que ver con el método de enfoque de momentos utilizado por aov. La idea lmeres maximizar una probabilidad marginal obtenida integrando los efectos aleatorios invisibles. En este caso, el nivel de fertilidad aleatorio de cada campo. El capítulo 2 de Pinheiro y Bates describe la fealdad de este proceso. La RXmatriz utilizada para obtener las sumas de cuadrados es su matriz de la ecuación 2.17, página 70 del texto. Esta matriz se obtiene de un grupo de descomposiciones QR sobre (entre otras cosas) la matriz de diseño de los efectos aleatorios y , donde es la varianza del efecto de campo . Este es ese factor faltante$R_{00}$ σ 2 $\sqrt{\sigma^2/\sigma_f^2}$$\sigma_f^2$ preguntaba, pero no entra en la solución de manera transparente o simple.

Asintóticamente, las estimaciones de los efectos fijos tienen distribución:

lo que significa que los componentes de se distribuyen independientemente. Si , los cuadrados de esos términos (divididos por ) siguen una distribución de chi-cuadrado. Dividir a través de la estimación de (otro chi-cuadrado cuando se divide por ) da una estadística F. No dividimos por el error cuadrático medio para el análisis entre grupos porque eso no tiene nada que ver con lo que está sucediendo aquí. Necesitamos comparar las sumas de cuadrados obtenidas a través de comparándolas con una estimación de . ß = 0 σ 2 σ 2 σ 2 R 00 σ $R_{00}\hat{\beta}$$\beta=0$$\sigma^2$$\sigma^2$$\sigma^2$$R_{00}$$\sigma^2$

Tenga en cuenta que no habría obtenido las mismas estadísticas F si los datos se hubieran desequilibrado. Tampoco habría obtenido las mismas estadísticas F si hubiera utilizado ML en lugar de REML.

La idea detrás aoves que los cuadrados medios esperados para el riego son una función de , y los efectos del riego. El cuadrado medio esperado para los residuos de campo es una función de y . Cuando los efectos de riego son 0, ambas cantidades estiman lo mismo y su relación sigue una distribución F.σ 2 f σ 2 σ 2 $\sigma^2$$\sigma_f^2$$\sigma^2$$\sigma_f^2$

Curiosamente, Bates y Pinheiro recomiendan usar el ANOVA sobre el ajuste de dos modelos y hacer una prueba de razón de probabilidad. Este último tiende a ser anti-conservador.

Por supuesto, si los datos no están equilibrados, ya no queda claro exactamente qué hipótesis está probando el ANOVA. Eliminé la primera observación de los datos de riego y la reinstalé. Aquí está la matriz nuevamente:$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 puede ver, las sumas de cuadrados para los parámetros de riego ahora también contienen parte del varietyefecto.