Pruebas de hipótesis para proporciones sobre datos geo-referenciados de El Salvador
DataIntelligence
09-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("datossalvador.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 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 336548## Obtenemos las proporciones de la variable nominal "ClaveCob1Cob2":
prop.table(table(data$ClaveCob1Cob2))##
## BosqueBosque BosqueNo Bosque No BosqueBosque No BosqueNo Bosque
## 0.30086528 0.02918448 0.03203526 0.63791499## Obtenemos los porcentajes de la variable nominal "ClaveCob1Cob2":
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 NO se distribuyen en la misma proporción que la población con un nivel de significación del 99%.
Cálculo de los límites del estadístico Z para distintos niveles de significación:
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.575829Valor del \(\alpha =\) 0.05 para dos colas:
alpha = .05
z.half.alpha = qnorm(1-alpha / 2)
c(-z.half.alpha , z.half.alpha)## [1] -1.959964 1.959964Cá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 población total:
data_m <- data[sample(NROW(data), NROW(data)*(1 - 0)),]
n = nrow(data_m)
n## [1] 527575#### Veamos los primeros 5 registros de la muestra:
head(data_m,5)## OBJECTID SAMPLE_ID PLOT_ID LON LAT ClaveCob1Cob2
## 286307 286307 7 19344 -87.97307 13.70159 BosqueNo Bosque
## 96007 96007 15 20983 -88.19639 13.99931 No BosqueNo Bosque
## 409865 409865 10 18010 -88.08313 13.50260 BosqueBosque
## 81681 81681 1 3980 -89.72409 14.03164 No BosqueNo Bosque
## 380683 380683 16 14777 -88.34185 13.55510 No BosqueNo Bosque# pbar es la proporción de la muestra
pbar = prop.table(table(data_m$ClaveCob1Cob2))
pbar ##
## BosqueBosque BosqueNo Bosque No BosqueBosque No BosqueNo Bosque
## 0.30086528 0.02918448 0.03203526 0.63791499# p0 es la proporcion en la poblacion:
p0_BB = 0.30086528
p0_BNoB = 0.02918448
p0_NoBB = 0.03203526
p0_NoBNoB = 0.63791499
# z_ es el valor p-value para la correspondiente categoría:
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 límite para \(\alpha\) = .01, z = [-2.575829; 2.575829].
Como hemos señalado, no puede ser alcanzado un tamaño muestral óptimo por la naturaleza aleatoria de la extracción, en cada extraccion muestral obtendremos un p-value distinto aunque aproximado a los límites z a medida que disminuya su tamaño.
Calculo para una muestra de 500:
Podemos creer que una muestra de 500 es aproximadamente siempre válida:
muestras_500 <- data[sample(NROW(data), NROW(data)*(1 - 0.999)),]
tamanio_subdata = nrow(muestras_500)
# Obtengamos el tamaño muestral:
tamanio_subdata## [1] 527# Listemos los primeros 5 datos de la muestra:
head(muestras_500,5)## OBJECTID SAMPLE_ID PLOT_ID LON LAT ClaveCob1Cob2
## 92031 92031 16 20998 -88.22408 14.00772 BosqueBosque
## 359702 359702 2 6390 -89.16474 13.58560 BosqueBosque
## 70867 70867 12 4138 -89.65940 14.05958 BosqueBosque
## 304968 304968 13 12138 -88.63845 13.67108 No BosqueNo Bosque
## 519281 519281 16 13677 -88.32150 13.21168 BosqueBosquepbar = prop.table(table(muestras_500$ClaveCob1Cob2))
# Obtengamos las proporciones de la muestra:
pbar ##
## BosqueBosque BosqueNo Bosque No BosqueBosque No BosqueNo Bosque
## 0.30360531 0.03036053 0.03795066 0.62808349# Calculemos los p-values de cada categoría:
z_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.1371496z_BNoB## BosqueNo Bosque
## 0.1603937z_NoBB## No BosqueBosque
## 0.7711624z_NoBNoB## No BosqueNo Bosque
## -0.4696111Grafica de los p-values por categoria versus el tamano muestral para 10 extracciones aleatorias:
Grafiquemos los p_value contra el tamaño muestral 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).
```r
# Declaramos vectores vacios:
z_Bos_Bos <- c()
z_Bos_no_Bos <- c()
z_No_Bos_Bos <- c()
z_No_Bos_no_Bos <- c()
muestras_decrecientes <- c()
n_row_vector <- c()
# Declaramos una matriz vacia:
la_matriz<-matrix()
# Declaramos una lista vacia que almacene nuestras matrices:
li<-list()
# Esta secuencia sirve para ir extrayendo muestras aleatorias cada vez que se recorre el ciclo, disminuyendo el tamano de estas en 5000 cada vez
my_seq <- seq(1, 527575, by=5000)
# Esta secuencia sirve para hacer el mismo procedimiento al que sirve la secuencia anterior
my_seq_2 <- seq(1, 2, by=1)
for (i in my_seq_2)
{
for (j in my_seq)
{
muestras_decrecientes <- data[sample(NROW(data), NROW(data)*(1)),]
muestras_decrecientes <- muestras_decrecientes[sample(NROW(muestras_decrecientes), NROW(muestras_decrecientes)- j),]
# muestras_decrecientes <- muestras_decrecientes[-sample(i)]
n_row <- nrow(muestras_decrecientes)
pbar = prop.table(table(muestras_decrecientes$ClaveCob1Cob2))
# pbar es la proporcion de la categoria asociada a la muestra:
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)
n_row_vector <- append(n_row_vector, n_row)
# Esta linea nos entrega el tamano de las muestras:
# print(n_row,10)
la_matriz <- z_Bos_Bos
la_matriz <- cbind(la_matriz, z_Bos_no_Bos)
la_matriz <- cbind(la_matriz, z_No_Bos_Bos)
la_matriz <- cbind(la_matriz, z_No_Bos_no_Bos)
# Agregamos una columna con el tamano muestral:
la_matriz <- cbind(la_matriz, n_row_vector)
}
li[[i]] <- la_matriz
}
# Esta linea nos entrega la primera matriz de la lista:
# li[[1]]
# Esta linea nos entrega la primera matriz de la lista:
# li[[2]]
# Transformamos la lista a un dataframe:
df <- as.data.frame(li)## Warning in (function (..., row.names = NULL, check.rows = FALSE, check.names =
## TRUE, : row names were found from a short variable and have been discarded# Con estas lineas imprimimos valores del dataframe y de una columna especifica:
# df
# df$z_No_Bos_no_Bos.1
# Con esta linea desplegamos el contenido de la lista:
# li
p = ggplot() +xlim(527574, 2574)+
geom_line(data = df, aes(x = n_row_vector, y = df$z_No_Bos_no_Bos), color = "blue") +
geom_line(data = df, aes(x = n_row_vector, y = df$z_No_Bos_no_Bos.1), color = "red") +
xlab('Tamaño muestral') +
ylab('p-value') + ggtitle('Gráfica del p-value para No Bosque no Bosque en dos extracciones aleatorias')
print(p)## Warning: Use of `df$z_No_Bos_no_Bos` is discouraged. Use `z_No_Bos_no_Bos`
## instead.## Warning: Use of `df$z_No_Bos_no_Bos.1` is discouraged. Use `z_No_Bos_no_Bos.1`
## instead.
########
head(z_Bos_Bos,10)## BosqueBosque BosqueBosque BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## -0.00209900 0.03808602 0.04228951 -0.01742022 0.22220742 -0.20934934
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## 0.11826064 0.09191316 0.68669222 -0.12793905plot(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.0002220545 -0.0086329213 0.0484899432 0.1393607705 0.1227207070
## BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque
## -0.1454878343 0.2068628810 0.1903235388 0.0119654309 -0.0142927409plot(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.0002325233 0.1194533561 0.0188050183 -0.2890496086 0.1812175471
## No BosqueBosque No BosqueBosque No BosqueBosque No BosqueBosque No BosqueBosque
## 0.2876646568 0.0973117352 0.0587741934 -0.0696149408 0.2175825813plot(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)