ANOVA script
ANOVA de una vía
Como lo prometido es deuda, te presento el scrip para realizar un análisis de varianza de una vía. Vale la pena aclarar que utilicé algunos paquetes extra, los cuales es necesario instalar, pero como todo en R la funciómn es muy sencilla además de que yo te la indicaré con un comentario dentro del texto. Sin más por agregar empezamos nuestro ejemplo.
En este caso voy a utilizar una base de datos almacenada en una url, esta base contiene información acerca de la comparación de rendimiento en distintas parcelas agrícolas de trigo sometidas a un tratamiento diferenciado de fertilizante. Nosotros esperaríamos encontrar diferencias significativas entre los tratamientos.
Carguémos los datos.
require(foreign) #esta función solicita el paquete determinado y lo carga, si no lo tienes lo instala.## Loading required package: foreignrequire(dplyr)## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, uniondata.agricola <- read.dta("http://www.stata-press.com/data/r12/yield.dta")#Cargo los datos en crudo
data.agricola <- select(data.agricola, yield,fertilizer) #selecciono sólo las columans que me interesan
summary(data.agricola)## yield fertilizer
## Min. :24.05 Min. :1
## 1st Qu.:37.35 1st Qu.:2
## Median :42.73 Median :3
## Mean :42.92 Mean :3
## 3rd Qu.:48.33 3rd Qu.:4
## Max. :60.43 Max. :5De la misma que manera que como te ocurrió en clase, el factor de tratamientos está en un formato equivocado, esto lo sabemos por que el resultado de la función summary() nos dice que la variable fertilizer es un número. Así que el siguiente paso será convertirlo en un factor y probar los supuestos de ANOVA.
Corrección de datos y supuesto de ANOVA
data.agricola$fertilizer <-factor(data.agricola$fertilizer) #convertimos en factor
summary(data.agricola)## yield fertilizer
## Min. :24.05 1:40
## 1st Qu.:37.35 2:40
## Median :42.73 3:40
## Mean :42.92 4:40
## 3rd Qu.:48.33 5:40
## Max. :60.43Ahora vemos que la variable de fertilizante es un factor con 5 tratamientos distintos. Continuemos a probar los supuestos.Primero, vamos a evaluar si el rendimiento en cada uno de los tratamientos es normal:
shapiro.test( data.agricola$yield[data.agricola$fertilizer == "1"] )##
## Shapiro-Wilk normality test
##
## data: data.agricola$yield[data.agricola$fertilizer == "1"]
## W = 0.98372, p-value = 0.8231shapiro.test( data.agricola$yield[data.agricola$fertilizer == "2"] )##
## Shapiro-Wilk normality test
##
## data: data.agricola$yield[data.agricola$fertilizer == "2"]
## W = 0.95954, p-value = 0.1616shapiro.test( data.agricola$yield[data.agricola$fertilizer == "3"] )##
## Shapiro-Wilk normality test
##
## data: data.agricola$yield[data.agricola$fertilizer == "3"]
## W = 0.96171, p-value = 0.1915shapiro.test( data.agricola$yield[data.agricola$fertilizer == "4"] )##
## Shapiro-Wilk normality test
##
## data: data.agricola$yield[data.agricola$fertilizer == "4"]
## W = 0.97247, p-value = 0.4294shapiro.test( data.agricola$yield[data.agricola$fertilizer == "5"] )##
## Shapiro-Wilk normality test
##
## data: data.agricola$yield[data.agricola$fertilizer == "5"]
## W = 0.98521, p-value = 0.8708Interpretación: Para todos los tratamientos de fertilización se obtiene un p-value mayor de 0.05 (0.8231, 0.1616, 0.1915, 0.4294 y 0.8708), por lo que no podemos rechazar la hipótesis nula (hipótesis de normalidad). Por lo tanto, podemos concluir que nuestros datos cumplen el supuesto de normalidad. Ahora evaluemos la homocedasticidad (homogeneidad de varianzas).
bartlett.test(yield~fertilizer, data =data.agricola)##
## Bartlett test of homogeneity of variances
##
## data: yield by fertilizer
## Bartlett's K-squared = 2.4365, df = 4, p-value = 0.656Interpretación: Con un p-value = 0.656, mayor de 0.05, por lo tanto, no podemos rechazar la hipótesis nula. Concluimos que existe la homogeneidad de varianzas.
Ya probamos nuestros supuestos pero vale la pena preguntarse si es pertinente la comparación, para esto vamos a hacer un boxplot().
boxplot(yield~fertilizer, data = data.agricola)
El gráfico parece indicar que existen diferencias -eso es bueno-, vamos a evaluar si éstas son significativas.
Estimación del modelo de ANOVA
En el ANOVA de una vía la hipótesis nula \(H0\) es que no hay diferencia entre las medias. Por otro lado, la hipótesis alternativa \(H1\) que al menos una de las medias difiere del resto. En nuestro caso sería que el rendimiento de cada parcela varía según el tratamiento de fertilización.
El modelo ANOVA contrastará las diferencias en las medias del rendimiento de las parcelas agrícolas según el tratamiento de fertilizante. Para estimar el modelo ANOVA de una vía se usa la función aov( ), que sigue la estructura aov(variable dependiente ~ factor, datos)
modelo <- aov(yield~fertilizer, data = data.agricola)
summary(modelo)## Df Sum Sq Mean Sq F value Pr(>F)
## fertilizer 4 1079 269.71 5.334 0.000426 ***
## Residuals 195 9860 50.56
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretación: Parece que nuestro modelo es significativo pues tiene un p-value de O.0004 (¡Yeah!), es decir, al menos una media es diferente, lo que sigue es hacer una prueba POST-HOC para evaluar cual de todas es la diferente y si existen grupos. Aquí la cosa se va a poner un poco compleja pues voy a usar algunos paquetes extra, si tiene dudas escribeme. Estoy casi seguro que le vas a entender eres bastante avispada.
Pruebas POST-HOC
require(multcomp) #este paquete que cargué tiene una herramientas para evaluar cuál tratamiento es distinto y formar grupos para saber cuales son iguales.## Loading required package: multcomp## Warning: package 'multcomp' was built under R version 3.4.4## Loading required package: mvtnorm## Loading required package: survival## Loading required package: TH.data## Warning: package 'TH.data' was built under R version 3.4.4## Loading required package: MASS## Warning: package 'MASS' was built under R version 3.4.4##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geysergrupos<- glht(modelo, linfct= mcp(fertilizer="Tukey")) #Generé un modelo lineal que comparaba un tratamiento contra otro
summary(grupos)##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = yield ~ fertilizer, data = data.agricola)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## 2 - 1 == 0 3.6227 1.5900 2.278 0.15629
## 3 - 1 == 0 0.4906 1.5900 0.309 0.99802
## 4 - 1 == 0 4.9228 1.5900 3.096 0.01899 *
## 5 - 1 == 0 -1.2383 1.5900 -0.779 0.93647
## 3 - 2 == 0 -3.1321 1.5900 -1.970 0.28469
## 4 - 2 == 0 1.3001 1.5900 0.818 0.92495
## 5 - 2 == 0 -4.8610 1.5900 -3.057 0.02117 *
## 4 - 3 == 0 4.4322 1.5900 2.788 0.04560 *
## 5 - 3 == 0 -1.7290 1.5900 -1.087 0.81288
## 5 - 4 == 0 -6.1611 1.5900 -3.875 0.00135 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)Interpretación: En el resultado vemos que al menos 4 medias son significativamente diferentes.
Ahora vamos a generar un boxplot bonito para representar grupos. Tratamiento con la misma letra serán iguales y aquellos con letras distintas serán significativamente diferentes.
Boxplot de grupos
tuk.cld <- cld(grupos) #esta función crea letras para representar los grupos en comparaciones pareadas como Tukey
old.par <- par(mai=c(1,1,1.25,1), no.readonly = TRUE)
plot(tuk.cld)
#Estos últimos dos comandos solo son para darle un margen más grande a la gráfica y se vean todas la letras
par(old.par)Conclusión
Aminta, como podrás ver realizamos todos los pasos que comprenden una ANOVA, si boservamos el último gráfico generado se nos muestran letras para saber que media es distinta a las demás. Así podemos saber que entre el tratamiento con letra “a” no existen diferencias, pero el “a” y “b” son distintos entre si. Por último si un tratamiento tiene dos letras quiere decir que pertenece a dos grupos.
Si tienes alguna duda no tengas miramientos en escribirme.