Ejemplo de análisis de varianza
Edimer David Jaramillo
2021-09-24
1 Datos de experimento
library(readxl)
library(tidyverse)
datos <- read_excel("crecimiento_pan.xlsx")
datos- Convertimos la variable tiempo a factor:
datos <- datos %>%
mutate(Tiempo = as.factor(Tiempo))
datos2 Descriptivos básicos
datos %>%
group_by(Tiempo) %>%
summarise(promedio = mean(Altura),
desviacion = sd(Altura))3 Boxplot
datos %>%
ggplot(mapping = aes(x = Tiempo, y = Altura)) +
geom_boxplot()
4 Análisis de varianza
4.1 Hipótesis
\[H_0: \mu_{35} = \mu_{40} = \mu_{45} \\ H_1: \mu_i \neq \mu_j\]
4.2 Nivel de significancia
En este caso vamos a utilizar un nivel de significancia del 5% (0.05).
4.3 Modelo
\[y_{ij} = \mu + \beta_j + \epsilon_{ij}\]
4.4 Ajuste del modelo
modelo <- aov(formula = Altura ~ Tiempo, data = datos)
summary(modelo)## Df Sum Sq Mean Sq F value Pr(>F)
## Tiempo 2 21.57 10.786 4.602 0.042 *
## Residuals 9 21.09 2.344
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Conclusión: como el valor p (0.042) es menor que el nivel de significancia (0.05) existe evidencia para rechazar la hipótesis nula, es decir, que al menos un par de tratamientos difieren estadísticamente entre sí.
4.5 Diagnósticos del modelo
4.5.1 Normalidad
\[H_0: \epsilon \sim N(\mu, \sigma) \\ H_1: \epsilon \nsim N(\mu, \sigma)\]
- Podemos obtener los residuales del modelo con la función residuals():
residuales <- residuals(modelo)- Gráfico cuantil-cuantil:
library(ggpubr)
ggqqplot(data = residuales)
- Prueba de shapiro wilk sobre los residuales:
shapiro.test(x = residuales)##
## Shapiro-Wilk normality test
##
## data: residuales
## W = 0.93225, p-value = 0.4046- Conclusión: como el valor p (0.4046) es mayor que el nivel de significancia (0.05), no existe evidencia para rechazar la hipótesis nula, es decir, que los residuales se distribuyen de forma normal
4.5.2 Homocedasticidad
\[H_0: Sí\ hay\ homocedasticidad \\ H_1: No\ hay\ homocedasticidad\]
- Prueba de Bartlett:
bartlett.test(residuales ~ datos$Tiempo)##
## Bartlett test of homogeneity of variances
##
## data: residuales by datos$Tiempo
## Bartlett's K-squared = 1.4732, df = 2, p-value = 0.4787Conclusión: como el valor p (0.4787) es mayor que el nivel de significancia (0.05) no existe evidencia para rechazar la hipótesis nula, es decir, que los residuales tienen varianza constante (son homocedásticos)
Prueba de Levene: en caso de que la normalidad de los residuales no se cumpla, se recomienda utilizar la prueba de Levene para validar la homocedasticidad de los residuales.
library(car)
leveneTest(residuales ~ datos$Tiempo)4.6 Estimativas del modelo
4.6.1 Estimación de medias
- Podemos usar la función model.tables() con el tipo “means” sobre el modelo que ejecutamos previamente.
model.tables(x = modelo, type = "means")## Tables of means
## Grand mean
##
## 7.333333
##
## Tiempo
## Tiempo
## 35 40 45
## 5.438 8.250 8.313- Estimación de medias con intervalos de confianza:
library(emmeans)
emmeans(object = modelo, specs = "Tiempo")## Tiempo emmean SE df lower.CL upper.CL
## 35 5.44 0.765 9 3.71 7.17
## 40 8.25 0.765 9 6.52 9.98
## 45 8.31 0.765 9 6.58 10.04
##
## Confidence level used: 0.954.6.2 Comparación de medias (Tukey)
- Tabla de comparación de medias:
TukeyHSD(x = modelo, which = "Tiempo", conf.level = 0.95)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = Altura ~ Tiempo, data = datos)
##
## $Tiempo
## diff lwr upr p adj
## 40-35 2.8125 -0.2099347 5.834935 0.0676257
## 45-35 2.8750 -0.1474347 5.897435 0.0618204
## 45-40 0.0625 -2.9599347 3.084935 0.9981643- Gráfico de comparación de medias:
library(broom)
TukeyHSD(x = modelo, which = "Tiempo", conf.level = 0.95) %>%
tidy() %>%
ggplot(mapping = aes(x = contrast, y = estimate,
ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_errorbar(width = 0.2) +
geom_hline(yintercept = 0, lty = 2, color = "red") +
coord_flip()