¿Cómo obtener Esfericidad en R para un diseño anidado dentro del sujeto?

Question

Gracias por leer. Estoy tratando de obtener valores de esfericidad para un diseño puramente dentro de sujetos. No he podido usar ezANOVA, o Anova(). Anova funciona si agrego un factor entre sujetos, pero no he podido obtener esfericidad para un diseño puramente dentro de sujetos. ¿Alguna recomendación?

Ya he leído el boletín de R, el apéndice del capítulo de fox, EZanova y cualquier cosa que haya podido encontrar en línea.

Mi ANOVA original

anova(aov(resp ~ sucrose*citral, random =~1 | subject, data = p12bl, subset = exps==1)) 
anova(aov(resp ~ sucrose*citral, random =~1 | subject/sucrose*citral, data = p12bl, subset = exps==1))

> str(subset(p12bl, exps==1))
'data.frame':   192 obs. of  12 variables:
 $ exps     : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Order    : int  1 1 1 1 1 1 1 1 1 1 ...
 $ threshold: Factor w/ 2 levels " Suprathreshold",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ SET      : Factor w/ 2 levels " A"," B": 1 1 1 1 1 1 1 1 1 1 ...
 $ subject  : Factor w/ 16 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ stim     : chr  "S1C1" "S1C1" "S1C1" "S1C1" ...
 $ resp     : num  6.01 5.63 0 2.57 6.81 ...
 $ id       : int  1 2 3 4 5 6 7 8 9 10 ...
 $ X1       : Factor w/ 1 level "S": 1 1 1 1 1 1 1 1 1 1 ...
 $ sucrose  : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...
 $ X3       : Factor w/ 1 level "C": 1 1 1 1 1 1 1 1 1 1 ...
 $ citral   : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 1 1 1 ...

subset(p12b,exps==1)
   exps Order       threshold SET observ S1C1 S1C2 S1C3 S1C4 S2C1 S2C2 S2C3 S2C4 S3C1 S3C2 S3C3 S3C4 S4C1 S4C2 S4C3 S4C4
1     1     1  Suprathreshold   A      1  6.0  7.1  7.5  8.6 15.0 15.4 15.0 13.1 16.9   13 13.1 16.5   24   16   21   20
2     1     1  Suprathreshold   A      2  5.6  0.8  4.0  5.6  5.6 11.3 12.9 14.5 18.5   15 12.9 14.5   24   26   29   28
3     1     1  Suprathreshold   A      3  0.0  0.0  1.7  0.0  5.0  8.4  8.4  5.0 11.7   20 18.5 16.8   29   37   37   30
4     1     1  Suprathreshold   A      4  2.6  3.3  9.1 16.3  5.4 10.0  9.6 16.8 13.5   12 22.2 23.1   19   20   22   23
5     1     1  Suprathreshold   A      5  6.8  5.3 15.4 14.5 11.5  8.3 14.5 14.2  8.9   17 11.2 15.1   24   23   19   19
6     1     1  Suprathreshold   A      6  2.6  2.8  2.6  5.2 13.4 15.6 13.7 13.0 13.7   15 16.0 18.9   22   24   25   25
7     1     1  Suprathreshold   A      7  1.3  5.8 10.2  9.8 11.9 12.3 17.7 16.7 11.4   19 19.2 21.1   16   19   18   19
8     1     1  Suprathreshold   A      8  2.0  5.6  3.9  2.0  4.9  5.2  7.5  4.9 20.2   21  8.2  9.5   30   26   32   45
9     1     1  Suprathreshold   A      9  9.4 11.3 11.7 12.1 14.7 13.8 12.6 14.9 15.2   15 15.9 13.9   17   18   15   18
10    1     1  Suprathreshold   A     10  4.5 17.8 18.5 21.6  5.8 10.9 17.0 20.2  6.6   10 17.8 18.7   12   12   16   19
11    1     1  Suprathreshold   A     11  9.8 13.0 16.1 18.0 10.5 11.6 15.4 17.3 10.1   14 15.2 16.7   13   15   15   17
12    1     1  Suprathreshold   A     12  9.6 10.4 13.3 11.3 12.1 12.6 13.6 13.6 14.9   16 15.1 16.3   16   18   18   17

Salida de muestra

ezANOVA( data = subset(p12bl, exps==1)  , dv= .(resp), sid = .(observ), within = .(sucrose,citral), between = NULL, collapse_within = FALSE)
Nota: el modelo tiene solo una intersección; se han sustituido pruebas equivalentes tipo III.
$ANOVA
          Efecto DFn DFd  SSn  SSd     F       p p<.05   pes
1        sucrose   3  33 4953 3263 16.70 9.0e-07     * 0.603
2         citral   3  33  410  553  8.16 3.3e-04     * 0.426
3 sucrose:citral   9  99   56  791  0.77 6.4e-01       0.066

Mensajes de advertencia: 1: Has eliminado uno o más sujetos del análisis. Refactorización de "observ" para ANOVA. 2: Demasiado pocos sujetos para Anova(), revirtiendo a aov(). Consulte la sección "Advertencia" de la ayuda en ezANOVA.

Answer 1

Prueba:

library(ez)
ezANOVA(data=subset(p12bl, exps==1),
  within=.(sucrose, citral),
  wid=.(subject),
  dv=.(resp)
  )

y el comando aov equivalente, sin esfericidad, es:

aov(resp ~ sucrose*citral + Error(subject/(sucrose*citral)), 
    data= subset(p12bl, exps==1))

Aquí está el equivalente usando Anova de car directamente:

library(car)
df1<-read.table("clipboard", header=T) #De copiar los datos en la pregunta anterior
sucrose<-factor(rep(c(1:4), each=4))
citral<-factor(rep(c(1:4), 4))
idata<-data.frame(sucrose,citral)

mod<-lm(cbind(S1C1, S1C2, S1C3, S1C4, S2C1, S2C2, S2C3, S2C4, 
        S3C1, S3C2, S3C3, S3C4, S4C1, S4C2, S4C3, S4C4)~1, data=df1)
av.mod<-Anova(mod, idata=idata, idesign=~sucrose*citral)
summary(av.mod)

Answer 2

¿Has probado el paquete car, de John Fox? Incluye la función Anova() que es muy útil cuando se trabaja con diseños experimentales. Debería darte el valor p corregido siguiendo la corrección de Greenhouse-Geisser y la corrección de Huynh-Feldt. Puedo publicar un ejemplo rápido en R si te preguntas cómo usarlo.

Además, hay un tutorial interesante sobre el uso de R con mediciones repetidas y modelos de efectos mixtos para experimentos y cuestionarios de psicología; ver la Sección 6.10 sobre esfericidad.

Como nota al margen, el Test de Esfericidad de Mauchly está disponible en mauchly.test(), pero no funciona con el objeto aov si mal no recuerdo. El Boletín de R de octubre de 2007 incluye una breve descripción de este tema.

Answer 3

Generalmente recomiendo evitar por completo este tipo de pruebas de esfericidad utilizando métodos modernos de modelado mixto. Si no estás trabajando con pocos sujetos, esto te dará una gran flexibilidad en el modelado de una estructura de covarianza apropiada, liberándote de la estricta suposición de esfericidad cuando sea necesario. Infiero a partir de la salida str que tienes 16 sujetos con 12 observaciones cada uno (asumo balance porque estás utilizando herramientas clásicas de método de momentos), lo cual debería ser suficiente para ajustar un modelo mixto con matrices de covarianza estructuradas a través de la (máxima) verosimilitud (restringida).

Sin estar cerca de tus datos no puedo ofrecer recomendaciones de modelos específicas, pero un buen lugar para comenzar en R sería reemplazar aov en tus especificaciones de modelo con lme (después de library(nlme)). La razón por la que funcionará es que has proporcionado erróneamente un argumento aleatorio de estilo nlme a aov (cuando, como señaló @Matt Albrecht, un término de Error hubiera sido apropiado). En nlme, con el argumento aleatorio establecido en ~ 1| y sin argumentos de correlation o weight, estás especificando una intercepción aleatoria para cada grupo, implicando que la covarianza de respuesta dentro de los grupos es ZGZ' + R = 1G1'+ \sigma^2 I ==> simetría compuesta con varianza entre grupos fuera de la diagonal y varianza entre grupos + varianza dentro de los grupos en la diagonal ==> una estructura esférica. A partir de ahí puedes comenzar a explorar (por ejemplo, utilizando los métodos gráficos incorporados), modelar y probar (por ejemplo, comparar criterios de información o usar LRTs para modelos anidados) las diversas formas de no esfericidad. Algunas de las herramientas para la componente de modelado son:

• Utilizando el argumento weights para modelar la varianza no constante (diagonales) dentro o entre grupos (por ejemplo, los cambios en la varianza del error entre niveles de sacarosa).
• Utilizando el argumento correlation para modelar la covarianza no constante (fuera de la diagonal) dentro de los grupos (por ejemplo, una estructura dentro de un grupo donde los errores residuales que están más cerca en el tiempo (por ejemplo, estructura AR1) o en el espacio (por ejemplo, estructura esférica) son más similares).
• Modelando pendientes aleatorias agregando términos al LHS de la | en la fórmula de random.

Aunque el proceso puede ser complejo con muchas posibles trampas, creo que te llevará a pensar más en el mecanismo generador de los datos, y cuando se combina con comprobaciones gráficas cuidadosas (recomiendo lattice porque nlme tiene excelentes métodos de trazado basados en lattice, pero ggplot también funciona bien), es probable que no solo tengas una mejor comprensión científica del proceso, sino también estimadores menos sesgados y más eficientes para hacer inferencias.

Answer 4

Ez ahora ha sido actualizado a la versión 2.0. Entre otras mejoras, se ha corregido el error que causaba que no funcionara para este ejemplo.