¿Cómo especificar contrastes específicos para ANOVA de medidas repetidas usando automóvil?

Estoy tratando de ejecutar medidas repetidas Anova en R seguido de algunos contrastes específicos en ese conjunto de datos. Creo que el enfoque correcto sería usar Anova()el paquete del automóvil.

Vamos a ilustrar mi pregunta con el ejemplo tomado del ?Anovauso de los OBrienKaiserdatos (Nota: omití el factor de género del ejemplo):
tenemos un diseño con un factor entre sujetos, tratamiento (3 niveles: control, A, B) y 2 repetidos -factores de medidas (dentro de los sujetos), fase (3 niveles: pretest, posttest, seguimiento) y hora (5 niveles: 1 a 5).

La tabla ANOVA estándar viene dada por (a diferencia del ejemplo (Anova) cambié a Sumas de cuadrados de tipo 3, eso es lo que mi campo quiere):

require(car)
phase <- factor(rep(c("pretest", "posttest", "followup"), c(5, 5, 5)),
levels=c("pretest", "posttest", "followup"))
hour <- ordered(rep(1:5, 3))
idata <- data.frame(phase, hour)
mod.ok <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser)
av.ok <- Anova(mod.ok, idata=idata, idesign=~phase*hour, type = 3)
summary(av.ok, multivariate=FALSE)

Ahora, imagine que la interacción de orden superior hubiera sido significativa (que no es el caso) y nos gustaría explorarla más a fondo con los siguientes contrastes:
¿Hay alguna diferencia entre las horas 1 y 2 y las horas 3 (contraste 1) y entre las horas 1 y 2 versus horas 4 y 5 (contraste 2) en las condiciones de tratamiento (A y B juntos)?
En otras palabras, ¿cómo especifico estos contrastes?

1. ((treatment %in% c("A", "B")) & (hour %in% 1:2)) versus ((treatment %in% c("A", "B")) & (hour %in% 3))
2. ((treatment %in% c("A", "B")) & (hour %in% 1:2)) versus ((treatment %in% c("A", "B")) & (hour %in% 4:5))

Mi idea sería ejecutar otro ANOVA omitiendo la condición de tratamiento no necesaria (control):

mod2 <- lm(cbind(pre.1, pre.2, pre.3, pre.4, pre.5, post.1, post.2, post.3, post.4, post.5, fup.1, fup.2, fup.3, fup.4, fup.5) ~ treatment, data=OBrienKaiser, subset = treatment != "control")
av2 <- Anova(mod2, idata=idata, idesign=~phase*hour, type = 3)
summary(av2, multivariate=FALSE)

Sin embargo, todavía no tengo idea de cómo configurar la matriz de contraste apropiada dentro del sujeto comparando las horas 1 y 2 con 3 y 1 y 2 con 4 y 5. Y no estoy seguro de si omitir el grupo de tratamiento no necesario es una buena idea, ya que cambia el término de error general.

Antes de ir Anova(), también estaba pensando en ir lme. Sin embargo, existen pequeñas diferencias en los valores de F y p entre el ANOVA del libro de texto y lo que se devuelve anove(lme) debido a posibles variaciones negativas en el ANOVA estándar (que no están permitidoslme ). Relacionado, alguien me señaló glsque permite ajustar medidas repetidas ANOVA, sin embargo, no tiene argumento de contraste.

Para aclarar: quiero una prueba F o t (usando sumas de cuadrados de tipo III) que responda si los contrastes deseados son significativos o no.

Actualizar:

Ya hice una pregunta muy similar sobre R-help, no hubo respuesta .

Una pregunta similar se planteó en R-help hace algún tiempo. Sin embargo, las respuestas tampoco resolvieron el problema.

Actualización (2015):

Como esta pregunta sigue generando algo de actividad, ahora se pueden especificar tesis y, básicamente, todos los demás contrastes con el afexpaquete en combinación con el lsmeanspaquete como se describe en la viñeta afex .

Este método generalmente se considera "anticuado", por lo que, aunque puede ser posible, la sintaxis es difícil y sospecho que menos personas saben cómo manipular los comandos de anova para obtener lo que desea. El método más común es usar glhtcon un modelo basado en la probabilidad de nlmeo lme4. (Sin embargo, estoy bienvenido a que otras respuestas me demuestren que están equivocadas).

Dicho esto, si tuviera que hacer esto, no me molestaría con los comandos de anova; Simplemente ajustaría el modelo equivalente usando lm, elegiría el término de error correcto para este contraste y calcularía la prueba F yo mismo (o equivalentemente, la prueba t ya que solo hay 1 df). Esto requiere que todo esté equilibrado y tenga esfericidad, pero si no lo tiene, probablemente debería estar utilizando un modelo basado en la probabilidad de todos modos. Es posible que pueda corregir algo la no esfericidad utilizando las correcciones Greenhouse-Geiser o Huynh-Feldt que (creo) usan la misma estadística F pero modifican la df del término de error.

Si realmente desea utilizar car, puede encontrar útiles las viñetas heplot ; describen cómo carse definen las matrices en el paquete.

Usando el método de caracal (para los contrastes 1 y 2 - 3 y 1 y 2 - 4 y 5), obtengo

      psiHat      tStat          F         pVal
1 -3.0208333 -7.2204644 52.1351067 2.202677e-09
2 -0.2083333 -0.6098777  0.3719508 5.445988e-01

Así es como obtendría esos mismos valores p:

Cambie la forma de los datos a formato largo y ejecútelos lmpara obtener todos los términos de SS.

library(reshape2)
d <- OBrienKaiser
d$id <- factor(1:nrow(d))
dd <- melt(d, id.vars=c(18,1:2), measure.vars=3:17)
dd$hour <- factor(as.numeric(gsub("[a-z.]*","",dd$variable)))
dd$phase <- factor(gsub("[0-9.]*","", dd$variable), 
                   levels=c("pre","post","fup"))
m <- lm(value ~ treatment*hour*phase + treatment*hour*phase*id, data=dd)
anova(m)

Haga una matriz de contraste alternativa para el término de la hora.

foo <- matrix(0, nrow=nrow(dd), ncol=4)
foo[dd$hour %in% c(1,2) ,1] <- 0.5
foo[dd$hour %in% c(3) ,1] <- -1
foo[dd$hour %in% c(1,2) ,2] <- 0.5
foo[dd$hour %in% c(4,5) ,2] <- -0.5
foo[dd$hour %in% 1 ,3] <- 1
foo[dd$hour %in% 2 ,3] <- 0
foo[dd$hour %in% 4 ,4] <- 1
foo[dd$hour %in% 5 ,4] <- 0

Comprueba que mis contrastes dan el mismo SS que los contrastes predeterminados (y el mismo que el del modelo completo).

anova(lm(value ~ hour, data=dd))
anova(lm(value ~ foo, data=dd))

Obtenga el SS y el df solo por los dos contrastes que quiero.

anova(lm(value ~ foo[,1], data=dd))
anova(lm(value ~ foo[,2], data=dd))

Obtenga los valores p.

> F <- 73.003/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 2.201148e-09
> F <- .5208/(72.81/52)
> pf(F, 1, 52, lower=FALSE)
[1] 0.5445999

Opcionalmente ajuste por esfericidad.

pf(F, 1*.48867, 52*.48867, lower=FALSE)
pf(F, 1*.57413, 52*.57413, lower=FALSE)

Si desea / tiene que usar contrastes con el término de error agrupado del ANOVA correspondiente, puede hacer lo siguiente. Desafortunadamente, esto será largo, y no sé cómo hacerlo más convenientemente. Aún así, creo que los resultados son correctos, ya que se verifican contra Maxwell y Delaney (ver más abajo).

Desea comparar grupos de su primer factor interno houren un diseño SPF-p.qr (notación de Kirk (1995): Diseño factorial de parcela dividida 1 entre factor treatmentcon grupos p, primero dentro de factor hourcon grupos q, segundo dentro de factor prePostFupcon r grupos). Lo siguiente supone treatmentgrupos y esfericidad de tamaño idéntico .

Nj    <- 10                                             # number of subjects per group
P     <- 3                                              # number of treatment groups
Q     <- 5                                              # number of hour groups
R     <- 3                                              # number of PrePostFup groups
id    <- factor(rep(1:(P*Nj), times=Q*R))                                  # subject
treat <- factor(rep(LETTERS[1:P], times=Q*R*Nj), labels=c("CG", "A", "B")) # treatment
hour  <- factor(rep(rep(1:Q, each=P*Nj), times=R))                         # hour
ppf   <- factor(rep(1:R, each=P*Q*Nj), labels=c("pre", "post", "fup"))     # prePostFup
DV    <- round(rnorm(Nj*P*Q*R, 15, 2), 2)               # some data with no effects
dfPQR <- data.frame(id, treat, hour, ppf, DV)           # data frame long format

summary(aov(DV ~ treat*hour*ppf + Error(id/(hour*ppf)), data=dfPQR)) # SPF-p.qr ANOVA

Primero tenga en cuenta que el efecto principal para houres el mismo después de promediar prePostFup, cambiando así al diseño SPF-pq más simple que solo contiene treatmenty hourcomo IV.

dfPQ <- aggregate(DV ~ id + treat + hour, FUN=mean, data=dfPQR)  # average over ppf
# SPF-p.q ANOVA, note effect for hour is the same as before
summary(aov(DV ~ treat*hour + Error(id/hour), data=dfPQ))

Ahora tenga en cuenta que en el ANOVA SPF-pq, hourse prueba el efecto para la interacción id:hour, es decir, esta interacción proporciona el término de error para la prueba. Ahora, los contrastes para los hourgrupos se pueden probar como en un ANOVA de una sola vía entre sujetos simplemente sustituyendo el término de error y los correspondientes grados de libertad. La manera fácil de obtener el SS y el df de esta interacción es ajustar el modelo lm().

(anRes <- anova(lm(DV ~ treat*hour*id, data=dfPQ)))
SSE    <- anRes["hour:id", "Sum Sq"]     # SS interaction hour:id -> will be error SS
dfSSE  <- anRes["hour:id", "Df"]         # corresponding df

Pero también calculemos todo manualmente aquí.

# substitute DV with its difference to cell / person / treatment group means
Mjk   <- ave(dfPQ$DV,           dfPQ$treat, dfPQ$hour, FUN=mean)  # cell means
Mi    <- ave(dfPQ$DV, dfPQ$id,                         FUN=mean)  # person means
Mj    <- ave(dfPQ$DV,           dfPQ$treat,            FUN=mean)  # treatment means
dfPQ$IDxIV <- dfPQ$DV - Mi - Mjk + Mj                             # interaction hour:id
(SSE  <- sum(dfPQ$IDxIV^2))               # SS interaction hour:id -> will be error SS
dfSSE <- (Nj*P - P) * (Q-1)               # corresponding df
(MSE  <- SSE / dfSSE)                     # mean square

$t = \frac{\hat{\psi} - 0}{||c|| \sqrt{MS_{E}}}$| El | c | El | Ψ = q Σ k = 1 c k M . k M S $c$$||c||$$\hat{\psi} = \sum\limits_{k=1}^{q} c_{k} M_{.k}$$MS_{E}$hour:id

Mj     <- tapply(dfPQ$DV, dfPQ$hour, FUN=mean)  # group means for hour
Nj     <- table(dfPQ$hour)                      # cell sizes for hour (here the same)
cntr   <- rbind(c(1, 1, -2,  0, 0),
                c(1, 1, -1, -1, 0))             # matrix of contrast vectors
psiHat <- cntr   %*% Mj                         # estimates psi-hat
lenSq  <- cntr^2 %*% (1/Nj)                     # squared lengths of contrast vectors
tStat  <- psiHat / sqrt(lenSq*MSE)              # t-statistics
pVal   <- 2*(1-pt(abs(tStat), dfSSE))           # p-values
data.frame(psiHat, tStat, pVal)

Para comparaciones múltiples, tendría que pensar en métodos de corrección , por ejemplo, Bonferroni.$\alpha$

Los cálculos correspondientes para el ejemplo de Maxwell y Delaney (2004) en la pág. 599f se puede encontrar aquí . Tenga en cuenta que M&D calcula el valor F, para ver que los resultados son idénticos, debe ajustar al cuadrado el valor de la estadística t. Ese código también incluye el análisis realizado con Anova()from car, así como un cálculo manual de las correcciones para el efecto principal del factor interno.$\hat{\epsilon}$