10.2 Intervalos de confianza bootstrap

En esta sección consideraremos el problema de construcción, mediante bootstrap, de un intervalo de confianza bilateral, con nivel de confianza \(1-\alpha\), para un parámetro \(\theta\) de la distribución \(F\). Una vez elegido el método bootstrap adecuado, teniendo en cuenta la información disponible en el contexto del que se trate, otro aspecto importante es el método para la construcción del intervalo de confianza bootstrap de forma que la probabilidad de cobertura sea lo más parecida posible al nivel nominal \(1-\alpha\).

Las diferencias entre los distintos métodos dependen del estadístico \(R\) empleado y de las suposiciones sobre su distribución. Como se comentó en la Sección 8.3.1, la función boot.ci() del paquete boot permite construir distintos tipos de intervalos de confianza dependiendo del parámetro type. En el Ejemplo 8.3 se ilustra la obtención de estimaciones por intervalo de confianza para la media empleando los distintos métodos bajo bootstrap uniforme (en el Capítulo ?? se incluyen ejemplos adicionales empleando bootstrap paramétrico y suavizado).

En esta sección se describirán brevemente los distintos métodos implementados en la función boot.ci(). Para un tratamiento más detallado, incluyendo los órdenes de los errores de cobertura, ver por ejemplo el Capítulo 4 de Cao y Fernández-Casal (2021) o el Capítulo 5 de Davison y Hinkley (1997).

10.2.1 Aproximación normal

Este método emplea las aproximaciones bootstrap del sesgo \(Sesgo^{\ast}\left( \hat{\theta}^{\ast} \right)\) y de la varianza \(Var^{\ast}\left( \hat{\theta}^{\ast} \right)\), y asume que la distribución del correspondiente estadístico studentizado es una normal estándar \[\frac{\hat{\theta} - Sesgo^{\ast}\left( \hat{\theta}^{\ast} \right) - \theta}{\sqrt{Var^{\ast}\left( \hat{\theta}^{\ast} \right)}} \underset{aprox}{\sim }\mathcal{N}\left( 0, 1 \right).\] De esta forma se obtiene la estimación por intervalo de confianza: \[\hat{I}_{norm}=\left( \hat{\theta} - Sesgo^{\ast}\left( \hat{\theta}^{\ast} \right) - z_{1-\alpha /2}\sqrt{Var^{\ast}\left( \hat{\theta}^{\ast} \right)},\hat{\theta} - Sesgo^{\ast}\left( \hat{\theta}^{\ast} \right) + z_{1 - \alpha /2}\sqrt{Var^{\ast}\left( \hat{\theta}^{\ast} \right)} \right).\]

Podemos obtener este intervalo de confianza estableciendo type = "norm" (o type = "all") en la llamada a la función boot.ci() (ver Ejemplo 8.3).

10.2.2 Método percentil directo

Este método se basa en la construcción del intervalo de confianza, mediante bootstrap, empleando como estadístico el estimador \[R = \hat{\theta}.\]

Una vez elegido el método de remuestreo, empleando un estimador, \(\hat{F}\,\), de la distribución poblacional, \(F\), la distribución en el muestreo de \(R = \hat{\theta}\) se aproxima directamente mediante la distribución bootstrap de \(R^{\ast}= \hat{\theta}^{\ast}\). A partir de las réplicas bootstrap del estimador aproximamos los cuantiles \(x_{\alpha /2}\) y \(x_{1-\alpha /2}\) (denotando por \(x_{\beta }\) el valor verificando \(P^{\ast}\left( R^{\ast }\leq x_{\beta } \right) =\beta\)), de forma que \[\begin{aligned} 1-\alpha &= 1-\frac{\alpha }{2}-\frac{\alpha }{2} = P^{\ast}\left( \hat{\theta}^{\ast}<x_{1-\alpha /2} \right) - P^{\ast}\left( \hat{\theta}^{\ast}\leq x_{\alpha /2} \right) \\ &= P^{\ast}\left( x_{\alpha /2}<\hat{\theta}^{\ast}<x_{1-\alpha /2} \right), \end{aligned}\] y asumimos que esto aproxima lo que ocurre con la distribución poblacional \[P\left( x_{\alpha /2} < \hat{\theta} < x_{1-\alpha /2} \right) \approx 1-\alpha.\] De donde se obtiene el intervalo de confianza bootstrap calculado por el método percentil directo \[\hat{I}_{perc}=\left( x_{\alpha /2}, x_{1-\alpha /2} \right).\]

Una ventaja de los intervalos construidos con este método es que son invariantes frente a transformaciones del estimador (en el caso de que fuese más adecuado trabajar en otra escala, no sería necesario conocer la transformación). Sin embargo, como se comentó en la Sección 8.1, la precisión puede verse seriamente afectada en el caso de estimadores sesgados.

Podemos obtener este intervalo de confianza estableciendo type = "perc" (o type = "all") en la llamada a la función boot.ci() (ver Ejemplo 8.3).

10.2.3 Método percentil básico

En este método se emplea como estadístico el estimador centrado (no estandarizado) \[R = \hat{\theta}-\theta.\] De forma análoga, la distribución en el muestreo de \(R\) se aproxima mediante la distribución bootstrap de \[R^{\ast}= \hat{\theta}^{\ast}-\theta \left( \hat{F} \right) = \hat{\theta}^{\ast}-\hat{\theta}.\] A partir de las réplicas bootstrap del estadístico se aproximan los cuantiles \(x_{\alpha /2}\) y \(x_{1-\alpha /2}\) tales que \[1-\alpha = P^{\ast}\left( x_{\alpha /2}<R^{\ast}<x_{1-\alpha /2} \right),\] tomándolo como aproximación de lo que ocurre con la distribución poblacional \[\begin{aligned} 1-\alpha &\approx P\left( x_{\alpha /2}<R<x_{1-\alpha /2} \right) \\ &= P\left( x_{\alpha /2} < \hat{\theta}-\theta < x_{1-\alpha /2} \right) \\ &= P\left( \hat{\theta} - x_{1-\alpha /2} < \theta <\hat{\theta} -x_{\alpha /2} \right). \end{aligned}\] De donde se obtiene el intervalo de confianza bootstrap calculado por el método percentil básico \[\hat{I}_{basic}=\left( \hat{\theta} - x_{1-\alpha /2},\hat{\theta} - x_{\alpha /2} \right).\]

Podemos obtener este intervalo de confianza estableciendo type = "basic" (o type = "all") en la llamada a la función boot.ci() (ver Ejemplo 8.3).

10.2.4 Método percentil-t

Este método bootstrap, construye un intervalo de confianza bootstrap a partir del estadístico studentizado: \[R = \frac{\hat \theta - \theta}{\sqrt{\widehat{Var}(\hat \theta)}}.\] Procediendo de modo análogo a los casos anteriores, se aproximan los cuantiles \(x_{\alpha /2}\) y \(x_{1-\alpha /2}\) tales que \[1-\alpha = P^{\ast}\left( x_{\alpha /2}<R^{\ast}<x_{1-\alpha /2} \right),\] a partir de los cuales se obtiene el intervalo de confianza bootstrap calculado por el método percentil-t (o percentil studentizado) \[\hat{I}_{stud}=\left( \hat{\theta} - x_{1-\alpha /2}\sqrt{\widehat{Var}(\hat \theta)},\hat{\theta} - x_{\alpha /2}\sqrt{\widehat{Var}(\hat \theta)} \right).\]

Si uno de los componentes del vector de estadísticos (por defecto el segundo) es una estimación de la varianza del estimador (por defecto el primer componente), podemos obtener este intervalo de confianza estableciendo type = "stud" (o type = "all") en la llamada a la función boot.ci() (ver Ejemplo 8.3). En caso de que el primer y segundo componente del vector de estadísticos no sean el estimador y su varianza estimada, respectivamente, habrá que emplear el argumento index, que debe ser un vector de longitud 2 con las posiciones correspondientes.

Hay una variante de este método, denominada percentil-t simetrizado, en la que se asume que la distribución del estadístico es simétrica (aunque no está implementado en boot.ci()). Asumiendo que esta suposición es correcta, podemos calcular los cuantiles de forma más eficiente (ya que dispondríamos del doble de información sobre las colas de la distribución). En lugar de tomar cuantiles que dejen colas iguales (\(\frac{\alpha }{2}\)) a la izquierda y a la derecha, respectivamente, se considera el valor \(x_{1-\alpha }\) que verifica \(P^{\ast}\left( \left\vert R^{\ast}\right\vert \leq x_{1-\alpha } \right) =1-\alpha\). Así se obtiene el intervalo de confianza bootstrap \[\hat{I}_{simstud}=\left( \hat{\theta} - x_{1-\alpha}\sqrt{\widehat{Var}(\hat \theta)},\hat{\theta} + x_{1-\alpha}\sqrt{\widehat{Var}(\hat \theta)} \right).\]

Ejemplo 10.2 (IC bootstrap para la media mediante el método percentil-*t* simetrizado)

Continuando con el Ejemplo 8.2 de inferencia sobre la media con varianza desconocida, podríamos obtener una estimación por intervalo de confianza del tiempo de vida medio de los microorganismos empleando el método bootstrap percentil-t simetrizado con el siguiente código:

muestra <- simres::lifetimes
n <- length(muestra)
alfa <- 0.05
x_barra <- mean(muestra)
cuasi_dt <- sd(muestra)

# Remuestreo
set.seed(1)
B <- 1000
estadistico_boot <- numeric(B)
for (k in 1:B) {
  remuestra <- sample(muestra, n, replace = TRUE)
  x_barra_boot <- mean(remuestra)
  cuasi_dt_boot <- sd(remuestra)
  estadistico_boot[k] <- sqrt(n) * abs(x_barra_boot - x_barra)/cuasi_dt_boot
}

# Aproximación bootstrap del pto crítico
pto_crit <- quantile(estadistico_boot, 1 - alfa)

# Construcción del IC
ic_inf_boot <- x_barra - pto_crit * cuasi_dt/sqrt(n)
ic_sup_boot <- x_barra + pto_crit * cuasi_dt/sqrt(n)
IC_boot <- c(ic_inf_boot, ic_sup_boot)
names(IC_boot) <- paste0(100*c(alfa/2, 1-alfa/2), "%")
IC_boot
##      2.5%     97.5% 
## 0.4334742 1.1771924

10.2.5 Método BCa

El método \(BCa\) (bias-corrected and accelerated) propuesto por Efron (1987) considera una transformación de forma que la distribución se aproxime a la normalidad, construye el intervalo en esa escala asumiendo normalidad y transforma el resultado a la escala original empleando la distribución bootstrap. El intervalo obtenido es de la forma: \[\hat{I}_{perc}=\left( x_{\alpha /2}, x_{1-\alpha /2} \right),\] donde
\[x_u = \hat G^{-1}\left(\Phi\left(z + \frac{z + z_u}{1-a(z+z_u)}\right) \right),\] siendo \(\hat G\) la distribución empírica de \(\hat{\theta}^{\ast}\), \(\Phi(z)\) la función de distribución de la normal estándar, \(z_u = \Phi^{-1}(u)\) el correspondiente cuantil y:

  • \(z = \Phi^{-1}(\hat G(\hat\theta))\) un factor de corrección de sesgo.
  • \(a\) la denominada constante aceleradora (o corrección de asimetría), que suele ser aproximada mediante jackknife.

Podemos obtener este intervalo de confianza estableciendo type = "bca" (o type = "all") en la llamada a la función boot.ci() (ver Ejemplo 8.3). Para más detalles ver Sección 5.3.2 de Davison y Hinkley (1997).

Ejercicio 10.2

Como continuación del ejemplo mostrado en la Sección 8.3.3, y de los Ejercicios 8.1 y 10.1, emplear el paquete boot para obtener estimaciones por intervalo de confianza del coeficiente de correlación lineal \(r\) entre prestige e income del conjunto de datos Prestige (mediante bootstrap uniforme multidimensional). En el caso del método percentil-t, como se indicó en el Ejercicio 8.1, considerar el estimador de la varianza: \[\widehat{Var}(r) = \frac{1 - r^2}{n - 2}.\]

Comparar los resultados con la aproximación paramétrica implementada en la función cor.test y descrita en la siguiente sección.

10.2.6 Ejemplo: IC bootstrap para el coeficiente de correlación

Supongamos de nuevo que queremos estudiar la correlación entre dos variables \(X\) e \(Y\) a partir del coeficiente de correlación lineal de Pearson: \[\rho =\frac{ Cov \left( X, Y \right) } { \sigma \left( X \right) \sigma \left( Y \right) },\] empleando como estimador el coeficiente de correlación muestral: \[r=\frac{\sum_{i=1}^{n}(x_i-\overline{x})(y_i-\overline{y})} {\sqrt{ \sum_{i=1}^{n}(x_i-\overline{x})^{2}} \sqrt{\sum_{i=1}^{n}(y_i-\overline{y})^{2}}}.\]

Para realizar inferencias sobre el coeficiente de correlación, como aproximación más simple, se puede considerar que la distribución muestral de \(r\) es aproximadamente normal (ver Ejercicio 8.1) y emplear el estadístico:

\[\begin{equation} \frac{r -\rho}{\sqrt{\frac{1 - r^2}{n - 2}}} \underset{aprox}{\sim } t_{n-2} \tag{10.1} \end{equation}\]

Pero esta aproximación solo sería válida en el caso de muestras grandes (o si la distribución bivariante de \((X, Y)\) es aproximadamente normal) cuando la correlación entre las variables es débil o moderada. En caso contrario la distribución muestral de \(r\) puede ser muy asimétrica y los resultados obtenidos con el estadístico anterior no ser muy adecuados (esto concuerda con lo observado en la Sección 8.3.3, al emplear bootstrap uniforme multidimensional para hacer inferencia sobre \(R = r -\rho\)). Para evitar este problema se suelen obtener intervalos de confianza para \(\rho\) empleando la transformación \(Z\) de Fisher (1915): \[Z = \frac{1}{2}\ln \left( \frac{1+r}{1-r} \right) = \operatorname{arctanh}(r),\] que es una transformación (aprox.) normalizadora y estabilizadora de la varianza. Suponiendo que \((X, Y)\) es normal bivariante y que hay independencia entre las observaciones: \[Z \sim \mathcal{N}\left( \frac{1}{2}\ln \left( \frac{1+\rho}{1-\rho} \right), \frac{1}{n-3} \right).\] El intervalo de confianza asintótico se obtiene empleando la aproximación normal tradicional en la escala \(Z\) y aplicando posteriormente la transformación inversa: \[r = \frac{\exp(2Z)-1}{\exp(2Z)+1} = \operatorname{tanh}(Z).\]

Esta aproximación está implementada en la función cor.test() del paquete base stat de R31, además de que también realiza el contraste \(H_0: \rho = 0\) empleando el estadístico (10.1).

Continuando con el ejemplo de la Sección 8.3.3 (y de los Ejercicios 8.1, 10.1 y 10.2), para obtener un intervalo de confianza para el coeficiente de correlación lineal entre las variables income y prestige del conjunto de datos Prestige, podríamos emplear el siguiente código:

data(Prestige, package="carData")
# with(Prestige, cor.test(income, prestige))
cor.test(Prestige$income, Prestige$prestige)
## 
##  Pearson's product-moment correlation
## 
## data:  Prestige$income and Prestige$prestige
## t = 10.224, df = 100, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6044711 0.7983807
## sample estimates:
##       cor 
## 0.7149057

También es de esperar que mejore la precisión de los intervalos de confianza bootstrap si se emplea una transformación que estabilice la varianza del estimador, especialmente en el caso del método basado en la aproximación normal y del bootstrap percentil básico. La función boot.ci() del paquete boot permite obtener intervalos de confianza calculados en una escala transformada del estadístico, mediante los parámetros:

  • h: función vectorial que define la transformación. Los intervalos se calculan en la escala de \(h(t)\) y se aplica la función inversa (si se especifica) para transformarlos a la escala original.

  • hinv: (opcional) función inversa de la transformación (si no se especifica solo se calculan los intervalos en la escala transformada).

  • hdot: (opcional en el método percentil o básico) función derivada de la transformación (empleada por algunos métodos para aproximar la varianza en la escala transformada mediante el método delta).

Por ejemplo, para considerar la transformación \(Z\) de Fisher en este caso, se podría emplear el siguiente código:

library(boot)

statistic <- function(data, i){
  remuestra <- data[i, ]
  cor(remuestra$income, remuestra$prestige)
}

set.seed(1)
res.boot <- boot(Prestige, statistic, R = 1000)

h <- function(t) atanh(t)
hdot <- function(t) 1/(1 - t^2)
hinv <- function(t) tanh(t)

# boot.ci(res.boot, type = "norm", h = h)
boot.ci(res.boot, type = "norm", h = h, hdot = hdot, hinv = hinv)
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = res.boot, type = "norm", h = h, hdot = hdot, 
##     hinv = hinv)
## 
## Intervals : 
## Level      Normal        
## 95%   ( 0.6016,  0.7858 )  
## Calculations on Transformed Scale;  Intervals on Original Scale

Esto sería en principio preferible a trabajar en la escala original, ya que la distribución bootstrap en la escala transformada se aproximaría más a la normalidad:

ht <- h(res.boot$t)
hist(ht, freq = FALSE, breaks = "FD", main = "")
curve(dnorm(x, mean=mean(ht), sd=sd(ht)), lty = 2, add = TRUE)
Distribución bootstrap en la escala transformada.

Figura 10.1: Distribución bootstrap en la escala transformada.

  1. Se puede obtener el código tecleando en la consola stats:::cor.test.default.↩︎