redondear = options(digits=3) 

read.table("plomo.txt", h=T)->plomo
head(plomo)

##   Plomo Laboratorio
## 1   2.3        LabA
## 2   4.1        LabA
## 3   4.9        LabA
## 4   2.5        LabA
## 5   3.1        LabA
## 6   3.7        LabA

aov (y~Factor)

o

aov (y~Tratamiento)

str(plomo)

## 'data.frame':    28 obs. of  2 variables:
##  $ Plomo      : num  2.3 4.1 4.9 2.5 3.1 3.7 6.5 4 4.2 6.3 ...
##  $ Laboratorio: Factor w/ 5 levels "LabA","LabB",..: 1 1 1 1 1 1 2 2 2 2 ...

Modelo del ANOVA

model1<-aov(Plomo~Laboratorio, data=plomo)
summary(model1)

##             Df Sum Sq Mean Sq F value  Pr(>F)    
## Laboratorio  4   53.1   13.28    10.3 6.2e-05 ***
## Residuals   23   29.6    1.29                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

library(ggplot2)
ggplot(plomo, aes(Laboratorio, Plomo,colour = Laboratorio)) + 
  geom_point() + geom_boxplot()

Resumen Estadistico

attach(plomo)
tapply(Plomo, Laboratorio, mean)

## LabA LabB LabC LabD LabE 
## 3.43 5.08 2.83 3.74 7.07

tapply(Plomo, Laboratorio, sd)

##  LabA  LabB  LabC  LabD  LabE 
## 0.993 1.215 1.098 1.003 1.497

tapply(Plomo, Laboratorio, length)

## LabA LabB LabC LabD LabE 
##    6    5    6    7    4

Supuestos (Normalidad)

lm(Plomo~Laboratorio, data=plomo)->lmplomo
shapiro.test(residuals(lmplomo))

## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(lmplomo)
## W = 0.9, p-value = 0.02

#No cumple con la normalidad de los residuos

Supuestos (Heterocedasticidad)

library(car)
ncvTest(lmplomo)

## Non-constant Variance Score Test 
## Variance formula: ~ fitted.values 
## Chisquare = 0.816    Df = 1     p = 0.366

par(mfrow=c(2,2))
plot(lmplomo)

Tenemos asimetria y homoedasticidad

Aplicaremos el ANOVA no parametrico

"Kruskall-Wallis"

Test No Parametrico

kruskal.test(Plomo~Laboratorio, data=plomo)

## 
##  Kruskal-Wallis rank sum test
## 
## data:  Plomo by Laboratorio
## Kruskal-Wallis chi-squared = 10, df = 4, p-value = 0.005

#Tenemos una resultado significativo
#Me demuestra que al menos una de las medias es diferente,
#pero no me dice cual?

pairwise.wilcox.test(plomo$Plomo,plomo$Laboratorio, 
                     p.adjust.method = "none", exact=F)

## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  plomo$Plomo and plomo$Laboratorio 
## 
##      LabA LabB LabC LabD
## LabB 0.06 -    -    -   
## LabC 0.47 0.02 -    -   
## LabD 0.67 0.12 0.15 -   
## LabE 0.01 0.18 0.01 0.01
## 
## P value adjustment method: none

El Kruskal Wallis test revela diferencias significativas entre los grupos (\(x^{2}\)=10; gl=4; P<0.05)

La prueba a post-hoc usando Mann-Whitney tests sin correccion muestra resultados significantes entre los grupo A y E (p < 0.05), B y C (p < 0.05), C y E (p < 0.05), y D y E (p < 0.05).

Post hoc

Util en situaciones donde haya realmente diferencias entre grupos, las pruebas post-hoc no lo detectan. Por lo que se reconocen como los metodos mas seguros, dado que detectan las verdaderas diferencias significativas.

Pone a prueba si las muestras proceden de distribuciones con la misma mediana. No asume la normalidad, sino como una prueba de la igualdad de las medianas.

Existen multiples tecnicas ("holm", "hochberg", "hommel", "bonferroni", "BH", "BY", "fdr", "none") y lo que intenta es reducir posibilidad de cometer errores tipo I (no acepta la hipotesis nula, siendo esta verdadera).

El efecto de esta correccion es controlar la probabilidad de cometer el error tipo I.

Bonferroni son extremadamente conservadores y reducen en gran medida el poder de la prueba, con el fin de controlar el numero de falsos positivos. Se requiere de muestras grandes >6.(esto es arbitrario, algunos autores mencionan >20)

pairwise.wilcox.test(plomo$Plomo,plomo$Laboratorio, 
                     p.adjust.method = "b", exact=F)

## 
##  Pairwise comparisons using Wilcoxon rank sum test 
## 
## data:  plomo$Plomo and plomo$Laboratorio 
## 
##      LabA LabB LabC LabD
## LabB 0.6  -    -    -   
## LabC 1.0  0.2  -    -   
## LabD 1.0  1.0  1.0  -   
## LabE 0.1  1.0  0.1  0.1 
## 
## P value adjustment method: bonferroni

Si nuestros residuos son parametrico (normalidad y homocedasticidad)

TukeyHSD(model1)

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Plomo ~ Laboratorio, data = plomo)
## 
## $Laboratorio
##            diff    lwr    upr p adj
## LabB-LabA  1.65 -0.385  3.679 0.152
## LabC-LabA -0.60 -2.537  1.337 0.888
## LabD-LabA  0.31 -1.557  2.176 0.988
## LabE-LabA  3.64  1.476  5.808 0.000
## LabC-LabB -2.25 -4.279 -0.215 0.025
## LabD-LabB -1.34 -3.302  0.628 0.292
## LabE-LabB  1.99 -0.256  4.246 0.100
## LabD-LabC  0.91 -0.957  2.776 0.609
## LabE-LabC  4.24  2.076  6.408 0.000
## LabE-LabD  3.33  1.229  5.435 0.001

plot(TukeyHSD(model1), cex.axis=.6)