Anova de una vía
Junior Pastor PÉREZ MOLINA
20 de octubre de 2018
Cargo los datos desde la dirección donde esta guardados
rm(list = ls()) #Remove all objects
graphics.off() #Remove all graphics
cat("\014") #Remove script in windows console
Lechuga_hojas=read.delim("https://raw.githubusercontent.com/JPASTORPM/Database/master/Lechuga_hojas.csv",header=T,sep=";",dec=".")
str(Lechuga_hojas)## 'data.frame': 70 obs. of 3 variables:
## $ X : int 1 2 3 4 5 6 7 8 9 10 ...
## $ Tratamiento: Factor w/ 5 levels "0%","100%","25%",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ c.hojas : int 7 5 6 7 5 7 6 7 7 7 ...names(Lechuga_hojas)## [1] "X" "Tratamiento" "c.hojas"ANOVA, ¿Cuál aplicamos, Fisher ó Kruskal-Wallis?
¿Fisher ó Kruskal-Wallis? Solo lo sabremos déspues de realizar la aov() para obtener los residuos y aplicar shapiro.test() sobre los residuos y ver si se acepta o rechaza el supuesto de normalidad (p>0.05: Fisher; p<0.05: Kruskal-Wallis).
El segundo paso es determinar si las repeticiones por tratamiento (n) son balanceado (mesmo número n) o no balanceados (diferente número n}) usando tapply() Balanceado se refiere a mismo número de repeticiones (n) por tratamiento.
Vamos a realizar un ejemplo, la variable de interés es c.hojas (crecimiento del número de hojas) y el factor tratamiento (0%, 25%, 50%, 75% y 100% de concentración de abono orgánico.
c.hojas <- aov(c.hojas ~ Tratamiento, data=Lechuga_hojas)
shapiro.test(c.hojas$residuals) #Opción 1, ver p-value##
## Shapiro-Wilk normality test
##
## data: c.hojas$residuals
## W = 0.97583, p-value = 0.1939ifelse((shapiro.test(c.hojas$residuals))$p.value>0.05,"SE ACEPTA EL SUPUESTO DE NORMALIDAD", "SE RECHAZA EL SUPUESTO DE NORMALIDAD, DATOS ASIMÉTRICOS") #Opción 2, la función te dice## [1] "SE ACEPTA EL SUPUESTO DE NORMALIDAD"#¿Cuando están balanceados? Cuando tiene el mismo número de repeticiónes por tratamiento
balanceado<-tapply(Lechuga_hojas$c.hojas, Lechuga_hojas$Tratamiento, length)
balanceado## 0% 100% 25% 50% 75%
## 14 14 14 14 14# Estos datos están balanceados SE ACEPTA EL SUPUESTO DE NORMALIDAD y REPETICIONES POR TRATAMIENTO SON BALANCEADOS, debo aplicar Fisher with post-hoc Tukey HSD, de acuerdo con los siguientes códigos:
c.hojas <- aov(c.hojas ~ Tratamiento, data=Lechuga_hojas) #Notar que aqui usamos " aov() "
summary.lm(c.hojas)##
## Call:
## aov(formula = c.hojas ~ Tratamiento, data = Lechuga_hojas)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.1429 -1.1071 -0.0714 0.8571 2.8571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.1429 0.3509 17.504 < 2e-16 ***
## Tratamiento100% -3.7857 0.4963 -7.628 1.34e-10 ***
## Tratamiento25% -1.0000 0.4963 -2.015 0.04805 *
## Tratamiento50% -1.5714 0.4963 -3.166 0.00235 **
## Tratamiento75% -3.1429 0.4963 -6.333 2.57e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.313 on 65 degrees of freedom
## Multiple R-squared: 0.546, Adjusted R-squared: 0.5181
## F-statistic: 19.55 on 4 and 65 DF, p-value: 1.337e-10Comparacion<-TukeyHSD(c.hojas)
plot(Comparacion)
SE ACEPTA EL SUPUESTO DE NORMALIDAD y DATOS POR TRATAMIENTO SON NO BALANCEADOS, debo aplicar Fisher with post-hoc LSD, de acuerdo con los siguientes códigos:
c.hojas<-oneway.test (c.hojas ~ Tratamiento, data=Lechuga_hojas) #Notar que aqui usamos " oneway.test() "
c.hojas##
## One-way analysis of means (not assuming equal variances)
##
## data: c.hojas and Tratamiento
## F = 22.434, num df = 4.000, denom df = 32.312, p-value = 6.098e-09Comparacion<-pairwise.t.test (Lechuga_hojas$c.hojas, Lechuga_hojas$Tratamiento, p.adj = "none", exact= F) #Notar que usamos en p.adj= "none" para indicar que es LSD
Comparacion<-ifelse(Comparacion$p.value<0.05, "Diferente", "Igual")
Comparacion## 0% 100% 25% 50%
## 100% "Diferente" NA NA NA
## 25% "Diferente" "Diferente" NA NA
## 50% "Diferente" "Diferente" "Igual" NA
## 75% "Diferente" "Igual" "Diferente" "Diferente"SE RECHAZA EL SUPUESTO DE NORMALIDAD, debo aplicar la opción no parametrica llamada " Kruskal Wallis " y post-hoc que va depender si los datos son balanceados (Bonferroni, muestras mayores a 6), o no balanceados (LSD y/o muestras menores a 6) de acuerdo con los siguientes códigos:
c.hojas<-kruskal.test (c.hojas ~ Tratamiento, data=Lechuga_hojas)
c.hojas##
## Kruskal-Wallis rank sum test
##
## data: c.hojas by Tratamiento
## Kruskal-Wallis chi-squared = 38.984, df = 4, p-value = 7.019e-08SE RECHAZA EL SUPUESTO DE NORMALIDAD y DATOS POR TRATAMIENTO SON BALANCEADOS cuyas muestras por tratamientos son mayores que 6: Kruskal Wallis con post-hoc Bonferroni
c.hojas<-kruskal.test (c.hojas ~ Tratamiento, data=Lechuga_hojas)
c.hojas##
## Kruskal-Wallis rank sum test
##
## data: c.hojas by Tratamiento
## Kruskal-Wallis chi-squared = 38.984, df = 4, p-value = 7.019e-08Comparacion<-pairwise.wilcox.test(Lechuga_hojas$c.hojas, Lechuga_hojas$Tratamiento, p.adj = "bonferroni", exact= F)
Comparacion<-ifelse(Comparacion$p.value<0.05, "Diferente", "Igual")
Comparacion## 0% 100% 25% 50%
## 100% "Diferente" NA NA NA
## 25% "Igual" "Diferente" NA NA
## 50% "Diferente" "Diferente" "Igual" NA
## 75% "Diferente" "Igual" "Diferente" "Igual"SE RECHAZA EL SUPUESTO DE NORMALIDAD y DATOS POR TRATAMIENTO SON NO BALANCEADOS cuyas muestras PUEDEN O NO SER mayores 6: Kruskal Wallis con post-hoc LSD
c.hojas<-kruskal.test (c.hojas ~ Tratamiento, data=Lechuga_hojas)
c.hojas##
## Kruskal-Wallis rank sum test
##
## data: c.hojas by Tratamiento
## Kruskal-Wallis chi-squared = 38.984, df = 4, p-value = 7.019e-08Comparacion<-pairwise.wilcox.test (Lechuga_hojas$c.hojas, Lechuga_hojas$Tratamiento, p.adj = "none", exact= F) #Notar que usamos en p.adj= "none" para indicar que es LSD
Comparacion<-ifelse(Comparacion$p.value<0.05, "Diferente", "Igual")
Comparacion## 0% 100% 25% 50%
## 100% "Diferente" NA NA NA
## 25% "Igual" "Diferente" NA NA
## 50% "Diferente" "Diferente" "Igual" NA
## 75% "Diferente" "Igual" "Diferente" "Diferente"Ahora construimos el gráfico
library(Rmisc)## Loading required package: lattice## Loading required package: plyrerror.bar.vertical<-function(x, y, se.y, col){arrows(x, y-se.y, x, y+se.y, code=3, angle=90, length=0.25, col=col)}
function.mean<-function(data, variable){
sum = summarySE(data, measurevar= variable, groupvars=c("Tratamiento_2"), na.rm=TRUE)
sum<-sum[c(1,3,5,6)]
sum<-data.frame(sum)
names(sum)<-c("Tratamiento","Mean","S.E.","C.I.95")
sum[order(sum[,1], decreasing = TRUE),]
sum$Tratamiento_2<-c("0%","25%","50%","75%","100%")
sum
}
#--------------
Lechuga_hojas$Tratamiento_2<-as.character(Lechuga_hojas$Tratamiento)
Lechuga_hojas$Tratamiento_2[Lechuga_hojas$Tratamiento=="0%"]<-"1"
Lechuga_hojas$Tratamiento_2[Lechuga_hojas$Tratamiento=="25%"]<-"2"
Lechuga_hojas$Tratamiento_2[Lechuga_hojas$Tratamiento=="50%"]<-"3"
Lechuga_hojas$Tratamiento_2[Lechuga_hojas$Tratamiento=="75%"]<-"4"
Lechuga_hojas$Tratamiento_2[Lechuga_hojas$Tratamiento=="100%"]<-"5"
data.c.hojas<-function.mean(data=Lechuga_hojas, variable="c.hojas")
#--------------
par(mfrow=c(1,1),mgp = c(1.75,0.5,0), mar = c(3,3,1,1))
barplot(data.c.hojas$Mean, beside=T, col = c("gray"),border=c("black"),names.arg = ((as.character(data.c.hojas$Tratamiento_2))),legend.text=FALSE,ylim=c(0,7),ylab="Crecimiento número hojas (± E.E.)",xlab="Concentración de abono orgánico",main="")
error.bar.vertical(c(0.7, 1.9, 3.15, 4.3, 5.5),data.c.hojas$Mean,data.c.hojas$S.E., col="red")
text(x= c(0.7, 1.9, 3.15, 4.3, 5.5), y= data.c.hojas$Mean*1.125, labels=c("a", "b", "b", "c", "c"), col="darkred", cex=1.25)
legend("topright", expression(paste("F= 19.55; g.l.= 4, 65; p<0.05")), col="darkred", merge = F, bg = NULL, bty='n', cex=1.25)
box()
Fig. 1. Crecimiento del número de hojas de L. sativa en cuatro tratamiento de concentración de abono orgánico (0%, 25%, 50% y 100%). Letras iguales muestran no diferencias estadísticamente significativas entre los tratamientos, Tukey HSD (p>0.05).
.