Para la resolución de este problema se siguió la siguiente aproximación:
Se parte de dos premisas interpretativas: en el QQ Plot se considera que entre los puntos estén más estrechamente distirbuidos sobre la diagonal teórica, entonces tenemos una distribución normal. Por el lado de la Prueba de Shapiro-Wilks, contamos con una H0 (hipótesis nula) donde la distribución es normal junto con un nivel de significancia estandar de 0.05; así, p-valores por encima de él no nos permitirían rechazar la hipótesis mencionada.
Teniendo esto en cuenta, lo que vemos para el grupo con el 50% de la población enferma es que solo desde la muestra n = 50 empezamos a ver normalidad en su distribución. Los valores empíricos en las gráficas QQ muestran que en las muestras pequeñas estos se agrupan bajo algunos pocos valores Y debido a la limitaciones mismas de los parametros en n reducidos. Para este caso, y teniendo en cuenta que es solo un escenario dentro de los posibles de aleatoriedad (seed), solo desde la muestra n = 60 encontramos evidencias para no rechanazar la hipotesis nula y, por tanto, vemos distribuciones normales donde la diagonal de los QQ es armónica respecto a los datos teóricos y empíricos.
library(gridExtra)
library(ggplot2)
# Probabilidad para el 50-50
set.seed(12)
n_poblacion <- 1000
poblacion_50_50 <- c(rep(1, n_poblacion / 2), rep(0, n_poblacion / 2))
# Función para calcular la proporción muestral
calcular_proporcion <- function(n_muestra, poblacion) {
muestra <- sample(poblacion, n_muestra, replace = TRUE)
return(mean(muestra))
}
# Simulación con diferentes tamaños de muestra
replicaciones <- 500
tamaños_muestra <- c(5, 10, 15, 20, 30, 50, 60, 100, 200, 500)
resultados_50_50 <- list()
for (n in tamaños_muestra) {
proporciones <- replicate(replicaciones, calcular_proporcion(n, poblacion_50_50))
resultados_50_50[[as.character(n)]] <- proporciones
}
# Crear una lista para almacenar los gráficos QQ
qq_plots <- list()
# Generar los gráficos QQ y almacenarlos en la lista
for (n in tamaños_muestra) {
# Obtener todas las proporciones generadas
proporciones <- resultados_50_50[[as.character(n)]]
# Crear un gráfico QQ para las proporciones usando ggplot
p <- ggplot(data.frame(Proporción = proporciones), aes(sample = Proporción)) +
stat_qq() +
stat_qq_line() +
ggtitle(paste("QQ plot para n =", n, "(50-50)")) +
theme_minimal()
# Añadir el gráfico a la lista
qq_plots[[as.character(n)]] <- p
# Prueba de Shapiro-Wilk sobre todas las proporciones
shapiro_test <- shapiro.test(proporciones)
cat("Prueba de Shapiro-Wilk para n =", n, "y 50% enfermas | p-valor =", shapiro_test$p.value, "\n")
}
## Prueba de Shapiro-Wilk para n = 5 y 50% enfermas | p-valor = 3.015778e-14
## Prueba de Shapiro-Wilk para n = 10 y 50% enfermas | p-valor = 1.178117e-09
## Prueba de Shapiro-Wilk para n = 15 y 50% enfermas | p-valor = 2.159901e-07
## Prueba de Shapiro-Wilk para n = 20 y 50% enfermas | p-valor = 5.419183e-06
## Prueba de Shapiro-Wilk para n = 30 y 50% enfermas | p-valor = 0.0001600657
## Prueba de Shapiro-Wilk para n = 50 y 50% enfermas | p-valor = 0.003857002
## Prueba de Shapiro-Wilk para n = 60 y 50% enfermas | p-valor = 0.01472343
## Prueba de Shapiro-Wilk para n = 100 y 50% enfermas | p-valor = 0.1243522
## Prueba de Shapiro-Wilk para n = 200 y 50% enfermas | p-valor = 0.2914177
## Prueba de Shapiro-Wilk para n = 500 y 50% enfermas | p-valor = 0.2411296
# Mostrar todos los gráficos QQ
grid.arrange(grobs = qq_plots, ncol = 3) # Por ejemplo, 3 gráficos por fila
En el caso con el 10% de la población enferma, tenemos un escenario disimil al anteior. Como se podrá observar, los p-valores bajos en todos los tamaños de muestra indican que las proporciones muestrales no siguen una distribución normal, incluso con tamaños de muestra grandes. Esto puede estar relacionado con la baja prevalencia (solo el 10% está enfermo), lo que genera una distribución altamente asimétrica. Esto NO debe asumirse como un rechazo del Teorema del Límite Central en cuanto este experimento de simulación esta sujeto a una de las infinitas posibilidades de aleatoriedad. Por ejemplo, si en el código base se cambia el seed de 12 a 10, se generarían en este caso dos muestras con distirbución normal (n = 200 y n = 500). Estos resultados no se dejan explícitos acá pero se hace la aclaración para la replicabilidad.
# Crear la población con 10% de plantas enfermas
poblacion_10_90 <- c(rep(1, n_poblacion * 0.1), rep(0, n_poblacion * 0.9))
# Simulación con diferentes tamaños de muestra
resultados_10_90 <- list()
for (n in tamaños_muestra) {
proporciones <- replicate(replicaciones, calcular_proporcion(n, poblacion_10_90))
resultados_10_90[[as.character(n)]] <- proporciones
}
# Crear una lista para almacenar los gráficos QQ
qq_plots <- list()
# Mostrar resultados y gráficos
for (n in tamaños_muestra) {
proporciones <- resultados_10_90[[as.character(n)]]
p <- ggplot(data.frame(Proporción = proporciones), aes(sample = Proporción)) +
stat_qq() +
stat_qq_line() +
ggtitle(paste("QQ plot para n =", n, "(10-90)")) +
theme_minimal()
# Añadir el gráfico a la lista
qq_plots[[as.character(n)]] <- p
# Prueba de Shapiro-Wilk
shapiro_test <- shapiro.test(proporciones)
cat("Prueba de Shapiro-Wilk para n =", n, "y 10% enfermas | p-valor =", shapiro_test$p.value, "\n")
}
## Prueba de Shapiro-Wilk para n = 5 y 10% enfermas | p-valor = 1.026133e-29
## Prueba de Shapiro-Wilk para n = 10 y 10% enfermas | p-valor = 4.916357e-22
## Prueba de Shapiro-Wilk para n = 15 y 10% enfermas | p-valor = 6.472866e-18
## Prueba de Shapiro-Wilk para n = 20 y 10% enfermas | p-valor = 1.205613e-14
## Prueba de Shapiro-Wilk para n = 30 y 10% enfermas | p-valor = 6.648186e-12
## Prueba de Shapiro-Wilk para n = 50 y 10% enfermas | p-valor = 2.499437e-09
## Prueba de Shapiro-Wilk para n = 60 y 10% enfermas | p-valor = 6.704392e-09
## Prueba de Shapiro-Wilk para n = 100 y 10% enfermas | p-valor = 1.129913e-05
## Prueba de Shapiro-Wilk para n = 200 y 10% enfermas | p-valor = 0.0005525855
## Prueba de Shapiro-Wilk para n = 500 y 10% enfermas | p-valor = 0.02318499
grid.arrange(grobs = qq_plots, ncol = 3) # Por ejemplo, 3 gráficos por fila
En el caso del 90% de la población enferma, sucede una situación similar a la de la proporción 10-90, y aunque a partir de n = 200 y n = 500, los p-valores son más altos (0.0116 y 0.0139, respectivamente), siguen estando en la zona de rechazo, en este caso para la normalidad. Acá valdría la pena dar la discusión sobre las pruebas de normalidad para estos casos teniendo en cuenta que según vemos en la gráficas QQ, desde la muestra n = 100 parece existir, por lo general, una concordancia de los valores empíricos y teóricos pero que debido a las repeticiones, en algunos casos se puede distorsionar.
# Crear la población con 90% de plantas enfermas
poblacion_90_10 <- c(rep(1, n_poblacion * 0.9), rep(0, n_poblacion * 0.1))
# Simulación con diferentes tamaños de muestra
resultados_90_10 <- list()
for (n in tamaños_muestra) {
proporciones <- replicate(replicaciones, calcular_proporcion(n, poblacion_90_10))
resultados_90_10[[as.character(n)]] <- proporciones
}
# Crear una lista para almacenar los gráficos QQ
qq_plots <- list()
# Mostrar resultados y gráficos
for (n in tamaños_muestra) {
proporciones <- resultados_90_10[[as.character(n)]]
p <- ggplot(data.frame(Proporción = proporciones), aes(sample = Proporción)) +
stat_qq() +
stat_qq_line() +
ggtitle(paste("QQ plot para n =", n, "(90-10)")) +
theme_minimal()
# Añadir el gráfico a la lista
qq_plots[[as.character(n)]] <- p
# Prueba de Shapiro-Wilk
shapiro_test <- shapiro.test(proporciones)
cat("Prueba de Shapiro-Wilk para n =", n, "y 90% enfermas | p-valor =", shapiro_test$p.value, "\n")
}
## Prueba de Shapiro-Wilk para n = 5 y 90% enfermas | p-valor = 2.378894e-27
## Prueba de Shapiro-Wilk para n = 10 y 90% enfermas | p-valor = 2.988038e-21
## Prueba de Shapiro-Wilk para n = 15 y 90% enfermas | p-valor = 9.225024e-19
## Prueba de Shapiro-Wilk para n = 20 y 90% enfermas | p-valor = 1.430499e-16
## Prueba de Shapiro-Wilk para n = 30 y 90% enfermas | p-valor = 2.582591e-12
## Prueba de Shapiro-Wilk para n = 50 y 90% enfermas | p-valor = 3.709934e-09
## Prueba de Shapiro-Wilk para n = 60 y 90% enfermas | p-valor = 1.787977e-07
## Prueba de Shapiro-Wilk para n = 100 y 90% enfermas | p-valor = 0.0003431724
## Prueba de Shapiro-Wilk para n = 200 y 90% enfermas | p-valor = 0.00124656
## Prueba de Shapiro-Wilk para n = 500 y 90% enfermas | p-valor = 0.1301138
grid.arrange(grobs = qq_plots, ncol = 3) # Por ejemplo, 3 gráficos por fila
Ante los resultados anteriores realicé una graficación de la distribución de la densidad de la frecuencia del estimador P para cuatro de las muestras en cada caso. Se observa así lo siguiente:
Definitivamente las muestras pequeñas, debajo incluso del umbral de n = 50, generan distribuciones alejadas de la normalidad.
Dependiendo de la proporción original de las tasas de enfermedad en la población, las muestras tienden a tener sesgos hacia la derecha o la izquierda, el cual si bien se reduce con el aumento del n, no desaparece totalmente.
Puede que la prueba Shapiro Wilks y su evaluación con base en el p-valor como estadístico lleve a conclusiones, en algunos casos, poco prácticas.
library(gridExtra)
library(ggplot2)
# Gráficos de densidad para tamaños 200 y 500 (50-50)
plots_50_50 <- list()
for (n in c(20,50,200, 500)) {
proporciones <- resultados_50_50[[as.character(n)]]
p <- ggplot(data.frame(Proporción = proporciones), aes(x = Proporción)) +
geom_density(fill = "blue", alpha = 0.5) +
ggtitle(paste("Distribución de estimadores p para n =", n, "(50-50)")) +
xlab("Proporción estimada") +
ylab("Densidad")
plots_50_50[[as.character(n)]] <- p
}
# Combinar los gráficos de 50-50 en una cuadrícula de 2 columnas y 1 fila
grid.arrange(plots_50_50[[as.character(20)]],plots_50_50[[as.character(50)]],plots_50_50[[as.character(200)]], plots_50_50[[as.character(500)]], ncol = 2, nrow=2)
# Gráficos de densidad para tamaños 200 y 500 (10-90)
plots_10_90 <- list()
for (n in c(20,50,200, 500)) {
proporciones <- resultados_10_90[[as.character(n)]]
p <- ggplot(data.frame(Proporción = proporciones), aes(x = Proporción)) +
geom_density(fill = "red", alpha = 0.5) +
ggtitle(paste("Distribución de estimadores p para n =", n, "(10-90)")) +
xlab("Proporción estimada") +
ylab("Densidad")
plots_10_90[[as.character(n)]] <- p
}
# Combinar los gráficos de 10-90 en una cuadrícula de 2 columnas y 1 fila
grid.arrange(plots_10_90[[as.character(20)]],plots_10_90[[as.character(50)]],plots_10_90[[as.character(200)]], plots_10_90[[as.character(500)]], ncol = 2, nrow=2)
# Gráficos de densidad para tamaños 200 y 500 (90-10)
plots_90_10 <- list()
for (n in c(20,50,200, 500)) {
proporciones <- resultados_90_10[[as.character(n)]]
p <- ggplot(data.frame(Proporción = proporciones), aes(x = Proporción)) +
geom_density(fill = "green", alpha = 0.5) +
ggtitle(paste("Distribución de estimadores p para n =", n, "(90-10)")) +
xlab("Proporción estimada") +
ylab("Densidad")
plots_90_10[[as.character(n)]] <- p
}
# Combinar los gráficos de 90-10 en una cuadrícula de 2 columnas y 1 fila
grid.arrange(plots_90_10[[as.character(20)]],plots_90_10[[as.character(50)]],plots_90_10[[as.character(200)]], plots_90_10[[as.character(500)]], ncol = 2, nrow=2)