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("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
## 453276 453276 1 11221 -88.80310 13.40759 No BosqueNo Bosque
## 108998 108998 17 8569 -89.29747 13.97294 No BosqueNo Bosque
## 301299 301299 24 6917 -89.00796 13.67752 No BosqueNo Bosque
## 251547 251547 7 7333 -88.93500 13.75016 BosqueBosque
## 136917 136917 17 3132 -89.53740 13.91635 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
## 36645 36645 10 9740 -89.31832 14.18092 No BosqueNo Bosque
## 54538 54538 18 4407 -89.48405 14.11578 No BosqueNo Bosque
## 359959 359959 9 6400 -89.07224 13.58663 No BosqueNo Bosque
## 207713 207713 20 20083 -88.26929 13.80910 No BosqueBosque
## 218946 218946 1 12560 -88.75055 13.79661 BosqueBosquepbar = prop.table(table(muestras_500$ClaveCob1Cob2))
# Obtengamos las proporciones de la muestra:
pbar ##
## BosqueBosque BosqueNo Bosque No BosqueBosque No BosqueNo Bosque
## 0.26565465 0.03984820 0.02087287 0.67362429# 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
## -1.762433z_BNoB## BosqueNo Bosque
## 1.454352z_NoBB## No BosqueBosque
## -1.455187z_NoBNoB## No BosqueNo Bosque
## 1.705689Grafica 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.002099000 0.013956452 -0.078939709 -0.002192837 -0.151166235 -0.037114288
## BosqueBosque BosqueBosque BosqueBosque BosqueBosque
## -0.261937336 -0.088274801 0.399415068 0.003886885plot(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.0406768532 0.1971320815 0.1725530722 0.3395307168
## BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque BosqueNo Bosque
## 0.2483827200 -0.0963372093 -0.1228763843 -0.0560998095 -0.1938872274plot(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.0298040214 0.0661660472 0.1392735712 0.1015087196
## No BosqueBosque No BosqueBosque No BosqueBosque No BosqueBosque No BosqueBosque
## 0.2155707426 0.1214634805 0.2610576523 -0.2160039059 0.0786119385plot(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)