Wie funktioniert das Bootstrapping in R?

22

Ich habe mir das Boot-Paket in R angeschaut und obwohl ich eine Reihe guter Grundlagen für die Verwendung gefunden habe, muss ich noch etwas finden, das genau beschreibt, was "hinter den Kulissen" passiert. In diesem Beispiel wird beispielsweise gezeigt , wie Standard-Regressionskoeffizienten als Ausgangspunkt für eine Bootstrap-Regression verwendet werden. Es wird jedoch nicht erläutert, wie die Bootstrap-Prozedur die Bootstrap-Regressionskoeffizienten ableitet. Es scheint, als ob es eine Art von iterativem Prozess gibt, aber ich kann nicht genau herausfinden, was los ist.

r regression bootstrap regression-coefficients zgall1
quelle
3
Es ist lange her, dass ich es das letzte Mal geöffnet habe, daher weiß ich nicht, ob es Ihre Frage beantworten wird, aber das Boot-Paket basiert insbesondere auf den in Davison, AC, & Hinkley, DV (1997) beschriebenen Methoden. Bootstrap-Methoden und ihre Anwendung. Cambridge: Cambridge University Press. (Der Verweis wird auch in der Hilfedatei für das Boot- Paket zitiert .)
Gala

Antworten:

34

Es gibt verschiedene "Varianten" oder Formen des Bootstraps (z. B. nicht parametrisch, parametrisch, Restresampling und viele mehr). Der Bootstrap im Beispiel wird als nicht parametrischer Bootstrap oder Resampling von Groß- und Kleinschreibung bezeichnet (siehe hier , hier , hier und hier für Anwendungen in der Regression). Die Grundidee ist, dass Sie Ihre Probe als Population behandeln und wiederholt neue Proben mit Ersatz daraus ziehen . Alle ursprünglichen Beobachtungen haben die gleiche Wahrscheinlichkeit, in die neue Stichprobe aufgenommen zu werden. Anschließend berechnen und speichern Sie die interessierende (n) Statistik (en). Dies kann der Mittelwert, der Median oder der Regressionskoeffizient unter Verwendung der neu gezogenen Stichprobe sein. Dies wird mal wiederholt . In jeder Iteration werden einige Beobachtungen aus Ihrer Originalprobe mehrmals gezogen, während einige Beobachtungen möglicherweise überhaupt nicht gezogen werden. Nach n Iterationen haben Sie n Bootstrap-Schätzungen der interessierenden Statistik (en) gespeichert (z. B. wenn n = 1000 und die interessierende Statistik der Mittelwert ist, haben Sie 1000 Bootstrap-Schätzungen des Mittelwerts). Zuletzt werden zusammenfassende Statistiken wie der Mittelwert, der Median und die Standardabweichung der n Bootstrap-Schätzungen berechnet.nnnn=1000n

Bootstrapping wird häufig verwendet für:

  1. Berechnung der Konfidenzintervalle (und Schätzung der Standardfehler)
  2. Schätzung der Verzerrung der Punktschätzungen

Es gibt verschiedene Methoden zum Berechnen von Konfidenzintervallen basierend auf den Bootstrap-Beispielen ( dieses Dokument enthält Erläuterungen und Anleitungen). Eine sehr einfache Methode zur Berechnung eines 95% -Konfidenzintervalls besteht darin, nur die empirischen 2,5- und 97,5-Perzentile der Bootstrap-Beispiele zu berechnen (dieses Intervall wird als Bootstrap-Perzentilintervall bezeichnet; siehe Code unten). Die einfache Perzentilintervallmethode wird in der Praxis nur selten verwendet, da es bessere Methoden gibt, wie zum Beispiel das vorspannungskorrigierte und beschleunigte Bootstrap (BCa). BCa-Intervalle passen sich sowohl der Vorspannung als auch der Neigung in der Bootstrap-Verteilung an.

n

Lassen Sie uns das Beispiel von der Website wiederholen, aber unsere eigene Schleife verwenden, die die oben skizzierten Ideen einbezieht (wiederholt mit Ersetzung zeichnen):

#-----------------------------------------------------------------------------
# Load packages
#-----------------------------------------------------------------------------

require(ggplot2)
require(pscl)
require(MASS)
require(boot)

#-----------------------------------------------------------------------------
# Load data
#-----------------------------------------------------------------------------

zinb <- read.csv("http://www.ats.ucla.edu/stat/data/fish.csv")
zinb <- within(zinb, {
  nofish <- factor(nofish)
  livebait <- factor(livebait)
  camper <- factor(camper)
})

#-----------------------------------------------------------------------------
# Calculate zero-inflated regression
#-----------------------------------------------------------------------------

m1 <- zeroinfl(count ~ child + camper | persons, data = zinb,
               dist = "negbin", EM = TRUE)

#-----------------------------------------------------------------------------
# Store the original regression coefficients
#-----------------------------------------------------------------------------

original.estimates <- as.vector(t(do.call(rbind, coef(summary(m1)))[, 1:2]))

#-----------------------------------------------------------------------------
# Set the number of replications
#-----------------------------------------------------------------------------

n.sim <- 2000

#-----------------------------------------------------------------------------
# Set up a matrix to store the results
#-----------------------------------------------------------------------------

store.matrix <- matrix(NA, nrow=n.sim, ncol=12)

#-----------------------------------------------------------------------------
# The loop
#-----------------------------------------------------------------------------

set.seed(123)

for(i in 1:n.sim) {

  #-----------------------------------------------------------------------------
  # Draw the observations WITH replacement
  #-----------------------------------------------------------------------------

  data.new <- zinb[sample(1:dim(zinb)[1], dim(zinb)[1], replace=TRUE),]

  #-----------------------------------------------------------------------------
  # Calculate the model with this "new" data
  #-----------------------------------------------------------------------------

  m <- zeroinfl(count ~ child + camper | persons,
                data = data.new, dist = "negbin",
                start = list(count = c(1.3711, -1.5152, 0.879),
                             zero = c(1.6028, -1.6663)))

  #-----------------------------------------------------------------------------
  # Store the results
  #-----------------------------------------------------------------------------

  store.matrix[i, ] <- as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))

}


#-----------------------------------------------------------------------------
# Save the means, medians and SDs of the bootstrapped statistics
#-----------------------------------------------------------------------------

boot.means <- colMeans(store.matrix, na.rm=T)

boot.medians <- apply(store.matrix,2,median, na.rm=T)

boot.sds <- apply(store.matrix,2,sd, na.rm=T)

#-----------------------------------------------------------------------------
# The bootstrap bias is the difference between the mean bootstrap estimates
# and the original estimates
#-----------------------------------------------------------------------------

boot.bias <- colMeans(store.matrix, na.rm=T) - original.estimates

#-----------------------------------------------------------------------------
# Basic bootstrap CIs based on the empirical quantiles
#-----------------------------------------------------------------------------

conf.mat <- matrix(apply(store.matrix, 2 ,quantile, c(0.025, 0.975), na.rm=T),
ncol=2, byrow=TRUE)
colnames(conf.mat) <- c("95%-CI Lower", "95%-CI Upper")

Und hier ist unsere Übersichtstabelle:

#-----------------------------------------------------------------------------
# Set up summary data frame
#-----------------------------------------------------------------------------

summary.frame <- data.frame(mean=boot.means, median=boot.medians,
sd=boot.sds, bias=boot.bias, "CI_lower"=conf.mat[,1], "CI_upper"=conf.mat[,2])

summary.frame

      mean  median       sd       bias CI_lower CI_upper
1   1.2998  1.3013  0.39674 -0.0712912  0.51960   2.0605
2   0.2527  0.2486  0.03208 -0.0034461  0.19898   0.3229
3  -1.5662 -1.5572  0.26220 -0.0509239 -2.12900  -1.0920
4   0.2005  0.1986  0.01949  0.0049019  0.16744   0.2418
5   0.9544  0.9252  0.48915  0.0753405  0.03493   1.9025
6   0.2702  0.2688  0.02043  0.0009583  0.23272   0.3137
7  -0.8997 -0.9082  0.22174  0.0856793 -1.30664  -0.4380
8   0.1789  0.1781  0.01667  0.0029513  0.14494   0.2140
9   2.0683  1.7719  1.59102  0.4654898  0.44150   8.0471
10  4.0209  0.8270 13.23434  3.1845710  0.58114  57.6417
11 -2.0969 -1.6717  1.56311 -0.4306844 -8.43440  -1.1156
12  3.8660  0.6435 13.27525  3.1870642  0.33631  57.6062

Einige Erklärungen

  • Die Differenz zwischen dem Mittelwert der Bootstrap-Schätzungen und den ursprünglichen Schätzungen wird als "Bias" in der Ausgabe von bezeichnet boot
  • Was die Ausgabe von bootAufrufen "Standardfehler" ist die Standardabweichung der Bootstrap-Schätzungen

Vergleichen Sie es mit der Ausgabe von boot:

#-----------------------------------------------------------------------------
# Compare with boot output and confidence intervals
#-----------------------------------------------------------------------------

set.seed(10)
res <- boot(zinb, f, R = 2000, parallel = "snow", ncpus = 4)

res

Bootstrap Statistics :
       original       bias    std. error
t1*   1.3710504 -0.076735010  0.39842905
t2*   0.2561136 -0.003127401  0.03172301
t3*  -1.5152609 -0.064110745  0.26554358
t4*   0.1955916  0.005819378  0.01933571
t5*   0.8790522  0.083866901  0.49476780
t6*   0.2692734  0.001475496  0.01957823
t7*  -0.9853566  0.083186595  0.22384444
t8*   0.1759504  0.002507872  0.01648298
t9*   1.6031354  0.482973831  1.58603356
t10*  0.8365225  3.240981223 13.86307093
t11* -1.6665917 -0.453059768  1.55143344
t12*  0.6793077  3.247826469 13.90167954

perc.cis <- matrix(NA, nrow=dim(res$t)[2], ncol=2)
    for( i in 1:dim(res$t)[2] ) {
  perc.cis[i,] <- boot.ci(res, conf=0.95, type="perc", index=i)$percent[4:5] 
}
colnames(perc.cis) <- c("95%-CI Lower", "95%-CI Upper")

perc.cis 

      95%-CI Lower 95%-CI Upper
 [1,]      0.52240       2.1035
 [2,]      0.19984       0.3220
 [3,]     -2.12820      -1.1012
 [4,]      0.16754       0.2430
 [5,]      0.04817       1.9084
 [6,]      0.23401       0.3124
 [7,]     -1.29964      -0.4314
 [8,]      0.14517       0.2149
 [9,]      0.29993       8.0463
[10,]      0.57248      56.6710
[11,]     -8.64798      -1.1088
[12,]      0.33048      56.6702

#-----------------------------------------------------------------------------
# Our summary table
#-----------------------------------------------------------------------------

summary.frame

      mean  median       sd       bias CI_lower CI_upper
1   1.2998  1.3013  0.39674 -0.0712912  0.51960   2.0605
2   0.2527  0.2486  0.03208 -0.0034461  0.19898   0.3229
3  -1.5662 -1.5572  0.26220 -0.0509239 -2.12900  -1.0920
4   0.2005  0.1986  0.01949  0.0049019  0.16744   0.2418
5   0.9544  0.9252  0.48915  0.0753405  0.03493   1.9025
6   0.2702  0.2688  0.02043  0.0009583  0.23272   0.3137
7  -0.8997 -0.9082  0.22174  0.0856793 -1.30664  -0.4380
8   0.1789  0.1781  0.01667  0.0029513  0.14494   0.2140
9   2.0683  1.7719  1.59102  0.4654898  0.44150   8.0471
10  4.0209  0.8270 13.23434  3.1845710  0.58114  57.6417
11 -2.0969 -1.6717  1.56311 -0.4306844 -8.43440  -1.1156
12  3.8660  0.6435 13.27525  3.1870642  0.33631  57.6062

Vergleichen Sie die Spalten "bias" und "std. Error" mit der Spalte "sd" unserer eigenen Übersichtstabelle. Unsere 95% -Konfidenzintervalle sind den Konfidenzintervallen, die boot.cimit der Perzentilmethode berechnet wurden , sehr ähnlich (jedoch nicht alle: Betrachten Sie die untere Grenze des Parameters mit Index 9).

COOLSerdash
quelle
Danke für die ausführliche Antwort. Wollen Sie damit sagen, dass es sich bei den Koeffizienten um den Durchschnitt der 2000 erzeugten Koeffizientensätze handelt?
zgall1
4
t
"Die Grundidee ist, dass Sie Ihre Stichprobe als Population behandeln und wiederholt neue Stichproben mit Ersatz daraus ziehen" - wie können Sie die Größe der neuen Stichproben bestimmen?
Sinusx
1
@ Sinusx Normalerweise zeichnen Sie Samples mit der gleichen Größe wie das Original-Sample. Die entscheidende Idee ist, die Probe mit Ersatz zu zeichnen . So werden einige Werte aus dem Original-Sample mehrfach und einige Werte gar nicht gezeichnet.
COOLSerdash
6

Sie sollten sich auf die Funktion konzentrieren, die bootals "Statistik" -Parameter übergeben wird, und deren Aufbau beachten.

f <- function(data, i) {
  require(pscl)
  m <- zeroinfl(count ~ child + camper | persons,
    data = data[i, ], dist = "negbin",
    start = list(count = c(1.3711, -1.5152, 0.879), zero = c(1.6028, -1.6663)))
  as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
}

Das Argument "data" empfängt einen gesamten Datenrahmen, das Argument "i" jedoch eine Stichprobe der vom "boot" generierten und von 1: NROW (data) übernommenen Zeilenindizes. Wie Sie aus diesem Code ersehen können, wird "i" verwendet, um ein Neo-Sample zu erstellen, das an übergeben wird, zeroinlund dann werden nur ausgewählte Teile der Ergebnisse zurückgegeben.

Stellen wir uns vor, dass "i" {1,2,3,3,3,6,7,7,10} ist. Die Funktion "[" gibt nur die Zeilen mit 3 Kopien von Zeile 3 und 2 Kopien von Zeile 7 zeroinl()zurück. bootDies wäre die Basis für eine einzelne Berechnung, und dann werden die Koeffizienten als Ergebnis dieser Replikation des Prozesses zurückgegeben. Die Anzahl solcher Replikate wird durch den "R" -Parameter gesteuert.

Da statisticin diesem Fall nur die Regressionskoeffizienten zurückgegeben werden, gibt die bootFunktion diese akkumulierten Koeffizienten als Wert von "t" zurück. Weitere Vergleiche können mit anderen Boot-Package-Funktionen durchgeführt werden.

DWin
quelle