Bootstrap in R: Tutorial

Was ist Bootstrap?

Bootstrap ist ein Verfahren, um aus Stichprobendaten auf eine Grundgesamtheit zu schließen. Bradley Efron hat es 1979 in diesem Paper vorgestellt. Bootstrap basiert auf Ziehen mit Zurücklegen aus der Stichprobe. Mit dieser Technik lässt sich der Standardfehler beliebiger Statistiken schätzen und ein Konfidenzintervall (CI) dafür bestimmen. Bootstrap ist besonders hilfreich, wenn es für das CI keine geschlossene Form gibt oder sie sehr kompliziert ist.

Angenommen, wir haben eine Stichprobe mit n Elementen: X = {x₁, x₂, …, x_n} und interessieren uns für ein CI einer Statistik T = t(X). Das Bootstrap-Vorgehen ist einfach: Wir wiederholen R-mal Folgendes: Für die i-te Wiederholung ziehen wir mit Zurücklegen n Elemente aus der vorhandenen Stichprobe (einige werden mehrfach gewählt). Diese neue Stichprobe nennen wir i-te Bootstrap-Stichprobe, Xi, und berechnen die gewünschte Statistik Ti = t(Xi).

Am Ende erhalten wir R Werte unserer Statistik: T₁, T₂, …, T_R. Das sind die Bootstrap-Realisationen von T bzw. die Bootstrap-Verteilung von T. Darauf basierend können wir ein CI für T berechnen. Es gibt mehrere Wege, der Perzentilansatz ist der einfachste.

Bootstrap in der Praxis

Wir nutzen erneut den bekannten iris-Datensatz. Schau dir die ersten Zeilen an:

 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

Angenommen, wir möchten CIs für den Median von Sepal.Length und Sepal.Width sowie für Spearmans Rangkorrelationskoeffizienten zwischen beiden bestimmen. Wir verwenden das boot-Paket in R und eine Funktion namens ... boot. Dafür brauchen wir eine Funktion, die aus den resampleten Daten unsere(n) Statistik(en) berechnet. Sie sollte mindestens zwei Argumente haben: ein dataset und einen Vektor mit den indices der Elemente aus dem Datensatz, die in die Bootstrap-Stichprobe gelangt sind.

Wenn wir CIs für mehrere Statistiken gleichzeitig berechnen möchten, muss unsere Funktion sie als einen einzelnen Vektor zurückgeben.

Für unser Beispiel könnte das so aussehen:

 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 wählt die gewünschten Elemente (deren Nummern in indices stehen) aus data und berechnet den Korrelationskoeffizienten der ersten beiden Spalten (method='s' wählt Spearman, method='p' ergibt Pearson) sowie deren Mediane.

Wir können auch zusätzliche Argumente erlauben, z. B. den Korrelations-Typ auswählbar machen:

 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}

Um das Tutorial allgemeiner zu halten, nutze ich die letztere Version von foo.

Jetzt können wir die Funktion boot verwenden. Wir übergeben den Namen des Datensatzes, die gerade erstellte Funktion, die Anzahl der Wiederholungen (R) und alle zusätzlichen Argumente unserer Funktion (wie cor.type). Unten setze ich set.seed für reproduzierbare Ergebnisse.

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

Die Funktion boot gibt ein Objekt der Klasse ... (genau!) boot zurück. Es hat zwei interessante Elemente. $t enthält die R Werte unserer Statistik(en), die durch das Bootstrap-Verfahren erzeugt wurden (Bootstrap-Realisationen von 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 enthält die Werte unserer Statistik(en) im ursprünglichen, vollständigen Datensatz:

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

Wenn wir das boot-Objekt ausgeben, erhalten wir weitere Informationen:

 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 entspricht $t0. bias ist die Differenz zwischen dem Mittelwert der Bootstrap-Realisationen (den Werten aus $t), also der Bootstrap-Schätzer für T, und dem Wert im Originaldatensatz (aus $t0).

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

std. error ist der Standardfehler des Bootstrap-Schätzers und entspricht der Standardabweichung der Bootstrap-Realisationen.

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

Verschiedene Arten von Bootstrap-CIs

Bevor wir CIs bestimmen, lohnt sich immer ein Blick auf die Verteilung der Bootstrap-Realisationen. Wir können die Funktion plot verwenden, mit index als Angabe, auf welche der in foo berechneten Statistiken wir schauen wollen. Hier ist index=1 Spearmans Korrelationskoeffizient zwischen Kelchlänge und -breite, index=2 ist der Median der Kelchlänge und index=3 der Median der Kelchbreite.

plot(myBootstrap, index=1)\n

Die Verteilung der Bootstrap-Korrelationskoeffizienten wirkt recht normalverteilt. Bestimmen wir dafür ein CI. Wir nutzen boot.ci. Standardmäßig werden 95%-CIs berechnet, das lässt sich mit dem Parameter conf ändern.

 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 liefert 5 Arten von Bootstrap-CIs. Eine davon, das studentisierte Intervall, ist speziell. Es benötigt eine Schätzung der Bootstrap-Varianz. Da wir diese nicht angegeben haben, gibt R die Warnung aus: bootstrap variances needed for studentized intervals. Varianzschätzungen lassen sich mit einem Bootstrap zweiter Stufe oder (einfacher) per Jackknife gewinnen. Das sprengt hier den Rahmen, daher konzentrieren wir uns auf die übrigen vier CI-Typen.

Wenn wir nicht alle sehen möchten, können wir im Argument type auswählen. Mögliche Werte sind norm, basic, stud, perc, bca oder ein Vektor daraus.

 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

Die Funktion boot.ci erzeugt ein Objekt der Klasse ... (genau!) bootci. Dessen Elemente heißen wie die im type-Argument gewählten CI-Typen. $norm ist ein Vektor mit drei Elementen: Konfidenzniveau sowie untere und obere Grenze.

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

$basic, $stud, $perc und $bca sind Vektoren mit fünf Elementen und enthalten zusätzlich die zur CI-Berechnung verwendeten Perzentile (dazu gleich mehr):

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

Ein bisschen Notation (sorry!)

Um die verschiedenen CI-Typen einzuordnen, brauchen wir etwas Notation. Sei also:

\n
• t^⋆ der Bootstrap-Schätzer (Mittelwert der Bootstrap-Realisationen),
\n
• t₀ der Wert unserer Statistik im Originaldatensatz,
\n
• se^⋆ der Standardfehler des Bootstrap-Schätzers,
\n
• b die Verzerrung (Bias) des Bootstrap-Schätzers, b = t^⋆ − t₀,
\n
• α das Konfidenzniveau, typischerweise α = 0.95,
\n
• z_α das $1-\\frac \\alpha 2$-Quantil der Standardnormalverteilung,
\n
• θ_α das α-Perzentil der Verteilung der Bootstrap-Realisationen.
\n

Perzentil-CI

Mit dieser Notation lautet das Perzentil-CI:

Das heißt: Wir nehmen einfach die entsprechenden Perzentile. Nicht mehr, nicht weniger.

Normal-CI

Ein klassisches Wald-CI wäre:

im Bootstrap-Fall sollten wir jedoch den Bias korrigieren. Es wird zu:

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

Basic-CI

Das Perzentil-CI wird generell nicht empfohlen, da es bei schiefen oder seltsam auslaufenden Verteilungen schwach performt. Das Basic-CI (auch pivotal oder empirical CI) ist robuster. Die Idee: Wir berechnen die Differenzen zwischen jeder Bootstrap-Replikation und t₀ und verwenden die Perzentile dieser Verteilung. Details findest du z. B. in L. Wassermans All of Statistics

Die endgültige Formel für das Basic-CI ist:

BC_α-CI

BC_α steht für bias-corrected, accelerated. Die Formel ist nicht besonders kompliziert, aber wenig intuitiv, daher spare ich sie hier aus. Wenn du Details möchtest, lies den Artikel von Thomas J. DiCiccio und Bradley Efron.

Die im Namen erwähnte Beschleunigung erfordert spezielle Perzentile der Bootstrap-Realisationen. Es kann passieren, dass dies extreme Perzentile sind, möglicherweise Ausreißer. BC_α kann in solchen Fällen instabil sein.

Schauen wir uns das BC_α-CI für den Median der Kelchblattbreite an. Im Originaldatensatz ist dieser Median genau 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

Wir erhalten ein BC_α-CI von (2.9, 2.9). Das ist ungewöhnlich, aber R hat uns gewarnt: extreme order statistics used as endpoints. Schauen wir uns an, was passiert ist:

 plot(myBootstrap, index=3)

Die Verteilung der Bootstrap-Realisationen ist ungewöhnlich. Die große Mehrheit (über 90%) sind 3er-Werte.

 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

In so einer Situation springt die „Beschleunigung\" im BC_α-Verfahren leicht zu Extremwerten. Hier ist das Perzentil-CI, das die Grundlage für das BC_α-CI bildet, (3, 3.1). Wenn wir es nun „nach links\" erweitern, müssen wir das 0,002-Perzentil verwenden. Die Erweiterung „nach rechts\" geht ebenfalls an die Extreme: das 0,997-Perzentil.

Replikation der Ergebnisse

Manchmal müssen Bootstrap-Replikationen wiederhergestellt werden. Wenn wir R nutzen können, ist das dank set.seed unproblematisch. Die Replikation in anderer Software wäre jedoch deutlich aufwendiger. Und bei großem R ist auch eine Neuberechnung in R nicht immer praktikabel (z. B. aus Zeitgründen).

Wir können das Problem lösen, indem wir die Indizes der Elemente des Originaldatensatzes speichern, die jede Bootstrap-Stichprobe bilden. Genau das macht die Funktion boot.array (mit dem Argument indices=T).

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

Jede Zeile entspricht einer Bootstrap-Stichprobe. Unsere erste Stichprobe enthält beispielsweise folgende Elemente:

 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

Mit indices=F (Standard) bekommen wir die Antwort auf die Frage: „Wie oft kam jedes Element des Originaldatensatzes in jeder Bootstrap-Stichprobe vor?\". Zum Beispiel in der ersten Stichprobe: Das 1. Element kam einmal vor, das 2. gar nicht, das 3. einmal, das 4. zweimal usw.

 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

Mit so einer Tabelle können wir die Bootstrap-Realisationen außerhalb von R rekonstruieren. Oder auch in R selbst, wenn wir set.seed nicht nutzen und die Rechnungen nicht erneut durchführen wollen.

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

Prüfen wir, ob die Ergebnisse identisch sind:

 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

Ja, sind sie!

Wenn du mehr über Machine Learning in R lernen möchtest, sieh dir den DataCamp-Kurs Machine Learning Toolbox an und lies das Tutorial Machine Learning in R für Einsteiger.