En este ejercicio programaremos y compararemos las pruebas de igualdad de varianza de Bartlett y la de Levne. La hipótesis para ambas pruebas es:
\[ H_0: \forall (i,j)\;|\;\sigma_i^2=\sigma_j^2 \]
\[ H_1:\exists (i,j) \;|\; \sigma_i^2\neq\sigma_j^2 \]
Para comparar las pruebas utilizaremos los siguientes datos:
val <- c(3, 1, 5, 4, 2, 6, 8, 7, 4, 5, 12, 6, 9, 3, 15, 20, 14, 11, 17, 8)
fac <- c(rep("A", 5), rep("B", 5), rep("C", 5), rep("D", 5))
La prueba de Bartlett utiliza el siguiente estadístico de prueba:
\[ T=\frac {(N-k)\ln s^2_p-\sum_{i=1}^k(N_i-1)\ln^2_i} {1+(1/(3(k-1)))((\sum_{i=1}^k 1/(N_i-1))-1/(N-k))} \]
Donde \( T\sim\chi^2_{k-1} \), utilizando esta ecuación se diseñó la siguiente función:
# La función requiere dos vectores, el primero conteniendo los valores de
# las observaciones, y el segundo los tratamientos (factores)
prueba.bartlett <- function(valores, factores) {
factores <- as.factor(factores)
N <- length(valores)
k <- length(levels(factores))
ni <- rep(NA, k)
for (i in 1:k) {
ni[i] <- sum(factores == levels(factores)[i])
}
s_i <- rep(NA, k)
for (i in 1:k) {
s_i[i] <- var(valores[factores == levels(factores)[i]])
}
sp <- sum((ni - 1) * s_i)/(N - k)
numerador <- (N - k) * log(sp) - sum((ni - 1) * log(s_i))
denominador <- 1 + (1/(3 * (k - 1))) * (sum(1/(ni - 1)) - 1/(sum(ni - 1)))
x2 <- numerador/denominador
df <- k - 1
p.val <- pchisq(x2, df, lower.tail = 0)
matrix(c(x2, df, p.val), nrow = 1, dimnames = list(c(""), c("ji-cuadrada",
"Grados de libertad", "p-valor")))
}
Utilizando la función diseñada a los datos de prueba obtenemos:
prueba.bartlett(val, fac)
## ji-cuadrada Grados de libertad p-valor
## 7.402 3 0.06013
Mientras que si utilizamos la función ya contenida en R, obtenemos los mismos resultados:
bartlett.test(val, fac)
##
## Bartlett test of homogeneity of variances
##
## data: val and fac
## Bartlett's K-squared = 7.402, df = 3, p-value = 0.06013
Por cierto, para los datos de prueba no rechazamos la hipótesis nula, con un 95% de confianza
La prueba de Levene utiliza el valor absoluto de los datos centrados dentro de cada grupo, dado el factor de la observación:
\[ z_{ij}=|y_{ij}-\bar{y}_i| \]
Y el estadístico de prueba es:
\[ L=\frac {\sum_{i=1}^kn_i(\bar{z}_{i.}-\bar{z}_{..})^2/(k-1)} {\sum_{i=1}^k\sum_{j=i}^{n_i}n_i(z_{ij}-\bar{z}_{i.})^2/(N-k)} \]
Y \( L\sim F_{k-1,n_i-k} \) ; pero este estadistico de prueba es en realidad la F de una tabla de análisis de varianza con \( z_{ij} \). Tomando esto en cuenta se diseño la función, aprovechando esta característica:
# Requiere como la Bartlett, un vector de valores y uno de factores (grupos)
prueba.levene <- function(valores, factores) {
factores <- as.factor(factores)
z <- rep(NA, length(valores))
for (i in 1:length(valores)) {
z[i] <- abs(valores[i] - mean(valores[factores == factores[i]]))
}
summary(aov(z ~ factores))
}
Utilizando los datos de prueba obtenemos:
prueba.levene(val, fac)
## Df Sum Sq Mean Sq F value Pr(>F)
## factores 3 28.8 9.6 2.74 0.077 .
## Residuals 16 56.0 3.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Si comparamos los resultados con los de la prueba incluida en el paquete “car”, obtenemos los mismos resultados:
library(car)
## Warning: package 'car' was built under R version 3.0.3
leveneTest(val, fac)
## Warning: fac coerced to factor.
## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.74 0.077 .
## 16
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Comparativamente entre las pruebas de Bartlett y Levene, llegamos para el caso de nuestros datos de prueba, a la misma conclusión. Sin embargo tomando en cuenta que existen diferencias entre los p-valores de las pruebas, concluimos que no siempre se llega a la misma conclusión mediante ambas pruebas. Podríamos simular algunos datos, pero no es el caso para este ejercicio.