Pruebas de hipótesis para proporciones sobre datos geo-referenciados de El Salvador
DataIntelligence
04-06-2020
Introducción:
Tenemos un conjunto de 527575 observaciones agrupadas en cuatro categorías:
- Bosque Bosque
- Bosque no Bosque
- No Bosque Bosque
- No Bosque no Bosque
Deseamos optimizar la información de la población extrayendo una muestra representativa que mantenga las proporciones de las categorías a cierto nivel de significación, digamos, con un 99% de certeza.
Para ello aplicaremos experimentos, que consisten en extraer muestras y aplicarles un estadístico para que al compararlos con con un valor límite de una distribución normal asociado al 0.05% decidamos el no rechazar una hipótesis nula o aceptar la alternativa.
Nuestra hipótesis nula será que la muestra que extraemos contendrá las mismas proporciones que la población con un nivel de significación estadístico.
Por el caracter experimental de las pruebas de hipótesis, no se puede determinar un tamaño muestral óptimo n, sin embargo, estimaremos su tamaño inicial graficando el momento en el que el tamaño de las muestras excede el valor límite asociado al nivel de significación (y por lo tanto rechazamos nuestra hipótesis nula).
Las pruebas de hipótesis son experimentos, por lo que la solución es un problema eminentemente empirico.
# Leemos los datos, construímos un dataframe y obtenemos el número de observaciones:
data <- read_excel("datos salvador.xlsx")
data <- data.frame(data)
nrow(data)## [1] 527575# Obtenemos el nombre de las columnas del dataset que cargamos:
colnames(data)## [1] "OBJECTID" "SAMPLE_ID" "PLOT_ID" "LON"
## [5] "LAT" "ClaveCob1Cob2"## Obtenemos las categorias y las frecuencias de la variable nominal "ClaveCob1Cob2":
count(data, 'ClaveCob1Cob2')## ClaveCob1Cob2 freq
## 1 BosqueBosque 158729
## 2 BosqueNo Bosque 15397
## 3 No BosqueBosque 16901
## 4 No BosqueNo Bosque 336548prop.table(table(data$ClaveCob1Cob2))##
## BosqueBosque BosqueNo Bosque No BosqueBosque No BosqueNo Bosque
## 0.30086528 0.02918448 0.03203526 0.63791499## Obtenemos los porcentajes
prop.table(table(data$ClaveCob1Cob2))*100##
## BosqueBosque BosqueNo Bosque No BosqueBosque No BosqueNo Bosque
## 30.086528 2.918448 3.203526 63.791499El problema de investigación es encontrar una muestra aleatoria que satisfaga estas cuatro proporciones simultáneamente utilizando un criterio estadístico con una significación del 99%.
Aplicaremos un test de hipótesis para proporciones poblacionales en el que nuestra hipótesis cero será el que la muestra que obtengamos satisface la proporción de la data original con un 99% de significación. Para ello utilizaremos el estadístico p-value y lo compararemos con el valor z asociado a un 0.01% de error (dos colas).
Hipótesis estadísticas:
\(H_0:\) Las categorías de nuestra muestra se distribuyen en la misma proporción que la población con un nivel de significación del 99%.
\(H_1:\) Las categorías de nuestra muestra se distribuyen en la misma proporción que la población con un nivel de significación del 99%.
Valor del \(\alpha =\) 0.01 para dos colas:
alpha = .01
z.half.alpha = qnorm(1-alpha / 2)
c(-z.half.alpha , z.half.alpha)## [1] -2.575829 2.575829Cálculo del estadístico p-value:
Como ejercicio consideremos una muestra que sea el total de la poblacion y apliquemos el estadistico p-value. Observaremos que éste tiende a 0 para todas las categorías (y no es cero exacto por el problema del redondeo computacional).
\(Z=\frac{{\hat{p}}-{p}_{0}}{\sqrt{\frac{{p}_{0}(1-{p}_{0})}{n}}}\)
Muestra igual a la poblacion total:
data_m <- data[sample(NROW(data), NROW(data)*(1 - 0)),]
n = nrow(data_m)
n## [1] 527575head(data_m,10)## OBJECTID SAMPLE_ID PLOT_ID LON LAT ClaveCob1Cob2
## 16033 16033 8 10284 -89.22688 14.29901 BosqueBosque
## 65833 65833 18 4238 -89.52068 14.07926 No BosqueNo Bosque
## 307836 307836 14 948 -89.37774 13.66522 No BosqueNo Bosque
## 490907 490907 2 16694 -88.10083 13.30326 No BosqueNo Bosque
## 282762 282762 7 1232 -90.00651 13.70330 BosqueBosque
## 294968 294968 8 1102 -89.84005 13.68733 BosqueBosque
## 496259 496259 20 16561 -88.24805 13.28486 BosqueNo Bosque
## 230690 230690 7 7486 -88.93522 13.77727 No BosqueNo Bosque
## 295935 295935 15 12175 -88.86966 13.68777 No BosqueNo Bosque
## 13149 13149 4 10343 -89.33837 14.31614 BosqueBosquepbar = prop.table(table(data_m$ClaveCob1Cob2))
pbar ##
## BosqueBosque BosqueNo Bosque No BosqueBosque No BosqueNo Bosque
## 0.30086528 0.02918448 0.03203526 0.63791499#pbar[1]
p0_BB = 0.30086528 # constante
p0_BNoB = 0.02918448
p0_NoBB = 0.03203526
p0_NoBNoB = 0.63791499
z_BB = (pbar[1] - p0_BB)/sqrt(p0_BB*(1-p0_BB)/n)
z_BNoB = (pbar[2] - p0_BNoB)/sqrt(p0_BNoB*(1-p0_BNoB)/n)
z_NoBB = (pbar[3] - p0_NoBB)/sqrt(p0_NoBB*(1-p0_NoBB)/n)
z_NoBNoB = (pbar[4] - p0_NoBNoB)/sqrt(p0_NoBNoB*(1-p0_NoBNoB)/n)Observamos que el estadistico p-value tiende a cero:
z_BB## BosqueBosque
## -2.881792e-07z_BNoB## BosqueNo Bosque
## -1.665296e-05z_NoBB## No BosqueBosque
## -1.793916e-05z_NoBNoB## No BosqueNo Bosque
## -2.432802e-06Mientras menor sea la muestra cada vez será mayor el estadístico p-value hasta que alcance el valor limite para alpha = .01, z = [-2.575829; 2.575829]
Como hemos señalado, no puede ser aslcanzado un tamaño muestral óptimo por la naturaleza aleatoria de la extraccion, cada vez dará un p-value distinto aunque aproximado a medida que aumenta el tamaño muestral.
Podemos creer que una muestra de 500 es aproximadamente siempre válida:
data_sub_m <- data[sample(NROW(data), NROW(data)*(1 - 0.99905)),]
tamanio_subdata = nrow(data_sub_m)
tamanio_subdata## [1] 501head(data_sub_m,10)## OBJECTID SAMPLE_ID PLOT_ID LON LAT ClaveCob1Cob2
## 218168 218168 14 7576 -89.04620 13.79469 No BosqueNo Bosque
## 384304 384304 14 190 -89.61664 13.54537 No BosqueNo Bosque
## 217858 217858 13 2068 -89.79502 13.79623 No BosqueNo Bosque
## 63076 63076 6 4276 -89.55804 14.08766 No BosqueNo Bosque
## 183720 183720 25 15706 -88.33412 13.84491 No BosqueNo Bosque
## 473858 473858 3 16998 -88.31343 13.35661 No BosqueNo Bosque
## 193254 193254 15 7760 -89.23145 13.82935 No BosqueNo Bosque
## 113040 113040 15 3512 -89.36227 13.96365 No BosqueNo Bosque
## 272518 272518 18 19491 -87.76944 13.72051 BosqueNo Bosque
## 88441 88441 21 21024 -88.13142 14.01721 No BosqueBosquepbar = prop.table(table(data_sub_m$ClaveCob1Cob2))
pbar ##
## BosqueBosque BosqueNo Bosque No BosqueBosque No BosqueNo Bosque
## 0.31536926 0.02395210 0.03992016 0.62075848z_BB = (pbar[1] - p0_BB)/sqrt(p0_BB*(1-p0_BB)/tamanio_subdata)
z_BNoB = (pbar[2] - p0_BNoB)/sqrt(p0_BNoB*(1-p0_BNoB)/tamanio_subdata)
z_NoBB = (pbar[3] - p0_NoBB)/sqrt(p0_NoBB*(1-p0_NoBB)/tamanio_subdata)
z_NoBNoB = (pbar[4] - p0_NoBNoB)/sqrt(p0_NoBNoB*(1-p0_NoBNoB)/tamanio_subdata)
z_BB## BosqueBosque
## 0.7078474z_BNoB## BosqueNo Bosque
## -0.6957835z_NoBB## No BosqueBosque
## 1.002239z_NoBNoB## No BosqueNo Bosque
## -0.7990263Grafiquemos los p_value para cada una de nuestras categorías para determinar cuando se acerca al valor límite en el que se rechaza la hipótesis nula (línea punteada roja).
# bosque a bosque
z_BB## BosqueBosque
## 0.7078474Grafica de los p_values:
my_seq <- seq(0.02, 0.5, by=0.01)
z_Bos_Bos <- c()
z_Bos_no_Bos <- c()
z_No_Bos_Bos <- c()
z_No_Bos_no_Bos <- c()
muestras_decrecientes <- data[sample(NROW(data), NROW(data)*(1 - 0.001)),]
for (i in my_seq)
{
muestras_decrecientes <- muestras_decrecientes[sample(NROW(muestras_decrecientes), NROW(muestras_decrecientes)*(1 - i)),]
n_row = nrow(muestras_decrecientes)
pbar = prop.table(table(muestras_decrecientes$ClaveCob1Cob2))
z_Bosque_Bosque = (pbar[1] - p0_BB)/sqrt(p0_BB*(1-p0_BB)/n_row)
z_Bos_Bos <- append(z_Bos_Bos, z_Bosque_Bosque)
z_Bosque_no_Bosque = (pbar[2] - p0_BNoB)/sqrt(p0_BNoB*(1-p0_BNoB)/n_row)
z_Bos_no_Bos <- append(z_Bos_no_Bos, z_Bosque_no_Bosque)
z_No_Bosque_Bosque = (pbar[3] - p0_NoBB)/sqrt(p0_NoBB*(1-p0_NoBB)/n_row)
z_No_Bos_Bos <- append(z_No_Bos_Bos, z_No_Bosque_Bosque)
z_No_Bosque_no_Bosque = (pbar[4] - p0_NoBNoB)/sqrt(p0_NoBNoB*(1-p0_NoBNoB)/n_row)
z_No_Bos_no_Bos <- append(z_No_Bos_no_Bos, z_No_Bosque_no_Bosque)
}
# Gráfica del p-value para Bosque Bosque
z_Bos_Bos## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## -0.035563878 -0.158684694 -0.364297926 -0.588236446
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## -0.565057997 -0.614174205 -0.534392950 -0.707929982
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## -0.590284190 0.139977732 0.252569951 0.380155810
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## 0.955833302 0.419486632 0.662024133 0.885991768
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## 0.562657656 0.486895005 0.703636038 1.475485849
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## 0.420686320 0.709262184 1.132112606 1.202839664
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## 1.313055657 0.929146362 0.626177725 0.291642352
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## 0.643102542 1.007149263 0.562793127 1.149454377
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## 0.638585016 -0.368482774 -0.126494232 -1.385215260
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## -0.419047889 -0.415186584 -0.085787584 -0.008437333
## BosqueBosque BosqueBosque No BosqueNo Bosque No BosqueNo Bosque
## -0.860892025 -0.716732639 2.640310642 1.524384060
##
## NA NA NA NA
##
## NAhead(z_Bos_Bos,10)## BosqueBosque BosqueBosque BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## -0.03556388 -0.15868469 -0.36429793 -0.58823645 -0.56505800 -0.61417420
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## -0.53439295 -0.70792998 -0.59028419 0.13997773plot(z_Bos_Bos, type="l", col="blue", lwd=1, xlab="Tamaño muestral")
abline(h=-2.575829, col="red", lwd=0.5, lty=6)
abline(h=2.575829, col="red", lwd=0.5, lty=6)
# Gráfica del p-value para Bosque no Bosque
head(z_Bos_no_Bos,10)## BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque
## -0.10712489 -0.17387822 -0.11847436 0.04418541 0.02854689
## BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque
## 0.21253250 -0.04675782 -0.37522414 0.22129038 -0.15329908plot(z_Bos_no_Bos, type="l", col="green", lwd=1, xlab="Tamaño muestral")
abline(h=-2.575829, col="red", lwd=0.5, lty=6)
abline(h=2.575829, col="red", lwd=0.5, lty=6)
# Gráfica del p-value para No Bosque Bosque
head(z_No_Bos_Bos,10)## No BosqueBosque No BosqueBosque No BosqueBosque No BosqueBosque No BosqueBosque
## -0.003192284 -0.048022140 0.073965835 0.356618191 0.543540345
## No BosqueBosque No BosqueBosque No BosqueBosque No BosqueBosque No BosqueBosque
## 0.619083855 0.708380526 1.111176327 0.584270160 0.976931664plot(z_No_Bos_Bos, type="l", col="magenta", lwd=1, xlab="Tamaño muestral")
abline(h=-2.575829, col="red", lwd=0.5, lty=6)
abline(h=2.575829, col="red", lwd=0.5, lty=6)
# Gráfica del p-value para No Bosque no Bosque
head(z_No_Bos_no_Bos,10)## No BosqueNo Bosque No BosqueNo Bosque No BosqueNo Bosque No BosqueNo Bosque
## 0.07261163 0.22990964 0.36202347 0.41519270
## No BosqueNo Bosque No BosqueNo Bosque No BosqueNo Bosque No BosqueNo Bosque
## 0.33006278 0.28481713 0.26677704 0.39983698
## No BosqueNo Bosque No BosqueNo Bosque
## 0.27170938 -0.43784868plot(z_No_Bos_no_Bos, type="l", col="purple", lwd=1, xlab="Tamaño muestral")
abline(h=-2.575829, col="red", lwd=0.5, lty=6)
abline(h=2.575829, col="red", lwd=0.5, lty=6)
Conclusión
Podemos extraer una muestra inicial con el 60% de la población y aplicarle el test estadístico p-value para corroborar que ésta muestra específica sea representativa de la población.
#
#
#