Los datos
Supongamos que tenemos un conjunto de datos d con dos factores entre-sujetos (es decir, grupos), group y condition, y dos en el de-factores de tema (es decir, de medidas repetidas de los factores), topic y problem (he subido los datos a pastebin, por lo que todo el mundo debería ser capaz de obtener de ella):
> d <- read.table(url("http://pastebin.com/raw.php?i=4hRFyaRj"), colClasses = c(rep("factor", 6), "numeric"))
> str(d)
'data.frame': 2928 obs. of 6 variables:
$ code : Factor w/ 183 levels "A03U","A08C",..: 1 1 1 1 1 1 1 1 1 1 ...
$ group : Factor w/ 2 levels "control","experimental": 2 2 2 2 2 2 2 2 2 2 ...
$ condition: Factor w/ 3 levels "alternatives",..: 3 3 3 3 3 3 3 3 3 3 ...
$ topic : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 2 2 2 2 3 3 ...
$ problem : Factor w/ 4 levels "AC","DA","MP",..: 3 4 1 2 3 4 1 2 3 4 ...
$ mean : num 94.5 94.5 86.5 84.5 80 46.5 73.5 43.5 51 39 ...
Los datos de un experimento conductual en el que los participantes en seis grupos (2 niveles de group veces 3 niveles de condition) ha trabajado en 16 de tareas (para cada uno de los 4 topics 4 diferentes problems). La asignación de los participantes a grupo/condición fue completamente al azar. Presentación de tareas fue al azar, en la medida en que el problema fue bloqueado dentro del tema (es decir, para cada tema todos los problemas donde se presentan de forma secuencial), pero el fin del problema y el tema fue al azar.
Actualización: El factor de identificación del participante (en el que el tema y el problema son anidados) code.
El Problema
Cómo puedo colocar este conjunto de datos utilizando la lme?
(Nota al margen: yo también considerar el uso de lme4, pero yo soy una especie de miedo a no tener los valores de p, si hay algo fácil de digerir como los valores de p, yo también consideraría lme4 una opción).
Hasta ahora me las arreglé para adaptarse a un lme modelo con sólo uno de sujeto factor, pero no dos (ver más abajo).
Lo he intentado
Me puede caber un lme modelo si tengo sólo uno de sujeto factor:
require(nlme)
m1 <- lme(mean ~ condition*group*problem, random = ~1|code/problem,
data = d, subset = topic == "1")
anova(m1)
numDF denDF F-value p-value
(Intercept) 1 531 12101 <.0001
condition 2 177 31 <.0001
group 1 177 2 0.2178
problem 3 531 35 <.0001
condition:group 2 177 1 0.3672
condition:problem 6 531 24 <.0001
group:problem 3 531 1 0.2180
condition:group:problem 6 531 2 0.0281
Este (especialmente el df) muy bien se corresponden con los resultados de un estándar de ANOVA (utilizando
ez):
require(ez)
ezANOVA(subset(d, topic == "1"), dv = .(mean), wid = .(code), between = .(condition, group), within = .(problem))$ANOVA
Warning: Data is unbalanced (unequal N per group). Make sure you specified a well-considered value for the type argument to ezANOVA().
Effect DFn DFd F p p<.05 ges
2 condition 2 177 30.69 0.000000000003611248905859672 * 0.13079
3 group 1 177 1.53 0.217821969825403999321267179 0.00374
5 problem 3 531 34.85 0.000000000000000000014254103 * 0.10028
4 condition:group 2 177 1.01 0.367225806638525886782531416 0.00492
6 condition:problem 6 531 24.40 0.000000000000000000000000142 * 0.13503
7 group:problem 3 531 1.48 0.217959293081550348203379031 0.00472
8 condition:group:problem 6 531 2.38 0.028119961573665430004664856 * 0.01499
Tratando de adaptarse a este tipo de datos con dos de sujeto factores en lme falla (ya sea por código o por dfs):
m2 <- lme(mean ~ condition*group*problem*topic, random = ~1|code/(problem*topic), data = d)
# fails: Error in getGroups.data.frame(dataMix, groups) :
# Invalid formula for groups
m3 <- lme(mean ~ condition*group*problem*topic, random = ~1|code/problem/topic, data = d)
# the next model takes some time (probably already an indicator, that it is the wrong model)
# and produces wrong denominator df!
# with both factors as ANOVA
m4 <- ezANOVA(d, dv = .(mean), wid = .(code), between = .(condition, group), within = .(problem, topic))
#effects are the same
all(row.names(anova(m3))[-1] == m4$ANOVA$Effect)
#denominator dfs are not:
anova(m3)$denDF[-1] == m4$ANOVA$DFd
# only for effects with topic:
row.names(anova(m3))[-1][!(anova(m3)$denDF[-1] == m4$ANOVA$DFd)]
ACTUALIZACIÓN: Como la precisión o error de anidación es poco claro yo aquí proporcionar el equivalente aov de llamadas (este es el "estándar" de la modelo a través de la aov), lo cual coincide con los resultados de ezANOVA. La crítica término de error es Error(code/(problem*topic)):
m5 <- aov(mean ~ (condition*group*problem*topic) + Error(code/(problem*topic)), d)
summary(m5)
