En un diseño de medidas repetidas (RM-ANOVA) se evalúa el efecto de un factor sobre una variable independiente medida en dos o más ocasiones para un mismo individuo o para varios individuos que pertenecen a un mismo grupo (e.g., miembros de una misma familia). Como los múltiples valores en la VD están correlacionados entre ellos ya no es posible usar ANOVA en su forma tradicional, donde los errores de estimación son independientes entre los niveles del factor.
Diseño de estudios intra sujetos
El uso de diseños intra-sujetos tiene algunas desventajas y algunas ventajas en comparación con el diseño entre-sujetos. Algunas de sus desventajas se pueden controlar mediante un diseño apropiado.
Desventajas +Efectos de orden: ocurre cuando el orden de los tratamientos es el mismo para todos los sujetos. Es imposible discernir si las diferencias observadas se deben a los distintos tratamientos o al “paso del tiempo”1 No olvidar que el tiempo no causa nada. Por ejemplo, el proceso de oxidación no se debe al paso del tiempo sino a ua reacción química que ocurre concurrentemente. Para controlar este efecto la mitad de los sujetos reciben los tratamientos en el orden A-B, y la otra mitad en el orden B-A.2 Esto se conoce como cross-over design porque todos los sujetos cruzan al otro tratamiento Cuando hay más de dos tratamientos aplicados a un mismo grupo de sujetos hay que implementar un diseño conocido como diseño de cuadrados latinos (del inglés latin squared designs). +Efecto de acarreo diferencial (del inglés carry-over effect): el efecto del tratamiento A tiene un efecto que persiste cuando se implementa el tratamiento B, que no es equivalente al efecto que tiene B sobre A cuando B se implementa primero. Si el efecto de A perdura en el tiempo, la única solución sería usar un diseño entre-sujetos.
Ventajas Requieren menos participantes Tienen más poder incluso con muestras del mismo tamaño que en un diseño entre-sujetos ya que controlan por el efecto de diferencias individuales (el error es más pequeño).
Análisis
El siguiente ejemplo corresponde a un grupo de 12 niños observados a los 30, 36, 42 y 48 meses. La VD en este caso es la inteligencia del grupo de 12 niños y el propósito del análisis es determinar si la inteligencia va cambiando con el tiempo.
Months30
Months36
Months42
Months48
108
96
110
122
103
117
127
133
96
107
106
107
84
85
92
99
118
125
125
116
110
107
96
91
129
128
123
128
90
84
101
113
84
104
100
88
96
100
103
105
105
114
105
112
113
117
132
130
Para RM-ANOVA cada sujeto se considera un factor. Si este fuese un ANOVA factorial le llamaríamos a este diseño un ANOVA de 4x12. La diferencia es que el factor “sujeto” se considera una factor aleatorio en vez de un factor fijo3 Un factor fijo es replicable ya que un investigador puede rehacer el estudio usando los mismos niveles del factor. Los sujetos no son replicables, en el sentido que son una muestra de los niveles que un factor puede asumir. Análogamente a lo que vimos en ANOVA factorial, la pregunta de investigación en RM-ANOVA es si existen diferencias entre las medias del factor A (i.e., meses) controlando por las diferencias asociadas a cada caso. Dicho de otra manera, cuánto peor es la predicción de \(Y_{ij}\) en base a la gran media (\(\mu\)) y el promedio de cada persona (\(\pi\)) en comparación con el mismo modelo pero más el efecto del factor investigado (\(\alpha\)).
El valor predicho en el modelo saturado sería: \[\hat{Y}_{ij}=\hat{\mu}+\hat{\alpha}_j+\hat{\pi}_i\]
Mientras que el valor predicho en el modelo restringido sería: \[\hat{Y}_{ij}=\hat{\mu}+\hat{\pi}_i\]
Si graficamos los valores esperados para el modelo saturado obtenemos el gráfico de más abajo. Cada punto en el gráfico representa el valor predicho para el sujeto i en el tiempo j, que es el factor de interés.
En comparación, el modelo restringido muestra que el puntaje esperado para cada sujeto i es el promedio del sujeto i ignorando el factor tiempo.
La diferencia entre los valores observados y predichos de cada modelo se usa para obtener el error del modelo saturado y restringido respectivamente. Específicamente, la suma de cuadrados del modelo saturado corresponde a la suma de cuadrados de la interacción \(A\times S\), que se obtiene mediante:
\[E_F=\sum_{i=1}^n\sum_{j=1}^a(Y_{ij}-\overline{Y}_{.j}+\overline{Y}_{i.}+\overline{Y}_{..})=SS_{A\times S} \] Mientras que la suma de cuadrados del modelo restringido es matemáticamente equivalente a error asociado a la interacción más el efecto principal y se obtiene mediante:
\[E_R=\sum_{i=1}^n\sum_{j=1}^a(Y_{ij}-\overline{Y}_{i.})=SS_{A\times S}+SS_A \] Por lo que \(E_R-E_F=(SS_{A\times S}+SS_A)-SS_{A\times S}=SS_A\). Dividiendo por los grados de libertad correspondiente, en definitiva, el valor para F se obtiene mediante:
\[F=\frac{SS_A/(a-1)}{SS_{A\times S}/(n-1)(a-1)} \] Si calculamos \(E_R\) y \(E_F\) obtendremos 2558 y 2006, respectivamente, por lo que el valor de F sería 3.03, que implica un valor p de .042: \[F=\frac{552/3}{2006/33}=\frac{184}{60.7879}=3.03\] En R el cálculo de F se obtiene mediante la función aov pero especificando que el Error es la interacción entre sujeto y condición como se muestra más abajo. Como se puede apreciar, los valores resultantes reproducen correctamente lo que obtuvimos en el paso anterior.En este link] de UCLA hay algunas estrategias adicionales para modelar cambios longitudinales.
##
## Error: ID
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 11 6624 602.2
##
## Error: ID:Month
## Df Sum Sq Mean Sq F value Pr(>F)
## Month 3 552 184.00 3.027 0.0432 *
## Residuals 33 2006 60.79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
El valor de F de 3.03 es estadísticamente significativo. Sin embargo, para interpretar correctamente F hay que chequear un supuesto de RM.ANOVA que se llama esfericidad. El valor para F obtenido en el paso anterior se llama F no ajustado, para subrayar que se asume esfericidad. Si no se cumple el supuesto de esfericidad hay que usar un valor de F ajustado (más adelante).
Supuesto de esfericidad
Un supuesto de RM-ANOVA es que la varianza de las diferencias entre cada par de niveles debe ser la misma. Si no se cumple este supuesto la probabilidad de cometer un error de Tipo I puede ser de .10 o incluso .15 cuando \(\alpha\) nominal es .05. En nuestrio ejemplo tenemos 4 niveles y queremos saber si las varianzas de \(Y_{l}-Y_{m}\) son todas iguales.
Para la data d1 la diferencia entre 30 meses y 36 meses es
(dif.30m36m<-d1[,1]-d1[,2])
## [1] 12 -14 -11 -1 -7 3 1 6 -20 -4 -9 -4
Y su varianza es 79.82
var(dif.30m36m)
## [1] 79.81818
Mientras que la varianza de la diferencia entre 36 meses y 42 meses es 91.27
(dif.36m42m<-d1[,2]-d1[,3])
## [1] -14 -10 1 -7 0 11 5 -17 4 -3 9 -15
var(dif.36m42m)
## [1] 91.27273
Formalmente, este supuesto se puede evaluar usando el Test de Esfericidad de Mauchly. En R el test de Mauchly requiere que un objeto de clase SSD que se puede obtener mediante el paquete car. Los resultados bajo el título Univariate Type III Repeated-Measures ANOVA Assuming Sphericity replican los obtenidos en el paso anterior. Sin embargo, dado que el test de Mauchly es estadísticamente significativo (ver Mauchly Tests for Sphericity) tenemos que mirar el valor p bajo el título Greenhouse-Geisser and Huynh-Feldt Corrections.
#require matrix form
month.levels<-c(1,2,3,4)
month.factor<-as.factor(month.levels)
month.frame<-data.frame(month.factor)
d1.matrix<-as.matrix(d1)
month.model<-lm(d1.matrix~1)
require(car)
## Loading required package: car
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
##
## recode
#El nombre de la función del paquete car es "Anova" con "A" mayúscula
rm.aov<-Anova(month.model, idata = month.frame, idesign = ~month.factor)
## Note: model has only an intercept; equivalent type-III tests substituted.
summary(rm.aov)
##
## Type III Repeated Measures MANOVA Tests:
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Response transformation matrix:
## (Intercept)
## Months30 1
## Months36 1
## Months42 1
## Months48 1
##
## Sum of squares and products for the hypothesis:
## (Intercept)
## (Intercept) 2239488
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.98831 929.7391 1 11 0.0000000000055863
## Wilks 1 0.01169 929.7391 1 11 0.0000000000055863
## Hotelling-Lawley 1 84.52174 929.7391 1 11 0.0000000000055863
## Roy 1 84.52174 929.7391 1 11 0.0000000000055863
##
## Pillai ***
## Wilks ***
## Hotelling-Lawley ***
## Roy ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: month.factor
##
## Response transformation matrix:
## month.factor1 month.factor2 month.factor3
## Months30 1 0 0
## Months36 0 1 0
## Months42 0 0 1
## Months48 -1 -1 -1
##
## Sum of squares and products for the hypothesis:
## month.factor1 month.factor2 month.factor3
## month.factor1 972 540 216
## month.factor2 540 300 120
## month.factor3 216 120 48
##
## Multivariate Tests: month.factor
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.4274916 2.240098 3 9 0.15284
## Wilks 1 0.5725084 2.240098 3 9 0.15284
## Hotelling-Lawley 1 0.7466993 2.240098 3 9 0.15284
## Roy 1 0.7466993 2.240098 3 9 0.15284
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value Pr(>F)
## (Intercept) 559872 1 6624 11 929.7391 0.000000000005586 ***
## month.factor 552 3 2006 33 3.0269 0.04322 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## month.factor 0.24265 0.017718
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## month.factor 0.60954 0.07479 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## month.factor 0.7248502 0.06353773
Tanto el test de Greenhouse-Geisser (GG) como el test de Huynh-Feldt (HF) se basan en el valor de \(\varepsilon\) que indica en qué medida la matriz de covarianza se aleja del supuesto de esfericidad. A medida que \(\varepsilon\) se aleja de 1 más se aleja la matriz de covaranza del supuesto de esfericidad. Como se observa en la tabla de resultados GG (\(\varepsilon=.61\)) es más conservador que HF (\(\varepsilon=.72\)). No obstante, en ambos casos el valor ajustado de F no es estadísticamente significativo, p > .05. El test de HF es preferible porque produce valores más ajustados a \(\alpha\) nominal (T&F, 2007, p. 287).
Tamaño del efecto como asociación entre variables
El cálculo de \(\omega^2\) sigue la misma lógica que vimos para el caso de ANOVA entre-sujetos usando la siguiente fórmula4 T&F (2007) usan \(\omega^2_p\) (omega cuadrado parcial) que no considera el efecto de los sujetos en el error, por lo que suele arrojar valores demasiado altos.
\[\omega^2=\frac{(a-1)(MS_A-MS_{A\times S})}{SS_{total}+MS_S} \] Usando como referencia los resultados del modelo anterior y reemplazando en la fórmula obtenemos:
Esto es, la edad de los niños explica aproximadamente un 4% de la varianza de sus puntajes de inteligencia.
Tamaño del efecto como diferencia de medias estandarizada
Al igual que vimos en ANOVA entre-sujetos podemos obtener una medida del tamaño del efecto ( d) que sea independiente de la escala de medición utilizada. Para calcular la varianza error que usaremos en el denominador del cálculo de d tenemos que tener en consideración que el error del factor A es ahora la sumatoria de \(SS_S\) y \(SS_{A\times S}\). En ANOVA entre-sujetos el error no distinguía entre el error atribuible a los sujetos y la interacción entre el factor y los sujetos. Luego, la estimación de la desviación estándar está dada por:
\[sd=\sqrt{\frac{SS_S+SS_{A\times S}}{a(n-1)}}\] Por lo que,
\[sd=\sqrt{\frac{6624+2006}{4(12-1)}}=14.00\] En nuestra data tenemos las siguientes medias:
#MARGIN = 1 se obtiene las medias de cada fila
#MARGIN = 2 se obtiene las medias de cada columna
apply(d1,MARGIN = 2, FUN = "mean")
Luego, la diferencia estandarizada ( d) entre la inteligencia al mes 48 y al mes 30 sería (112-103)/14=0.64
Análisis de tendencia
Tal como vimos en relación a ANOVA entre sujetos, es posible usar contrastes para evaluar diferencias entre medias o grupos de medias. Sin embargo, en RM-ANOVA es usualmente interesante evaluar si existe una tendencia. En el ejemplo, tenemos que las medias de los grupos van subiendo progresivamente, desde 103 hasta 112. Imaginemos que este fuese un diseño entre-sujetos (niños distintos en cada grupo)
## Df Sum Sq Mean Sq F value Pr(>F)
## Month 3 552 184.0 0.938 0.43
## Residuals 44 8630 196.1
Podemos hacer un contraste para evaluar la diferencia entre el grupo de niños entre 30 y 36 meses y el grupo de niños entre 42 y 48 meses.5 Recordemos que la suma de los contrastes debe ser 0 y que idealmente los contrastes deben sumar |1|.
Sin embargo, en este caso es más interesante evaluar si hay una progresión lineal entre el mes 30 y el mes 42. En este caso tenemos 4 niveles (1,2,3, y 4). Construimos los contrastes como la diferencia entre cada nivel y la media de los 4 niveles (i.e., 2.5).
\[c_1=(1-2.5)=-1.5\]\[c_2=(2-2.5)=-0.5\]
\[c_3=(3-2.5)=0.5\]\[c_4=(2-2.5)=-0.5\] Luego, el valor predicho para el contraste es:
Sin embargo, el valor p para el contraste es > .05, lo que es consistente con los resultados iniciales de ANOVA. Sin embargo, este es el resultado ignorando el efecto de los sujetos. El suma de cuadrados para un contraste en el contexto de RM-ANOVA se obtiene mediante:
\[F=\frac{SS_\psi}{MS_{A\times S}}\] donde \(SS_\psi\) se obtiene mediante:
Podemos intentar resolver el mismo problema con el paquete emmeans usando los contrastes especificados pero teniendo cuidado que ahora el modelo incorpore el efecto de los sujetos.
En este caso, el valor t resultante es estadísticamente significativo. Si elevamos el valor t al cuadrado obtenemos el valor F calculado en el paso anterior.6 ¡BAM!
\[t^2=F=8.88\] Habiendo dicho todo lo anterior, en el caso de RM-ANOVA casi siempre se viola el supuesto de esfericidad, en cuyo caso los resultados obtenidos a través de este procedimiento no son correctos.
Diseños con más de un factor intrasujeto
Es posible analizar modelos con más de un factor intra-sujetos. En el siguiente ejemplo, se le presentó la letra T o I a 10 participantes. En algunas ocasiones la letra objetivo (T o I) aparece en un grupo con otras letras (noise present) y en otras aparece sola (noise absent). Adicionalmente, la letra objetivo se muestra al centro de la pantalla (0°) o desplazada 4° y 8°. Por lo tanto, cada persona tiene 6 puntajes, uno para cada combinación de los factores. Por lo tanto en este diseño hay dos efectos principales (ruido y grado de rotación, con 2 y 3 niveles respectivamente) más un efecto de interacción (ruido \(\times\) grado de rotación)
kable(d)
Absent0
Absent4
Absent8
Present0
Present4
Present8
420
420
480
480
600
780
420
480
480
360
480
600
480
480
540
660
780
780
420
540
540
480
780
900
540
660
540
480
660
720
360
420
360
360
480
540
480
480
600
540
720
840
480
600
660
540
720
900
540
600
540
480
720
780
480
420
540
540
660
780
Para poder analizar la data vamos a llevarla del formato ancho (wide) a un formato largo (long) donde una columna representa los niveles de la variables independiente intra-sujeto. Además, tendremos cuidado de decirle a R que las variables independientes son factores.
Preparación de la data
require(AMCP)
data("chapter_12_table_1")
# Add an ID variable (by row number)
C12_T1_PP <- chapter_12_table_1
C12_T1_PP$ID <- as.factor(1:nrow(C12_T1_PP))
# Note that the "sep=-2" partitions between the 1st and 2nd value from the end.
# The "-ID" here means to leave ID as it is (i.e., without "gathering" it). We
# use "Noise.Angle" variable before seperating into the group and the angle.
C12_T1_PP = C12_T1_PP %>%
gather(key=Noise.Angle, value=Reaction_Time, -ID) %>%
separate(col=Noise.Angle, into=c("Noise", "Angle"), sep=-1) %>%
arrange(ID, Noise, Angle)
# Display the data
C12_T1_PP
# Examining the structure of the data we see that "Angle" and "Noise" are not factors,
# thus the following code.
C12_T1_PP$Angle <- as.factor(C12_T1_PP$Angle)
C12_T1_PP$Noise <- as.factor(C12_T1_PP$Noise)
str(C12_T1_PP)
##
## Error: ID
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 9 292140 32460
##
## Error: Noise
## Df Sum Sq Mean Sq
## Noise 1 285660 285660
##
## Error: ID:Angle
## Df Sum Sq Mean Sq F value Pr(>F)
## Angle 2 289920 144960 40.72 0.000000209 ***
## Residuals 18 64080 3560
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Error: ID:Noise
## Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 9 76140 8460
##
## Error: ID:Angle:Noise
## Df Sum Sq Mean Sq F value Pr(>F)
## Angle:Noise 2 105120 52560 45.31 0.0000000942 ***
## Residuals 18 20880 1160
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
La complejidad del análisis reside en determinar para cada efecto cuál es el error apropiado. En general, el denominador de F es la interacción de sujetos con el efecto investigado
\[F=\frac{MS_{effect}}{MS_{effect \times S}}\]
Por ejemplo, para el efecto del ruido, tenemos que la media cuadrática para este efecto es 285.660 con 1 df. El denominador de este efecto sería la media cuadrática de ruido \(\times\) sujeto, que es 8.460 con 9 df. Por lo tanto el efecto de este factor es:
\[F=\frac{285660}{8460}=33.77\] El valor crítico para F en este caso es 5.12 por lo que F observado es estadísticamente significativo
qf(.95,1,9)
## [1] 5.117355
Diseños con un factor intra-sujeto y un factor inter-sujeto
En este caso tenemos un factor inter-sujeto además del factor intra-sujeto. Usaremos la data de Tabachnick y Fidell (2007, p.324) para ejemplificar el análisis
d_wide <- read.table(header=TRUE, text="
book month1 month2 month3
1 science 1 3 6
2 science 1 4 8
3 science 3 3 6
4 science 5 5 7
5 science 2 4 5
6 mystery 3 1 0
7 mystery 4 4 2
8 mystery 5 3 2
9 mystery 4 2 0
10 mystery 4 5 3
11 romance 4 2 0
12 romance 2 6 1
13 romance 3 3 3
14 romance 6 2 1
15 romance 3 3 2
")
###################
## Esta sintaxis permite tener una versión de la data
## donde el factor está representado en una columna
## Para más detalles ver
## http://www.cookbook-r.com/Manipulating_data/Converting_data_between_wide_and_long_format/
d_wide$Subject<-as.factor(1:nrow(d_wide))
library(tidyr)
d_long<-gather(d_wide,month,n_novels,month1:month3,factor_key = TRUE)
str(d_long)
El análisis en este caso consiste en agregar los factores inter- e inter-sujetos y su interacción especificando como denominador los factores intrasujeto (en este caso Subject y month)
Lo importante en este caso es entender cuál es denominador para la media cuadrática de cada efecto. Como se ve a continuación, el error depende del efecto que se está evaluando. Si el factor inter-sujeto es A y el factor intrasujeto es B, F se calcula de la siguiente manera:
