Analisis de Varianza 2
Alfredo Cárdenas-Rivera
22/10/2020
CLASE 12
ANÁLISIS DE VARIANZA
dos vías
1. OBJETIVO
El objetivo de este capítulo es evaluar la igualdad de medias entre 3 o más grupos independientes
================================================================================
2. PUNTOS CLAVE
- La prueba de ANOVA nos permite evaluar la diferencia estadística entre medias de tres a más grupos.
| Parámetro | Muestra | Población |
|---|---|---|
| Proporción | \(\hat{p}\) | \(p\) |
| Media | \(\bar{x}\) | \(\mu\) |
| Desviación estandar | \(s\) | \(\sigma\) |
| Varianza | \(s^2\) | \(\sigma^2\) |
| \(H_0\) es verdadera | \(H_0\) es falsa | |
|---|---|---|
| Rechaza \(H_0\) | ERROR TIPO I (\(\alpha\)) | Decisión correcta |
| Acepta \(H_0\) | Decisión correcta | ERROR TIPO II (\(\beta\)) |
================================================================================
3. TEMARIO
- ANOVA de un factor
- ANOVA de dos factores
================================================================================
4. ANOVA DE DOS FACTORRES
La intención de ese análisis es evaluar la interacción de dos factores en los grupos muestrales. Para demostrar este análisis vamos a utilizar la base de datos de crecimiento bacteriano. Primero seleccionaremos la muestra que vamos a analizar: Lecturas de Chloridium en T:3. Para ello usaremos la función subset()
growth <- read.delim("C:/Users/acard/CloudStation/R-project/bioestadistica_2020/chp4/data/curva de crecimiento final.txt", stringsAsFactors=TRUE)
db <- subset(x = growth, subset = tiempo1 == "T:3" & sp == "Chloridium", select = c(tiempo2, Medio.BBM, Medio.BG.11, Medio.F.2))
head(db)## tiempo2 Medio.BBM Medio.BG.11 Medio.F.2
## 4 12h 2585.889 1912.556 1965.889
## 11 12h 2997.000 5834.333 3595.889
## 18 12h 576.100 1681.444 9908.889
## 25 24h 4505.111 2479.222 3489.222
## 32 24h 6618.444 5020.667 8435.000
## 39 24h 2853.667 3151.444 11729.444Una vez generada la base de datos, reordenareemos esta para que podamos hacer el análisis de ANOVA de dos vías, con la función gather()
library(tidyr)
db <- tidyr::gather(data = db, key = "medio", value = "lectura", -tiempo2)
db## tiempo2 medio lectura
## 1 12h Medio.BBM 2585.889
## 2 12h Medio.BBM 2997.000
## 3 12h Medio.BBM 576.100
## 4 24h Medio.BBM 4505.111
## 5 24h Medio.BBM 6618.444
## 6 24h Medio.BBM 2853.667
## 7 12h Medio.BG.11 1912.556
## 8 12h Medio.BG.11 5834.333
## 9 12h Medio.BG.11 1681.444
## 10 24h Medio.BG.11 2479.222
## 11 24h Medio.BG.11 5020.667
## 12 24h Medio.BG.11 3151.444
## 13 12h Medio.F.2 1965.889
## 14 12h Medio.F.2 3595.889
## 15 12h Medio.F.2 9908.889
## 16 24h Medio.F.2 3489.222
## 17 24h Medio.F.2 8435.000
## 18 24h Medio.F.2 11729.444La base de datos así ordenada, tiene 3 columnas, las dos primeras corresponden a factores (“tiempo2”, “medio”) y la última corresponde a la variable numérica continua independiente (“lectura”). Ahora podemos realizar el análisis de dos vías.
res <- aov(formula = lectura ~ tiempo2*medio, data = db)
summary.aov(res)## Df Sum Sq Mean Sq F value Pr(>F)
## tiempo2 1 16481898 16481898 2.094 0.173
## medio 2 40180629 20090314 2.553 0.119
## tiempo2:medio 2 5115097 2557548 0.325 0.729
## Residuals 12 94446174 7870515print("---------------")## [1] "---------------"summary(res)## Df Sum Sq Mean Sq F value Pr(>F)
## tiempo2 1 16481898 16481898 2.094 0.173
## medio 2 40180629 20090314 2.553 0.119
## tiempo2:medio 2 5115097 2557548 0.325 0.729
## Residuals 12 94446174 7870515El valor de F para cada factor se obtine diviendo el valor de la media de la suma de cuadrados del factor entre el del residual. Por ejemplo, si queremos calcular el valor de F para el factor “medio” tendríamos que dividir 20090314 entre 7111519
20090314/7111519## [1] 2.8250384.1 Prueba de hipótesis
En este análisis la hipotesis nula se defifne como No hay efecto entre los factores “tiempo2” y “medio”. La hipótesis alternantica quedaría definida como: “Existe una interacción o relación entre los factores”
table(db$medio, db$tiempo2)##
## 12h 24h
## Medio.BBM 3 3
## Medio.BG.11 3 3
## Medio.F.2 3 3#plot(res)Antes de continuar con la interepretración de los resultados, revisemos los requisitos que se necesitan para realizar este análisis
4.2 Requisitos
- Normalidad: las muestras provienen de poblaciones con distribución normal
- Variación: las muestras presenta desviacicones similiares
- Muestreo: las muestras son aletorias cuantitativas continuas independientes
- Bifactorial: los valores se clasifican de dos formas independientes
- Diseño equilibrado: todas las muestras tienen la misma n
4.3 Análisis visual de los datos
Para una análisis visual de los datos vamos utilizar las funciones definidas en el paquete “ggplot2”, para ello primero debemos cargar la librería:
library(ggplot2)Ahora proyectaremos nuestros datos en una gráfica de caja y bigotes divididos por las dos variables
ggplot(data = db, aes(x=medio, y=lectura)) +
geom_boxplot(aes(fill=tiempo2))
Si deseamos dividdir cada uno de los grupos, lo podemos hacer con la siguiente línea de código
p <- ggplot(data = db, aes(x = medio, y = lectura)) +
geom_boxplot(aes(fill = tiempo2))
p + facet_wrap( ~ tiempo2, scales="free")
Si deseamos agregar los puntos correspondientes a cada dato, y agregar títulos, lo podemos hacer de la siguiente manera
p <- ggplot(data = db, aes(x = medio, y = lectura)) + geom_boxplot(aes(fill = tiempo2))
p <- p + geom_point(aes(y=lectura, group=tiempo2), position = position_dodge(width=0.75))
p <- p + facet_wrap( ~ tiempo2, scales="free")
p <- p + xlab("Medio") + ylab("Lectura") + ggtitle("CRECIMINTO BACTERIANO")
p <- p + guides(fill=guide_legend(title="Ciclo de luz"))
p 
4.4 Evaluación gráfica del modelo de ANOVA
Para la evaluación gráfica de nuestro modelo de ANOVA lo que tenemos que hacer es imprimir el resultado del análisis con la función ""
#layout(matrix(c(1,2,3,4),2,2)) # optional layout
#plot(res) # diagnostic plots4.5 Evaluación de la desviación estandar
library(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, uniondb %>%
group_by(tiempo2) %>%
summarise(
count_tiempo2 = n(),
mean_time = mean(lectura, na.rm = TRUE),
sd_time = sd(lectura, na.rm = TRUE)
)## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 2 x 4
## tiempo2 count_tiempo2 mean_time sd_time
## <fct> <int> <dbl> <dbl>
## 1 12h 9 3451. 2836.
## 2 24h 9 5365. 3070.db %>%
group_by(medio) %>%
summarise(
count = n(),
mean = mean(lectura, na.rm = TRUE),
sd = sd(lectura, na.rm = TRUE)
)## `summarise()` ungrouping output (override with `.groups` argument)## # A tibble: 3 x 4
## medio count mean sd
## <chr> <int> <dbl> <dbl>
## 1 Medio.BBM 6 3356. 2033.
## 2 Medio.BG.11 6 3347. 1709.
## 3 Medio.F.2 6 6521. 4019.4.6 Análisis de homogenicidad
Para analizar la homogenicidad de las varianzas podemos aplicar la prueba de levene, sin embargo este es poco sensible para distribuciones normales. Esta prueba está definida en el paquete “car”, y la función es leveneTest()
library(car)## Loading required package: carData##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodecar::leveneTest(lectura~tiempo2, data = db)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 1 0.1459 0.7075
## 16car::leveneTest(lectura~medio, data = db)## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 2 5.5728 0.0155 *
## 15
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Otra alternativa para hacer este análisis es el análisis de medias que asume que las varianzas no son iguale. Esta función está definida en oneway.test()
oneway.test(lectura~tiempo2, data = db)##
## One-way analysis of means (not assuming equal variances)
##
## data: lectura and tiempo2
## F = 1.8871, num df = 1.0, denom df = 15.9, p-value = 0.1886oneway.test(lectura~medio, data = db)##
## One-way analysis of means (not assuming equal variances)
##
## data: lectura and medio
## F = 1.5751, num df = 2.0000, denom df = 9.3241, p-value = 0.25744.7 Analisis de la normalidad
La normalidad de la distribución de los datos lo podemos graficar con la función plot(), definido en el segundo gráfico, que corresponde a un gráfica de normalidad de QQ.
plot(res, 2)
Otra forma para evaluar graficamente la distribución de los residuales está defindo en el paquete “car” con la función residualPlot()
car::residualPlot(res)
Para evaluar la normalidad de la distribución de los residuales de esta gráfica, podemos hacer un test de Shapiro a los residuales
residuos <- residuals(object = res)
cbind(residuos)## residuos
## 1 532.8927
## 2 944.0037
## 3 -1476.8963
## 4 -153.9630
## 5 1959.3700
## 6 -1805.4070
## 7 -1230.2217
## 8 2691.5553
## 9 -1461.3337
## 10 -1071.2223
## 11 1470.2227
## 12 -399.0003
## 13 -3191.0000
## 14 -1561.0000
## 15 4752.0000
## 16 -4395.3333
## 17 550.4447
## 18 3844.8887shapiro.test(x = residuos )##
## Shapiro-Wilk normality test
##
## data: residuos
## W = 0.97749, p-value = 0.9205Análisis pareado de ANOVA (poshoc)
Para hacer este análisis podemos usar la función TukeyHSD()
TukeyHSD(res)## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = lectura ~ tiempo2 * medio, data = db)
##
## $tiempo2
## diff lwr upr p adj
## 24h-12h 1913.804 -967.6739 4795.281 0.17348
##
## $medio
## diff lwr upr p adj
## Medio.BG.11-Medio.BBM -9.424167 -4330.626 4311.778 0.9999813
## Medio.F.2-Medio.BBM 3164.687000 -1156.515 7485.889 0.1663375
## Medio.F.2-Medio.BG.11 3174.111167 -1147.091 7495.313 0.1648147
##
## $`tiempo2:medio`
## diff lwr upr p adj
## 24h:Medio.BBM-12h:Medio.BBM 2606.0777 -5087.985 10300.140 0.8566268
## 12h:Medio.BG.11-12h:Medio.BBM 1089.7813 -6604.281 8783.844 0.9961971
## 24h:Medio.BG.11-12h:Medio.BBM 1497.4480 -6196.614 9191.510 0.9839404
## 12h:Medio.F.2-12h:Medio.BBM 3103.8927 -4590.170 10797.955 0.7510724
## 24h:Medio.F.2-12h:Medio.BBM 5831.5590 -1862.503 13525.621 0.1852720
## 12h:Medio.BG.11-24h:Medio.BBM -1516.2963 -9210.359 6177.766 0.9830373
## 24h:Medio.BG.11-24h:Medio.BBM -1108.6297 -8802.692 6585.433 0.9958812
## 12h:Medio.F.2-24h:Medio.BBM 497.8150 -7196.247 8191.877 0.9999132
## 24h:Medio.F.2-24h:Medio.BBM 3225.4813 -4468.581 10919.544 0.7222704
## 24h:Medio.BG.11-12h:Medio.BG.11 407.6667 -7286.396 8101.729 0.9999677
## 12h:Medio.F.2-12h:Medio.BG.11 2014.1113 -5679.951 9708.174 0.9445305
## 24h:Medio.F.2-12h:Medio.BG.11 4741.7777 -2952.285 12435.840 0.3619722
## 12h:Medio.F.2-24h:Medio.BG.11 1606.4447 -6087.618 9300.507 0.9782153
## 24h:Medio.F.2-24h:Medio.BG.11 4334.1110 -3359.951 12028.173 0.4504479
## 24h:Medio.F.2-12h:Medio.F.2 2727.6663 -4966.396 10421.729 0.8330973summary(res) # display Type I ANOVA table## Df Sum Sq Mean Sq F value Pr(>F)
## tiempo2 1 16481898 16481898 2.094 0.173
## medio 2 40180629 20090314 2.553 0.119
## tiempo2:medio 2 5115097 2557548 0.325 0.729
## Residuals 12 94446174 7870515drop1(res,~.,test="F") # type III SS and F Tests## Single term deletions
##
## Model:
## lectura ~ tiempo2 * medio
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 94446174 290.52
## tiempo2 1 10187461 104633635 290.36 1.2944 0.2774
## medio 2 14878418 109324592 289.15 0.9452 0.4157
## tiempo2:medio 2 5115097 99561271 287.47 0.3250 0.7287