1 Introducción

El modelo ANOVA de dos vías es la puerta a la extensión de los diseños factoriales con todos los criterios de clasificación posibles. En teoría es posible hacer un experimento y su análisis con un número muy grande de factores con muchos niveles. Sin embargo la interpretación de los análisis de diseños experimentales con muchos factores con muchos niveles puede ser muy complicado. Es por ello que se recomienda que los experimentos sean tan sencillos como sea posible o sensible desde un punto de vista estadístico y biológico.

Como es de esperarse el primer modelo más complejo que el ANOVA de dos vías es aquel que utiliza tres criterios de clasificación. Estadísticamente, es complicado hacer el análisis de dichos diseños, máxime si existen factores de agrupación como los bloques experimentales. De hecho, los paquetes estadísticos de R especializados en la modelación con efectos mixtos únicamente permiten la especificación de efectos de segundo nivel, es decir, interacciones de dos factores, independientemente de la existencia de más factores experimentales.

Para esta clase haremos el análisis completo de la base npk con efectos mixtos utilizando el paquete lme4.

2 Ejemplos de análisis

2.1 ANOVA de tres vías

2.1.1 Análisis de npk con interacción de \(3^{er}\) grado

Para una descripción de la base npk pueden dirigirse a la clase Diseños con efectos mixtos. Aquí nos limitaremos a analizar la base de datos comenzando por el modelo:

yield ~ N * P * K

lo que equivale a la fórmula expandida:

yield ~ N + P + K + N:P + N:K + P:K + N:P:K

Si nos referimos a la clase de interpretación de las interacciones, recordaremos que para una interacción de tres factores la hipótesis nula que estamos probando es la de no diferencia entre los grupos que recibieron nitrógeno, fósforo y potasio con aquellos que no los recibieron todos. La misma interpretación aplica para los términos de segundo orden (N:P, N:K y K:P).

La \(SS-II\) para este modelo es más complicada que para el ANOVA de dos vías. Para comprender la manera en que los efectos se están estimando veamos cuáles son las sumas de cuadrados que se ajustan y los modelos con los que se evalúa la significancia de la suma de cuadrados para cada término del modelo:

Término Modelo nulo Modelo alternativo
N P + K + P:K N + P + K + P:K
P N + K + N:K N + P + K + N:K
K N + P + N:P N + P + K + N:P
N:P N + P + K + N:K + P:K N + P + K + N:P + N:K + P:K
N:K N + P + K + N:P + P:K N + P + K + N:P + N:K + P:K
P:K N + P + K + N:P + N:K N + P + K + N:P + N:K + P:K
N:P:K N + P + K + N:P + N:K + P:K N + P + K + N:P + N:K + P:K + N:P:K

En esta tabla, modelo nulo es el modelo con que el modelo alternativo se compara. Nota que el término para el cual se estiman los efectos siempre está presente en el modelo alternativo. Además, nunca hay términos de mayor grado que incluyan al término cuyo efecto se está estimando, lo cual se conoce como principio de marginalidad.

Continuando con el análisis:

m1 <- lme4::lmer(yield ~ N * P * K + (1|block), npk)
car::Anova(m1, type = 2)
## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: yield
##         Chisq Df Pr(>Chisq)    
## N     12.2587  1  0.0004631 ***
## P      0.5441  1  0.4607262    
## K      6.1657  1  0.0130252 *  
## N:P    1.3783  1  0.2403915    
## N:K    2.1460  1  0.1429445    
## P:K    0.0312  1  0.8598063    
## N:P:K  0.4832  1  0.4869680    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Podemos ver que ninguno de los términos de grado mayor a primero son significativos y que tienen sumas de cuadrados más bajas. Entonces es posible que un modelo más sencillo sea el que mejor explique los datos observados. El mismo proceso exhaustivo de eliminación de términos no significativos nos lleva a la misma conclusión que vimos en el tutorial de intepretación de las interacciones.

2.2 Análisis de heartattack

El paquete datasets contiene las bases más populares para enseñanza en R, y datarium está especializado en datos para su análisis con ANOVA. La base que analizaremos hoy está contenida en datarium y se llama heartattack. En breve, esta base de datos fue generada en un estudio para medir la efectividad de tres medicamentos (A, B y C) para disminuir los niveles de colesterol sanguíneo en hombres y mujeres (male y female) en dos categorías de riesgo (low y high) de infarto al miocardio.

colesterol <- datarium::heartattack
knitr::kable(head(colesterol))
gender risk drug cholesterol id
male low A 5.238730 1
male low A 5.075693 2
male low A 4.678654 3
male low A 5.361076 4
male low A 4.957381 5
male low A 4.828897 6 
knitr::kable(tail(colesterol))
gender risk drug cholesterol id
female high C 4.953143 67
female high C 5.468236 68
female high C 5.367793 69
female high C 5.311907 70
female high C 5.755958 71
female high C 5.258370 72

Tenemos entonces que la variable de respuesta es cholesterol, y que las variables independientes son gender, risk y drug. La columna id es la identidad del voluntario en el estudio. Los niveles con que cuenta cada factor experimental son 2, 2, y 3, por lo que tenemos 12 combinaciones totales de tratamientos:

tabla.combn <- expand.grid( #Creando la tabla de combinaciones
      gender = c("male", "female"),
      risk = c("low", "high"),
      drug = c("A", "B", "C")
      )
knitr::kable(tabla.combn, caption = "Tabla de combinaciones de tratamientos de `heartattack`")
Table 2.1: Tabla de combinaciones de tratamientos de heartattack
gender risk drug
male low A
female low A
male high A
female high A
male low B
female low B
male high B
female high B
male low C
female low C
male high C
female high C

2.2.1 Análisis previo de normalidad

Comenzaremos por verificar visualmente si cholesterol cumple con el supuesto de normalidad:

boxplot(cholesterol ~ gender + risk + drug, colesterol)

Parece que en este caso también tenemos razones para sospechar que colesterol, dentro de cada grupo experimental no es normal. Entonces revisaremos con shapiro.test dentro de cada grupo. Vamos a utilizar un método que nos permitirá hacer la prueba con sólo un comando:

library(tidyverse); library(rstatix)
shap.test <- colesterol %>%
  group_by(gender, risk, drug) %>%
  shapiro_test(cholesterol)
knitr::kable(shap.test, caption = "Resultados de la prueba de Shapiro dentro de cada grupo experimental.")
Table 2.2: Resultados de la prueba de Shapiro dentro de cada grupo experimental.
gender risk drug variable statistic p
male high A cholesterol 0.9584342 0.8075952
male high B cholesterol 0.9017300 0.3842567
male high C cholesterol 0.9554783 0.7843029
male low A cholesterol 0.9821936 0.9619276
male low B cholesterol 0.9202722 0.5073236
male low C cholesterol 0.9240513 0.5350067
female high A cholesterol 0.7141795 0.0086901
female high B cholesterol 0.9393324 0.6538245
female high C cholesterol 0.9712317 0.9006042
female low A cholesterol 0.9325506 0.5999377
female low B cholesterol 0.9267738 0.5554278
female low C cholesterol 0.9575847 0.8009463

Podemos ver que sólo hay un grupo (female, high, A) que no cumple con el supuesto de normalidad. Veamos entonces si este grupo solo puede afectar la validez de toda la prueba. Haremos ANOVA con lm ya que no hay ningún factor de agrupación de tratamientos experimentales (bloques).

m1 <- lm(cholesterol ~ gender * risk * drug, data = colesterol)
car::Anova(m1, type = 2)
## Anova Table (Type II tests)
## 
## Response: cholesterol
##                  Sum Sq Df F value    Pr(>F)    
## gender           1.3672  1 16.1957 0.0001625 ***
## risk             7.8251  1 92.6988   8.8e-14 ***
## drug             1.2354  2  7.3177 0.0014328 ** 
## gender:risk      0.0119  1  0.1411 0.7084867    
## gender:drug      0.5636  2  3.3384 0.0422001 *  
## risk:drug        0.1204  2  0.7131 0.4942214    
## gender:risk:drug 1.2504  2  7.4063 0.0013345 ** 
## Residuals        5.0649 60                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Como podemos ver, la triple interacción gender:risk:drug es bastante significativa, y de hecho explica una buena parte de la variación, pero también tenemos dos términos de segundo orden, gender:risk y risk:drug que no son significativos, y de hecho sus sumas de cuadrados son bastante pequeñas. Es importante notar esto, pues si recordamos el principio de marginalidad, la base para \(SS-II\), indica que si un término de mayor orden incluye otro término de orden menor ambos deberán estar en el modelo. Este principio es sólo eso, un principio y se puede ignorar, porque de manera práctica podríamos ajustar un modelo que no incluya esos términos que no fueron significativos. Antes de proceder, podemos verificar que el modelo sea correcto:

par(mfrow = c(2,2))
plot(m1)

Este diagnóstico sugiere que el modelo es adecuado para los datos. Entonces sólo vamos a hacer una comparación con un modelo que no incluya al término con la menor suma de cuadrados (gender:risk), y compararemos ambos modelos usando anova:

m2 <- lm(cholesterol ~ gender + risk + drug + gender:drug + risk:drug + gender : risk : drug, data = colesterol)
anova(m1, m2)
## Analysis of Variance Table
## 
## Model 1: cholesterol ~ gender * risk * drug
## Model 2: cholesterol ~ gender + risk + drug + gender:drug + risk:drug + 
##     gender:risk:drug
##   Res.Df    RSS Df   Sum of Sq F Pr(>F)
## 1     60 5.0649                        
## 2     60 5.0649  0 -4.4409e-15

Y la diferencia en la suma de cuadrados explicada por ambos modelos es mínima, por lo que es probable que eliminar ese término no beneficie ni perjudique. Por lo tanto concluimos que el efecto de la droga depende tanto de si el que la recibe es hombre o mujer y se est@ está en un grupo de alto riesgo. Ahora, para ver si el nivel de colesterol fue significativamente diferente de la media global podemos volver a ajustar un modelo con las medias para male y female como intercepto y ver el coeficiente para las triples interacciones:

m3 <- lm(cholesterol ~ 0 + gender * risk * drug, colesterol)
summary(m3)
## 
## Call:
## lm(formula = cholesterol ~ 0 + gender * risk * drug, data = colesterol)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.69387 -0.20807  0.01759  0.16801  0.58370 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## gendermale                   6.1256     0.1186  51.644  < 2e-16 ***
## genderfemale                 5.2092     0.1186  43.918  < 2e-16 ***
## risklow                     -1.1022     0.1677  -6.571 1.34e-08 ***
## drugB                       -0.6868     0.1677  -4.094 0.000129 ***
## drugC                       -0.8625     0.1677  -5.142 3.14e-06 ***
## genderfemale:risklow         0.7912     0.2372   3.335 0.001466 ** 
## genderfemale:drugB           0.8394     0.2372   3.538 0.000785 ***
## genderfemale:drugC           1.0059     0.2372   4.240 7.84e-05 ***
## risklow:drugB                0.4944     0.2372   2.084 0.041436 *  
## risklow:drugC                0.7571     0.2372   3.191 0.002253 ** 
## genderfemale:risklow:drugB  -1.0297     0.3355  -3.069 0.003221 ** 
## genderfemale:risklow:drugC  -1.1895     0.3355  -3.546 0.000766 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2905 on 60 degrees of freedom
## Multiple R-squared:  0.9974, Adjusted R-squared:  0.9968 
## F-statistic:  1882 on 12 and 60 DF,  p-value: < 2.2e-16

Los dos efectos estimados para las triples interacciones son significativamente diferentes de su respectivo intercepto, y efectivamente disminuyeron la concentracion de colesterol.

3 Aplicación del conocimiento

Dirígete a classroom para completar la actividad correspondiente a esta sección.

Regresar al índice del curso