Cuando queremos analizar el efecto de más de un factor sobre una variable dependiente debemos usar ANOVA factorial. La principal diferencia entre ANOVA unifactorial y ANOVA factorial es que este último evalúa la presencia de un efecto de interacción entre los factores estudiados. Veamos un ejemplo de un ANOVA con dos factores, cada uno con dos niveles (que se designa como ANOVA de 2$$2). El propósito de este estudio es evaluar el efecto del biofeedback (Biofeedback) y una droga (DrugTherapy) para tratar la hipertensión. Un tratamiento efectivo en este caso debiese tener el efecto de bajar la presión arterial, que en este caso es la variable dependiente (Score). La combinación de los niveles de la variable independiente da lugar a cuatro condiciones experimentales: sólo biofeedback, sólo droga, biofeedback+droga (de aquí en adelante B+D), y sin tratamiento. Los resultados de cada condición se presentan más abajo.
d1<-readRDS(file="biofeedback_2x2")
#Asignar nombres a los niveles de la VI
levels(d1$Group) <- c('B+D', 'Biofeedback','Droga','Ninguno')
En este caso se aprecia que la combinación de biofeedback y drogas tiene un resultado más favorable que cualquier tratamiento por separado o la falta de tratamiento.
Hay dos tipos de efectos a analizar en este caso. Por un lado, los efectos principales se refieren al efecto de cada factor ignorando el nivel del segundo factor. Esto se puede obtener mediante un contraste. Para el caso de biofeedback podemos combinar la información del grupo B+D y Biofeedback y compararlo con la media de Droga y Ninguno.1 En general es preferible usar contrastes que sumen 1 y -1. en este caso usaremos contrastes que suman 2 y -2, porque me parece más fácil entender el problema usando estos valores
aov1 = aov(Score ~ as.factor(Group),
data = d1)
summary(aov1)## Df Sum Sq Mean Sq F value Pr(>F)
## as.factor(Group) 3 1540 513.3 8.213 0.00155 **
## Residuals 16 1000 62.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Requiere el paquete emmeans
#install.packages("emmeans")
require(emmeans)## Loading required package: emmeansaov1.emm = emmeans(aov1, specs =~Group)
# Agreguemos un contraste para la diferencia entre la
# condición sin música y las otras dos condiciones
# combinadas
Contrasts.1 = list(Biofeedback_vs_todos = c(1,1,-1,-1))
contrast(aov1.emm, Contrasts.1)## contrast estimate SE df t.ratio p.value
## Biofeedback_vs_todos -20 7.071068 16 -2.828 0.0121Para evaluar el efecto de la droga (ignorando el efecto del biofeedback) podemos combinar la información del grupo B+D y Droga y compararlo con Biofeedback y Ninguno.
Contrasts.2 = list(Droga_vs_todos = c(1,-1,1,-1))
contrast(aov1.emm, Contrasts.2)## contrast estimate SE df t.ratio p.value
## Droga_vs_todos -24 7.071068 16 -3.394 0.0037Evaluar cualquiera de estos efectos principales implica ignorar el efecto del otro nivel. Pero también podemos evaluar el efecto de interacción. Para entender la lógica de este paso podemos hacer un tercer contraste (ortogonal a los anteriores) para evaluar la siguiente hipótesis nula.
\[Ho:\mu_{B+D}-\mu_{D}=\mu_{B}-\mu_{ninguno}\]
Si la diferencia entre B+D y D es igual a la diferencia entre B y ningún tratamiento, esto indicaría que el efecto de B+D no es más que la suma de los efectos de cada tratamiento por separado. En el caso que se rechace Ho el efecto de ambos tratamientos combinados es mayor que el efecto de cualquiera de ellos por separado.
Una hipótesis equivalente se puede plantear de la siguiente manera: \[Ho:\mu_{B+D}-\mu_{B}=\mu_{D}-\mu_{ninguno}\] En general, el contraste para evaluar el efecto interacción se podría especificar como \[Ho:\mu_{B+D}-\mu_{B}-\mu_{D}+\mu_{ninguno}=0\] \[Ho:(1)\mu_{B+D}-(1)\mu_{B}-(1)\mu_{D}+(1)\mu_{ninguno}=0\]
Contrasts.3 = list(interaccion_BxD = c(1,-1,-1,1))
contrast(aov1.emm, Contrasts.3)## contrast estimate SE df t.ratio p.value
## interaccion_BxD -16 7.071068 16 -2.263 0.0379El valor de -16 quiere decir que el efecto de B+D es distinto que el efecto combinado (o la suma) de B y D. Por eso los efectos de interacción se llaman también efectos no-aditivos. Para visualizar el efecto de interacción es conveniente graficar las medias.
# La data tiene dos variables dummy para indexar a qué celda pertenece cada
# caso. Así que vamos a hacer algunas modificaciones para facilitar ver el gráfico
# de interacción. Puedes comparar el antes y después corriendo sólo el código
# para generar el gráfico de interacción primero (al final de este chunk)
d1$DrugTherapy<-as.factor(d1$DrugTherapy)
levels(d1$DrugTherapy)<-c("Ausencia","Presencia")
d1$DrugTherapy <- factor(d1$DrugTherapy, levels = c("Presencia", "Ausencia"))
d1$Biofeedback<-as.factor(d1$Biofeedback)
levels(d1$Biofeedback)<-c("Ausencia","Presencia")
d1$Biofeedback <- factor(d1$Biofeedback, levels = c("Presencia", "Ausencia"))Aquí voy a usar la función summarySE que está disponible en www.cookbook-r.com
Gráfico de Interacción con interaction.plot
d1.means <- summarySE(d1, measurevar="Score", groupvars=c("Biofeedback","DrugTherapy"))
kable(d1.means)| Biofeedback | DrugTherapy | N | Score | sd | se | ci |
|---|---|---|---|---|---|---|
| Presencia | Presencia | 5 | 168 | 7.905694 | 3.535534 | 9.816216 |
| Presencia | Ausencia | 5 | 188 | 7.905694 | 3.535534 | 9.816216 |
| Ausencia | Presencia | 5 | 186 | 7.905694 | 3.535534 | 9.816216 |
| Ausencia | Ausencia | 5 | 190 | 7.905694 | 3.535534 | 9.816216 |
Gráfico de interacción con ggplot2
Gráfico de interacción con ggplot2
Podemos clarificar qué evalúa el efecto de interacción reemplazando por las medias en la Ho (sólo con fines ilustrativos).
\[Ho:(1)\mu_{B+D}-(1)\mu_{B}-(1)\mu_{D}+(1)\mu_{ninguno}=0\] \[Ho:(1)\mu_{B+D}-(1)\mu_{B}=(1)\mu_{D}-(1)\mu_{ninguno}\] \[Ho:(1)(168)-(1)(188)=(1)(186)-(1)(190)\] \[Ho:-20=-4\]
La diferencia entre B+D y B es de -20 puntos, mientras que la diferencia entre D y ningún tratamiento es de sólo -4 puntos. Si el efecto de combinar B+D fuese nada más que la suma de B y D, entonces la diferencia entre B+D y B debiese ser aproximadamente la misma que la diferencia entre D y la ausencia de tratamiento.
Para determinar la significancia estadística de los efectos principales y el efecto de interacción es necesario hacer una comparación de modelos. Para calcular el error del modelo saturado usaremos como referencia la media de cada celda (que en nuestro caso corresponde a la intersección de los dos niveles). El modelo saturado se puede formular como: \[Y_{ijk}=\mu_{jk}+\varepsilon_{ijk}\] donde i corresponde a cada sujeto; j y k representan el nivel al que pertenece en cada factor; y \(\mu\) es la media de cada celda. Por lo tanto, el error de este modelo se puede formular como:
\[E_F=SS_W=\sum(Y_{ijk}-\hat{Y}_{jk}(F))^2\]
donde \(\hat{Y}_{jk}(F)=\overline{Y}_{jk}\), es decir, el valor predicho para cada sujeto es el promedio de su celda definida por una columna (j) y una fila (k). En este caso el error del modelo saturado sería 1000.
\[E_F=SS_W=(158-168)^2+(163-168)^2...(200-190)^2+(180-190)^2=1000\]
El asunto es que ahora tenemos 3 modelos restringidos: uno para cada efecto principal2 Siempre que n es el mismo para cada combinación de los factores los resultados de la prueba F son idénticos al que obtendríamos usualmente más uno adicional para el efecto de interacción.
\[E_R-E_F=SS_A=nb\sum_{j=1}^a(\overline{Y}_{j.}-\overline{Y}_{..})^2\]
donde n es el número de observaciones por celda y b es el número de niveles del factor b. Así, tenemos que el efecto del biofeedback es:
\[SS_A=(5)(2)[(177-183)^2+(189-183)^2]\] \[SS_A=(10)[(-6)^2+(6)^2]=720\] Mientras que el efecto de la droga es:
\[SS_B=(10)[(178-183)^2+(188-183)^2]\] \[SS_B=(10)[(-5)^2+(5)^2]=500\]
Por último, para el efecto de interacción, el valor predicho para cada celda es el efecto aditivo de los factores A y B. Por ejemplo, el valor esperado para la condición B+D es igual a la gran media (183) más el efecto del biofeedback (177-183=-6) y más el efecto de la droga (178-183=-5). Es decir, el valor predicho para la celda B+D=183-6-5=172. En general, este modelo se puede parametrizar de la siguiente manera:3 Esto es así porque si no hay interacción entonces el valor de \(\hat{\alpha\beta}=0\)
\[Y_{ijk}=\mu+\alpha_j+\beta_j+\varepsilon_{ijk}\] Por lo que el valor predicho se puede estimar usando: \[\hat{Y}_{ijk}(R)=\overline{Y}_{..}+(\overline{Y}_j-\overline{Y}_{..})+(\overline{Y}_k-\overline{Y}_{..})=\overline{Y}_j+\overline{Y}_k-\overline{Y}_{..}\] Por ejemplo, para cualquier sujeto del grupo B+D el valor predicho es
\[\hat{Y}_{ijk}(R)=\overline{Y}_j+\overline{Y}_k-\overline{Y}_{..}=178+177-183=172\]
Finalmente, se obtiene el error para este modelo como la diferencia entre cada valor observado y su correspondiente valor predicho, o bien, mediante la fórmula:
\[E_R-E_F=SS_{AB}=n\sum_{j=1}^a\sum_{k=1}^b(\overline{Y}_{jk}-\overline{Y}_{j.}-\overline{Y}_{k.}+\overline{Y}_{..})\] Reemplazando, tenemos: \[SS_{AB}=(5)(168-177-178+183)^2+(5)(186-177-188+183)^2+(5)(188-189-178+183)^2+(5)(190-189-188+183)^2\] \[SS_{AB}=(5)(-4)^2+(5)(4)^2+(5)(4)^2+(5)(-4)^2=320\] Finalmente, para obtener F debemos dividir cada término por sus correspondientes grados de libertad. Para calcular \(MS_W\) debemos tener en cuenta que se estimó una media para cada combinación de niveles (ab) y que cada ab tiene el mismo número de individuos.
\[df_F=nab-ab\]
\[df_F=(5)(2\times2)-(2\times2)=16\] Mientras que para \(SS_A\) o \(SS_B\) el número de grados de libertad corresponde es \(a-1\) o \(b-1\). En este caso, los grados de libertad de cada efecto principal es 1.
Por último, para el efecto de interacción, los grados de libertad se obtienen como \(df_R-df_F=(a-1)(b-1)\) por lo que en nuestro caso los grados de libertad para \(SS_{AB} = df_R-df_F=(2-1)(2-1)=1\)
Podemos comprobar que hemos obtenido bien los componentes para calcular F si hacemos ANOVA factorial en R.
aov2<-aov(Score~Biofeedback*DrugTherapy,data=d1)
summary(aov2)## Df Sum Sq Mean Sq F value Pr(>F)
## Biofeedback 1 500 500.0 8.00 0.01211 *
## DrugTherapy 1 720 720.0 11.52 0.00371 **
## Biofeedback:DrugTherapy 1 320 320.0 5.12 0.03792 *
## Residuals 16 1000 62.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Es interesante notar que la suma de cuadrados asociada al factor Group es la suma de las sumas de cuadrado de los efectos principales más el efecto de interacción (i.e., 500+720+320=1540).
aov1<-aov(Score~Group,data=d1)
summary(aov1)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 3 1540 513.3 8.213 0.00155 **
## Residuals 16 1000 62.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Una última observación sobre este tema, es que aunque se realizan 3 test, en la práctica no se controla la inflación del error Tipo I porque se asume que cada contraste proviene de una familia de tests (preguntas de investigación) distintas.
Para evaluar la magnitud del efecto de los efectos principales y del efecto de interacción usamos \(R^2\) o \(\omega^2\). Cualquiera de estas medidas indica qué proporción de la variabilidad de la de la data se asocia a un factor. Una forma de calcular \(R^2\) es comparar el efecto con la varianza total del modelo.
\[R^2=\frac{E_R-E_F}{E_F}=\frac{SS_{efecto}}{SS_{total}}=\frac{SS_{efecto}}{SS_W+SS_A+SS_B+SS_{AB}}\]
donde \(SS_{efecto}\) corresponde al efecto del factor A, B, o su interacción (AB). Sin embargo, en un diseño factorial es más conveniente usar una medida que se llama \(R^2\) parcial que sólo considera en el denominador \(SS_W\) y \(SS_{efecto}\)
\[R^2_p=\frac{E_R-E_F}{E_F}=\frac{SS_{efecto}}{SS_W+SS_{efecto}}\]
aov2<-aov(Score~Biofeedback*DrugTherapy,data=d1)
summary(aov2)## Df Sum Sq Mean Sq F value Pr(>F)
## Biofeedback 1 500 500.0 8.00 0.01211 *
## DrugTherapy 1 720 720.0 11.52 0.00371 **
## Biofeedback:DrugTherapy 1 320 320.0 5.12 0.03792 *
## Residuals 16 1000 62.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Tomando como referencia la tabla de resultados de ANOVA de 2$$2 obtenemos los siguientes valores para \(R^2\) parcial:
\[R^2_p=\frac{SS_{efecto}}{SS_W+SS_{efecto}}\] Para el el factor biofeedback:
\[R^2_p=\frac{500}{1000+500}=.33\] \[R^2_p=\frac{720}{1000+720}=.42\] \[R^2_p=\frac{320}{1000+320}=.24\] El nombre “parcial” tiene que ver con la forma en que se calcula la varianza error (o el modelo saturado) en un diseño factorial. Dado que en el cálculo de \(MS_W\) se compara cada valor observado con la media de su respectiva celda, este valor representa el error del modelo controlando por todos los factores. Entonces, cuando se evalúa el efecto de cualquier factor, el denominador refleja la varianza error que queda después de evaluar el efecto del otro factor.4 Compare \(MS_W\) del diseño factorial con \(MS_W\) de un modelo usando sólo el facto A o sólo el factor B. En consecuencia, \(R^2\) parcial va a ser siempre mayor que \(R^2\) (porque el denominador es más pequeño).
Para obtener \(\omega^2\) parcial usamos la siguiente fórmula:
\[\hat{\omega^2}_{p}=\frac{SS_{efecto}-df_{efecto}MS_W}{SS_{efecto}+(N-df_{efecto})MS_W}\]
Por ejemplo, para el factor biofeedback,\(\hat{\omega^2}_{p}\) sería:
\[\hat{\omega^2}_{p}=\frac{500-(1)(62.5)}{500+(20-1)62.5}=\frac{437.5}{1687.5}=.26\]
O, usando la fórmula de T&F (2007) \[\hat{\omega^2}_{p}=\frac{df_{efecto}(MS_{efecto}-MS_{error})/N}{(df_{efecto}(MS_{efecto}-MS_{error}))/N+MS_{error}}\] \[\hat{\omega^2}_{p}=\frac{(1)(500-62.5)/20}{((1)(500-62.5))/20+62.5}=\frac{437.5/20}{437.5/20+62.5}=.26\] #Comparaciones Cuando el efecto de la interacción no es estadísticamente significativo, podemos evaluar las diferencias entre medias marginales usando el mismo enfoque que vimos en la clase de comparaciones. Cuando hay un efecto de interacción la interpretación de las medias marginales es ambigua. Sin embargo, podemos evaluar lo que se conoce como efectos simples para los efectos principales (seguido de comparaciones simples) o contrastes de interacción para el(o los) efectos de interacción. ##Efectos simples Dado que en nuestro ejemplo observamos un efecto de interacción significativo, podemos hacer un test para evaluar el efecto del factor A eliminando el efecto del factor B. Esto se puede conseguir generando un subset de la data para cada nivel de B.
aov_A<-aov(Score~Biofeedback,data=subset(d1,DrugTherapy=="Presencia"))
summary(aov_A)## Df Sum Sq Mean Sq F value Pr(>F)
## Biofeedback 1 810 810.0 12.96 0.00698 **
## Residuals 8 500 62.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Y para el segundo nivel de B
aov_A<-aov(Score~Biofeedback,data=subset(d1,DrugTherapy=="Ausencia"))
summary(aov_A)## Df Sum Sq Mean Sq F value Pr(>F)
## Biofeedback 1 10 10.0 0.16 0.7
## Residuals 8 500 62.5Una dificultad con este procedimiento es el cálculo de \(MS_W\) que se basa en una parte de la data. Asumiendo que las varianzas de los grupos son homogéneas es preferible usar \(MS_W\) de todas las celdas. Para ello podemos usar un contraste.
aov3<-aov(Score~Group,data=d1)
summary(aov3)## Df Sum Sq Mean Sq F value Pr(>F)
## Group 3 1540 513.3 8.213 0.00155 **
## Residuals 16 1000 62.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1aov3.emm = emmeans(aov3, specs =~Group)
Contrasts.4 = list(BD_vs_D = c(1,0,-1,0))
contrast(aov3.emm, Contrasts.4)## contrast estimate SE df t.ratio p.value
## BD_vs_D -18 5 16 -3.6 0.0024La diferencia entre la condición B+D y D es de -18 puntos. El valor de la prueba t al cuadrado es igual al valor F de 12.96 obtenido anteriormente. Esto es así porque en este ejemplo todos los grupos tienen la misma varianza y, por lo tanto, el valor de \(MS_W\) va a ser el mismo independientemente de como divida la muestra para hacer las distintas comparaciones.
Como ilustración, el valor de F en este caso de puede obtener como:
mean(c(168,186))## [1] 1775*((168-177)^2+(186-177)^2)## [1] 810810/62.5## [1] 12.96El contraste para el efecto del biofeedback para el otro nivel de B se obtiene mediante:
Contrasts.5 = list(BD_vs_D = c(0,1,0,-1))
contrast(aov3.emm, Contrasts.5)## contrast estimate SE df t.ratio p.value
## BD_vs_D -2 5 16 -0.4 0.6944Cuando no hay una droga presente, la diferencia entre la presencia o ausencia de biofeedback no es significativa.
Una función del paquete emmeans es muy útil para resolver en un paso todas estas preguntas.
aov4<-aov(Score~Biofeedback*DrugTherapy,data=d1)
aov4.emm = emmeans(aov4,specs=~Biofeedback*DrugTherapy)
#Este modelo es correcto cuando NO hay interacción
emmeans(aov4.emm, pairwise~Biofeedback)## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## Biofeedback emmean SE df lower.CL upper.CL
## Presencia 178 2.5 16 172.7002 183.2998
## Ausencia 188 2.5 16 182.7002 193.2998
##
## Results are averaged over the levels of: DrugTherapy
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Presencia - Ausencia -10 3.535534 16 -2.828 0.0121
##
## Results are averaged over the levels of: DrugTherapyemmeans(aov4.emm, pairwise~DrugTherapy)## NOTE: Results may be misleading due to involvement in interactions## $emmeans
## DrugTherapy emmean SE df lower.CL upper.CL
## Presencia 177 2.5 16 171.7002 182.2998
## Ausencia 189 2.5 16 183.7002 194.2998
##
## Results are averaged over the levels of: Biofeedback
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## Presencia - Ausencia -12 3.535534 16 -3.394 0.0037
##
## Results are averaged over the levels of: Biofeedback#Este es el modelo correcto cuando hay una interacción significativa
emmeans(aov4.emm, pairwise~Biofeedback|DrugTherapy)## $emmeans
## DrugTherapy = Presencia:
## Biofeedback emmean SE df lower.CL upper.CL
## Presencia 168 3.535534 16 160.505 175.495
## Ausencia 186 3.535534 16 178.505 193.495
##
## DrugTherapy = Ausencia:
## Biofeedback emmean SE df lower.CL upper.CL
## Presencia 188 3.535534 16 180.505 195.495
## Ausencia 190 3.535534 16 182.505 197.495
##
## Confidence level used: 0.95
##
## $contrasts
## DrugTherapy = Presencia:
## contrast estimate SE df t.ratio p.value
## Presencia - Ausencia -18 5 16 -3.6 0.0024
##
## DrugTherapy = Ausencia:
## contrast estimate SE df t.ratio p.value
## Presencia - Ausencia -2 5 16 -0.4 0.6944