Tutorial: Bootstrap in R

Che cos'è il Bootstrap?

Il bootstrap è un metodo d'inferenza su una popolazione utilizzando dati campionari. Bradley Efron lo ha introdotto per la prima volta in questo articolo nel 1979. Il bootstrap si basa sul campionamento con reinserimento dai dati campionari. Questa tecnica può essere usata per stimare l'errore standard di qualsiasi statistica e per ottenere un intervallo di confidenza (CI) per essa. Il bootstrap è particolarmente utile quando il CI non ha una forma chiusa, o ne ha una molto complicata.

Supponiamo di avere un campione di n elementi: X = {x₁, x₂, …, x_n} e di essere interessati al CI per qualche statistica T = t(X). Il framework bootstrap è lineare. Ripetiamo semplicemente per R volte il seguente schema: per la i-esima ripetizione, campioniamo con reinserimento n elementi dal campione disponibile (alcuni verranno scelti più di una volta). Chiamiamo questo nuovo campione i-esimo campione bootstrap, Xi, e calcoliamo la statistica desiderata Ti = t(Xi).

Come risultato, otterremo R valori della nostra statistica: T₁, T₂, …, T_R. Li chiamiamo realizzazioni bootstrap di T o una distribuzione bootstrap di T. Sulla base di ciò, possiamo calcolare un CI per T. Ci sono diversi modi per farlo. Prendere i percentili sembra essere il più semplice.

Bootstrap in pratica

Usiamo (di nuovo) il ben noto dataset iris. Dai un'occhiata alle prime righe:

 head(iris)

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

Supponiamo di voler trovare i CI per la mediana di Sepal.Length, la mediana di Sepal.Width e il coefficiente di correlazione di rango di Spearman tra le due. Useremo il pacchetto boot di R e una funzione chiamata... boot. Per sfruttarne la potenza dobbiamo creare una funzione che calcoli la/le nostra/e statistica/he dai dati ricampionati. Dovrebbe avere almeno due argomenti: un dataset e un vettore contenente gli indices degli elementi del dataset che sono stati scelti per creare un campione bootstrap.

Se vogliamo calcolare i CI per più di una statistica alla volta, la nostra funzione deve restituirle come un unico vettore.

Per il nostro esempio, potrebbe essere così:

 library(boot)

foo <- function(data, indices){
  dt<-data[indices,]
  c(
    cor(dt[,1], dt[,2], method='s'),
    median(dt[,1]),
    median(dt[,2])
  )
}Offerto da

foo sceglie gli elementi desiderati (i cui numeri sono memorizzati in indices) da data e calcola il coefficiente di correlazione delle prime due colonne (method='s' serve per scegliere il coefficiente di rango di Spearman, method='p' restituirà il coefficiente di Pearson) e le loro mediane.

Possiamo anche aggiungere argomenti extra, ad es. lasciare all'utente la scelta del tipo di coefficiente di correlazione:

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

Per rendere questo tutorial più generale, userò la seconda versione di foo.

Ora possiamo usare la funzione boot. Dobbiamo fornirle il nome del dataset, la funzione che abbiamo appena creato, il numero di ripetizioni (R) e gli eventuali argomenti aggiuntivi della nostra funzione (come cor.type). Sotto, uso set.seed per rendere riproducibile questo esempio.

 set.seed(12345)
myBootstrap <- boot(iris, foo, R=1000, cor.type='s')Offerto da

La funzione boot restituisce un oggetto di classe chiamata... (sì, esatto!) boot. Ha due elementi interessanti. $t contiene R valori della/e nostra/e statistica/he generati dalla procedura bootstrap (realizzazioni bootstrap di T):

 head(myBootstrap$t)

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3Offerto da

$t0 contiene i valori della/e nostra/e statistica/he nel dataset originale, completo:

 myBootstrap$t0

## [1] -0.1667777  5.8000000  3.0000000Offerto da

Stampare l'oggetto boot in console fornisce ulteriori informazioni:

 myBootstrap

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

original è lo stesso di $t0. bias è la differenza tra la media delle realizzazioni bootstrap (quelle in $t), chiamata stima bootstrap di T, e il valore nel dataset originale (quello in $t0).

 colMeans(myBootstrap$t)-myBootstrap$t0

## [1]  0.002546391 -0.013350000  0.007900000Offerto da

std. error è l'errore standard della stima bootstrap, che è uguale alla deviazione standard delle realizzazioni bootstrap.

 apply(myBootstrap$t,2,sd)

## [1] 0.07573983 0.10295571 0.02726414Offerto da

Tipi diversi di CI Bootstrap

Prima di passare ai CI, vale sempre la pena dare un'occhiata alla distribuzione delle realizzazioni bootstrap. Possiamo usare la funzione plot, con index che indica a quale delle statistiche calcolate in foo vogliamo guardare. Qui index=1 è il coefficiente di correlazione di Spearman tra lunghezza e larghezza del sepalo, index=2 è la mediana della lunghezza del sepalo e index=3 è la mediana della larghezza del sepalo.

plot(myBootstrap, index=1)
Offerto da

La distribuzione dei coefficienti di correlazione bootstrap sembra piuttosto normale. Troviamone il CI. Possiamo usare boot.ci. Il valore predefinito è un CI al 95%, ma può essere cambiato con il parametro conf.

 boot.ci(myBootstrap, index=1)

## Warning in boot.ci(myBootstrap, index = 1): bootstrap variances needed for
## studentized intervals

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = myBootstrap, index = 1)
##
## Intervals :
## Level      Normal              Basic         
## 95%   (-0.3178, -0.0209 )   (-0.3212, -0.0329 )  
##
## Level     Percentile            BCa          
## 95%   (-0.3007, -0.0124 )   (-0.3005, -0.0075 )  
## Calculations and Intervals on Original ScaleOfferto da

boot.ci fornisce 5 tipi di CI bootstrap. Uno di questi, l'intervallo studentizzato, è particolare. Ha bisogno di una stima della varianza bootstrap. Non l'abbiamo fornita, quindi R stampa un avviso: bootstrap variances needed for studentized intervals. Le stime di varianza possono essere ottenute con un bootstrap di secondo livello o (più facilmente) con la tecnica jackknife. Questo esula un po' dallo scopo di questo tutorial, quindi concentriamoci sui restanti quattro tipi di CI bootstrap.

Se non vogliamo vederli tutti, possiamo selezionare quelli rilevanti nell'argomento type. I valori possibili sono norm, basic, stud, perc, bca o un vettore di questi.

 boot.ci(myBootstrap, index=1, type=c('basic','perc'))

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

La funzione boot.ci crea un oggetto di classe... (sì, l'hai indovinato!) bootci. I suoi elementi si chiamano esattamente come i tipi di CI usati nell'argomento type. $norm è un vettore di 3 elementi che contiene il livello di confidenza e gli estremi del CI.

 boot.ci(myBootstrap, index=1, type='norm')$norm

##      conf                       
## [1,] 0.95 -0.3177714 -0.02087672Offerto da

$basic, $stud, $perc e $bca sono vettori di 5 elementi che contengono anche i percentili usati per calcolare il CI (ci torneremo tra poco):

 boot.ci(myBootstrap, index=1, type='basic')$basic

##      conf                                    
## [1,] 0.95 975.98 25.03 -0.3211981 -0.03285178Offerto da

Un po' di notazione (scusa!)

Per capire cosa sono i vari tipi di CI, dobbiamo introdurre un po' di notazione. Quindi, sia:

• t^⋆ la stima bootstrap (media delle realizzazioni bootstrap),
• t₀ il valore della nostra statistica nel dataset originale,
• se^⋆ l'errore standard della stima bootstrap,
• b il bias della stima bootstrap, b = t^⋆ − t₀
• α il livello di confidenza, tipicamente α = 0.95,
• z_α il quantile $1-\frac \alpha 2$ della distribuzione normale standard,
• θ_α un α-percentile della distribuzione delle realizzazioni bootstrap.

CI a percentili

Con la notazione sopra, il CI a percentili è:

quindi prende semplicemente i percentili rilevanti. Niente di più.

CI normale

Un tipico CI di Wald sarebbe:

ma nel caso bootstrap, dovremmo correggerlo per il bias. Quindi diventa:

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

CI di base

Il CI a percentili in generale non è raccomandato perché si comporta male in presenza di distribuzioni con code strane. Il CI di base (chiamato anche pivotal o empirical) è molto più robusto a questo. L'idea è calcolare le differenze tra ciascuna replica bootstrap e t₀ e usare i percentili della loro distribuzione. I dettagli completi si trovano, ad esempio, in All of Statistics di L. Wasserman

La formula finale per il CI di base è:

BC_α CI

BC_α sta per bias-corrected, accelerated (corretto per bias, accelerato). La formula non è molto complicata ma è poco intuitiva, quindi la tralascio. Vedi l'articolo di Thomas J. DiCiccio e Bradley Efron se ti interessano i dettagli.

L'“accelerazione” menzionata nel nome del metodo richiede l'uso di specifici percentili delle realizzazioni bootstrap. A volte può capitare che questi siano percentili estremi, possibilmente outlier. In tali casi, il BC_α può risultare instabile.

Vediamo il CI BC_α per la mediana della larghezza del petalo. Nel dataset originale, questa mediana è esattamente 3.

 boot.ci(myBootstrap, index=3)

## Warning in boot.ci(myBootstrap, index = 3): bootstrap variances needed for
## studentized intervals

## 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 = myBootstrap, index = 3)
##
## Intervals :
## Level      Normal              Basic         
## 95%   ( 2.939,  3.046 )   ( 2.900,  3.000 )  
##
## Level     Percentile            BCa          
## 95%   ( 3.0,  3.1 )   ( 2.9,  2.9 )  
## Calculations and Intervals on Original Scale
## Warning : BCa Intervals used Extreme Quantiles
## Some BCa intervals may be unstableOfferto da

Otteniamo il CI BC_α (2.9, 2.9). È strano, ma per fortuna R ci ha avvisato che extreme order statistics used as endpoints. Vediamo esattamente cosa è successo qui:

 plot(myBootstrap, index=3)Offerto da

La distribuzione delle realizzazioni bootstrap è insolita. Una grande maggioranza di esse (oltre il 90%) vale 3.

 table(myBootstrap$t[,3])

##
##  2.9 2.95    3 3.05  3.1 3.15  3.2
##    1    1  908   24   63    1    2Offerto da

In una situazione del genere, l'“accelerazione” nel metodo BC_α “salta” facilmente a valori estremi. Qui, il CI a percentili, che è alla base del CI BC_α, è (3, 3.1). Ora, quando lo espandiamo “a sinistra” dobbiamo usare il percentile 0.002. L'espansione “a destra” va anch'essa verso gli estremi: percentile 0.997.

Replicazione dei risultati

A volte dobbiamo ricreare le repliche bootstrap. Se possiamo usare R per farlo, non è un problema. La funzione set.seed risolve la questione. Ricreare le repliche bootstrap in altro software sarebbe molto più complicato. Inoltre, per R grande, ricalcolare in R potrebbe non essere un'opzione (per mancanza di tempo, ad esempio).

Possiamo affrontare il problema salvando gli indici degli elementi del dataset originale che hanno formato ciascun campione bootstrap. È ciò che fa la funzione boot.array (con argomento indices=T).

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

Ogni riga è un campione bootstrap. Ad esempio, il nostro primo campione contiene i seguenti elementi:

 tableOfIndices[1,]

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

Impostando indices=F (predefinito), otterremo una risposta alla domanda "Quante volte ciascun elemento del dataset originale è apparso in ciascun campione bootstrap?". Per esempio, nel primo campione: il primo elemento del dataset è apparso una volta, il secondo elemento non è apparso, il terzo elemento è apparso una volta, il quarto due volte e così via.

 tableOfAppearances<-boot.array(myBootstrap)
tableOfAppearances[1,]

##   [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
##  [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
##  [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
## [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
## [141] 2 0 1 1 1 1 0 0 1 0Offerto da

Una tabella come questa ci permette di ricreare le realizzazioni bootstrap al di fuori di R. O in R stesso quando non vogliamo usare set.seed e rifare tutti i calcoli da capo.

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

Verifichiamo se i risultati coincidono:

 head(t(onceAgain))

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3

head(myBootstrap$t)

##             [,1] [,2] [,3]
## [1,] -0.26405188 5.70    3
## [2,] -0.12973299 5.80    3
## [3,] -0.07972066 5.75    3
## [4,] -0.16122705 6.00    3
## [5,] -0.20664808 6.00    3
## [6,] -0.12221170 5.80    3

all(t(onceAgain)==myBootstrap$t)

## [1] TRUEOfferto da

Sì, coincidono!

Se vuoi imparare di più sul Machine Learning in R, segui il corso di DataCamp Machine Learning Toolbox e dai un'occhiata al tutorial Machine Learning in R for beginners.