ANOVA de medidas repetidas con lme/lmer en R para dos factores intra-sujetos

Question

Estoy intentando utilizar lme del nlme para reproducir los resultados de aov para ANOVAs de medidas repetidas. He hecho esto para un experimento de medidas repetidas de un solo factor y para un experimento de dos factores con un factor entre sujetos y un factor dentro de los sujetos, pero estoy teniendo problemas para hacerlo para un experimento de dos factores con dos factores dentro de los sujetos.

A continuación se muestra un ejemplo. A y B son factores de efecto fijo y subject es un factor de efecto aleatorio.

set.seed(1)
d <- data.frame(
    Y = rnorm(48),
    subject = factor(rep(1:12, 4)),
    A = factor(rep(1:2, each=24)),
    B = factor(rep(rep(1:2, each=12), 2)))

summary(aov(Y ~ A*B + Error(subject/(A*B)), data=d))  # Standard repeated measures ANOVA

library(nlme)
# Attempts:
anova(lme(Y ~ A*B, data=d, random = ~ 1 | subject))  # not same as above
anova(lme(Y ~ A*B, data=d, random = ~ 1 | subject/(A+B)))  # gives error

No pude ver una explicación de esto en el libro de Pinheiro y Bates, pero puede que lo haya pasado por alto.

Answer 1

Con qué encaja aov se llama strip plot, y es difícil de encajar con lme porque el subject:A y subject:B se cruzan los efectos aleatorios.

Su primer intento equivale a aov(Y ~ A*B + Error(subject), data=d) que no incluye todos los efectos aleatorios; tu segundo intento es la idea correcta, pero la sintaxis para efectos aleatorios cruzados usando lme es muy complicada.

Utilizando lme del nlme el código sería

lme(Y ~ A*B, random=list(subject=pdBlocked(list(~1, pdIdent(~A-1), pdIdent(~B-1)))), data=d)

Utilizando lmer del lme4 el código sería algo como

lmer(Y ~ A*B + (1|subject) + (1|A:subject) + (1|B:subject), data=d)    

Estos hilos de R-ayuda puede ser útil (y para dar crédito, que es de donde saqué el nlme código de).

http://www.biostat.wustl.edu/archives/html/s-news/2005-01/msg00091.html

http://permalink.gmane.org/gmane.comp.lang.r.lme4.devel/3328

http://www.mail-archive.com/r-help@stat.math.ethz.ch/msg10843.html

Este último enlace remite a la página 165 de Pinheiro/Bates; también puede ser útil.

EDIT: También tenga en cuenta que en el conjunto de datos que tiene, algunos de los componentes de la varianza son negativos, que no se permite el uso de efectos aleatorios con lme, por lo que los resultados difieren. Se puede crear un conjunto de datos con todos los componentes de varianza positivos utilizando una semilla de 8. Los resultados coinciden. Véase esta respuesta para más detalles.

Tenga en cuenta también que lme de nlme no calcula correctamente los grados de libertad del denominador, por lo que los estadísticos F coinciden pero no los valores p, y lmer de lme4 no lo intenta porque es muy complicado en presencia de efectos aleatorios cruzados desequilibrados, y puede que ni siquiera sea sensato hacerlo. Pero eso es más de lo que quiero tratar aquí.

Answer 2

Su primer intento es la respuesta correcta si eso es todo lo que está intentando hacer. nlme() calcula los componentes between y within, no necesita especificarlos.

El problema que te encuentras no es porque no sepas especificar el modelo, es porque ANOVA de medidas repetidas y efectos mixtos no son lo mismo. A veces los resultados del ANOVA y del modelo de efectos mixtos coinciden. Este es especialmente el caso cuando usted agrega sus datos como lo haría para un ANOVA y calcula ambos a partir de eso. Pero en general, cuando se hace correctamente, aunque las conclusiones pueden ser similares, los resultados casi nunca son los mismos. Los datos de su ejemplo no son como las medidas repetidas reales, en las que a menudo se tienen réplicas de cada medida dentro de S. Cuando se hace un ANOVA, normalmente se agregan esas réplicas para obtener una estimación del efecto para cada sujeto. En la modelización de efectos mixtos no se hace tal cosa. Se trabaja con los datos brutos. Cuando lo haga, verá que los resultados nunca son los mismos entre ANOVA y lme().

[Como nota al margen, el uso de lmer() (del paquete lme4) en lugar de lme() me da valores de SS y MS que coinciden exactamente con el ANOVA para los efectos de tu ejemplo, sólo que las F son diferentes].