R Bootstrap-statistiek & betrouwbaarheidsintervallen (CI) Tutorial

Wat is bootstrap?

Bootstrap is een methode om uitspraken te doen over een populatie op basis van steekproefdata. Bradley Efron introduceerde het voor het eerst in dit artikel in 1979. Bootstrap berust op steekproeven met teruglegging uit de steekproefdata. Met deze techniek kun je de standaardfout van om het even welke statistiek schatten en er een betrouwbaarheidsinterval (CI) voor verkrijgen. Bootstrap is vooral nuttig wanneer een CI geen gesloten vorm heeft, of wanneer die erg ingewikkeld is.

Stel, we hebben een steekproef van n elementen: X = {x₁, x₂, …, x_n} en we zijn geïnteresseerd in een CI voor een of andere statistiek T = t(X). Het bootstrap-kader is rechttoe rechtaan. We herhalen R keer het volgende schema: Neem bij de i-de herhaling een steekproef met teruglegging van n elementen uit de beschikbare steekproef (sommige worden meer dan eens gekozen). Noem deze nieuwe steekproef de i-de bootstrapsteekproef, Xi, en bereken de gewenste statistiek Ti = t(Xi).

Als resultaat krijgen we R waarden van onze statistiek: T₁, T₂, …, T_R. We noemen dit bootstraprealisaties van T of een bootstrapverdeling van T. Op basis daarvan kunnen we een CI voor T berekenen. Er zijn verschillende manieren om dit te doen. Percentielen nemen lijkt de eenvoudigste.

Bootstrap in de praktijk

Laten we (weer eens) de bekende iris-dataset gebruiken. Kijk even naar de eerste paar rijen:

 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  setosaMogelijk gemaakt door

Stel dat we CIs willen vinden voor de mediaan van Sepal.Length, de mediaan van Sepal.Width en Spearmans rangcorrelatiecoëfficiënt tussen deze twee. We gebruiken het boot-pakket van R en een functie genaamd... boot. Om die kracht te benutten, moeten we een functie maken die onze statistiek(en) uit de geresamplede data berekent. Die moet minstens twee argumenten hebben: een dataset en een vector met indices van de elementen uit de dataset die gekozen zijn om een bootstrapsteekproef te vormen.

Als we CIs voor meer dan één statistiek tegelijk willen berekenen, moet onze functie ze als één enkele vector teruggeven.

Voor ons voorbeeld kan dat er zo uitzien:

 library(boot)

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

foo kiest de gewenste elementen (welke nummers in indices staan) uit data en berekent de correlatiecoëfficiënt van de eerste twee kolommen (method='s' kiest Spearmans rangcoëfficiënt, method='p' geeft Pearsons coëfficiënt) en hun medianen.

We kunnen ook extra argumenten toevoegen, bv. de gebruiker het type correlatiecoëfficiënt laten kiezen:

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

Om deze tutorial algemener te maken, gebruik ik de laatste versie van foo.

Nu kunnen we de functie boot gebruiken. We moeten de naam van de dataset meegeven, de zojuist gemaakte functie, het aantal herhalingen (R) en extra argumenten van onze functie (zoals cor.type). Hieronder gebruik ik set.seed om dit voorbeeld reproduceerbaar te maken.

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

De functie boot geeft een object terug van de klasse... (ja, je raadt het al!) boot. Het heeft twee interessante elementen. $t bevat R waarden van onze statistiek(en) die door de bootstrapprocedure zijn gegenereerd (bootstraprealisaties van 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    3Mogelijk gemaakt door

$t0 bevat de waarden van onze statistiek(en) in de oorspronkelijke, volledige dataset:

 myBootstrap$t0

## [1] -0.1667777  5.8000000  3.0000000Mogelijk gemaakt door

Het printen van het boot-object naar de console geeft nog wat extra informatie:

 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.02726414Mogelijk gemaakt door

original is hetzelfde als $t0. bias is het verschil tussen het gemiddelde van de bootstraprealisaties (die uit $t), de zogeheten bootstrapschatting van T, en de waarde in de oorspronkelijke dataset (die uit $t0).

 colMeans(myBootstrap$t)-myBootstrap$t0

## [1]  0.002546391 -0.013350000  0.007900000Mogelijk gemaakt door

std. error is de standaardfout van de bootstrapschatting, gelijk aan de standaarddeviatie van de bootstraprealisaties.

 apply(myBootstrap$t,2,sd)

## [1] 0.07573983 0.10295571 0.02726414Mogelijk gemaakt door

Verschillende types bootstrap-CI's

Voordat we met CI aan de slag gaan, is het altijd de moeite waard om naar de verdeling van de bootstraprealisaties te kijken. We kunnen de functie plot gebruiken, met index die aangeeft naar welke van de in foo berekende statistieken we willen kijken. Hier is index=1 Spearmans correlatiecoëfficiënt tussen kelklengte en -breedte, index=2 is de mediaan van kelklengte en index=3 is de mediaan van kelkbreedte.

plot(myBootstrap, index=1)
Mogelijk gemaakt door

De verdeling van de bootstrapcorrelatiecoëfficiënten lijkt vrij normaal. Laten we er een CI voor vinden. We kunnen boot.ci gebruiken. Standaard geeft het 95%-CI's, maar dat kun je aanpassen met de parameter 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 ScaleMogelijk gemaakt door

boot.ci levert 5 types bootstrap-CI's. Eén daarvan, het studentized interval, is uniek. Het heeft een schatting van de bootstrapvariantie nodig. Die hebben we niet aangeleverd, dus R geeft een waarschuwing: bootstrap variances needed for studentized intervals. Variantieschattingen kunnen verkregen worden met een bootstrap op tweede niveau of (makkelijker) met de jackknife-techniek. Dat valt enigszins buiten de scope van deze tutorial, dus we focussen op de overige vier types bootstrap-CI's.

Als we ze niet allemaal willen zien, kunnen we relevante kiezen in het argument type. Mogelijke waarden zijn norm, basic, stud, perc, bca of een vector hiervan.

 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 ScaleMogelijk gemaakt door

De functie boot.ci maakt een object van klasse... (ja, je hebt het geraden!) bootci. De elementen heten precies zoals de types CI die in het argument type zijn gebruikt. $norm is een vector met 3 elementen die het betrouwbaarheidsniveau en de CI-grenzen bevat.

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

##      conf                       
## [1,] 0.95 -0.3177714 -0.02087672Mogelijk gemaakt door

$basic, $stud, $perc en $bca zijn vectoren met 5 elementen die ook de percentielen bevatten die zijn gebruikt om het CI te berekenen (daar komen we zo op terug):

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

##      conf                                    
## [1,] 0.95 975.98 25.03 -0.3211981 -0.03285178Mogelijk gemaakt door

Een beetje notatie (sorry!)

Om gevoel te krijgen bij wat de verschillende types CI zijn, moeten we wat notatie introduceren. Laat dus:

t^⋆ de bootstrapschatting zijn (het gemiddelde van de bootstraprealisaties),
t₀ de waarde van onze statistiek in de oorspronkelijke dataset zijn,
se^⋆ de standaardfout van de bootstrapschatting zijn,
b de bias van de bootstrapschatting zijn b = t^⋆ − t₀
α het betrouwbaarheidsniveau zijn, typisch α = 0.95,
z_α de $1-\frac \alpha 2$-kwantiel van de standaardnormale verdeling zijn,
θ_α het α-percentiel van de verdeling van bootstraprealisaties zijn.

Percentiel-CI

Met bovenstaande notatie is het percentiel-CI:

(θ_{(1 − α)/2}, θ_{1 − (1 − α)/2})

dus dit neemt simpelweg de relevante percentielen. Niet meer dan dat.

Normaal CI

Een typisch Wald-CI zou zijn:

t₀ ± z_α ⋅ se^⋆

maar in het bootstrapgeval moeten we het corrigeren voor bias. Het wordt dan:

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

Basic CI

Het percentiel-CI wordt over het algemeen niet aanbevolen omdat het slecht presteert bij vreemd-gedaante staartverdelingen. Het basic CI (ook wel pivotal of empirical CI genoemd) is hier veel robuuster tegen. Het idee erachter is om de verschillen te berekenen tussen elke bootstraprealisatie en t₀ en de percentielen van hun verdeling te gebruiken. Volledige details vind je bijvoorbeeld in L. Wassermans All of Statistics

De uiteindelijke formule voor het basic CI is:

(2t₀ − θ_{1 − (1 − α)/2}, 2t₀ − θ_{(1 − α)/2})

BC_α-CI

BC_α staat voor bias-corrected, accelerated. De formule is niet heel ingewikkeld maar wel enigszins onintuïtief, dus ik sla die over. Zie het artikel van Thomas J. DiCiccio en Bradley Efron als je de details wilt.

De acceleratie in de naam van de methode vereist het gebruik van specifieke percentielen van de bootstraprealisaties. Soms kan het gebeuren dat dit extreme percentielen zijn, mogelijk uitschieters. BC_α kan in zulke gevallen instabiel zijn.

Laten we naar het BC_α-CI voor de mediaan van kelkbreedte kijken. In de oorspronkelijke dataset is deze mediaan precies 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 unstableMogelijk gemaakt door

We krijgen het BC_α-CI (2,9, 2,9). Dat is vreemd, maar gelukkig waarschuwde R ons dat extreme order statistics used as endpoints. Laten we precies zien wat hier gebeurde:

 plot(myBootstrap, index=3)Mogelijk gemaakt door

De verdeling van de bootstraprealisaties is ongebruikelijk. Een overgrote meerderheid (meer dan 90%) is 3.

 table(myBootstrap$t[,3])

##
##  2.9 2.95    3 3.05  3.1 3.15  3.2
##    1    1  908   24   63    1    2Mogelijk gemaakt door

In zo'n situatie springt de "acceleratie" in de BC_α-methode gemakkelijk naar extreme waarden. Hier is het percentiel-CI, dat de basis vormt voor het BC_α-CI, (3, 3,1). Wanneer we het nu "naar links" uitbreiden, moeten we het 0,002-percentiel gebruiken. De uitbreiding "naar rechts" gaat ook naar extremen: het 0,997-percentiel.

Reproductie van resultaten

Soms moeten we bootstraprealisaties opnieuw creëren. Als we daarvoor R kunnen gebruiken, is dat geen probleem. De functie set.seed lost het op. Het reproduceren van bootstraprealisaties in andere software zou een stuk ingewikkelder zijn. Bovendien kan voor grote R herberekenen in R ook geen optie zijn (bijvoorbeeld wegens tijdgebrek).

We kunnen dit aanpakken door de indices op te slaan van de elementen uit de oorspronkelijke dataset die elke bootstrapsteekproef vormden. Dat is wat de functie boot.array (met argument indices=T) doet.

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

Elke rij is één bootstrapsteekproef. Bijv., onze eerste steekproef bevat de volgende elementen:

 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  69Mogelijk gemaakt door

Door indices=F te zetten (de standaard) krijgen we antwoord op de vraag: "Hoe vaak kwam elk element van de oorspronkelijke dataset voor in elke bootstrapsteekproef?". Bijvoorbeeld, in de eerste steekproef: het 1e element van de dataset kwam één keer voor, het 2e element kwam helemaal niet voor, het 3e element één keer, het 4e twee keer, enzovoort.

 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 0Mogelijk gemaakt door

Een dergelijke tabel stelt ons in staat om een bootstraprealisatie buiten R te reconstrueren. Of in R zelf wanneer we set.seed niet willen gebruiken en alles opnieuw willen uitrekenen.

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

Laten we controleren of de resultaten hetzelfde zijn:

 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] TRUEMogelijk gemaakt door

Ja, dat zijn ze!

Wil je meer leren over Machine Learning in R? Volg dan de DataCamp-cursus Machine Learning Toolbox en bekijk de tutorial Machine Learning in R voor beginners.

Onderwerpen

R Data Science

Onderwerpen

Technologieën

Bootstrap in R-tutorial

Verken met AI

Wat is bootstrap?

Bootstrap in de praktijk

Verschillende types bootstrap-CI's

Een beetje notatie (sorry!)

Percentiel-CI

Normaal CI

Basic CI

BC_α-CI

Reproductie van resultaten

Onderwerpen

Leer meer over R en Machine Learning

Machine Learning met caret in R

Machine Learning for Marketing Analytics in R

Unsupervised learning in R

Gerelateerd

AI vanaf nul leren in 2026: een complete gids van de experts

Onderwerpen

Technologieën

Leer meer over R en Machine Learning

Machine Learning met caret in R

Machine Learning for Marketing Analytics in R

Unsupervised learning in R

Wat is bootstrap?

Bootstrap in de praktijk

Verschillende types bootstrap-CI's

Een beetje notatie (sorry!)

Percentiel-CI

Normaal CI

Basic CI

BCα-CI

Reproductie van resultaten

Onderwerpen

Leer meer over R en Machine Learning

Gerelateerd

Leer meer over R en Machine Learning

BC_α-CI