Capitolo 23 Introduzione alle Tecniche Bootstrap

23.1 Introduzione

Il bootstrap è uno strumento estremamente potente e applicabile in molte aree, che può essere usato per quantificare l’incertezza associata con uno stimatore dato. La potenza del bootstrap si fonda sul fatto che può essere applicato facilmente ad un’ampia classe di metodi statistici, non solo quelli per le previsioni numeriche, e includere alcuni di quelli per cui è difficile ottenere una misura di variabilità e non è fornita automaticamente dai software statistici. Supponiamo per esempio di essere interessati a fare inferenza sul rapporto tra due misure di medie della popolazione. Non ci sono semplici metodi tradizionali per fare inferenza in questi casi, e per valutarne l’incertezza. Il bootstrap può essere usato efficacemente in tali situazioni non standard.

Il bootstrap (non parametrico) funziona indipendentemente dalla distribuzione dei dati nella popolazione. Rispetto a questo aspetto, il bootstrap fornisce un approccio non parametrico generale per fare inferenza sulle caratteristiche di una popolazione di interesse.

Si suppnga di avere un set di \(n\) osservazioni da cui voler trarre stime per un parametro (o un insieme di parametri) \(T\), la cui stima la indichiamo con \(t\). I passaggi per applicare il bootstrap non parametrico sono:

Genera \(B\) nuovi campioni (usualmente alcune centinaia sono sufficienti), chiamati campioni bootstrap, o ricampionamenti, campionando con reinserimento dal campione casuale originario. Ogni ricampionamento avrà la stessa dimensione dell’originale campione casuale, \(n\). Campionare con reinserimento significa che, dopo aver estratto a caso una osservazione del campione originario, la dobbiamo reimmettere nel campione originario stesso prima di estrarre l’osservazione successiva. Di conseguenza, avremo che ogni osservazione potrà essere campionata una volta, pià di una volta, o anche mai, all’interno dello stesso ricampionamento.
Per ciascun ricampionamento, calcola la statistica di interesse, indicata come \(t(b)\), con \(b = 1, \cdots ,B;\) la distribuzione di queste statistiche di ricampionamento è chiamata “distribuzione bootstrap”.
La distribuzione bootstrap fornisce informazioni sulle caratteristiche (forma, centralità, dispersione, ecc.) della distribuzione campionaria della statistica d’interesse. In altre parole, la distrbuzione bootstrap fornisce una approssimazione della distribuzione campionaria della statistica.

Di conseguenza, l’idea alla base della tecnica è che il campione originario rappresenta la popolazione da cui campionare. Quindi, il ricampionamento da questo campione rappresenta una “stima” di ciò che uno otterrebbe se si estraessero diversi campioni dalla popolazione. La distribuzione bootstrap di una statistica, basata su molti ricampionamenti, rappresenta una “stima” della distribuzione campionaia della statistica \(t\), basata su molti campioni.

23.2 Esempio: semplice applicazione del bootstrap

L’esecuzione di una analisi di bootstrap in R richiede solo due passaggi: 1. Per prima cosa dobbiamo creare una funzione che calcola la statistica d’interesse. 2. Poi possiamo usare la funzione boot() del package boot per eseguire il bootstrap campionando ripetutamente con reinserimento osservazioni dal dataset originario.

A livello illustrativo, usiamo il dataset Auto del package ISLR per determinare la variabilità dei termini di intercetta e coefficiente angolare per il modello di regressione lineare che usa horsepower per “prevedere” mpg. Confeonteremo quindi le stime ottenute usando il bootstrap con quelle ottenute con il metodo OLS classico.

summary(Auto)

##       mpg          cylinders      displacement     horsepower        weight    
##  Min.   : 9.00   Min.   :3.000   Min.   : 68.0   Min.   : 46.0   Min.   :1613  
##  1st Qu.:17.00   1st Qu.:4.000   1st Qu.:105.0   1st Qu.: 75.0   1st Qu.:2225  
##  Median :22.75   Median :4.000   Median :151.0   Median : 93.5   Median :2804  
##  Mean   :23.45   Mean   :5.472   Mean   :194.4   Mean   :104.5   Mean   :2978  
##  3rd Qu.:29.00   3rd Qu.:8.000   3rd Qu.:275.8   3rd Qu.:126.0   3rd Qu.:3615  
##  Max.   :46.60   Max.   :8.000   Max.   :455.0   Max.   :230.0   Max.   :5140  
##                                                                                
##   acceleration        year           origin                      name    
##  Min.   : 8.00   Min.   :70.00   Min.   :1.000   amc matador       :  5  
##  1st Qu.:13.78   1st Qu.:73.00   1st Qu.:1.000   ford pinto        :  5  
##  Median :15.50   Median :76.00   Median :1.000   toyota corolla    :  5  
##  Mean   :15.54   Mean   :75.98   Mean   :1.577   amc gremlin       :  4  
##  3rd Qu.:17.02   3rd Qu.:79.00   3rd Qu.:2.000   amc hornet        :  4  
##  Max.   :24.80   Max.   :82.00   Max.   :3.000   chevrolet chevette:  4  
##                                                  (Other)           :365

Creiamo quindi una semplice funzione, boot.fn(), che prende in input il dataset Auto ed un insieme di indici per le osservazioni, e ritorna le stime del coefficiente angolare e dell’intercetta. Possiamo quindi applicare questa funzione per calcolare le stime dell’intero dataset usando l’usuale formula per la stime dei coefficienti della regressione lineare:

boot.fn <- function(data, index) {
    return(coef(lm(mpg ~ horsepower, data = data, subset = index)))
}
boot.fn(Auto, 1:392)

## (Intercept)  horsepower 
##  39.9358610  -0.1578447

La funzione boot.fn() può essere usata per creare stime bootstrap per l’intercetta e il coefficiente angolare campionando casualmente con reinserimento le osservazioni:

set.seed(10)
boot.fn(Auto, sample(1:392, 392, replace = TRUE))

## (Intercept)  horsepower 
##  39.5905422  -0.1584437

Le stime risultanti sono una stima campionaria bootstrap dei parametri intercetta e coefficiente angolare.

Ora possiamo usare la funzione boot() del package boot per calcolare gli errori standard su 1000 stime bootstrap per i termini di intercetta e coefficiente angolare:

library(boot)

set.seed(20)
res <- boot(data = Auto, statistic = boot.fn, R = 1000)
print(res)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Auto, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##       original       bias    std. error
## t1* 39.9358610  0.022492354 0.850524153
## t2* -0.1578447 -0.000351767 0.007262426

Notte l’uso del parametro R = 1000, che indica quanti ricampionamenti devono essere fatti (a tutti gli effetti corrisponde al valore del parametro \(B\) indicato nella parte di introduzione all’algoritmo).

Questo indica che la stima bootstrap dell’errore standard dell’intercetta è 0.8505242, e che la stima bootstrap per l’errore standard dell’intercetta è 0.0072624. I corrispondenti errori standard OLS sono dati da

summary(lm(mpg ~ horsepower, data = Auto))$coef

##               Estimate  Std. Error   t value      Pr(>|t|)
## (Intercept) 39.9358610 0.717498656  55.65984 1.220362e-187
## horsepower  -0.1578447 0.006445501 -24.48914  7.031989e-81

La cosa interessante è che queste stime sono in qualche modo differenti da quelle ottenute usando il bootstrap. Questo dovrebbe indicare qualche problema con il bootstrap? Nei fatti, questo dovrebbe suggerire il contrario, poiché l’usuale inferenza in un modello di regressione lineare richiede determinati assunti (linearità, omoscedasticità, normalità, ecc.) che raramente sono soddisfatti nella pratica.

Come ulteriore esempio, calcoliamo qui di seguito la stime dell’errore standard bootstrap e della regressione lineare standard che risultano dall’adattamento del modello quadratico sugli stessi dati:

boot.fn <- function(data, index) {
    return(coefficients(lm(mpg ~ horsepower + I(horsepower^2), data = data,
    subset = index)))
}
set.seed(10)
res <- boot(data = Auto, statistic = boot.fn, R = 1000)
print(res)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Auto, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##         original        bias     std. error
## t1* 56.900099702  1.364490e-01 2.0493438886
## t2* -0.466189630 -2.244414e-03 0.0325800279
## t3*  0.001230536  7.912760e-06 0.0001179559

summary(lm(mpg ~ horsepower + I(horsepower^2), data = Auto))$coef

##                     Estimate   Std. Error   t value      Pr(>|t|)
## (Intercept)     56.900099702 1.8004268063  31.60367 1.740911e-109
## horsepower      -0.466189630 0.0311246171 -14.97816  2.289429e-40
## I(horsepower^2)  0.001230536 0.0001220759  10.08009  2.196340e-21

Infine, per testare il bootstrap su dati correttamente specificati “per costruzione”, potremo provare col seguente:

set.seed(10)

x <- runif(n = 500, max = 200)
y <- 2 + 3*x +rnorm(500)

dsim <- data.frame(x=x,y=y)
boot.fn <- function(data, index) {
    return(coefficients(lm(y ~ x, data = data,
    subset = index)))
}

res <- boot(data = dsim, statistic = boot.fn, R = 1000)
print(res)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = dsim, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original        bias     std. error
## t1* 2.055483 -2.687525e-03 0.0946359049
## t2* 2.999537  1.937389e-05 0.0008330455

summary(lm(y ~ x, data = dsim))$coef

##             Estimate  Std. Error    t value     Pr(>|t|)
## (Intercept) 2.055483 0.094825118   21.67657 7.085419e-74
## x           2.999537 0.000827011 3626.96147 0.000000e+00

In questo caso, dove il modello è correttamente specificato “per costruzione”, si può vedere come i risultati bootstrap e “standard” siano piuttosto simili.