10 Inferencia por intervalos

En esta parte consideramos hacer inferencia a una población. Quiséramos decir algo acerca de la mediana de una población cuando solo tenemos información parcial.

En este ejemplo, queremos aprender acerca de los precios de casas en una región determinada (los precios están en miles de dólares).

Supongamos que tenemos la siguiente datos, que sólo son parte de la población:

Código

set.seed(6622)
muestra_casas <- datos_casas |> slice_sample(n = 20)
muestra_casas

# A tibble: 20 × 3
      id precio_miles orden
   <dbl>        <dbl> <int>
 1   764        337     580
 2  1037        316.    562
 3  1108        275.    519
 4  1400        137.    176
 5   803        189     353
 6   403        108      84
 7   217        210     401
 8   352        190     355
 9  1094        146     210
10   763        215.    412
11  1155        202.    382
12   154        235     450
13   590         79.5    25
14   439         90.4    46
15   594        140     189
16   951        129     144
17  1376        239     457
18  1418        340     583
19   517        158     250
20   550        263     501

Código

ggplot(muestra_casas, aes(sample = precio_miles ) ) +
  geom_qq(distribution = stats::qunif)

Supongamos que queremos estimar el valor \(m\) que es la mediana poblacional. Este valor tiene la propiedad de que separa a los datos en dos grupos de tamaño igual.

Tip

Estimar en este caso significa dar un rango donde es creíble que está el valor poblacional. No tiene sentido dar un sólo valor pues no conocemos la población total, así que tendremos incertidumbre de dónde está la mediana verdadera

Para construir estos rangos haremos algunos cálculos básicos.

En primer lugar, si observáramos un valor \(y\) en la muestra, ¿cuál es la probabilidad de que caiga por arriba de la mediana? Si no sabemos cómo se extrajo la muestra, o se extrajo bajo un mecanismo complicado, esta pregunta es difícil de contestar.

Sin embargo, si este valor se extrajo tomando un valor al azar entre todos los de la población, entonces sabemos que la probabilidad es

\[P(y <= m) = 0.5\]

Aleatorización al rescate

Supondremos que cada elemento de la muestra se extrajo de la población total al azar, cada uno de manera independiente de los otros.

Para poblaciones finitas, esto equivale a:

Seleccionamos un elemento al azar de la población, donde todos tienen la misma probabildad de ser seleccionado, por ejemplo seleccionando un número de una bolsa que corresponde a cada elemento de la población
Repetimos 1, considerando que es posible sacar al mismo elemento más de una vez.

A este proceso le llamamos muestreo aleatorio simple con reemplazo, que cuando la población es grande y la muestra relativamente chica, prácticamente es igual que hacer el muestreo sin reemplazo. Más adelante discutiremos cómo generalizar a procesos aleatorios diferentes.

Con esta idea básica, podríamos por ejemplo calcular: ¿cuál es la probabilidad de que la mediana esté entre el mínimo y el máximo de la muestra?

Denotemos por \(y_1,y_2,\ldots, y_n\) la muestra observada. Denotaremos a la muestra ordenada del valor más chico al valor más grande como

\[y_{(1)}, y_{(2)}, \ldots, y_{(n)}\]

Nuestra primera pregunta entonces es:

¿Cuál es la probabilidad de el valor de la mediana verdadera esté en el intervalo \([y_{(1)}, y_{(20)}]\)?

Respuesta: La probabilidad de que la mediana esté por debajo del mínimo \(y_{(1)}\) de la muestra es \((1/2)^n\), pues todos los valores de la muestra tienen que caer por arriba de la mediana. La probabilidad de que la mediana esté por arriba del máximo \(y_{(1)}\) es también \((1/2)^n\). Así que la probabilidad de el verdadero valor de la mediana esté en este intervalo es de \(2(1/2)^n\).

Esta es la probabilidad de tirar solamente soles o solamente águilas en \(n\) tiradas de una moneda, Si \(X\) denota el número de soles obtenidos en \(n\) tiradas de moneda, buscamos calcular la probabilidad de que la mediana esté en intervalo como

\[1 - P(X > n - 1) - P(X \leq 0) = 1 - 2(1/2)^n\]

Código

#! code-fold: false
n <- 20
bajo <- 1
alto <- n - bajo + 1
1 - pbinom(alto - 1, 20, 0.5, lower.tail = FALSE) - pbinom(bajo - 1, 20, 0.5)

[1] 0.9999981

Código

1- 2*(1/2)^20

[1] 0.9999981

Prácticamente es seguro que este intervalo contenga a la media verdadero. Esto no es muy informativo. Probemos ahora con \([y_{(2)}, y_{(n-1)}]\). Por un argumento similar, la probabilidad de la mediana verdadera esté en este intervalo es:

Código

#! code-fold: false
n <- 20
bajo <- 2
alto <- n - bajo + 1
1 - pbinom(alto - 1, 20, 0.5, lower.tail = FALSE) - pbinom(bajo - 1, 20, 0.5)

[1] 0.9999599

Esto también es altamente probable. Sin embargo, la probabilidad de la mediana poblacional esté entre \([y_{(7)}, y_{(14)}]\) es al menos

Código

#! code-fold: false
n <- 20
bajo <- 7
alto <- n - bajo + 1
1 - pbinom(alto - 1, 20, 0.5, lower.tail = FALSE) - pbinom(bajo - 1, 20, 0.5)

[1] 0.8846817

Así, este intervalo tiene probabilidad de casi 90% de probabilidad de contener al verdadero valor:

Esta es una probabilidad que se refiere al conjunto de muestras posibles de los datos. Es altamente probable que este intervalo cubra al verdadero valor.

Ahora extraemos el intervalo correspondiente:

Código

intervalo <- muestra_casas |> mutate(pos = dense_rank(precio_miles)) |> 
  arrange(pos) |> 
  filter(pos %in% c(7, 14))
intervalo |> select(pos, precio_miles)

# A tibble: 2 × 2
    pos precio_miles
  <int>        <dbl>
1     7          146
2    14          235

Y este es un intervalo de confianza de cerca de 90% para la mediana de la población. Cuando tomamos muestras más grandes, podemos obtener mejores precisiones. Por ejemplo, para una muestra de 151 casas, usaríamos

Código

n <- 151
bajo <- 65
alto <- n - bajo - 1
cobertura <- 1 - pbinom(alto - 1, n, 0.5, lower.tail = FALSE) - pbinom(bajo - 1, n, 0.5)
cobertura

[1] 0.8921164

Tomamos una muestra y extraemos el intervalo:

Código

set.seed(3555881)
muestra_2 <- datos_casas |> slice_sample(n = 151, replace = TRUE) |> 
  mutate(pos = row_number(precio_miles)) 
muestra_2 |> 
  arrange(pos) |> 
  filter(pos %in% c(bajo, alto)) |> 
  select(precio_miles, pos)

# A tibble: 2 × 2
  precio_miles   pos
         <dbl> <int>
1          141    65
2          165    85

Esto intervalos tienen la garantía inferencial, y no es necesario hacer ningún supuesto, excepto en que la muestra se escoge al azar.

10.1 Motivación de remuestreo

Podemos pensar en este procedimiento de una manera diferente, sin tener que hacer los cálculos de probabilidades para las estadísticas de orden.

Supongamos que escogemos la muestra de casas al azar de manera independiente. Nuestra pregunta es:

Dada la información de la muestra, qué podemos decir de los valores que no observamos?, y finalmente, ¿qué información nos da sobre la mediana poblacional?

Nótese en primer lugar, igual que antes, que no podemos decir algo definitivo acerca de la muestra particular que obtuvimos. Sin embargo, podemos decir qué puede pasar con alta probabilidad, quitando casos excepcionalmente raros, y basar nuestra inferencia en los casos con la mayor parte de la probabilidad.

Nuestra muestro ordenada la escribimos como \(y_{(1)}\leq y_{(2)}\leq \cdots \leq y_{(n)}\). Podemos considerar los intervalos aleatorios:

\[I_1 = [0, y_{(1)}), I_2 = [y_{(1)}, y_{(2)}), I_3 = [y_{(2)}, y_{(3)}),\ldots, I_{n+1} = [y_{(n)}, \infty ]\]

que cubren todos los valores posibles que puede haber en la población.

Ahora consideramos un dato \(y\) extraído también al azar de la población. Es igualmente probable que \(y\) caiga en cualquiera de los intervalos \(I\). Para demostrar esto, considera que pasaría si tomaras una muestra de tamaño \(n+1\). Todos los ordenamientos de \(y_1, y_2, \ldots, y_n, y_{n+1}\) son igualmente probables: esto implica que es igualmente probable que \(y_{n+1}\) caiga en cualquiera de estos intervalos (probabilidad \(1/{(n+1)}\)).

Podemos imputar nuestra población entera tomando escogiendo para cada valor no observado un intervalo al azar de manera equiprobable.
Calculamos dónde puede estar la mediana para esta población.
Repetimos los dos anteriores una gran cantidad de veces para ver todos los posibles valores que podría tomar la mediana dada la información que tenemos de la muestra.

Veamos un ejemplo chico ilustrado para entender la idea. Supongamos que tenemos una población de tamaño 25, y que tomamos una muestra de 8 elementos al azar. Los valores obtenidos, ordenados, son 12, 15, 16, 17, 18, 20, 23, 30, 32. Tenemos entonces 9 intervalos que considerar.

Hacemos ahora el paso el paso 1, y la población imputada ordenada de chica a grande es:

Código

set.seed(7212)
muestra <- c(12, 15, 16, 17, 18, 20, 23, 30, 32)
graficar_pob <- function(muestra, n_pob){
  n_muestra <- length(muestra)
  intervalos <- sample(1:(n_muestra + 1), n_pob - n_muestra, replace = TRUE) |> sort()
  muestra_ord <- sort(muestra)
  tab_muestra <- tibble(pos = 1:length(muestra), valor = as.character(muestra_ord))
  tab_graf <- bind_rows(tab_muestra, tibble(pos = intervalos - 0.5, valor = "*")) |> 
    arrange(pos)
  tab_graf$valor 
}
graficar_pob(muestra, 25) |> paste(collapse = " ")

[1] "* 12 * * * 15 * * 16 * * 17 * * 18 * 20 * * 23 30 * 32 * *"

En este caso, contando casos encontramos que la mediana está entre 17 y 18, por ejemplo podríamos aproximar con 17.5. Hagamos otra repetición:

Código

graficar_pob(muestra, 25) |> paste(collapse = " ")

[1] "* 12 15 16 * * * 17 * * * 18 * 20 * 23 * * * 30 * * * 32 *"

Y en este ejemplo, la mediana está entre 18 y 20. Podríamos aproximar a 19.5. Repetimos este proceso varias veces y examinamos los resultados:

Código

set.seed(12)
pob_imputadas <- tibble(rep = 1:10) |> 
  mutate(poblacion = map_chr(rep, ~ graficar_pob(muestra, 20) |> paste(collapse = " ")))
pob_imputadas

# A tibble: 10 × 2
     rep poblacion                                       
   <int> <chr>                                           
 1     1 12 * * * * 15 16 17 * * 18 * 20 * 23 * 30 * 32 *
 2     2 12 15 16 * 17 * 18 20 * * 23 * * * 30 * 32 * * *
 3     3 12 * 15 * 16 * * * * 17 * 18 20 * * 23 30 * 32 *
 4     4 12 * 15 * 16 * 17 * * 18 * 20 * * * 23 * 30 32 *
 5     5 12 15 16 * 17 * 18 20 * 23 * * * * * * 30 * 32 *
 6     6 12 * 15 * 16 * 17 * * * 18 * 20 * 23 * 30 32 * *
 7     7 12 * * * 15 * * 16 17 * 18 * * 20 23 * 30 32 * *
 8     8 * 12 * * 15 16 * 17 18 * 20 * * * * 23 30 * * 32
 9     9 12 15 * 16 * 17 * * 18 20 * * 23 * * * 30 * 32 *
10    10 * * * 12 * * 15 * * 16 17 * 18 20 * 23 * 30 * 32

Nótese que como al final vamos a dar un intervalo con extremos en los datos, en realidad no importa mucho qué valor esté representado por las estrellas (excluyendo el caso extremo de que la todos más de la mitad de las estrellas estén por abajo/arriba de todos los datos, lo cual es poco probable). Podemos contar las estrellas en cada posición como un valor en los datos, y ver entre qué valores cae la mediana. De esta forma podríamos encontrar un intervalo que contenga con probabilidad alta a la mediana poblacional.

Aunque arriba imputamos la población completa, una manera más simple de hacer esto es la siguiente, de forma aproximada:

Tomamos una muestra con reemplazo de tamaño 8 de la muestra \(12,15,16,17,20,23,30,32\).
Calculamos la mediana de esta nueva remuestra.
Repetimos los dos anteriores una gran cantidad de veces para ver todos los posibles valores que podría tomar la mediana dada la información que tenemos de la muestra.

Código

remuestreo_med <- map_df(1:1000, function(rep){
  remuestra <- sample(muestra, length(muestra), replace = TRUE)
  tibble(rep = rep, mediana = median(remuestra))
})
ggplot(remuestreo_med, aes(x = mediana)) + geom_histogram(breaks = muestra) +
  geom_rug(data = tibble(muestra = muestra), aes(x= muestra),colour = "red")

Y podemos calcular por ejemplo intervalo de 80% tomando cuantiles de las medias remuestreadas:

Código

quantile(remuestreo_med$mediana, c(0.1, 0.9))

10% 90% 
 16  23

10.2 Ejemplo: precios de casas

Este proceso se llama remuestreo, y consiste en que una vez que tenemos un estimador para una cantidad poblacional, podemos remuestrar la muestra original, calcular nuestro estimador en cada caso, y así obtener una distribución de posibles valores que puede tener la mediana poblacional. Lo hacemos a continuación para nuestro ejemplo de la mediana de los precios de casas:

Código

remuestrear <- function(muestra){
  slice_sample(muestra_casas, prop = 1, replace = TRUE) |> 
  summarise(mediana = median(precio_miles))
}
remuestras <- map_df(1:10000, ~ remuestrear(muestra))

Código

ggplot(remuestras, aes(x = mediana)) +
  geom_histogram(breaks = muestra_casas |> pull(precio_miles)) +
  geom_rug(data = muestra_casas, aes(x = precio_miles), colour = "red")

Código

quantile(remuestras$mediana, c(0.051, 0.948), type = 8)

 5.1% 94.8% 
  146   235

y obtenemos un resultado idéntico al primer argumento de arriba. La cobertura es cercana ser exacta, y no afecta al resultado lo que sucede en las colas de la distribución, gracias a que estamos usando la mediana.

10.3 Remuestreo para otras estadísticas

Podemos repetir este método para otras estadísticas como la media. Para el ejemplo de casas:

remuestrear_media <- function(muestra){
  slice_sample(muestra, prop = 1, replace = TRUE) |> 
  summarise(media = mean(precio_miles))
}
remuestras_media <- map_df(1:10000, ~ remuestrear_media(muestra_casas))

Código

ggplot(remuestras_media) +
  geom_histogram(aes(x = media), binwidth = 1.5)

Obtenemos un intervalo del 90% para la media obteniendo los cuantiles de estas simulaciones:

Código

quantile(remuestras_media$media, c(0.05, 0.95)) |> round(1)

   5%   95% 
171.9 228.5

Intervalos de bootstrap

Motivamos el remuestreo (bootstrap) usando el ejemplo de estimación de la mediana población, donde es posible encontrar un intervalo exacto. El remuestreo se puede usar para muchas otras estadísticas difíciles de tratar teóricamente, y en general su comportamiento es bueno en relación a la cobertura.

Sin embargo, para estadísticas arbitrarias es difícil cumplir de manera exacta las garantías inferenciales (por ejemplo, cuando la distribución tiene una cola muy larga). Veremos como tratar este problema más adelante.