Tutorial de bootstrap en R

¿Qué es el bootstrap?

El bootstrap es un método de inferencia sobre una población a partir de datos muestrales. Bradley Efron lo presentó por primera vez en este artículo en 1979. El bootstrap se basa en remuestrear con reemplazo a partir de los datos de la muestra. Esta técnica permite estimar el error estándar de cualquier estadístico y obtener su intervalo de confianza (IC). Es especialmente útil cuando el IC no tiene una forma cerrada o esta es muy compleja.

Supongamos que tenemos una muestra de n elementos: X = {x₁, x₂, …, x_n} y nos interesa el IC de algún estadístico T = t(X). El esquema del bootstrap es directo. Repetimos R veces lo siguiente: en la repetición i-ésima, muestreamos con reemplazo n elementos de la muestra disponible (algunos se elegirán más de una vez). Llamamos a esta nueva muestra la muestra bootstrap i-ésima, Xi, y calculamos el estadístico deseado Ti = t(Xi).

Como resultado, obtenemos R valores de nuestro estadístico: T₁, T₂, …, T_R. Los llamamos realizaciones bootstrap de T o la distribución bootstrap de T. A partir de ella, podemos calcular un IC para T. Hay varias formas de hacerlo; usar percentiles suele ser la más sencilla.

Bootstrap en acción

Usemos (otra vez) el conocido conjunto de datos iris. Echa un vistazo a las primeras filas:

 head(iris)\n\n##   Sepal.Length Sepal.Width Petal.Length Petal.Width Species\n## 1          5.1         3.5          1.4         0.2  setosa\n## 2          4.9         3.0          1.4         0.2  setosa\n## 3          4.7         3.2          1.3         0.2  setosa\n## 4          4.6         3.1          1.5         0.2  setosa\n## 5          5.0         3.6          1.4         0.2  setosa\n## 6          5.4         3.9          1.7         0.4  setosa

Imagina que queremos obtener IC para la mediana de Sepal.Length, la mediana de Sepal.Width y el coeficiente de correlación por rangos de Spearman entre ambas. Usaremos el paquete boot de R y una función llamada… boot. Para aprovechar su potencial debemos crear una función que calcule nuestro(s) estadístico(s) a partir de los datos remuestreados. Debe tener al menos dos argumentos: un dataset y un vector con los indices de los elementos del conjunto de datos que se eligieron para crear una muestra bootstrap.

Si queremos calcular IC para más de un estadístico a la vez, nuestra función debe devolverlos en un único vector.

En nuestro ejemplo podría verse así:

 library(boot)\n\nfoo <- function(data, indices){\n  dt<-data[indices,]\n  c(\n    cor(dt[,1], dt[,2], method='s'),\n    median(dt[,1]),\n    median(dt[,2])\n  )\n}

foo elige los elementos deseados (cuyos números están en indices) de data y calcula el coeficiente de correlación de las dos primeras columnas (method='s' selecciona el coeficiente de Spearman; con method='p' obtendrías el de Pearson) y sus medianas.

También podemos añadir argumentos extra; por ejemplo, dejar que el usuario elija el tipo de coeficiente de correlación:

 foo <- function(data, indices, cor.type){\n  dt<-data[indices,]\n  c(\n    cor(dt[,1], dt[,2], method=cor.type),\n    median(dt[,1]),\n    median(dt[,2])\n  )\n}

Para hacer este tutorial más general, usaré la última versión de foo.

Ahora podemos usar la función boot. Tenemos que pasarle el nombre del conjunto de datos, la función que acabamos de crear, el número de repeticiones (R) y cualquier argumento adicional de nuestra función (como cor.type). Abajo uso set.seed para que el ejemplo sea reproducible.

 set.seed(12345)\nmyBootstrap <- boot(iris, foo, R=1000, cor.type='s')

La función boot devuelve un objeto de clase… (sí, ¡lo has adivinado!) boot. Tiene dos elementos interesantes. $t contiene los R valores de nuestro(s) estadístico(s) generados por el procedimiento bootstrap (las realizaciones bootstrap de T):

 head(myBootstrap$t)\n\n##             [,1] [,2] [,3]\n## [1,] -0.26405188 5.70    3\n## [2,] -0.12973299 5.80    3\n## [3,] -0.07972066 5.75    3\n## [4,] -0.16122705 6.00    3\n## [5,] -0.20664808 6.00    3\n## [6,] -0.12221170 5.80    3

$t0 contiene los valores de nuestro(s) estadístico(s) en el conjunto de datos original y completo:

 myBootstrap$t0\n\n## [1] -0.1667777  5.8000000  3.0000000

Imprimir el objeto boot en la consola proporciona más información:

 myBootstrap\n\n##\n## ORDINARY NONPARAMETRIC BOOTSTRAP\n##\n##\n## Call:\n## boot(data = iris, statistic = foo, R = 1000, cor.type = \"s\")\n##\n##\n## Bootstrap Statistics :\n##       original       bias    std. error\n## t1* -0.1667777  0.002546391  0.07573983\n## t2*  5.8000000 -0.013350000  0.10295571\n## t3*  3.0000000  0.007900000  0.02726414

original es lo mismo que $t0. bias es la diferencia entre la media de las realizaciones bootstrap (las de $t), llamada estimación bootstrap de T, y el valor en el conjunto de datos original (el de $t0).

 colMeans(myBootstrap$t)-myBootstrap$t0\n\n## [1]  0.002546391 -0.013350000  0.007900000

std. error es el error estándar de la estimación bootstrap, que equivale a la desviación estándar de las realizaciones bootstrap.

 apply(myBootstrap$t,2,sd)\n\n## [1] 0.07573983 0.10295571 0.02726414

Tipos de IC bootstrap

Antes de ponernos con los IC, siempre conviene mirar la distribución de las realizaciones bootstrap. Podemos usar la función plot, con index indicando a cuál de los estadísticos calculados en foo queremos mirar. Aquí index=1 es el coeficiente de correlación de Spearman entre la longitud y la anchura del sépalo, index=2 es la mediana de la longitud del sépalo e index=3 es la mediana de la anchura del sépalo.

plot(myBootstrap, index=1)\n

La distribución de los coeficientes de correlación bootstrap parece bastante normal. Busquemos su IC. Podemos usar boot.ci. Por defecto calcula IC al 95%, aunque se puede cambiar con el parámetro conf.

 boot.ci(myBootstrap, index=1)\n\n## Warning in boot.ci(myBootstrap, index = 1): bootstrap variances needed for\n## studentized intervals\n\n## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS\n## Based on 1000 bootstrap replicates\n##\n## CALL :\n## boot.ci(boot.out = myBootstrap, index = 1)\n##\n## Intervals :\n## Level      Normal              Basic         \n## 95%   (-0.3178, -0.0209 )   (-0.3212, -0.0329 )  \n##\n## Level     Percentile            BCa          \n## 95%   (-0.3007, -0.0124 )   (-0.3005, -0.0075 )  \n## Calculations and Intervals on Original Scale

boot.ci proporciona 5 tipos de IC bootstrap. Uno de ellos, el intervalo studentized, es especial: necesita una estimación de la varianza bootstrap. No la proporcionamos, así que R muestra el aviso: bootstrap variances needed for studentized intervals. Las estimaciones de varianza pueden obtenerse con un bootstrap de segundo nivel o (más fácil) con la técnica jackknife. Esto queda fuera del alcance de este tutorial, así que centrémonos en los otros cuatro tipos de IC bootstrap.

Si no queremos verlos todos, podemos elegir los relevantes en el argumento type. Los valores posibles son norm, basic, stud, perc, bca o un vector con varios de ellos.

 boot.ci(myBootstrap, index=1, type=c('basic','perc'))\n\n## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS\n## Based on 1000 bootstrap replicates\n##\n## CALL :\n## boot.ci(boot.out = myBootstrap, type = c(\"basic\", \"perc\"), index = 1)\n##\n## Intervals :\n## Level      Basic              Percentile     \n## 95%   (-0.3212, -0.0329 )   (-0.3007, -0.0124 )  \n## Calculations and Intervals on Original Scale

La función boot.ci crea un objeto de clase… (sí, ¡lo has adivinado!) bootci. Sus elementos se llaman igual que los tipos de IC usados en el argumento type. $norm es un vector de 3 elementos que contiene el nivel de confianza y los límites del IC.

 boot.ci(myBootstrap, index=1, type='norm')$norm\n\n##      conf                       \n## [1,] 0.95 -0.3177714 -0.02087672

$basic, $stud, $perc y $bca son vectores de 5 elementos que también incluyen los percentiles usados para calcular el IC (volveremos a esto más adelante):

 boot.ci(myBootstrap, index=1, type='basic')$basic\n\n##      conf                                    \n## [1,] 0.95 975.98 25.03 -0.3211981 -0.03285178

Un poco de notación (¡perdón!)

Para entender qué son los distintos tipos de IC, necesitamos introducir algo de notación. Sea:

\n
• t^⋆ la estimación bootstrap (media de las realizaciones bootstrap),
\n
• t₀ el valor de nuestro estadístico en el conjunto de datos original,
\n
• se^⋆ el error estándar de la estimación bootstrap,
\n
• b el sesgo de la estimación bootstrap, b = t^⋆ − t₀
\n
• α el nivel de confianza, típicamente α = 0.95,
\n
• z_α el cuantil $1-\\frac \\alpha 2$ de la distribución normal estándar,
\n
• θ_α el α-percentil de la distribución de las realizaciones bootstrap.
\n

IC por percentiles

Con la notación anterior, el IC por percentiles es:

es decir, tomamos los percentiles correspondientes. Nada más.

IC normal

Un IC de Wald típico sería:

pero en el caso bootstrap conviene corregirlo por sesgo. Así queda:

$$ t_0 - b \\pm z_\\alpha \\cdot se^\\star \\\\ 2t_0 - t^\\star \\pm z_\\alpha \\cdot se^\\star$$

IC básico

El IC por percentiles en general no se recomienda porque se comporta mal con distribuciones de colas raras. El IC básico (también llamado pivotal o empirical) es mucho más robusto. La idea es calcular las diferencias entre cada réplica bootstrap y t₀ y usar los percentiles de su distribución. Los detalles completos pueden verse, por ejemplo, en All of Statistics de L. Wasserman.

La fórmula final del IC básico es:

BC_α (bias-corrected and accelerated)

BC_α viene de bias-corrected, accelerated. Su fórmula no es muy complicada pero sí poco intuitiva, así que la omitiremos. Consulta el artículo de Thomas J. DiCiccio y Bradley Efron si te interesan los detalles.

La aceleración que menciona el método requiere usar percentiles específicos de las realizaciones bootstrap. A veces puede ocurrir que estos sean percentiles extremos, quizá valores atípicos. En esos casos, BC_α puede resultar inestable.

Veamos el IC BC_α para la mediana de la anchura del pétalo. En el conjunto de datos original, esta mediana es exactamente 3.

 boot.ci(myBootstrap, index=3)\n\n## Warning in boot.ci(myBootstrap, index = 3): bootstrap variances needed for\n## studentized intervals\n\n## Warning in norm.inter(t, adj.alpha): extreme order statistics used as\n## endpoints\n\n## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS\n## Based on 1000 bootstrap replicates\n##\n## CALL :\n## boot.ci(boot.out = myBootstrap, index = 3)\n##\n## Intervals :\n## Level      Normal              Basic         \n## 95%   ( 2.939,  3.046 )   ( 2.900,  3.000 )  \n##\n## Level     Percentile            BCa          \n## 95%   ( 3.0,  3.1 )   ( 2.9,  2.9 )  \n## Calculations and Intervals on Original Scale\n## Warning : BCa Intervals used Extreme Quantiles\n## Some BCa intervals may be unstable

Obtenemos un IC BC_α (2.9, 2.9). Es raro, pero por suerte R nos avisó con extreme order statistics used as endpoints. Veamos qué ha pasado exactamente:

 plot(myBootstrap, index=3)

La distribución de las realizaciones bootstrap es inusual. Una gran mayoría (más del 90%) son 3.

 table(myBootstrap$t[,3])\n\n##\n##  2.9 2.95    3 3.05  3.1 3.15  3.2\n##    1    1  908   24   63    1    2

En una situación así, la \"aceleración\" del método BC_α salta fácilmente a valores extremos. Aquí, el IC por percentiles, que sirve de base para el IC BC_α, es (3, 3.1). Al expandirlo \"hacia la izquierda\" tenemos que usar el percentil 0.002. La expansión \"hacia la derecha\" también se va a los extremos: percentil 0.997.

Reproducir los resultados

A veces necesitamos recrear las réplicas bootstrap. Si podemos usar R, no hay problema: la función set.seed lo soluciona. Recrearlas en otro software sería mucho más complicado. Además, para R grandes, recalcular en R puede no ser una opción (por falta de tiempo, por ejemplo).

Podemos afrontarlo guardando los índices de los elementos del conjunto original que formaron cada muestra bootstrap. Esto es lo que hace la función boot.array (con el argumento indices=T).

 tableOfIndices<-boot.array(myBootstrap, indices=T)

Cada fila es una muestra bootstrap. Por ejemplo, nuestra primera muestra contiene los siguientes elementos:

 tableOfIndices[1,]\n\n##   [1] 109  12 143  65  41  28 105  62  23 102  37  34  55  53  38   3  11\n##  [18] 146  31 122  85 141 118 116 117  78 100  13  98  49  59  88  24  56\n##  [35] 136   4  59  67 126 118 104 101  70  13  12  19  21 149  73 133  67\n##  [52]  86   6  88 131 105 121 145 121  83  70  68  71 111  76  73 122 116\n##  [69]  85 144 114 139  89 124  64  27  81  78  58  39  17  50  37 117  76\n##  [86]  97 114  64  17  93 141  65 124  80 137  54  37  57  70 128  10  55\n## [103]  13  53  59  45 116  14  77 118 108 138  50  78  49 104  49   1  85\n## [120]  28  43  82  74  64  33  55  32  59  62 112   8  11  96  30  14  30\n## [137]  38  85  66  85  97 107   4  18  76  35  31 133  27  69

Si establecemos indices=F (valor por defecto), obtendremos la respuesta a la pregunta: \"¿Cuántas veces apareció cada elemento del conjunto original en cada muestra bootstrap?\". Por ejemplo, en la primera muestra: el primer elemento apareció una vez, el segundo no apareció, el tercero apareció una vez, el cuarto dos veces, y así sucesivamente.

 tableOfAppearances<-boot.array(myBootstrap)\ntableOfAppearances[1,]\n\n##   [1] 1 0 1 2 0 1 0 1 0 1 2 2 3 2 0 0 2 1 1 0 1 0 1 1 0 0 2 2 0 2 2 1 1 1 1\n##  [36] 0 3 2 1 0 1 0 1 0 1 0 0 0 3 2 0 0 2 1 3 1 1 1 4 0 0 2 0 3 2 1 2 1 1 3\n##  [71] 1 0 2 1 0 3 1 3 0 1 1 1 1 0 5 1 0 2 1 0 0 0 1 0 0 1 2 1 0 1 1 1 0 2 2\n## [106] 0 1 1 1 0 1 1 0 2 0 3 2 3 0 0 2 2 0 2 0 1 0 1 0 0 1 0 2 0 0 1 1 1 1 0\n## [141] 2 0 1 1 1 1 0 0 1 0

Una tabla así nos permite recrear las realizaciones bootstrap fuera de R. O en el propio R cuando no queremos usar set.seed y repetir todos los cálculos.

 onceAgain<-apply(tableOfIndices, 1, foo, data=iris, cor.type='s')

Comprobemos si los resultados coinciden:

 head(t(onceAgain))\n\n##             [,1] [,2] [,3]\n## [1,] -0.26405188 5.70    3\n## [2,] -0.12973299 5.80    3\n## [3,] -0.07972066 5.75    3\n## [4,] -0.16122705 6.00    3\n## [5,] -0.20664808 6.00    3\n## [6,] -0.12221170 5.80    3\n\nhead(myBootstrap$t)\n\n##             [,1] [,2] [,3]\n## [1,] -0.26405188 5.70    3\n## [2,] -0.12973299 5.80    3\n## [3,] -0.07972066 5.75    3\n## [4,] -0.16122705 6.00    3\n## [5,] -0.20664808 6.00    3\n## [6,] -0.12221170 5.80    3\n\nall(t(onceAgain)==myBootstrap$t)\n\n## [1] TRUE

¡Sí, coinciden!

Si quieres aprender más sobre Machine Learning en R, haz el curso Machine Learning Toolbox de DataCamp y echa un vistazo al tutorial Machine Learning en R para principiantes.