ANOVA de parcela dividida: pruebas de comparación de modelos en R

¿Cómo puedo probar los efectos en un ANOVA de parcela dividida usando comparaciones de modelos adecuados para usar con los argumentos Xy Mde anova.mlm()en R? Estoy familiarizado con ?anova.mlmDalgaard (2007) [1]. Desafortunadamente, solo cepilla los diseños de parcelas divididas. Hacer esto en un diseño completamente al azar con dos factores dentro de los sujetos:

N  <- 20  # 20 subjects total
P  <- 3   # levels within-factor 1
Q  <- 3   # levels within-factor 2
DV <- matrix(rnorm(N* P*Q), ncol=P*Q)           # random data in wide format
id <- expand.grid(IVw1=gl(P, 1), IVw2=gl(Q, 1)) # intra-subjects layout of data matrix

library(car)        # for Anova()
fitA <- lm(DV ~ 1)  # between-subjects design: here no between factor
resA <- Anova(fitA, idata=id, idesign=~IVw1*IVw2)
summary(resA, multivariate=FALSE, univariate=TRUE)  # all tests ...

Las siguientes comparaciones de modelos conducen a los mismos resultados. El modelo restringido no incluye el efecto en cuestión, pero todos los demás efectos del mismo orden o menos, el modelo completo agrega el efecto en cuestión.

anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1
anova(fitA, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2
anova(fitA, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,
                      X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2

Un diseño Split-Splot con un factor interno y otro entre sujetos:

idB  <- subset(id, IVw2==1, select="IVw1")          # use only first within factor
IVb  <- gl(2, 10, labels=c("A", "B"))               # between-subjects factor
fitB <- lm(DV[ , 1:P] ~ IVb)                        # between-subjects design
resB <- Anova(fitB, idata=idB, idesign=~IVw1)
summary(resB, multivariate=FALSE, univariate=TRUE)  # all tests ...

Estos son los anova()comandos para replicar las pruebas, pero no sé por qué funcionan. ¿Por qué las pruebas de las siguientes comparaciones de modelos conducen a los mismos resultados?

anova(fitB, idata=idB, X=~1, test="Spherical") # IVw1, IVw1:IVb
anova(fitB, idata=idB, M=~1, test="Spherical") # IVb

Dos factores dentro de los sujetos y un factor entre sujetos:

fitC <- lm(DV ~ IVb)  # between-subjects design
resC <- Anova(fitC, idata=id, idesign=~IVw1*IVw2)
summary(resC, multivariate=FALSE, univariate=TRUE)  # all tests ...

¿Cómo puedo replicar los resultados dados anteriormente con las comparaciones de modelos correspondientes para usar con los argumentos Xy ? ¿Cuál es la lógica detrás de estas comparaciones modelo?Manova.mlm()

EDITAR: suncoolsu señaló que, a todos los efectos prácticos, los datos de estos diseños deben analizarse utilizando modelos mixtos. Sin embargo, todavía me gustaría entender cómo replicar los resultados de summary(Anova())with anova.mlm(..., X=?, M=?).

[1]: Dalgaard, P. 2007. Nuevas funciones para el análisis multivariante. R News, 7 (2), 2-7.

El Xy Mbásicamente especificar los dos modelos que desea comparar, pero sólo en cuanto a los efectos en un mismo sujeto; luego muestra resultados para la interacción de todos los efectos entre sujetos (incluida la intercepción) con los efectos dentro del sujeto que cambiaron entre Xy M.

Sus ejemplos fitBson más fáciles de entender si agregamos los valores predeterminados para Xy M:

anova(fitB, idata=idB, M=~1, X=~0, test="Spherical") # IVb
anova(fitB, idata=idB, M=diag(3), X=~1, test="Spherical") # IVw1, IVw1:IVb

El primer modelo es el cambio de los efectos no dentro del sujeto (todos tienen la misma media) a una media diferente para cada uno, por lo que hemos agregado el idefecto aleatorio, que es lo correcto para probar la intercepción general y el efecto global entre sujetos en.

El segundo modelo anuncia la id:IVw1interacción, que es lo correcto para probar IVw1y los IVw1:IVbtérminos en contra. Como solo hay un efecto dentro del sujeto (con tres niveles), el valor predeterminado del diag(3)segundo modelo lo explicará; sería equivalente a correr

anova(fitB, idata=idB, M=~IVw1, X=~1, test="Spherical") # IVw1, IVw1:IVb

Para usted fitC, creo que estos comandos recrearán el Anovaresumen.

anova(fitC, idata=id, M=~1, X=~0, test="Spherical") #IVb
anova(fitC, idata=id, M=~IVw1 + IVw2, X=~IVw2, test="Spherical") # IVw1
anova(fitC, idata=id, M=~IVw1 + IVw2, X=~IVw1, test="Spherical") # IVw2
anova(fitC, idata=id, M=~IVw1 + IVw2 + IVw1:IVw2,
                  X=~IVw1 + IVw2, test="Spherical")          # IVw1:IVw2

Ahora, como descubriste, estos comandos son realmente complicados. Afortunadamente, ya no hay muchas razones para usarlos. Si está dispuesto a asumir la esfericidad, solo debe usar aov, o para una sintaxis aún más fácil, simplemente use lmy calcule las pruebas F correctas usted mismo. Si no está dispuesto a asumir la esfericidad, el uso lmees realmente el camino a seguir, ya que obtiene mucha más flexibilidad que con las correcciones GG y HF.

Por ejemplo, aquí está el código aovy lmpara su fitA. Primero debe tener los datos en formato largo; Aquí hay una forma de hacerlo:

library(reshape)
d0 <- data.frame(id=1:nrow(DV), DV)
d0$IVb <- IVb
d0 <- melt(d0, id.vars=c(1,11), measure.vars=2:10)
id0 <- id
id0$variable <- factor(levels(d0$variable), levels=levels(d0$variable))
d <- merge(d0, id0)
d$id <- factor(d$id)

Y aquí está el lm andcódigo aov`:

anova(lm(value ~ IVw1*IVw2*id, data=d))
summary(aov(value ~ IVw1*IVw2 + Error(id/(IVw1*IVw2)), data=d))

Los diseños de parcelas divididas se originaron en la agricultura, de ahí el nombre. Pero con frecuencia ocurren y yo diría: el caballo de batalla de la mayoría de los ensayos clínicos. La trama principal se trata con un nivel de un factor, mientras que los niveles de algún otro factor pueden variar con las subtramas. El diseño surge como resultado de restricciones en una aleatorización completa. Por ejemplo: un campo puede dividirse en cuatro subtramas. Puede ser posible plantar diferentes variedades en subparcelas, pero solo se puede usar un tipo de riego para todo el campo. No es la distinción entre divisiones y bloques.. Los bloques son características de las unidades experimentales que tenemos la opción de aprovechar en el diseño experimental, porque sabemos que están allí. Las divisiones, por otro lado, imponen restricciones sobre qué asignaciones de factores son posibles. Imponen requisitos en el diseño que impiden una aleatorización completa.

Se usan mucho en ensayos clínicos en los que un factor es fácil de cambiar y otro tarda mucho más en cambiar. Si el experimentador debe hacer todas las corridas para cada nivel del factor difícil de cambiar consecutivamente, se obtiene un diseño de diagrama dividido con el factor difícil de cambiar que representa el factor de diagrama completo.

Aquí hay un ejemplo: en una prueba de campo agrícola, el objetivo era determinar los efectos de dos variedades de cultivos y cuatro métodos de riego diferentes. Ocho campos estaban disponibles, pero solo se puede aplicar un tipo de riego a cada campo. Los campos se pueden dividir en dos partes con una variedad diferente en cada parte. El factor de parcela completo es el riego, que debe asignarse aleatoriamente a los campos. Dentro de cada campo, se asigna la variedad.

Así es como haces esto en R:

install.packages("faraway")
data(irrigation)
summary(irrigation)

library(lme4)

R> (lmer(yield ~ irrigation * variety + (1|field), data = irrigation))
Linear mixed model fit by REML 
Formula: yield ~ irrigation * variety + (1 | field) 
   Data: irrigation 
  AIC  BIC logLik deviance REMLdev
 65.4 73.1  -22.7     68.6    45.4
Random effects:
 Groups   Name        Variance Std.Dev.
 field    (Intercept) 16.20    4.02    
 Residual              2.11    1.45    
Number of obs: 16, groups: field, 8

Fixed effects:
                       Estimate Std. Error t value
(Intercept)               38.50       3.02   12.73
irrigationi2               1.20       4.28    0.28
irrigationi3               0.70       4.28    0.16
irrigationi4               3.50       4.28    0.82
varietyv2                  0.60       1.45    0.41
irrigationi2:varietyv2    -0.40       2.05   -0.19
irrigationi3:varietyv2    -0.20       2.05   -0.10
irrigationi4:varietyv2     1.20       2.05    0.58

Correlation of Fixed Effects:
            (Intr) irrgt2 irrgt3 irrgt4 vrtyv2 irr2:2 irr3:2
irrigation2 -0.707                                          
irrigation3 -0.707  0.500                                   
irrigation4 -0.707  0.500  0.500                            
varietyv2   -0.240  0.170  0.170  0.170                     
irrgtn2:vr2  0.170 -0.240 -0.120 -0.120 -0.707              
irrgtn3:vr2  0.170 -0.120 -0.240 -0.120 -0.707  0.500       
irrgtn4:vr2  0.170 -0.120 -0.120 -0.240 -0.707  0.500  0.500

Básicamente, lo que dice este modelo es que el riego y la variedad son efectos fijos y la variedad está anidada dentro del riego. Los campos son los efectos aleatorios y gráficamente serán algo así como

I_1 | I_2 | I_3 | I_4

V_1 V_2 | V_1 V_2 | V_1 V_2 | V_1 V_2

Pero esta era una variante especial con un efecto de trama completo fijo y un efecto de subtrama. Puede haber variantes en las que una o más son aleatorias. Puede haber diseños más complicados como split-split ... diseños de trama. Básicamente, puedes volverte loco y loco. Pero dada la estructura subyacente y la distribución (es decir, fija o aleatoria, anidada o cruzada, ...) se entiende claramente, lmer-Ninjano tendrá problemas en el modelado. Puede ser interpretación será un desastre.

Con respecto a las comparaciones, diga que tiene lmer1y lmer2:

anova(lmer1, lmer2)

le dará la prueba apropiada basada en la estadística de prueba chi-sq con grados de libertad iguales a la diferencia de parámetros.

cf: Faraway, J., Extendiendo modelos lineales con R.

Casella, G., diseño estadístico