ANOVA
Miguel Tripp
1. Generalidades
El análisis de variance (ANOVA o ANDEVA) se emplea cuando se desea comparar las medias de dos o más grupos.
La hipótesis nula de la que parten los diferentes tipos de ANOVA es que la media de la variable estudiada es la misma en los diferentes grupos, en contraposición de la hipótesis alternativa de que al menos dos medias difieren de forma significativa.
\(H_0\): No hay diferencias entre las medias de los diferentes grupos: \(\mu_1 = \mu_2 ... \mu_k = \mu\)
\(H_1\): Al menos un par de medias son significativamente diferentes
El funcionamiento básico de la ANOVA consiste en calcular la media de cada uno de los grupos para despúes comparar la varianza de estas medias (varianza explicada por grupo) frente a la varianza promedio dentro de los grupos.
Bajo la hipótesis nula de que las observaciones de los distintos grupos proceden de la misma población (misma media y varianza), la varianza ponderara entre grupos será la misma que la varianza promedio dentro de los grupos.
Conforme las medias de los grupos estén más alejadas entre ellas, la varianza entre medias se incrementará y dejará de ser igual a la varianza promedio dentro de los grupos.
El estadístico estudiado en ANOVA, conocido como \(F_{ratio}\) es el ratio entre la varianza de las medias de los grupos y el promedio de la varianza dentro de los grupos.
De esta manera, si \(S_1^{2}\) es la varianza de una muestra de tamaño \(N_1\) extraida de una población normal de varianza \(\sigma_1^{2}\) y \(S_2^{2}\) es la varianza de una muestra tamaño \(N_2\) extraida de una poblacion normal\(N_2\) con varianza \(\sigma_2^{2}\), el cociente:
\(F = \frac{S_1^{2}/\sigma_1^{2}} {S_2^{2}/\sigma_2^{2}}\)
se distribuye como una variable F con (\(N_1\) y \(N_2\)) grados de libertad.
2: Diferencias entre Prueba-T y ANOVA
Prueba T de Student
- La prueba T examina si las medias de dos poblaciones difieren una de la otra.
- Se usa cuando la desviación estandar no se conoce o cuando N es pequeño (30).
- La prueba se basa en el estadístico T, el cual asume que la variable tiene una distribución normal
- La hipótesis nula toma la forma \(H_0: \mu_1 = \mu_2\) en contraparte con la hipótesis alternativa \(H_A: \mu_1 \ne \mu_2\). Los grados de libertad son \(n_1 + n_2 -2\).
ANOVA
- Se usa para comparar las medias de dos o mas poblaciones.
- Se asume que las muestras se obtienen de una población ocn dsitribución normal y con varianzas iguales.
- El principio básico es probar la varianza entre las medias poblacionales analizando la cantidad de variación dentro de grupos, en proporción a la variación entre grupos. *La \(H_0\) establece que las medias de todas las poblaciones son iguales, mientras que la \(H_A\) dice que al menos una pobación es diferente.
3: Entendiendo ANOVA
Esta sección esta tomada del libro The new statistics with R. An introduction for biologist.
Para esta sección, usaremos los datos de altura de 30 plantas de maíz tras polinización cruzada o propia. La pregutna que trataremos de resolver es si la endogamia reduce la condición de la planta.
Empecemos por cargar los paquetes necesarios y abrir la tabla Darwin.csv
library(tidyverse)
library(rstatix)
library(arm)darwin = read_csv("Data/Darwin.csv")## Parsed with column specification:
## cols(
## pot = col_character(),
## pair = col_double(),
## type = col_character(),
## height = col_double()
## )ggplot(darwin, aes(x = type, y = height, col = type, shape = type))+
geom_point()
Esta primera inspección visual sugiere que la altura promedio es mayor en la planta con cruza.
3.1: Suma de cuadrados
La estrategia general de la ANOVA es cuantificar la variabilidad en los datos y dividirla entre la variabilidad entre y dentro de los grupos, para calcular la proporcion o signal to noise ratio.
La suma de cuadrados representa una medida de variación o desviación con respecto a la media. El cálculo de la suma total de los cuadrados considera tanto la suma de los cuadrados de los factores como la de aleatoriedad o error.
3.2 ANOVA como modelo lineal
En el fondo, tanto la ANOVA, el análisis de covarianza (ANCOVA) y la regresión lineal son modelos lineales muy relacionados.
Para hacer una ANOVA usamos la función lm()
ls0 <- lm(height ~ 1, data = darwin)En este caso, usamos 1 para especificar que no tenemos una variable explanatoria. Por lo que solo nos da el promedio total (gran media).
que es igual:
mean(darwin$height)## [1] 18.88333Ahora utilizaremos la función display() del paquete arm para extraer la información del modelo
display(ls0)## lm(formula = height ~ 1, data = darwin)
## coef.est coef.se
## (Intercept) 18.88 0.58
## ---
## n = 30, k = 1
## residual sd = 3.18, R-Squared = 0.00La primera fila corresponde al intercepto, que en este caso es simplemente la media total.
El tamaño de muestra n nos dice el número de puntos (y de filas en la tabla) y el número de parámetros estimados -media total- indicado como k.
Ahora ajustaremos el modelo lineal utilizando la altura de la planta como variable de respuesta en función del tipo de polinización
ls1 <- lm(height ~ 1 + type, data = darwin)En este modelo, el intercepto ya no es la gran media, tal como se observa con la función display:
display(ls1)## lm(formula = height ~ 1 + type, data = darwin)
## coef.est coef.se
## (Intercept) 20.19 0.76
## typeSelf -2.62 1.07
## ---
## n = 30, k = 2
## residual sd = 2.94, R-Squared = 0.18En este caso, la segunda fila “typeself” se produce con el nombre la variable explanatoria (type), seguido sin espacio del nombre de uno de los niveles de este factor (self). Por eliminación, la fila “intercepto” corresponde al primer factor (Cross).
Es importante notar que el valor del segundo factor (typeself) corresponde a la diferencía entre este valor y la media del otro nivel.
La parte inferior de los resultados con display(ls1) indica el número de paramtros calculados k que es de 2, que corresponde a los dos tratamientos de nuestro factor.
R square es la proporción de la suma de cuadrados explicado por el modelo estadistico.
Residual SD corresponde a la desviación de los residuales (no explicados por el modelo).
3.3: Tablas ANOVA
Una tabla ANOVA resume los calculos de las sumas de cuadrados
La tabla de anova resultante de modelo con la función lm()se puede obtener con la función anova()
anova(ls1)## Analysis of Variance Table
##
## Response: height
## Df Sum Sq Mean Sq F value Pr(>F)
## type 1 51.352 51.352 5.9395 0.02141 *
## Residuals 28 242.083 8.646
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En esta tabla, la primera columna indica la fuente de variación (tratamiento) y la variación residual (no explicada por el modelo)
El valor F se calcula dividiendo la varianza de tratamiento entre la varianza residual del error. En este caso, el valor de 5.9 significa que la señal es seis veces mayor que el ruido
El último argumento proporciona el valor de probabilidad el cual se calcula a partir del valor F.
Finalmente, el resutlado se reportaria como:
“La altura de las plantas con polinización cruzada fue significativamente mayor que la altura de las plantas auto-polinizadas (\(F_{1,28}\) = 5.9; P = 0.02).”
4: Ejercicio. ¿Hay diferencias en el peso corporal de los pingüinos Adelie machos entre islas?
Para este ejercicio, vamos a utilziar los datos de penguins_size.csv y filtramos solamente los machos de la especie Adelie
pin <- read_csv("Data/penguins_size.csv") ## Parsed with column specification:
## cols(
## species = col_character(),
## island = col_character(),
## culmen_length_mm = col_double(),
## culmen_depth_mm = col_double(),
## flipper_length_mm = col_double(),
## body_mass_g = col_double(),
## sex = col_character()
## )pin_Ade <- pin %>%
filter(species == "Adelie", sex == "MALE")- Descripciónd de los datos utilizando histogramas de frecuencia.
ggplot(pin_Ade, aes(x = body_mass_g, col = island, fill = island))+
geom_histogram(alpha = 0.3)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
- Probar si los datos siguen una distribución normal utilizando qqplots
ggplot(pin_Ade, aes(sample = body_mass_g, col = island))+
stat_qq()+
stat_qq_line()+
facet_grid(.~ island)
- Prueba de Normalidad con Shapiro-Wilks
pin_Ade_shap <- pin_Ade %>%
group_by(island) %>%
shapiro_test(body_mass_g)
pin_Ade_shap## # A tibble: 3 x 4
## island variable statistic p
## <chr> <chr> <dbl> <dbl>
## 1 Biscoe body_mass_g 0.948 0.292
## 2 Dream body_mass_g 0.978 0.796
## 3 Torgersen body_mass_g 0.981 0.919- Probar si las varianzas de los datos son homogeneas entre las islas mediante una prueba de Levene
pin_Ade_lev <- pin_Ade %>%
levene_test(body_mass_g ~ island)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.pin_Ade_lev## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 70 0.113 0.893- Generamos el modelo para realizar el análisis de varianza
pin_modelo <- lm(body_mass_g ~ island, data = pin_Ade)
display(pin_modelo)## lm(formula = body_mass_g ~ island, data = pin_Ade)
## coef.est coef.se
## (Intercept) 4050.00 74.98
## islandDream -4.46 100.19
## islandTorgersen -15.22 104.87
## ---
## n = 73, k = 3
## residual sd = 351.67, R-Squared = 0.00Graficas de los residuales
plot(pin_modelo)
- Finalmente usamos la función
anovapara generar la tabla ANOVA
anova(pin_modelo)## Analysis of Variance Table
##
## Response: body_mass_g
## Df Sum Sq Mean Sq F value Pr(>F)
## island 2 2793 1397 0.0113 0.9888
## Residuals 70 8657241 1236756.1. Anova con rstatixs
pin_anova <- pin_Ade %>%
anova_test(body_mass_g ~ island)## Coefficient covariances computed by hccm()pin_anova## ANOVA Table (type II tests)
##
## Effect DFn DFd F p p<.05 ges
## 1 island 2 70 0.011 0.989 0.000323- Visualizacióin de los datos con boxplots
ggplot(pin_Ade, aes(x = island, y = body_mass_g, fill = island))+
geom_boxplot(outlier.shape = "")+
geom_point(position = position_jitter(0.1))+
labs(x = "Isla", y = "Peso corporal (g)",
caption = "No hay diferencias en el peso corporal entre las tres islas \n F(2,70) = 0.01; P = 0.98")+
theme_light()
5. Ejercicio. ¿Hay diferencias en el peso corporal entre especies, sin importar lugar o sexo?
ggplot(pin, aes(x = species, y = body_mass_g, col = species))+
geom_boxplot(outlier.shape = "")## Warning: Removed 2 rows containing non-finite values (stat_boxplot).
Histograma de frecuencias
ggplot(pin, aes(x = body_mass_g, col = species, fill = species)) +
geom_histogram(alpha = 0.3)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.## Warning: Removed 2 rows containing non-finite values (stat_bin).
Summary
pin_spe_sum <- pin %>%
group_by(species) %>%
get_summary_stats(body_mass_g)
pin_spe_sum## # A tibble: 3 x 14
## species variable n min max median q1 q3 iqr mad mean sd
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Adelie body_ma~ 151 2850 4775 3700 3350 4000 650 445. 3701. 459.
## 2 Chinst~ body_ma~ 68 2700 4800 3700 3488. 3950 462. 371. 3733. 384.
## 3 Gentoo body_ma~ 123 3950 6300 5000 4700 5500 800 556. 5076. 504.
## # ... with 2 more variables: se <dbl>, ci <dbl>QQplots
ggplot(pin, aes(sample = body_mass_g, col = species))+
stat_qq()+
stat_qq_line()+
facet_grid(.~species)## Warning: Removed 2 rows containing non-finite values (stat_qq).## Warning: Removed 2 rows containing non-finite values (stat_qq_line).
Distribucion de normalidad
pin_spe_shap <- pin %>%
group_by(species) %>%
shapiro_test(body_mass_g)
pin_spe_shap## # A tibble: 3 x 4
## species variable statistic p
## <chr> <chr> <dbl> <dbl>
## 1 Adelie body_mass_g 0.981 0.0324
## 2 Chinstrap body_mass_g 0.984 0.561
## 3 Gentoo body_mass_g 0.986 0.234convertir datos a raiz cuadrada
pin_sqrt <- pin %>%
mutate(body_mass_g = sqrt(body_mass_g))Re-analizar la normalidad tras la transformación
pin_spe_shap_sqrt <- pin_sqrt %>%
group_by(species) %>%
shapiro_test(body_mass_g)
pin_spe_shap_sqrt## # A tibble: 3 x 4
## species variable statistic p
## <chr> <chr> <dbl> <dbl>
## 1 Adelie body_mass_g 0.985 0.0940
## 2 Chinstrap body_mass_g 0.987 0.676
## 3 Gentoo body_mass_g 0.987 0.264Ahora se verifica la homogeneidad de varianzas con una prueba de levene
pin_spe_leve <- pin_sqrt %>%
levene_test(body_mass_g ~ species)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.pin_spe_leve## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 2 339 2.73 0.0666En este caso, nuestros datos muestran varianzas homogeneas y distribución normal, por lo que podemos proceder con el ajusto del modelo lineal:
modelo_pin_sqrt<- lm(body_mass_g ~ species, data = pin_sqrt)
display(modelo_pin_sqrt)## lm(formula = body_mass_g ~ species, data = pin_sqrt)
## coef.est coef.se
## (Intercept) 60.72 0.29
## speciesChinstrap 0.30 0.52
## speciesGentoo 10.44 0.43
## ---
## n = 342, k = 3
## residual sd = 3.56, R-Squared = 0.66Resultados del ANOVA
anova(modelo_pin_sqrt)## Analysis of Variance Table
##
## Response: body_mass_g
## Df Sum Sq Mean Sq F value Pr(>F)
## species 2 8437.2 4218.6 331.93 < 2.2e-16 ***
## Residuals 339 4308.4 12.7
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(modelo_pin_sqrt)
Cuando tenemos datos con distribuciones diferentes, es posible realizar la prueba relajando el supuesto de homogeneidad de varianzas utilizando la función oneway.test() lo cual nos arroja el resultado con la prueba de Welch
oneway.test(body_mass_g ~ species, data = pin_sqrt)##
## One-way analysis of means (not assuming equal variances)
##
## data: body_mass_g and species
## F = 334.32, num df = 2.00, denom df = 187.54, p-value < 2.2e-16Equivalente con rstatixs
welch_anova <- pin_sqrt %>%
welch_anova_test(body_mass_g ~ species)
welch_anova## # A tibble: 1 x 7
## .y. n statistic DFn DFd p method
## * <chr> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 body_mass_g 344 334. 2 188. 1.45e-62 Welch ANOVAEntonces podemos graficar nuestros datos con mas información
ggplot(pin, aes(x = species, y = body_mass_g, fill = species))+
geom_boxplot(outlier.shape = "")+
geom_point(position = position_jitter(0.1), col = "grey45")+
labs(caption = "Existen diferencias en el peso corporal entre especies F(2,34.57) = 345.57; P < 0.001")## Warning: Removed 2 rows containing non-finite values (stat_boxplot).## Warning: Removed 2 rows containing missing values (geom_point).
6: Comparación multiple de medias. Contrastes post-hoc
Si un análisis de varianza resulta significativo, implica que al menos dos de las medias comparadas son significativamente distintas entre si, pero no se indica cuales.
Para identificarlas hay que comparar dos a dos las medias de todos los grupos introducidos en el análisis mediate una prueba t u otro test que compare 2 grupos.
Debido a la inflación del error de tipo I, cuantas más comparaciones se hagan más aumenta la probabilidad de encontrar diferencias significativas (para \(\alpha\) = 0.05, de cada 100 comparaciones se esperan 5 significativas solo por azar).
6.1: multiples t-test
Para hacer comparaciones entre cada grupo, podemos utilizar multiples t-test usando la función pairwise.t.test
pairwise.t.test(x = pin_sqrt$body_mass_g, g = pin_sqrt$species, p.adjust.method = "bonf", paired = FALSE, pool.sd = FALSE)##
## Pairwise comparisons using t tests with non-pooled SD
##
## data: pin_sqrt$body_mass_g and pin_sqrt$species
##
## Adelie Chinstrap
## Chinstrap 1 -
## Gentoo <2e-16 <2e-16
##
## P value adjustment method: bonferroniComparaciones multiples usando rstatix
pin_mul_t <- pin_sqrt %>%
pairwise_t_test(body_mass_g ~ species,
paired = FALSE, p.adjust.method = "bonf",pool.sd = FALSE)
pin_mul_t## # A tibble: 3 x 10
## .y. group1 group2 n1 n2 statistic df p p.adj p.adj.signif
## * <chr> <chr> <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 body~ Adelie Chins~ 152 68 -0.617 153. 5.38e- 1 1.00e+ 0 ns
## 2 body~ Adelie Gentoo 152 124 -23.6 266. 2.91e-67 8.73e-67 ****
## 3 body~ Chins~ Gentoo 68 124 -20.4 153. 1.90e-45 5.70e-45 ****7: Anotación de la gráfica usando ggpubr
El paquete ggpubr contiene diversas funciones para realizar pruebas estadisticas y anotar los resultados dentro de la gráfica
install.packages("ggpubr")library(ggpubr)ggplot(pin, aes(x = species, y = body_mass_g, fill = species))+
geom_boxplot(outlier.shape = "")+
geom_point(position = position_jitter(0.1), col = "grey45")+
labs(caption = "Existen diferencias en el peso corporal entre especies F(2,34.57) = 345.57; P < 0.001")+
stat_compare_means(data = pin_sqrt,method = "anova") #agregar resultados de anova## Warning: Removed 2 rows containing non-finite values (stat_boxplot).## Warning: Removed 2 rows containing non-finite values (stat_compare_means).## Warning: Removed 2 rows containing missing values (geom_point).
my_comparisons <- list( c("Adelie", "Chinstrap"), c("Adelie", "Gentoo"), c("Chinstrap", "Gentoo") )
ggplot(pin, aes(x = species, y = body_mass_g, fill = species))+
geom_boxplot(outlier.shape = "")+
geom_point(position = position_jitter(0.1), col = "grey45")+
labs(caption = "Existen diferencias en el peso corporal entre especies F(2,34.57) = 345.57; P < 0.001")+
stat_compare_means(data = pin_sqrt, method = "t.test", na.rm = TRUE, comparisons = my_comparisons)## Warning: Removed 2 rows containing non-finite values (stat_boxplot).## Warning: Removed 2 rows containing non-finite values (stat_signif).## Warning: Removed 2 rows containing missing values (geom_point).
8: Anova de dos vias
El análisis de variaza de dos vias, (tambien conocido como factorial), sirve para estudiar la relación entre una variable dependiente cuantitativa y dos variables independientes cualitativas (factores), cada uno con varios niveles.
El ANOVA de dos vias permite estudiar como influyen por si solos cada uno de los factores sobre la variable dependiente (modelo aditivo) asi como la influencia de las combinaciones que se puedan dar entre ellas (modelo con interacción)
8.1: Jugo de naranja y conejillos de indias
Para este ejemplo, utilizaremos la base de datos toothgrowth. Esta corresponde a un experimento realizado en conejillos de indias para demostrar el efecto de adminstrar jugo de naranja (OJ) y vitamina c (VC) a diferentes dosis (0.5, 1.0 y 2.0 mg).
Para abrir la tabla ejecutamos
data("ToothGrowth")
dientes <- ToothGrowth
summary(dientes)## len supp dose
## Min. : 4.20 OJ:30 Min. :0.500
## 1st Qu.:13.07 VC:30 1st Qu.:0.500
## Median :19.25 Median :1.000
## Mean :18.81 Mean :1.167
## 3rd Qu.:25.27 3rd Qu.:2.000
## Max. :33.90 Max. :2.000dientes <- dientes %>%
mutate(dose = factor(dose))Obten el resumen de los datos
dientes_sum <- dientes %>%
group_by(supp, dose) %>%
get_summary_stats(len, type = "mean_sd")
dientes_sum## # A tibble: 6 x 6
## supp dose variable n mean sd
## <fct> <fct> <chr> <dbl> <dbl> <dbl>
## 1 OJ 0.5 len 10 13.2 4.46
## 2 OJ 1 len 10 22.7 3.91
## 3 OJ 2 len 10 26.1 2.66
## 4 VC 0.5 len 10 7.98 2.75
## 5 VC 1 len 10 16.8 2.52
## 6 VC 2 len 10 26.1 4.80Visualización de los datos
ggplot(dientes, aes(x = supp, y = len, col = dose))+
geom_boxplot()
A partir de la representación gráfica y el calculo de las medias se puede intuir que existe una diferncia en el crecimiento del diente con la dosis.
A priori, parece que se satisfacen las condiciones necesarios para realizar un ANOVA, aunque se requiere hacer las pruebas correspondientes.
dientes_shap <- dientes %>%
group_by(dose, supp) %>%
shapiro_test(len)
dientes_shap## # A tibble: 6 x 5
## supp dose variable statistic p
## <fct> <fct> <chr> <dbl> <dbl>
## 1 OJ 0.5 len 0.893 0.182
## 2 VC 0.5 len 0.890 0.170
## 3 OJ 1 len 0.927 0.415
## 4 VC 1 len 0.908 0.270
## 5 OJ 2 len 0.963 0.815
## 6 VC 2 len 0.973 0.919dientes_leven <- dientes %>%
levene_test(len ~ supp * dose)
dientes_leven## # A tibble: 1 x 4
## df1 df2 statistic p
## <int> <int> <dbl> <dbl>
## 1 5 54 1.71 0.148Vamos a empezar nuestro análisis revisando cada uno de los factores principales (main effects)supp y dose
me_supp <- lm(len ~ supp, data = dientes)
anova(me_supp)## Analysis of Variance Table
##
## Response: len
## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.4 205.35 3.6683 0.06039 .
## Residuals 58 3246.9 55.98
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Esto nos indica que no hay evidencia para rechazar la hipotesis nula que las medias por el suplemento son diferentes. Por si solo, pareciera que el tipo de suppemento (OJ y VC) no tuvieran efecto en la longitud del diente, lo cual se puede visualizar en la siguiente gráfica:
ggplot(dientes, aes(x = supp, y = len))+
geom_boxplot()
Ahora, ajustemos un modelo con el siguiente efecto principal, dose.
me_dose <- lm(len ~ dose, data = dientes)
anova(me_dose)## Analysis of Variance Table
##
## Response: len
## Df Sum Sq Mean Sq F value Pr(>F)
## dose 2 2426.4 1213.2 67.416 9.533e-16 ***
## Residuals 57 1025.8 18.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1En este caso, hay evidencia para rechazar la hipótesis nula de manera que la dosis tene un efecto en la longitud del diente:
ggplot(dientes, aes(x = dose, y = len))+
geom_boxplot()
ahora ajustemos un modelo con interacción
dientes_mod_int <- lm(len ~ supp * dose, data = dientes)
anova(dientes_mod_int)## Analysis of Variance Table
##
## Response: len
## Df Sum Sq Mean Sq F value Pr(>F)
## supp 1 205.35 205.35 15.572 0.0002312 ***
## dose 2 2426.43 1213.22 92.000 < 2.2e-16 ***
## supp:dose 2 108.32 54.16 4.107 0.0218603 *
## Residuals 54 712.11 13.19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Dado que el modelo len ~ dose mostró diferencias significativas en las medias, esperariamos que el modelo factorial diera los mismos resultados.
Por otro lado, observamos que el efecto principal supp ahora es significativo. Esto indica que al controlar por el nivel de la dosis y la interacción dose*supp hay un efecto independiente por el tipo de suplemento.
ggplot(dientes, aes(x = supp, y = len, col = dose, group = dose))+
stat_summary(fun = mean, geom = "point")+
stat_summary(fun = mean, geom = "line")
Podemos usar la herramienta del paquete HH para visualizar esta interacción
install.packages("HH")library(HH)
interaction2wt(len ~ supp * dose, data = dientes)
8.2: Luz y fertilizantes
El siguiente ejmplo es tomado del libro The new statistics with R. An introduction for biologist. La base de datos consiste en la estimación de biomasa de una comunidad de plantas tras una exposición a fertilizante (F) y luz (L) en un diseño factorial.
biomasa <- read_csv("Data/fertilizantes_luz.csv")## Parsed with column specification:
## cols(
## Fert = col_character(),
## Light = col_character(),
## FL = col_character(),
## Biomass_m2 = col_double()
## )Visualización de los datos
ggplot(biomasa, aes(x = FL, y = Biomass_m2))+
geom_point(position = position_jitter(0.1))
De manera demostrativa, ajustaremos un modelo lineal con los cuatro niveles.
mod1 <- lm(Biomass_m2 ~ FL, data = biomasa)
display(mod1)## lm(formula = Biomass_m2 ~ FL, data = biomasa)
## coef.est coef.se
## (Intercept) 355.79 23.14
## FLF-L+ 30.12 32.72
## FLF+L- 93.69 32.72
## FLF+L+ 219.23 32.72
## ---
## n = 64, k = 4
## residual sd = 92.56, R-Squared = 0.47En este caso, vemos que el valor del último coeficiente es superior a lo que se esperario bajo un modelo completamente aditivo (30 + 94).
A pesar de que este análisis sugiere que hay una interacción, no lo hace explicito. Esto es porque trata a cada uno de los cuatro factores como si fueran completamente independiente el uno del otro, cuando en realidad es una combinación de dos.
Ahora, ajustamos un modelo que considera la interacción:
mod2 <- lm(Biomass_m2 ~ Light * Fert, data = biomasa)
display(mod2)## lm(formula = Biomass_m2 ~ Light * Fert, data = biomasa)
## coef.est coef.se
## (Intercept) 355.79 23.14
## LightL+ 30.13 32.72
## FertF+ 93.69 32.72
## LightL+:FertF+ 95.41 46.28
## ---
## n = 64, k = 4
## residual sd = 92.56, R-Squared = 0.47Tras usar la combinación factorial de los dos tratamientos, la última linea de la tabla de los coeficientes estimados resulta en un tamaño de la interacción de 95 g. Podemos relacionar esta estimación al modelo anterior al restar la prediccón aditiva (124) del promedio estimado con del modelo combinado (219.23)
219.23 - (30.12+93.69)## [1] 95.42Esto puede resumirse con summary()
summary(mod2)##
## Call:
## lm(formula = Biomass_m2 ~ Light * Fert, data = biomasa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -233.619 -42.842 1.356 67.961 175.381
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 355.79 23.14 15.376 < 2e-16 ***
## LightL+ 30.13 32.72 0.921 0.36095
## FertF+ 93.69 32.72 2.863 0.00577 **
## LightL+:FertF+ 95.41 46.28 2.062 0.04359 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 92.56 on 60 degrees of freedom
## Multiple R-squared: 0.4686, Adjusted R-squared: 0.442
## F-statistic: 17.63 on 3 and 60 DF, p-value: 2.528e-08En este caso, es importante no concluir que no hay evidencia de que el effecto principal main effect de luz tenga efecto debido a su valor no significativo de P = 0.36.
El que la luz tenga un efecto depende del nivel de fertilizante
En otras palabras, al interpretar los coeficientes del modelo hay que empezar de abajo hacia arriba.
anova(mod2)## Analysis of Variance Table
##
## Response: Biomass_m2
## Df Sum Sq Mean Sq F value Pr(>F)
## Light 1 96915 96915 11.3131 0.001345 **
## Fert 1 319889 319889 37.3413 8.019e-08 ***
## Light:Fert 1 36409 36409 4.2501 0.043587 *
## Residuals 60 513998 8567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Podemos utilizar la función anova() para comparar un par de modelos e identificar el efecto de la interacción.
#ajuste de modelo sin interaccion
mod3 <- lm(Biomass_m2 ~ Light+Fert, data = biomasa)
anova(mod2, mod3)## Analysis of Variance Table
##
## Model 1: Biomass_m2 ~ Light * Fert
## Model 2: Biomass_m2 ~ Light + Fert
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 60 513998
## 2 61 550407 -1 -36409 4.2501 0.04359 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1De esta manera podemos conlcluir que existe un efectivo interactivo significativo de la luz junto con el fertilizante (\(F_{1,60} = 4.3\), P = 0.04 )
HH::interaction2wt(Biomass_m2 ~ Light * Fert, data = biomasa)
8.3: Post-hoc
Al igual que como vimos en la sección 8.2, cuando tenemos una ANOVA significativa sabemos que al menos una media es diferente al resto, pero no sabemos cual.
En una ANOVA factoral donde no hay un interacción significativa es aconsejable realizar una ANOVA para cada uno de los factores principales significativos con sus respectivas pruebas posteriori
Cuando tenemos una interacción significativa, hay que ser mas cautelosos ya que la interacción puede opacar el effecto principal
Cuando hay una interacción significativa no es apropiado confiar en los ratios F para los efectos individuales A y B
Para evaluar cada una de los contrastes, podemos (Tukey´s Honestly Significant Difference), pero para ello necesitamos un objetvo aov
mod2_aov <- aov(mod2)
anova(mod2_aov)## Analysis of Variance Table
##
## Response: Biomass_m2
## Df Sum Sq Mean Sq F value Pr(>F)
## Light 1 96915 96915 11.3131 0.001345 **
## Fert 1 319889 319889 37.3413 8.019e-08 ***
## Light:Fert 1 36409 36409 4.2501 0.043587 *
## Residuals 60 513998 8567
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Y posteriorment utilizar la función base TukeyHSD()
tukey <- TukeyHSD(mod2_aov)plot(tukey)
Alternativamente, se puede utilizar la prueba de Tukey implementada en el paquete agricolae
tukey_agri <- agricolae::HSD.test(mod2_aov, trt = c("Fert", "Light"))
tukey_agri## $statistics
## MSerror Df Mean CV MSD
## 8566.629 60 441.5547 20.96141 86.47259
##
## $parameters
## test name.t ntr StudentizedRange alpha
## Tukey Fert:Light 4 3.737089 0.05
##
## $means
## Biomass_m2 std r Min Max Q25 Q50 Q75
## F-:L- 355.7937 67.82697 16 202.0 455.3 320.725 366.15 389.925
## F-:L+ 385.9187 98.19830 16 152.3 514.5 352.675 386.80 453.600
## F+:L- 449.4875 88.04292 16 257.9 598.0 402.925 445.90 520.050
## F+:L+ 575.0187 110.77706 16 383.0 750.4 509.650 549.85 654.725
##
## $comparison
## NULL
##
## $groups
## Biomass_m2 groups
## F+:L+ 575.0187 a
## F+:L- 449.4875 b
## F-:L+ 385.9187 bc
## F-:L- 355.7937 c
##
## attr(,"class")
## [1] "group"o con rstatixs
aov(Biomass_m2 ~ Light*Fert, data = biomasa) %>%
tukey_hsd()## Warning: Expected 2 pieces. Additional pieces discarded in 8 rows [1, 2, 3, 4,
## 5, 6, 7, 8].## # A tibble: 8 x 8
## term group1 group2 estimate conf.low conf.high p.adj p.adj.signif
## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 Light "L" L+ 77.8 31.5 124. 0.00135 **
## 2 Fert "F" F+ 141. 95.1 188. 0.0000000802 ****
## 3 Light:Fert "" L+:F 30.1 -56.3 117. 0.794 ns
## 4 Light:Fert ":F+" L 93.7 7.22 180. 0.0287 *
## 5 Light:Fert "L" L+:F+ 219. 133. 306. 0.0000000484 ****
## 6 Light:Fert ":F+" L 63.6 -22.9 150. 0.222 ns
## 7 Light:Fert "L+:F" L+:F+ 189. 103. 276. 0.0000017 ****
## 8 Light:Fert "L" L+:F+ 126. 39.1 212. 0.00168 **ggplot(biomasa, aes(x = FL, y = Biomass_m2))+
geom_point(position = position_jitter(0.1))+
stat_compare_means(method = "t.test", ref.group = "F-L-", aes(label = ..p.signif..))