Ejercicio homogeneidad de varianzas
2023-08-07
INTRODUCCIÓN.
El supuesto de homogeneidad de varianzas, también conocido como supuesto de homocedasticidad, considera que la varianza es constante (no varía) en los diferentes niveles de un factor, es decir, entre diferentes grupos.
A la hora de realizar contrastes de hipótesis o intervalos de confianza, cuando los tamaños de cada grupo son muy distintos ocurre que:
Si los grupos con tamaños muestrales pequeños son los que tienen mayor varianza, la probabilidad real de cometer un error de tipo I en los contrastes de hipótesis será menor de lo que se obtiene al hacer el test. En los intervalos, los límites superior e inferior reales son menores que los que se obtienen. La inferencia será por lo general más conservadora.
Si por el contrario, son los grupos con tamaños muestrales grandes los que tienen mayor varianza, entonces se tendrá el efecto contrario y las pruebas serán más liberales. Es decir, la probabilidad real de cometer un error de tipo I es mayor que la devuelta por el test y los intervalos de confianza verdaderos serán más amplios que los calculados.
Existen diferentes test que permiten evaluar la distribución de la varianza. Todos ellos consideran como hipótesis nula que la varianza es igual entre los grupos y como hipótesis alternativa que no lo es. La diferencia entre ellos es el estadístico de centralidad que utilizan:
Los test que trabajan con la media de la varianza son los más potentes cuando las poblaciones que se comparan se distribuyen de forma normal.
Utilizar la media truncada mejora el test cuando los datos siguen una distribución de Cauchy (colas grandes).
La mediana consigue mejorarlo cuando los datos siguen una distribución asimétrica.
# Librerías necesarias
#===================================
library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.2.3library(dplyr)## Warning: package 'dplyr' was built under R version 4.2.3##
## 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("iris")
iris<-filter(.data = iris, Species %in% c("versicolor", "virginica"))
iris## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 7.0 3.2 4.7 1.4 versicolor
## 2 6.4 3.2 4.5 1.5 versicolor
## 3 6.9 3.1 4.9 1.5 versicolor
## 4 5.5 2.3 4.0 1.3 versicolor
## 5 6.5 2.8 4.6 1.5 versicolor
## 6 5.7 2.8 4.5 1.3 versicolor
## 7 6.3 3.3 4.7 1.6 versicolor
## 8 4.9 2.4 3.3 1.0 versicolor
## 9 6.6 2.9 4.6 1.3 versicolor
## 10 5.2 2.7 3.9 1.4 versicolor
## 11 5.0 2.0 3.5 1.0 versicolor
## 12 5.9 3.0 4.2 1.5 versicolor
## 13 6.0 2.2 4.0 1.0 versicolor
## 14 6.1 2.9 4.7 1.4 versicolor
## 15 5.6 2.9 3.6 1.3 versicolor
## 16 6.7 3.1 4.4 1.4 versicolor
## 17 5.6 3.0 4.5 1.5 versicolor
## 18 5.8 2.7 4.1 1.0 versicolor
## 19 6.2 2.2 4.5 1.5 versicolor
## 20 5.6 2.5 3.9 1.1 versicolor
## 21 5.9 3.2 4.8 1.8 versicolor
## 22 6.1 2.8 4.0 1.3 versicolor
## 23 6.3 2.5 4.9 1.5 versicolor
## 24 6.1 2.8 4.7 1.2 versicolor
## 25 6.4 2.9 4.3 1.3 versicolor
## 26 6.6 3.0 4.4 1.4 versicolor
## 27 6.8 2.8 4.8 1.4 versicolor
## 28 6.7 3.0 5.0 1.7 versicolor
## 29 6.0 2.9 4.5 1.5 versicolor
## 30 5.7 2.6 3.5 1.0 versicolor
## 31 5.5 2.4 3.8 1.1 versicolor
## 32 5.5 2.4 3.7 1.0 versicolor
## 33 5.8 2.7 3.9 1.2 versicolor
## 34 6.0 2.7 5.1 1.6 versicolor
## 35 5.4 3.0 4.5 1.5 versicolor
## 36 6.0 3.4 4.5 1.6 versicolor
## 37 6.7 3.1 4.7 1.5 versicolor
## 38 6.3 2.3 4.4 1.3 versicolor
## 39 5.6 3.0 4.1 1.3 versicolor
## 40 5.5 2.5 4.0 1.3 versicolor
## 41 5.5 2.6 4.4 1.2 versicolor
## 42 6.1 3.0 4.6 1.4 versicolor
## 43 5.8 2.6 4.0 1.2 versicolor
## 44 5.0 2.3 3.3 1.0 versicolor
## 45 5.6 2.7 4.2 1.3 versicolor
## 46 5.7 3.0 4.2 1.2 versicolor
## 47 5.7 2.9 4.2 1.3 versicolor
## 48 6.2 2.9 4.3 1.3 versicolor
## 49 5.1 2.5 3.0 1.1 versicolor
## 50 5.7 2.8 4.1 1.3 versicolor
## 51 6.3 3.3 6.0 2.5 virginica
## 52 5.8 2.7 5.1 1.9 virginica
## 53 7.1 3.0 5.9 2.1 virginica
## 54 6.3 2.9 5.6 1.8 virginica
## 55 6.5 3.0 5.8 2.2 virginica
## 56 7.6 3.0 6.6 2.1 virginica
## 57 4.9 2.5 4.5 1.7 virginica
## 58 7.3 2.9 6.3 1.8 virginica
## 59 6.7 2.5 5.8 1.8 virginica
## 60 7.2 3.6 6.1 2.5 virginica
## 61 6.5 3.2 5.1 2.0 virginica
## 62 6.4 2.7 5.3 1.9 virginica
## 63 6.8 3.0 5.5 2.1 virginica
## 64 5.7 2.5 5.0 2.0 virginica
## 65 5.8 2.8 5.1 2.4 virginica
## 66 6.4 3.2 5.3 2.3 virginica
## 67 6.5 3.0 5.5 1.8 virginica
## 68 7.7 3.8 6.7 2.2 virginica
## 69 7.7 2.6 6.9 2.3 virginica
## 70 6.0 2.2 5.0 1.5 virginica
## 71 6.9 3.2 5.7 2.3 virginica
## 72 5.6 2.8 4.9 2.0 virginica
## 73 7.7 2.8 6.7 2.0 virginica
## 74 6.3 2.7 4.9 1.8 virginica
## 75 6.7 3.3 5.7 2.1 virginica
## 76 7.2 3.2 6.0 1.8 virginica
## 77 6.2 2.8 4.8 1.8 virginica
## 78 6.1 3.0 4.9 1.8 virginica
## 79 6.4 2.8 5.6 2.1 virginica
## 80 7.2 3.0 5.8 1.6 virginica
## 81 7.4 2.8 6.1 1.9 virginica
## 82 7.9 3.8 6.4 2.0 virginica
## 83 6.4 2.8 5.6 2.2 virginica
## 84 6.3 2.8 5.1 1.5 virginica
## 85 6.1 2.6 5.6 1.4 virginica
## 86 7.7 3.0 6.1 2.3 virginica
## 87 6.3 3.4 5.6 2.4 virginica
## 88 6.4 3.1 5.5 1.8 virginica
## 89 6.0 3.0 4.8 1.8 virginica
## 90 6.9 3.1 5.4 2.1 virginica
## 91 6.7 3.1 5.6 2.4 virginica
## 92 6.9 3.1 5.1 2.3 virginica
## 93 5.8 2.7 5.1 1.9 virginica
## 94 6.8 3.2 5.9 2.3 virginica
## 95 6.7 3.3 5.7 2.5 virginica
## 96 6.7 3.0 5.2 2.3 virginica
## 97 6.3 2.5 5.0 1.9 virginica
## 98 6.5 3.0 5.2 2.0 virginica
## 99 6.2 3.4 5.4 2.3 virginica
## 100 5.9 3.0 5.1 1.8 virginicaggplot(data=iris, aes(x=Species, y= Sepal.Length, colour = Species)) +
geom_boxplot() +
geom_point() +
theme_bw() +
theme(legend.position = "none")
aggregate(Sepal.Length~Species, data=iris, FUN=var)## Species Sepal.Length
## 1 versicolor 0.2664327
## 2 virginica 0.4043429Test de Bartlett de homogeneidad de las varianzas
var.test(x= iris[iris$Species == "versicolor", "Sepal.Length"],
y = iris[iris$Species == "virginica", "Sepal.Length"])##
## F test to compare two variances
##
## data: iris[iris$Species == "versicolor", "Sepal.Length"] and iris[iris$Species == "virginica", "Sepal.Length"]
## F = 0.65893, num df = 49, denom df = 49, p-value = 0.1478
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
## 0.3739257 1.1611546
## sample estimates:
## ratio of variances
## 0.6589276H0 = las varianzas son similares (homogeneidad u homocedasticidad) H1 = las varianzas son diferentes (heterogeneidad o heterocedasticidad)
El estadístico del contraste F 0 0.65 con p-valor = 0.14 > 0.05. Tenemos evidencia empírica suficiente para aceptar la hipótesis de homocedasticidad.
versicolor_sepal_length <- iris[iris$Species == "versicolor", "Sepal.Length"]
versicolor_sepal_length## [1] 7.0 6.4 6.9 5.5 6.5 5.7 6.3 4.9 6.6 5.2 5.0 5.9 6.0 6.1 5.6 6.7 5.6 5.8 6.2
## [20] 5.6 5.9 6.1 6.3 6.1 6.4 6.6 6.8 6.7 6.0 5.7 5.5 5.5 5.8 6.0 5.4 6.0 6.7 6.3
## [39] 5.6 5.5 5.5 6.1 5.8 5.0 5.6 5.7 5.7 6.2 5.1 5.7virginica_sepal_length <- iris[iris$Species == "virginica", "Sepal.Length"]
virginica_sepal_length## [1] 6.3 5.8 7.1 6.3 6.5 7.6 4.9 7.3 6.7 7.2 6.5 6.4 6.8 5.7 5.8 6.4 6.5 7.7 7.7
## [20] 6.0 6.9 5.6 7.7 6.3 6.7 7.2 6.2 6.1 6.4 7.2 7.4 7.9 6.4 6.3 6.1 7.7 6.3 6.4
## [39] 6.0 6.9 6.7 6.9 5.8 6.8 6.7 6.7 6.3 6.5 6.2 5.9shapiro.test(versicolor_sepal_length)##
## Shapiro-Wilk normality test
##
## data: versicolor_sepal_length
## W = 0.97784, p-value = 0.4647shapiro.test(virginica_sepal_length)##
## Shapiro-Wilk normality test
##
## data: virginica_sepal_length
## W = 0.97118, p-value = 0.2583Para la variable en sus dos categorías analizadas, tenemos evidencia empírica suficiente para aceptar la hipótesis de normalidad. El test de varianza es concluyente.
hist(versicolor_sepal_length, probability = TRUE, main = "Versicolor Sepal Length", xlab = "Length", col = "lightblue")
curve(dnorm(x, mean = mean(versicolor_sepal_length), sd = sd(versicolor_sepal_length)),
col = "red", lwd = 2, add = TRUE)
hist(virginica_sepal_length, probability = TRUE, main = "Virginica Sepal Length", xlab = "Length", col = "lightgreen")
curve(dnorm(x, mean = mean(virginica_sepal_length), sd = sd(virginica_sepal_length)),
col = "blue", lwd = 2, add = TRUE)
Test de Levene
H0 = Las varianzas son similares (homogeneidad u homocedasticidad) H1 = las varianzas son diferentes (heterogeneidad o heterocedasticidad)
iris## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 1 7.0 3.2 4.7 1.4 versicolor
## 2 6.4 3.2 4.5 1.5 versicolor
## 3 6.9 3.1 4.9 1.5 versicolor
## 4 5.5 2.3 4.0 1.3 versicolor
## 5 6.5 2.8 4.6 1.5 versicolor
## 6 5.7 2.8 4.5 1.3 versicolor
## 7 6.3 3.3 4.7 1.6 versicolor
## 8 4.9 2.4 3.3 1.0 versicolor
## 9 6.6 2.9 4.6 1.3 versicolor
## 10 5.2 2.7 3.9 1.4 versicolor
## 11 5.0 2.0 3.5 1.0 versicolor
## 12 5.9 3.0 4.2 1.5 versicolor
## 13 6.0 2.2 4.0 1.0 versicolor
## 14 6.1 2.9 4.7 1.4 versicolor
## 15 5.6 2.9 3.6 1.3 versicolor
## 16 6.7 3.1 4.4 1.4 versicolor
## 17 5.6 3.0 4.5 1.5 versicolor
## 18 5.8 2.7 4.1 1.0 versicolor
## 19 6.2 2.2 4.5 1.5 versicolor
## 20 5.6 2.5 3.9 1.1 versicolor
## 21 5.9 3.2 4.8 1.8 versicolor
## 22 6.1 2.8 4.0 1.3 versicolor
## 23 6.3 2.5 4.9 1.5 versicolor
## 24 6.1 2.8 4.7 1.2 versicolor
## 25 6.4 2.9 4.3 1.3 versicolor
## 26 6.6 3.0 4.4 1.4 versicolor
## 27 6.8 2.8 4.8 1.4 versicolor
## 28 6.7 3.0 5.0 1.7 versicolor
## 29 6.0 2.9 4.5 1.5 versicolor
## 30 5.7 2.6 3.5 1.0 versicolor
## 31 5.5 2.4 3.8 1.1 versicolor
## 32 5.5 2.4 3.7 1.0 versicolor
## 33 5.8 2.7 3.9 1.2 versicolor
## 34 6.0 2.7 5.1 1.6 versicolor
## 35 5.4 3.0 4.5 1.5 versicolor
## 36 6.0 3.4 4.5 1.6 versicolor
## 37 6.7 3.1 4.7 1.5 versicolor
## 38 6.3 2.3 4.4 1.3 versicolor
## 39 5.6 3.0 4.1 1.3 versicolor
## 40 5.5 2.5 4.0 1.3 versicolor
## 41 5.5 2.6 4.4 1.2 versicolor
## 42 6.1 3.0 4.6 1.4 versicolor
## 43 5.8 2.6 4.0 1.2 versicolor
## 44 5.0 2.3 3.3 1.0 versicolor
## 45 5.6 2.7 4.2 1.3 versicolor
## 46 5.7 3.0 4.2 1.2 versicolor
## 47 5.7 2.9 4.2 1.3 versicolor
## 48 6.2 2.9 4.3 1.3 versicolor
## 49 5.1 2.5 3.0 1.1 versicolor
## 50 5.7 2.8 4.1 1.3 versicolor
## 51 6.3 3.3 6.0 2.5 virginica
## 52 5.8 2.7 5.1 1.9 virginica
## 53 7.1 3.0 5.9 2.1 virginica
## 54 6.3 2.9 5.6 1.8 virginica
## 55 6.5 3.0 5.8 2.2 virginica
## 56 7.6 3.0 6.6 2.1 virginica
## 57 4.9 2.5 4.5 1.7 virginica
## 58 7.3 2.9 6.3 1.8 virginica
## 59 6.7 2.5 5.8 1.8 virginica
## 60 7.2 3.6 6.1 2.5 virginica
## 61 6.5 3.2 5.1 2.0 virginica
## 62 6.4 2.7 5.3 1.9 virginica
## 63 6.8 3.0 5.5 2.1 virginica
## 64 5.7 2.5 5.0 2.0 virginica
## 65 5.8 2.8 5.1 2.4 virginica
## 66 6.4 3.2 5.3 2.3 virginica
## 67 6.5 3.0 5.5 1.8 virginica
## 68 7.7 3.8 6.7 2.2 virginica
## 69 7.7 2.6 6.9 2.3 virginica
## 70 6.0 2.2 5.0 1.5 virginica
## 71 6.9 3.2 5.7 2.3 virginica
## 72 5.6 2.8 4.9 2.0 virginica
## 73 7.7 2.8 6.7 2.0 virginica
## 74 6.3 2.7 4.9 1.8 virginica
## 75 6.7 3.3 5.7 2.1 virginica
## 76 7.2 3.2 6.0 1.8 virginica
## 77 6.2 2.8 4.8 1.8 virginica
## 78 6.1 3.0 4.9 1.8 virginica
## 79 6.4 2.8 5.6 2.1 virginica
## 80 7.2 3.0 5.8 1.6 virginica
## 81 7.4 2.8 6.1 1.9 virginica
## 82 7.9 3.8 6.4 2.0 virginica
## 83 6.4 2.8 5.6 2.2 virginica
## 84 6.3 2.8 5.1 1.5 virginica
## 85 6.1 2.6 5.6 1.4 virginica
## 86 7.7 3.0 6.1 2.3 virginica
## 87 6.3 3.4 5.6 2.4 virginica
## 88 6.4 3.1 5.5 1.8 virginica
## 89 6.0 3.0 4.8 1.8 virginica
## 90 6.9 3.1 5.4 2.1 virginica
## 91 6.7 3.1 5.6 2.4 virginica
## 92 6.9 3.1 5.1 2.3 virginica
## 93 5.8 2.7 5.1 1.9 virginica
## 94 6.8 3.2 5.9 2.3 virginica
## 95 6.7 3.3 5.7 2.5 virginica
## 96 6.7 3.0 5.2 2.3 virginica
## 97 6.3 2.5 5.0 1.9 virginica
## 98 6.5 3.0 5.2 2.0 virginica
## 99 6.2 3.4 5.4 2.3 virginica
## 100 5.9 3.0 5.1 1.8 virginicalibrary(car)## Warning: package 'car' was built under R version 4.2.3## Loading required package: carData## Warning: package 'carData' was built under R version 4.2.3##
## Attaching package: 'car'## The following object is masked from 'package:dplyr':
##
## recodeleveneTest(y=iris$Sepal.Length, group = iris$Species, center = "median")## Levene's Test for Homogeneity of Variance (center = "median")
## Df F value Pr(>F)
## group 1 1.0245 0.3139
## 98El estadístico F asciende a 1.02 con p-valor de 0.31. Tenemos evidencia empírica suficiente para aceptar la hipótesis de igualdad de varianzas. Tenemos homocedasticidad.
Test de Bartlett
Para hacer el test de Bartlett previamente hay que tener filtrada la variable en sus categorías.
versicolor_sepal_length virginica_sepal_length
bartlett.test(list(versicolor_sepal_length, virginica_sepal_length))##
## Bartlett test of homogeneity of variances
##
## data: list(versicolor_sepal_length, virginica_sepal_length)
## Bartlett's K-squared = 2.0949, df = 1, p-value = 0.1478El p-valor asciende a 0.14 > 0.05. Tenemos evidencia empírica suficiente para aceptarla hipótesis de igualdad de vaianzas. Es decir, se verificca la homocedasticidad.
Test de Brown-Forsyth
library(HH)## Warning: package 'HH' was built under R version 4.2.3## Loading required package: lattice## Loading required package: grid## Loading required package: latticeExtra## Warning: package 'latticeExtra' was built under R version 4.2.3##
## Attaching package: 'latticeExtra'## The following object is masked from 'package:ggplot2':
##
## layer## Loading required package: multcomp## Warning: package 'multcomp' was built under R version 4.2.3## Loading required package: mvtnorm## Warning: package 'mvtnorm' was built under R version 4.2.3## Loading required package: survival## Loading required package: TH.data## Warning: package 'TH.data' was built under R version 4.2.3## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyser## Loading required package: gridExtra## Warning: package 'gridExtra' was built under R version 4.2.3##
## Attaching package: 'gridExtra'## The following object is masked from 'package:dplyr':
##
## combine##
## Attaching package: 'HH'## The following objects are masked from 'package:car':
##
## logit, vifhov(iris$Sepal.Length~iris$Species)##
## hov: Brown-Forsyth
##
## data: iris$Sepal.Length
## F = 1.0245, df:iris$Species = 1, df:Residuals = 98, p-value = 0.3139
## alternative hypothesis: variances are not identicalEl p-valor es 0.3139 > 0.05. Tenemos evidencia empírica suficiente para aceptar la hipótesis de homocedasticidad a la significación de 5%.
Test de Fligner-Killeen
Para hacer el test de Fligner-Killeen previamente hay que tener filtrada la variable en sus categorías.
versicolor_sepal_length virginica_sepal_length
fligner.test(versicolor_sepal_length, virginica_sepal_length)##
## Fligner-Killeen test of homogeneity of variances
##
## data: versicolor_sepal_length and virginica_sepal_length
## Fligner-Killeen:med chi-squared = 23.203, df = 20, p-value = 0.279El p-valor asociado es 0.279 > 0.05. Tenemos evidencia empírica suficiente para aceptar la homocedasticidad.
CONCLUSIÓN
Todos los tests anteriores se basan en distribuciones F de Fligner y necesitan la hipótesis de normalidad para poder efectuarse correctamente.
La batería de testS realizados confirma que tenemos en este caso homocedasticidad para la variable en cada categoría.
Esto sería el primer paso para hacer el ANOVA y hacer comparaciones de medias de Tukey.
- Normalidad — ANOVA
- Homocedasticidad — ANOVA — comparación de medias con TUKEY.
Información de la sesión
sesion_info <- devtools::session_info()
dplyr::select(
tibble::as_tibble(sesion_info$packages),
c(package, loadedversion, source)
)## # A tibble: 105 × 3
## package loadedversion source
## <chr> <chr> <chr>
## 1 abind 1.4-5 CRAN (R 4.2.0)
## 2 backports 1.4.1 CRAN (R 4.2.0)
## 3 base64enc 0.1-3 CRAN (R 4.2.0)
## 4 bslib 0.4.2 CRAN (R 4.2.2)
## 5 cachem 1.0.6 CRAN (R 4.2.2)
## 6 callr 3.7.3 CRAN (R 4.2.3)
## 7 car 3.1-2 CRAN (R 4.2.3)
## 8 carData 3.0-5 CRAN (R 4.2.3)
## 9 checkmate 2.2.0 CRAN (R 4.2.3)
## 10 cli 3.6.0 CRAN (R 4.2.2)
## # ℹ 95 more rows