Evaluación de los supuestos del análisis de varianza
Francisco Pec
30 de mayo de 2017
Los siguientes datos fueron obtenidos del estudio denominado “Evaluación del parasitismo de cepas de Metarhizium anisopliae (Metsch.) Sorokin, Beauveria bassiana (Bals.) Vuillemin y Paecilomyces lilacinus (Thom.) Samson en ninfas de chinche salivosa Aeneolamia postica (Walk.), bajo condiciones semicontroladas en la finca Belén, Santa Lucía Cotzumalguapa, Escuintla, Guatemala, C.A.”
Los tratamientos evaluados fueron los siguientes:
T1: Beauveria bassiana Bisa 0114
T2: Metarhizium anisopliae ILU 0114
T3: Metarhizium anisopliae BISA 01-2000
T4: Metarhizium anisopliae Bisa 0214
T5: Paecilomyces lilacinus BISA 0113
El experimento fue conducido bajo el diseño completamente al azar con 5 repeticiones, entre las variables medidas está el porcentaje de parasitismo.
Los datos obtenidos de mortalidad confirmada (parasitismo) de cada una de las cepas empleadas en cada unidad experimental se presentan a continuación
Importando y visualizando los datos
El objeto llamado ecepas recibirá el contenido de la base de datos cepasP.txt, con header=TRUE se indica que las variables tienen encabezado.
Para poder importar una base de datos con la función read.table es necesario cambiar el directorio. En RStudio el proceso es el siguiente: Session > Set Working Directory > Choose Directory y Seleccionar la carpeta (donde se encuentra la base de datos en formaro .txt delimitado por tabulaciones)
ecepas <-read.table('cepasP.txt', header =TRUE)
dim(ecepas) #Dimensión de la base de datos, dos variables y 25 observaciones## [1] 25 2head(ecepas) #Resumen de la base de datos## Tratamiento parasitismo
## 1 T1 4
## 2 T2 20
## 3 T3 16
## 4 T4 20
## 5 T5 36
## 6 T1 12ecepas #Para visualizar la base de datos completa## Tratamiento parasitismo
## 1 T1 4
## 2 T2 20
## 3 T3 16
## 4 T4 20
## 5 T5 36
## 6 T1 12
## 7 T2 24
## 8 T3 12
## 9 T4 28
## 10 T5 24
## 11 T1 4
## 12 T2 32
## 13 T3 12
## 14 T4 24
## 15 T5 40
## 16 T1 16
## 17 T2 16
## 18 T3 24
## 19 T4 24
## 20 T5 32
## 21 T1 0
## 22 T2 16
## 23 T3 20
## 24 T4 16
## 25 T5 44attach(ecepas) #Cargar las variables a RCargando los paquetes necesarios
require(agricolae) #Coeficiente de variación## Loading required package: agricolaerequire(laercio) #Para la prueba de Tukey## Loading required package: laerciorequire(car) #Durbin-Watson Test, Levene's Test## Loading required package: carrequire(ggplot2) #Gráficos avanzados## Loading required package: ggplot2require(outliers) #Cochran Test## Loading required package: outliersDefinición del modelo para el Análisis de la Varianza
tratamiento <- factor(Tratamiento) #Definiendo la variable tratamiento como factor
parasitismo<-as.numeric(parasitismo) #Definiendo la variable parasitismo como variable numérica
modp <- aov(parasitismo ~ tratamiento) #Modelo correspondiente a un experimento de un factor, bajo el diseño completamente al azarVerificación del cumplimiento de los supuestos del análisis de varianza para la variable porcentaje de parasitismo
La verificación de los supuestos del análisis de varianza garantiza que el modelo empleado sea adecuado para el análisis, lo que conlleva a que las conclusiones del ANDEVA estén acotadas en sus errores Tipo I y Tipo II.
Esto se realiza con los residuos del modelo, un residuo es la diferencia entre el valor real () y el valor que un modelo predice.
modp es el objeto que contiene el modelo creado anteriormente
residuos<-residuals(modp) #Creando un objeto llamado residuos que contiene los residuos el modelo
residuos #Visualizando los residuos## 1 2 3 4 5 6 7 8 9 10 11 12
## -3.2 -1.6 -0.8 -2.4 0.8 4.8 2.4 -4.8 5.6 -11.2 -3.2 10.4
## 13 14 15 16 17 18 19 20 21 22 23 24
## -4.8 1.6 4.8 8.8 -5.6 7.2 1.6 -3.2 -7.2 -5.6 3.2 -6.4
## 25
## 8.8Independencia de los residuos
Este supuesto requiere que la probabilidad de que el residuo de una observación cualquiera tenga un determinado valor, no debe depender de los valores de los otros residuos, en teoría esto se cumplió al realizar la aleatorización de los tratamientos, pero se realizó la prueba de Durbin Watson (DW) para tener certeza del cumplimiento de este supuesto, para ello se emplea la ecuación:
alt text
Donde n es el número de observaciones, ei se refiere a los residuos. Puesto que el valor de DW es aproximadamente igual a 2(1-Rho) donde Rho es la autocorrelación de la muestra de los residuos.
Si el valor del estadístico Durbin Watson (DW) está próximo a 2 entonces los residuos no están autocorrelacionados. Si su valor es 0 hay autocorrelación perfecta positiva. Si tiene un valor de 4 existe autocorrelación perfecta negativa.
Al realizar la prueba de independencia de residuos para la variable porcentaje de parasitismo se determinó que los residuos no están correlacionados, debido a que el DW está próximo a 2 y el valor de p (p-value) es superior al nivel de significancia de 5% (alfa=0.05) por lo que se concluye que existe independencia de los residuos
durbinWatsonTest(modp) ## lag Autocorrelation D-W Statistic p-value
## 1 -0.254918 2.397541 0.206
## Alternative hypothesis: rho != 0Distribución normal de los residuos
Este supuesto establece que los residuos deben estar normalmente distribuidos ei~N(0,1). A continuación, se presenta una curva de densidad observada de las frecuencias de los residuos y el gráfico de probabilidad normal, considerando que la curva de densidad observada se asemeja a la campana de Gauss, y que, en el gráfico de probabilidad normal, los cuantiles observados de los residuos de la variable porcentaje de parasitismo (representados por los puntos negros) se aproximan a la línea central que representa los cuantiles de una distribución normal teórica por lo que no existen indicios de incumplimiento del supuesto de normalidad, pero para tener la certeza del cumplimiento de este supuesto se realizó la prueba de Shapiro Wilk.
par(mfrow=c(1,2)) #Para dividir el área del gráfico en dos partes (una fila y dos columnas)
dplot<-density(residuos) #Creando un objeto llamado dplot que recibe un Density_Plot de los residuos
plot(dplot, #Graficando el objeto dplot
main="Curva de densidad observada", #Título principal de la gráfica
xlab = "Residuos", #Etiqueta del eje x
ylab = "Densidad") #Etiqueta del eje y
polygon(dplot, #Añadiendo el poligono
col = "green", #Definiendo el color del poligono
border = "black") #Color del borde del poligono
qqPlot(residuos, #Un gráfico Cuantil-Cuantil de los residuos
pch =20, #Forma de los puntos
main="QQ-Plot de los residuos", #Título principal
xlab = "Cuantiles teóricos", #Etiqueta eje x
ylab="Cuantiles observados de los residuos") #Etiqueta eje y
mtext(side=3, at=par("usr")[1], adj=0.7, cex=0.6, col="gray40", line=-21, #Posición del texto
text=paste("Francisco Pec --", #Texto
format(Sys.time(),
"%d/%m/%Y %H:%M:%S --"), #Fecha y Hora
R.version.string)) #Versión de R
Prueba de Shapiro Wilk
Se realizó una prueba Shapiro-Wilk para ratificar el cumplimiento del supuesto de normalidad de los residuos, evaluando las hipótesis:
Ho: Los residuos de la variable porcentaje de parasitismo se distribuyen normalmente con media cero y varianza constante ei~N(0,1).
Ha: los residuos de la variable porcentaje de parasitismo no siguen la distribución normal.
shapiro.test(residuals(modp)) #Prueba Shapiro-Wilk para los residuos de la variable##
## Shapiro-Wilk normality test
##
## data: residuals(modp)
## W = 0.97079, p-value = 0.6652Se acepta la hipótesis nula, debido a que el valor de p (p-value) es mayor al valor del nivel de significancia (alfa=0.05), por lo que se concluye que los residuos de la variable porcentaje de parasitismo están normalmente distribuidos con media cero y varianza constante.
Homogeneidad de varianzas.
Este es el supuesto más importante que los residuos deben cumplir para que el modelo empleado sea válido.
Con un box-plot podemos apreciar gráficamente que la dispersión de los residuos para cada tratamiento son muy distintas entre si, en el caso de esta inspección no se encontraron indicios de un incumplimiento de este supuesto.
plot(residuos ~ Tratamiento, #Box-plot de los residuos para cada tratamiento
xlab="Tratamiento",
ylab="Residuos del porcentaje de parasitismo",
col="green",
main="Boxplot para cada tratamiento evaluado")
mtext(side=3, at=par("usr")[2], adj=1,
cex=0.6, col="gray40", line=-21,
text=paste("Francisco Pec --",
format(Sys.time(), "%d/%m/%Y %H:%M:%S --"),
R.version.string))
En siguiente gráfica (plot1) se puede observar los valores predichos por el modelo para la variable porcentaje de parasitismo contra la raíz cuadrada de los residuos estandarizados, donde no se observa ninguna tendencia en la distribución de dichos valores lo cual implica en que no existen indicios de un incumpliendo del supuesto de homogeneidad de varianzas, pero para corroborar esto se realizó la prueba de C de Cochran.
plot1<-ggplot(modp, aes(.fitted, sqrt(abs(.stdresid))))+geom_point(na.rm=TRUE)
plot1<-plot1+stat_smooth(method="loess", na.rm = TRUE)+xlab("Valores predichos")
plot1<-plot1+ylab(expression(sqrt("|Residuos estandarizados|")))
plot1<-plot1+theme_bw()
plot(plot1)
La prueba de C de Cocharn, es una prueba unilateral del límite superior atípico de la varianza. La prueba C se utiliza para determinar si una sola estimación de una varianza es significativamente más grande que un grupo de varianzas, se presenta el resultado de la prueba de C de Cochran para la variable porcentaje de parasitismo.
cochran.test(parasitismo ~ tratamiento) #Esta prueba requiere que el diseño se balanceado (el número de repeticiones debe ser igual para cada tratamiento)##
## Cochran test for outlying variance
##
## data: parasitismo ~ tratamiento
## C = 0.30328, df = 5, k = 5, p-value = 0.9512
## alternative hypothesis: Group T5 has outlying variance
## sample estimates:
## T1 T2 T3 T4 T5
## 43.2 44.8 27.2 20.8 59.2Debido a que el valor de p es mayor al nivel de significancia (alfa=0.05) se acepta la hipótesis nula que plantea que las varianzas de los residuos de tratamientos son iguales, con lo anterior se determinó el cumplimiento de los supuestos del análisis de varianza para la variable porcentaje de parasitismo de ninfas de chinche salivosa, concluyendo que el modelo es adecuado para el análisis y no existe necesidad de realizar transformación de datos.
Una alternativa a la prueba ** C de Cochran está la prueba de Levene**
La prueba de Levene consiste en realizar un análisis de la varianza usando como variable dependiente el valor absoluto de los residuos.
leveneTest(parasitismo ~ tratamiento) #Esta prueba es útil aún cuando no se cumple con el supuesto de normalidad## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 4 0.2407 0.9119
## 20Debido a que el valor de p (p-value) es mayor a 0.05 se acepta la hipótesis nula, por lo que existe homogeneidad de varianzas, llegando a la misma conclusión.
Análisis de la varianza
Se realizó un análisis de varianza para determinar si al menos alguna de las cepas de hongos entomopatógenos empleadas produjo un incremento significativo sobre el porcentaje de parasitismo de ninfas de chinche salivosa
Paral ver los resultados es necesario que solicitar un resumen del modelo que se generó al início.
summary(modp)## Df Sum Sq Mean Sq F value Pr(>F)
## tratamiento 4 2057.0 514.2 13.17 2.05e-05 ***
## Residuals 20 780.8 39.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Se concluye que los tratamientos tienen un efecto diferenciado sobre el porcentaje de parasitismo en las ninfas de chinche salivosa.
Para solicitar el coeficiente de variación
cv.model(modp)## [1] 30.27229Para solicitar la prueba de Tukey
LTukey(modp, which="tratamiento", conf.level=0.95)##
## TUKEY TEST TO COMPARE MEANS
##
## Confidence level: 0.95
## Dependent variable: parasitismo
## Variation Coefficient: 30.27229 %
##
## Independent variable: tratamiento
## Factors Means
## T5 35.2 a
## T4 22.4 b
## T2 21.6 b
## T3 16.8 bc
## T1 7.2 c
##
## Se puede apreciar que el tratamiento 5 causa el porcentaje de mortalidad más alto y se diferencia significativamente de los demás tratamientos.
Finalmente un gráfico de barras
gr<- ggplot(ecepas,
aes(x=tratamiento, y=parasitismo))+ geom_bar(stat="identity")
gr<-gr+ ylab("Porcentaje de parasitismo")
gr<-gr+ xlab("Tratamientos")
gr