6.3 Intervalos de confianza
Hasta ahora hemos discutido la idea detrás del bootstrap y como se puede usar para estimar errores estándar. Comenzamos con el error estándar pues es la manera más común para describir la precisión de una estadística.
En términos generales, esperamos que \(\bar{x}\) este a una distancia de \(\mu_P\) menor a un error estándar el 68% del tiempo, y a menos de 2 errores estándar el 95% del tiempo.
Estos porcentajes están basados el teorema central del límite que nos dice que bajo ciertas condiciones (bastante generales) de \(P\) la distribución de \(\bar{x}\) se aproximará a una distribución normal:
\[\bar{x} \overset{\cdot}{\sim} N(\mu_P,\sigma_P^2/n)\]
Veamos algunos ejemplos de como funciona el Teorema del Límite Central, buscamos ver como se aproxima la distribución muestral de la media (cuando las observaciones provienen de distintas distribuciones) a una Normal conforme aumenta el tamaño de muestra. Para esto, aproximamos la distribución muestral de la media usando simulación de la población.
Vale la pena observar que hay distribuciones que requieren un mayor tamaño de muestra \(n\) para lograr una buena aproximación (por ejemplo la log-normal), ¿a qué se debe esto?
Para la opción de Elecciones tenemos una población de tamaño \(N=143,437\) y el objetivo es estimar la media del tamaño de la lista nominal de las casillas (datos de las elecciones presidenciales de 2012). Podemos ver como mejora la aproximación Normal de la distribución muestral conforme aumenta el tamaño de muestra \(n\); sin embargo, también sobresale que no es necesario tomar una muestra demasiado grande (\(n = 60\) ya es razonable).
En lo que sigue veremos distintas maneras de construir intervalos de confianza usando bootstrap.
Un intervalo de confianza \((1-2\alpha)\)% para un parámetro \(\theta\) es un intervalo \((a,b)\) tal que \(P(a \le \theta \le b) = 1-2\alpha\) para todo \(\theta \in \Theta\).
Y comenzamos con la versión bootstrap del intervalo más popular.
- Intervalo Normal con error estándar bootstrap. El intervalo para \(\hat{\theta}\) con un nivel de confianza de \(100\cdot(1-2\alpha)\%\) se define como:
\[(\hat{\theta}-z^{(1-\alpha)}\cdot \hat{se}_B, \hat{\theta}+z^{(1-\alpha)}\cdot \hat{se})\].
donde \(z^{(\alpha)}\) denota el percentil \(100\cdot \alpha\) de una distribución \(N(0,1)\).
este intervalo está soportado por el Teorema Central del Límite, sin embargo, no es adecuado cuando \(\hat{\theta}\) no se distribuye aproximadamente Normal.
Ejemplo: kurtosis
Supongamos que queremos estimar la kurtosis de una base de datos que consta de 799 tiempos de espera entre pulsasiones de un nervio (Cox, Lewis 1976).
\[\hat{\theta} = t(P_n) =\frac{1/n \sum_{i=1}^n(x_i-\hat{\mu})^3}{\hat{\sigma}^3}\]
library(ACSWR)
data("nerve")
head(nerve)
#> [1] 0.21 0.03 0.05 0.11 0.59 0.06
kurtosis <- function(x){
n <- length(x)
1 / n * sum((x - mean(x)) ^ 3) / sd(x) ^ 3
}
theta_hat <- kurtosis(nerve)
theta_hat
#> [1] 1.757943
kurtosis_boot <- function(x){
x_boot <- sample(x, replace = TRUE)
kurtosis(x_boot)
}
B <- 10000
kurtosis <- rerun(B, kurtosis_boot(nerve)) %>%
flatten_dbl()Usando el intervalo normal tenemos:
li_normal <- round(theta_hat - 1.96 * sd(kurtosis), 2)
ls_normal <- round(theta_hat + 1.96 * sd(kurtosis), 2)
c(li_normal, ls_normal)
#> [1] 1.44 2.08Una modificación común del intervalo normal es el intervalo t, estos intervalos son mejores en caso de muestras pequeñas (\(n\) chica).
- Intervalo \(t\) con error estándar bootstrap. Para una muestra de tamaño \(n\) el intervalo \(t\) con un nivel de confianza de \(100\cdot(1-2\alpha)\%\) se define como:
\[(\hat{\theta}-t^{(1-\alpha)}_{n-1}\cdot \hat{se}_B, \hat{\theta}+t^{(1-\alpha)}_{n-1}\cdot \hat{se}_B)\].
donde \(t^{(\alpha)}_{n-1}\) denota denota el percentil \(100\cdot \alpha\) de una distribución \(t\) con \(n-1\) grados de libertad.
n_nerve <- length(nerve)
li_t <- round(theta_hat + qt(0.025, n_nerve - 1) * sd(kurtosis), 2)
ls_t <- round(theta_hat - qt(0.025, n_nerve - 1) * sd(kurtosis), 2)
c(li_t, ls_t)
#> [1] 1.44 2.08Los intervalos normales y \(t\) se valen de la estimación bootstrap del error estándar; sin embargo, el bootstrap se puede usar para estimar la función de distribución de \(\hat{\theta}\) por lo que no es necesario hacer supuestos distribucionales para \(\hat{\theta}\) sino que podemos estimarla como parte del proceso de construir intervalos de confianza.
Veamos un histograma de las replicaciones bootstrap de \(\hat{\theta}^*\)
library(gridExtra)
nerve_kurtosis <- tibble(kurtosis)
hist_nerve <- ggplot(nerve_kurtosis, aes(x = kurtosis)) +
geom_histogram(binwidth = 0.05, fill = "gray30") +
geom_vline(xintercept = c(li_normal, ls_normal, theta_hat),
color = c("black", "black", "red"), alpha = 0.5)
qq_nerve <- ggplot(nerve_kurtosis) +
geom_abline(color = "red", alpha = 0.5) +
stat_qq(aes(sample = kurtosis), dparams = list(mean = mean(kurtosis), sd = sd(kurtosis)))
grid.arrange(hist_nerve, qq_nerve, ncol = 2, newpage = FALSE)
En el ejemplo anterior el supuesto de normalidad parece razonable, veamos como se comparan los cuantiles de la estimación de la distribución de \(\hat{\theta}\) con los cuantiles de una normal:
comma(q_kurt <- quantile(kurtosis,
probs = c(0.025, 0.05, 0.1, 0.9, 0.95, 0.975)))
comma(qnorm(p = c(0.025, 0.05, 0.1, 0.9, 0.95, 0.975), mean = theta_hat,
sd = sd(kurtosis)))
#> 2.5% 5% 10% 90% 95% 97.5%
#> "1.4" "1.5" "1.5" "2.0" "2.0" "2.1"
#> [1] "1.4" "1.5" "1.5" "2.0" "2.0" "2.1"Esto sugiere usar los cuantiles del histograma bootstrap para definir los límites de los intervalos de confianza:
- Percentiles. Denotemos por \(G\) la función de distribución acumulada de \(\hat{\theta}^*\) el intervalo percentil de \(1-2\alpha\) se define por los percentiles \(\alpha\) y \(1-\alpha\) de \(G\) \[(\theta^*_{\%,inf}, \theta^*_{\%,sup}) = (G^{-1}(\alpha), G^{-1}(1-\alpha))\] Por definición \(G^{-1}(\alpha)=\hat{\theta}^*(\alpha)\), esto es, el percentil \(100\cdot \alpha\) de la distribución bootstrap, por lo que podemos escribir el intervalo bootstrap como \[(\theta^*_{\%,inf}, \theta^*_{\%,sup})=(\hat{\theta}^*(\alpha),\hat{\theta}^*(1-\alpha))\]
ggplot(arrange(nerve_kurtosis, kurtosis)) +
stat_ecdf(aes(x = kurtosis)) +
geom_segment(data = data_frame(x = c(-Inf, -Inf, q_kurt[c(1, 6)]),
xend = q_kurt[c(1, 6, 1, 6)], y = c(0.025, 0.975, 0, 0),
yend = c(0.025, 0.975, 0.025, 0.975)), aes(x = x, xend = xend, y = y,
yend = yend), color = "red", size = 0.4, alpha = 0.5) +
labs(x = "Cuantiles muestrales", y = "ecdf")
Las expresiones de arriba hacen referencia a la situación bootstrap ideal donde el número de replicaciones bootstrap es infinito, en la práctica usamos aproximaciones. Y se procede como sigue:
Intervalo percentil:
Generamos B muestras bootstrap independientes \(\textbf{x}^{*1},..., \textbf{x}^{*B}\) y calculamos las replicaciones \(\hat{\theta}^{*b}=s(x^{*b}).\)
- Sea \(\hat{\theta}^{*}_B(\alpha)\) el percentil \(100\cdot\alpha\) de la distribución empírica de \(\hat{\theta}^{*}\), y \(\hat{\theta}^{*}_B(1-\alpha)\) el correspondiente al percentil \(100\cdot (1-\alpha)\), escribimos el intervalo de percentil \(1-2\alpha\) como \[(\theta^*_{\%,inf}, \theta^*_{\%,sup})\approx(\hat{\theta}^*_B(\alpha),\hat{\theta}^*_B(1-\alpha))\]
ls_per <- round(quantile(kurtosis, probs = 0.975), 2)
li_per <- round(quantile(kurtosis, probs = 0.025), 2)
stringr::str_c(li_normal, ls_normal, sep = ",")
stringr::str_c(li_per, ls_per, sep = ",")
#> [1] "1.44,2.08"
#> [1] "1.43,2.07"Si la distribución de \(\hat{\theta}^*\) es aproximadamente normal, entonces los intervalos normales y de percentiles serán similares.
Con el fin de comparar los intervalos creamos un ejemplo de simulación (ejemplo tomado de Efron and Tibshirani (1993)), generamos una muestra de tamaño 10 de una distribución normal estándar, supongamos que el parámetro de interés es \(e^{\mu}\) donde \(\mu\) es la media poblacional.
set.seed(137612)
x <- rnorm(10)
boot_sim_exp <- function(){
x_boot <- sample(x, size = 10, replace = TRUE)
exp(mean(x_boot))
}
theta_boot <- rerun(1000, boot_sim_exp()) %>% flatten_dbl()
theta_boot_df <- data_frame(theta_boot)
hist_emu <- ggplot(theta_boot_df, aes(x = theta_boot)) +
geom_histogram(fill = "gray30", binwidth = 0.08)
qq_emu <- ggplot(theta_boot_df) +
geom_abline(color = "red", alpha = 0.5) +
stat_qq(aes(sample = theta_boot),
dparams = list(mean = mean(theta_boot), sd = sd(theta_boot)))
grid.arrange(hist_emu, qq_emu, ncol = 2, newpage = FALSE)
La distribución empírica de \(\hat{\theta}^*\) es asimétrica, por lo que no esperamos que coincidan los intervalos.
# Normal
round(exp(mean(x)) - 1.96 * sd(theta_boot), 2)
#> [1] 0.36
round(exp(mean(x)) + 1.96 * sd(theta_boot), 2)
#> [1] 1.6
#Percentil
round(quantile(theta_boot, prob = 0.025), 2)
#> 2.5%
#> 0.53
round(quantile(theta_boot, prob = 0.975), 2)
#> 97.5%
#> 1.79La inspección del histograma deja claro que la aproximación normal no es conveniente en este caso, veamos que ocurre cuando aplicamos la transformación logarítmica.
hist_log <- ggplot(data_frame(theta_boot), aes(x = log(theta_boot))) +
geom_histogram(fill = "gray30", binwidth = 0.08)
qq_log <- ggplot(data_frame(theta_boot)) +
geom_abline(color = "red", alpha = 0.5) +
stat_qq(aes(sample = log(theta_boot)),
dparams = list(mean = mean(log(theta_boot)), sd = sd(log(theta_boot))))
grid.arrange(hist_log, qq_log, ncol = 2, newpage = FALSE)
Y los intervalos se comparan:
# Normal
round(mean(x) - 1.96 * sd(log(theta_boot)), 2)
#> [1] -0.63
round(mean(x) + 1.96 * sd(log(theta_boot)), 2)
#> [1] 0.58
#Percentil
round(quantile(log(theta_boot), prob = 0.025), 2)
#> 2.5%
#> -0.63
round(quantile(log(theta_boot), prob = 0.975), 2)
#> 97.5%
#> 0.58La transformación logarítmica convierte la distribución de \(\hat{\theta}\) en normal y por tanto los intervalos de \(\hat{\phi}^*=log(\hat{\theta}^*)\) son similares. La forma normal no es sorprendente pues \(\hat{\phi}^*=\bar{x}^*\).
Si mapeamos los intervalos normales calculados para \(log(\hat{\theta}^*)\) de regreso a la escala de \(\theta\) obtenemos intervalos similares a los calculados para \(\hat{\theta}^*\) usando percentiles:
exp(round(mean(x) - 1.96 * sd(log(theta_boot)), 2))
#> [1] 0.5325918
exp(round(mean(x) + 1.96 * sd(log(theta_boot)), 2))
#> [1] 1.786038Podemos ver que el método de aplicar una transformación, calcular intervalos usando la normal y aplicar la transformación inversa para volver a la escala original genera intervalos de confianza atractivos, el problema con este método es que requiere que conozcamos la transformación adecuada para cada parámetro.
Por otra parte, podemos pensar en el método del percentil como un algoritmo que incorpora la transformación de manera automática.
Lema. Supongamos que la transformación \(\hat{\phi}=m(\hat{\theta})\) normaliza la distribución de \(\hat{\theta}\) de manera perfecta, \[\hat{\phi} \approx N(\phi, c^2)\] para alguna desviación estándar \(c\). Entonces el intervalo de percentil basado en \(\hat{\theta}\) es igual a \[(m^{-1} (\hat{\phi}-z^{(1-\alpha)}c), m^{-1}(\hat{\phi}-z^{(\alpha)}c))\]
Se dice que el intervalo de confianza de percentiles es invariante a transformaciones.
Existen otras alternativas al método del percentil y cubren otras fallas del intervalo normal. Por ejemplo, hay ocasiones en que \(\hat{\theta}\) tiene una distribución normal sesgada: \[\hat{\theta} \approx N(\theta + sesgo, \hat{se}^2)\]
en este caso no existe una transformación \(m(\theta)\) que arregle el intervalo.
- Intervalos acelerados y corregidos por sesgo. Esta es una versión mejorada del intervalo de percentil, la denotamos \(BC_{a}\) (bias-corrected and accelerated).
Usaremos un ejemplo de Efron and Tibshirani (1993), los datos constan de los resultados en dos pruebas espaciales de 26 niños con discapacidad neurológico. Supongamos que queremos calcular un intervalo de confianza de 90% para \(\theta=var(A)\).
El estimador plugin es: \[\hat{\theta}=\sum_{i=1}^n(A_i-\bar{A})^2/n\] notemos que el estimador plug-in es ligeramente menor que el estimador usual insesgado: \[\hat{\theta}=\sum_{i=1}^n(A_i-\bar{A})^2/(n-1)\]
library(bootstrap)
spatial
#> A B
#> V1 48 42
#> V2 36 33
#> V3 20 16
#> V4 29 39
#> V5 42 38
#> V6 42 36
#> V7 20 15
#> V8 42 33
#> V9 22 20
#> V10 41 43
#> V11 45 34
#> V12 14 22
#> V13 6 7
#> V14 0 15
#> V15 33 34
#> V16 28 29
#> V17 34 41
#> V18 4 13
#> V19 32 38
#> V20 24 25
#> V21 47 27
#> V22 41 41
#> V23 24 28
#> V24 26 14
#> V25 30 28
#> V26 41 40
ggplot(spatial) +
geom_point(aes(A, B))
El estimador plug-in de \(\theta\) es
Notemos que es ligeramente menor que el estimador insesgado:
El método \(BC_{a}\) corrige el sesgo de manera automática, lo cuál es una de sus prinicipales ventajas comparado con el método del percentil.
Los extremos en los intervalos \(BC_{a}\) están dados por percentiles de la distribución bootstrap, los percentiles usados dependen de dos números \(\hat{a}\) y \(\hat{z}_0\), que se denominan la aceleración y la corrección del sesgo: \[BC_a : (\hat{\theta}_{inf}, \hat{\theta}_{sup})=(\hat{\theta}^*(\alpha_1), \hat{\theta}^*(\alpha_2))\] donde \[\alpha_1= \Phi\bigg(\hat{z}_0 + \frac{\hat{z}_0 + z^{(\alpha)}}{1- \hat{a}(\hat{z}_0 + z^{(\alpha)})}\bigg)\] \[\alpha_2= \Phi\bigg(\hat{z}_0 + \frac{\hat{z}_0 + z^{(1-\alpha)}}{1- \hat{a}(\hat{z}_0 + z^{(1-\alpha)})}\bigg)\] y \(\Phi\) es la función de distribución acumulada de la distribución normal estándar y \(z^{\alpha}\) es el percentil \(100 \cdot \alpha\) de una distribución normal estándar.
Notemos que si \(\hat{a}\) y \(\hat{z}_0\) son cero entonces \(\alpha_1=\alpha\)
y \(\alpha_2=1-\alpha\), obteniendo así los intervalos de percentiles.
El valor de la corrección por sesgo \(\hat{z}_0\) se obtiene de la
propoción de de replicaciones bootstrap menores a la estimación original
\(\hat{\theta}\),
\[z_0=\Phi^{-1}\bigg(\frac{\#\{\hat{\theta}^*(b) < \hat{\theta} \} }{B} \bigg)\]
a grandes razgos \(\hat{z}_0\) mide la mediana del sesgo de \(\hat{\theta}^*\), esto es, la discrepancia entre la mediana de \(\hat{\theta}^*\) y \(\hat{\theta}\) en unidades normales.
Por su parte la aceleración \(\hat{a}\) se refiere a la tasa de cambio del error estándar de \(\hat{\theta}\) respecto al verdadero valor del parámetro \(\theta\). La aproximación estándar usual \(\hat{\theta} \approx N(\theta, se^2)\) supone que el error estándar de \(\hat{\theta}\) es el mismo para toda \(\hat{\theta}\), esto puede ser poco realista, en nuestro ejemplo, donde \(\hat{\theta}\) es la varianza si los datos provienen de una normal \(se(\hat{\theta})\) depende de \(\theta\). Una manera de calcular \(\hat{a}\) es
\[\hat{a}=\frac{\sum_{i=1}^n (\hat{\theta}(\cdot) - \hat{\theta}(i))^3}{6\{\sum_{i=1}^n (\hat{\theta}(\cdot) - \hat{\theta}(i))^2\}^{3/2}}\]
Los intervalos \(BC_{a}\) tienen 2 ventajas teóricas:
Respetan transformaciones, esto nos dice que los extremos del intervalo se transforman de manera adecuada si cambiamos el parámetro de interés por una función del mismo.
Su exactitud, los intervalos \(BC_{a}\) tienen precisión de segundo orden, esto es, los errores de cobertura se van a cero a una tasa de 1/n.
Los intervalos \(BC_{a}\) están implementados en el paquete boot (boot.ci()) y
en el paquete bootstrap (bcanon()). La desventaja de los intervalos \(BC_{a}\) es
que requieren intenso cómputo estadístico, de acuerdo a Efron and Tibshirani (1993) al
menos \(B= 1000\) replicaciones son necesarias para reducir el error de muestreo.
Ante esto surgen los intervalos ABC (approximate bootstrap confidence intervals), que es un método para aproximar \(BC_{a}\) analíticamente (usando expansiones de Taylor), estos intervalos requieren que la estadística \(\hat{\theta} = s(x)\) este definida de manera suave sobre x (la mediana, por ejemplo, no es suave).
Usando la implementación del paquete bootstrap:
var_sesgada <- function(x) sum((x - mean(x)) ^ 2) / length(x)
bcanon(x = spatial[, 1], nboot = 2000, theta = var_sesgada, alpha = c(0.025, 0.975))
#> $confpoints
#> alpha bca point
#> [1,] 0.025 103.8402
#> [2,] 0.975 274.0533
#>
#> $z0
#> [1] 0.1383042
#>
#> $acc
#> [1] 0.06124012
#>
#> $u
#> [1] 164.3936 176.7200 174.5184 178.3776 172.0544 172.0544 174.5184 172.0544
#> [9] 175.9584 173.0400 168.5984 168.2016 155.1200 141.8144 177.9296 178.2816
#> [17] 177.6096 151.0176 178.1664 177.0656 165.8784 173.0400 177.0656 177.8400
#> [25] 178.3904 173.0400
#>
#> $call
#> bcanon(x = spatial[, 1], nboot = 2000, theta = var_sesgada, alpha = c(0.025,
#> 0.975))
Comapara el intervalo anterior con los intervalos
normal y de percentiles.
Otros intervalos basados en bootstrap incluyen los intervalos pivotales y los intervalos bootstrap-t. Sin embargo, BC y ABC son mejores alternativas.
- Intervalos pivotales. Sea \(\theta=s(P)\) y \(\hat{\theta}=s(P_n)\) definimos el pivote \(R=\hat{\theta}-\theta\). Sea \(H(r)\) la función de distribución acumulada del pivote: \[H(r) = P(R<r)\]
Definimos \(C_n^*=(a,b)\) donde: \[a=\hat{\theta}-H^{-1}(1-\alpha), b=\hat{\theta}-H^{-1}(\alpha)\] \(C_n^*\) es un intervalo de confianza de \(1-2\alpha\) para \(\theta\); sin embargo, \(a\) y \(b\) dependen de la distribución desconocida \(H\), la podemos estimar usando bootstrap: \[\hat{H}(r)=\frac{1}{B}\sum_{b=1}^B I(R^*_b \le r)\] y obtenemos \[C_n=(2\hat{\theta} - \hat{\theta}^*_{1-\alpha}, 2\hat{\theta} + \hat{\theta}^*_{1-\alpha})\]
Exactitud en intervalos de confianza. Un intervalo de \(95%\) de confianza exacto no captura el verdadero valor \(2.5%\) de las veces, en cada lado.
Un intervalo que sub-cubre un lado y sobre-cubre el otro es sesgado.
Los intervalos estándar y de percentiles tienen exactitud de primer orden: los errores de cobertura se van a cero a una tasa de \(1/\sqrt{n}\). Suelen ser demasido estrechos resultando en cobertura real menor a la nominal, sobretodo con muestras chicas.
Los intervalos \(BC_a\) tienen exactitud de segundo orden: los errores de cobertura se van a cero a una tasa de \(1/n\).
A pesar de que los intervalos \(BC_a\) pueden ser superiores a los intervalos normales y de percentiles, en la práctica es más común utilizar intervalos normales o de percentiles pues su implementación es más sencilla y son adecuados para un gran número de casos.
Referencias
Efron, Bradley, and Robert J. Tibshirani. 1993. An Introduction to the Bootstrap. Monographs on Statistics and Applied Probability 57. Boca Raton, Florida, USA: Chapman & Hall/CRC.