Bootstrap

Introducción

• El bootstrap es un método de simulación mediante remuestreo que se puede usar para estimar desviaciones estándar de un estadístico, intervalos de confianza de ese estadístico o incluso en contraste de hipótesis.

• El término bootstrap deriva de la frase to pull oneself up by one’s bootstrap.

• Se basa en el libro del siglo XVIII “Las aventuras del barón Munchausen” de Rudolph E. Raspe.

• El barón había caído en el fondo de un profundo lago y se le ocurrió escapar tirando de los cordones de sus propias botas…

Esquema del bootstrap

• Supongamos que se observa una muestra aleatoria simple \(x_1,x_2,\dots,x_n\) que se denota como el vector \(\textbf{x}\) en el que se calcula un estadístico dado \(s(\textbf{x})\).

• Por ejemplo, \(\textbf{x}\) es un grupo de datos observados y \(s(\textbf{x})\) es la mediana muestral.

• Para hacer bootstrap, se define una muestra bootstrap \(\textbf{x}^{\ast}=(x_1^{\ast},x_2^{\ast},\dots,x_n^{\ast})\) que se obtiene muestreando \(n\) veces con reemplazamiento entre los datos originales \(x_1,x_2,\dots,x_n\).

• Por ejemplo si \(n=7\) una posible muestra bootstrap podría ser \[ \textbf{x}^{\ast}=(x_5,x_4,x_7,x_2,x_4,x_1,x_5) \]

Algoritmo del bootstrap

• El algoritmo empieza generando un número elevado de muestras bootstrap \(\textbf{x}^{\ast 1},\textbf{x}^{\ast 2},\dots,\textbf{x}^{\ast B}\) cada una de tamaño \(n\).

• Los valores típicos de \(B\) para calcular los errores estándar suelen estar en un rango entre 500 y 1000.

• Para cada muestra bootstrap \(b=1,2,\dots,B\) se calcula el estadístico \(s(\textbf{x}^{\ast b})\).

• Por ejemplo si el estadístico \(s\) es la mediana, se calcula la mediana para cada muestra bootstrap.

• El estimador bootstrap del error estándar es la desviación estándar de las \(B\) muestras bootstrap \[ \widehat{se}_{Boot} = \sqrt{ \frac{1}{B-1} \sum_{b=1}^B (s(\textbf{x}^{\ast b}) - s(\cdot))^2 } \] donde \[ s(\cdot) = \frac{1}{B} \sum_{b=1}^B s(\textbf{x}^{\ast b}) \]

Desviación estándar de la mediana

Apliquemos el procedimiento al cálculo de la desviación estándar de la mediana de unos datos.

x <- faithful$waiting
n <- length(x)
B <- 1000
set.seed(123)

matriz_aux <- matrix(0, ncol=B, nrow=n)
for (i in 1:B) matriz_aux[ ,i] <- sample(x, rep=TRUE) 
bootmed <- apply(matriz_aux, 2, median)

o también

bootmed <- apply(matrix(sample(x, rep=TRUE, B*n), nrow=n), 2, median)

o mejor

bootmed <- replicate(B, median(sample(x, rep = TRUE)))

Dos gráficos:

hist(bootmed)

plot(density(bootmed))

La desviación estándar de la mediana es

sd(bootmed)

## [1] 0.9921368

y el intervalo de confianza

quantile(bootmed, c(0.025, 0.975))

##  2.5% 97.5% 
##  73.5  77.0

hist(bootmed)
abline(v=quantile(bootmed, 0.025), lty=2, lwd=2, col="blue")
abline(v=quantile(bootmed, 0.975), lty=2, lwd=2, col="blue")

Comparación de medianas

Supongamos que se ha observado a dos conjuntos de ratas que han recibido un tratamiento o no (grupo control) para prolongar su supervivencia después de una cirugía invasiva.

Los datos son:

tratadas <- c(94 , 197 , 16 , 38 , 99 , 141 , 23)
control <- c(52 , 104 , 146 , 10 , 51 , 30 , 40 , 27 , 46)

ecdf1 <- ecdf(tratadas)
ecdf2 <- ecdf(control)
plot(ecdf2, verticals = TRUE, do.points = FALSE , col ="blue", xlab = "dias",
ylab = " Distribución Empírica", main = " Controles (azul) / Tratados (naranja)")
plot(ecdf1, verticals = TRUE , col = "orange", do.points = FALSE , add = TRUE)

La variable observada no es, a priori, normal y la escasez de datos sugiere no utilizar el clásico test \(t\) de Student.

El test de Wilcoxon es

wilcox.test(tratadas,control)

## 
##  Wilcoxon rank sum exact test
## 
## data:  tratadas and control
## W = 36, p-value = 0.6806
## alternative hypothesis: true location shift is not equal to 0

Vamos a utilizar el bootstrap con la diferencia de medianas.

median(tratadas) - median(control)

## [1] 48

median_dif <- replicate(B, median(sample(tratadas, replace = TRUE)) - median(sample(control, replace = TRUE )))
sd(median_dif)

## [1] 40.92844

y podemos calcular el intervalo de confianza para la diferencia de medianas

quantile(median_dif, c(0.025,0.975))

##  2.5% 97.5% 
##   -29   101

o gráficamente

hist(median_dif, breaks = 20)
abline(v=quantile(median_dif, 0.025), lty=2, lwd=2, col="blue")
abline(v=quantile(median_dif, 0.975), lty=2, lwd=2, col="blue")

Vemos que el intervalo de confianza contiene al cero. No podemos rechazar la hipótesis de igualdad de medianas.

Intervalo de confianza para el coeficiente de determinación

Tomemos una regresión múltiple con los datos del data.frame prostate del paquete faraway.

library(faraway)
data("prostate")
fit <- lm(lpsa ~ ., data=prostate)
summary(fit)$r.squared

## [1] 0.6547541

En primer lugar definimos una función con algunos parámetros generales que permite calcular el coeficiente de determinación \(R^2\) para un remuestreo de los datos.

rsq <- function(formula, data, indices) {
  d <- data[indices,] # allows boot to select sample 
  fit <- lm(formula, data=d)
  summary(fit)$r.squared
} 

Y ahora podemos aplicar el método bootstrap

coef_det <- replicate(B, rsq(lpsa ~ ., data=prostate, indices=sample(1:dim(prostate)[1], rep=TRUE)))
quantile(coef_det, c(0.025,0.975))

##      2.5%     97.5% 
## 0.5638416 0.7844732

hist(coef_det, breaks = 20)
abline(v=quantile(coef_det, 0.025), lty=2, lwd=2, col="blue")
abline(v=quantile(coef_det, 0.975), lty=2, lwd=2, col="blue")

Vamos a utilizar el paquete boot tal como se explica en Quick R

# Bootstrap 95% CI for R-Squared
library(boot)

## 
## Attaching package: 'boot'

## The following objects are masked from 'package:faraway':
## 
##     logit, melanoma

# bootstrapping with 1000 replications 
results <- boot(data=prostate, statistic=rsq, R=1000, formula=lpsa ~ .)

Veamos el resultado

results

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = prostate, statistic = rsq, R = 1000, formula = lpsa ~ 
##     .)
## 
## 
## Bootstrap Statistics :
##      original     bias    std. error
## t1* 0.6547541 0.02741787  0.05743678

plot(results)

Y el intervalo de confianza al 95% es

boot.ci(results, type="bca")

## Warning in norm.inter(t, adj.alpha): extreme order statistics used as endpoints

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = results, type = "bca")
## 
## Intervals : 
## Level       BCa          
## 95%   ( 0.5050,  0.7333 )  
## Calculations and Intervals on Original Scale
## Warning : BCa Intervals used Extreme Quantiles
## Some BCa intervals may be unstable

Bootstrap con los residuos de un modelo lineal

Tomemos un caso de regresión lineal simple

library(MASS)
data(cats)
str(cats)

## 'data.frame':    144 obs. of  3 variables:
##  $ Sex: Factor w/ 2 levels "F","M": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Bwt: num  2 2 2 2.1 2.1 2.1 2.1 2.1 2.1 2.1 ...
##  $ Hwt: num  7 7.4 9.5 7.2 7.3 7.6 8.1 8.2 8.3 8.5 ...

n <- dim(cats)[1]
fit <- lm(Hwt ~ Bwt, data=cats)
summary(fit)$r.squared

## [1] 0.6466209

with(cats, plot(Bwt,Hwt))
abline(fit)

Ahora consideremos un bootstrap sobre los residuos del ajuste anterior

set.seed(123)
ynew <- predict(fit) + residuals(fit)[sample(n, rep=TRUE)]
fit2 <- lm(ynew ~ Bwt, data=cats)
with(cats, plot(Bwt,Hwt,col="red"))
abline(fit, col="red")
with(cats, points(Bwt,ynew, col="blue"))
abline(fit2, col="blue")

Si lo repetimos \(B=1000\) veces obtenemos

rsq2 <- numeric(B)
for(i in 1:B){
  ynew <- predict(fit) + residuals(fit)[sample(n, rep=TRUE)]
  fit2 <- lm(ynew ~ Bwt, data=cats)
  rsq2[i] <- summary(fit2)$r.squared
}  
quantile(rsq2, c(0.025,0.975))

##      2.5%     97.5% 
## 0.5701689 0.7217259

hist(rsq2, breaks = 20)
abline(v=summary(fit)$r.squared, lwd=2, col="red")
abline(v=quantile(rsq2, 0.025), lty=2, lwd=2, col="blue")
abline(v=quantile(rsq2, 0.975), lty=2, lwd=2, col="blue")

Ahora podemos utilizar este método de bootstrap con los residuos para calcular intervalos de confianza para los coeficientes del modelo.

bcoef <- matrix(0, nrow=B, ncol=2)
for(i in 1:B){
  ynew <- predict(fit) + residuals(fit)[sample(n, rep=TRUE)]
  fit2 <- lm(ynew ~ Bwt, data=cats)
  bcoef[i, ] <- coef(fit2)
}  
colnames(bcoef) <- names(coef(fit))
t(apply(bcoef,2,function(x) quantile(x, c(0.025,0.975))))

##                  2.5%     97.5%
## (Intercept) -1.733736 0.9487219
## Bwt          3.576271 4.5503629

Que se acerca mucho al resultado del intervalo de confianza bajo las condiciones de Gauss-Markov y suponiendo la normalidad de los errores:

confint(fit)

##                 2.5 %   97.5 %
## (Intercept) -1.725163 1.011838
## Bwt          3.539343 4.528782

Intervalos de confianza para los coeficientes de la regresión LTS

Utilizando el algoritmo que hemos visto en el apartado anterior podemos calcular intervalos de confianza para la regresión Least Trimmed Squares (LTS).

La regresión LTS minimiza la suma de cuadrados de los \(q\) residuos más pequeños \(\sum_{i=1}^q \hat{\epsilon}_{(i)}^2\), donde \(q\) es un número menor que \(n\) (tamaño muestral) y \((i)\) indica que los valores están ordenados. Este método es resistente ya que tiene un punto de colapso (breakdwon point) alto y puede tolerar un gran número de puntos atípicos (outliers) en función del valor de \(q\).

Para calcular este método de regresión en R se utiliza la función ltsreg() del paquete MASS.

ltsmod <- ltsreg(Hwt ~ Bwt, data=cats, nsamp="exact")
coef(ltsmod)

## (Intercept)         Bwt 
##  -0.4683219   3.8125000

Por defecto, el algoritmo del método no es determinista de forma que en los casos con pocos datos y pocas variables podemos forzarlo para que haga una búsqueda exhaustiva con la opción nsamp="exact".

Los intervalos de confianza para los coeficientes de esta regresión se pueden obtener con el método bootstrap explicado antes.

bcoef <- matrix(0, nrow=B, ncol=2)
for(i in 1:B){
  ynew <- predict(ltsmod) + residuals(ltsmod)[sample(n, rep=TRUE)]
  fit2 <- ltsreg(ynew ~ Bwt, data=cats, nsamp="best")
  bcoef[i, ] <- coef(fit2)
}  
colnames(bcoef) <- names(coef(ltsmod))
t(apply(bcoef,2,function(x) quantile(x, c(0.025,0.975))))

##                  2.5%    97.5%
## (Intercept) -3.964902 3.761229
## Bwt          2.527509 5.272349

Un gráfico para la pendiente:

plot(density(bcoef[,2]), main="Coeficiente de Bwt")
abline(v=coef(ltsmod)[2], lwd=2, col="red")
abline(v=quantile(bcoef[,2], 0.025), lty=2, lwd=2, col="blue")
abline(v=quantile(bcoef[,2], 0.975), lty=2, lwd=2, col="blue")

Observemos que esta distribución es mucho más puntiaguda que una distribución normal.

En particular, queda claro que este coeficiente no es cero ya que el intervalo de confianza no contiene ese valor.