Jak faktycznie działa ładowanie w R?

Patrzyłem na pakiet rozruchowy w R i chociaż znalazłem kilka dobrych starterów, jak go używać, to jeszcze nie znalazłem niczego, co dokładnie opisuje to, co dzieje się „za kulisami”. Na przykład w tym przykładzie przewodnik pokazuje, jak używać standardowych współczynników regresji jako punktu wyjścia do regresji bootstrapu, ale nie wyjaśnia, co właściwie robi procedura bootstrap, aby uzyskać współczynniki regresji bootstrap. Wygląda na to, że zachodzi pewien proces iteracyjny, ale nie wydaje mi się, żeby dokładnie rozgryźć, co się dzieje.

Istnieje kilka „smaków” lub form bootstrapu (np. Nieparametryczny, parametryczny, resampling resztkowy i wiele innych). Bootstrap w tym przykładzie nazywany jest nieparametrycznym bootstrapem lub resamplingiem wielkości liter (patrz tutaj , tutaj , tutaj i tutaj dla aplikacji w regresji). Podstawową ideą jest to, że traktujesz próbkę jako populację i wielokrotnie pobierasz z niej nowe próbki z wymianą . Wszystkie oryginalne obserwacje mają jednakowe prawdopodobieństwo, że zostaną wciągnięte do nowej próbki. Następnie obliczasz i przechowujesz statystyki będące przedmiotem zainteresowania, może to być średnia, mediana lub współczynniki regresji przy użyciu nowo narysowanej próbki. Jest to powtarzane $n$ razy. W każdej iteracji niektóre obserwacje z oryginalnej próbki są narysowane wiele razy, a niektóre obserwacje mogą w ogóle nie być narysowane. Po $n$ iteracjach masz $n$ zapisanych estymacji ładowania początkowego dla interesujących statystyk (np. Jeśli $n=1000$ a statystyki odsetek są średnią, masz 1000 szacunkowych wartości początkowej średniej). Na koniec obliczane są statystyki podsumowujące, takie jak średnia, mediana i odchylenie standardowe $n$ -oszacowań bootstrap.

Bootstrapping jest często używany do:

1. Obliczanie przedziałów ufności (i szacowanie błędów standardowych)
2. Oszacowanie odchylenia szacunków punktowych

Istnieje kilka metod obliczania przedziałów ufności na podstawie próbek bootstrap ( ten dokument zawiera wyjaśnienia i wskazówki). Jedną z bardzo prostych metod obliczania przedziału ufności 95% jest po prostu obliczenie empirycznych 2,5 i 97,5 percentyla próbek bootstrapu (ten przedział nazywa się interwałem percentyla bootstrap; patrz kod poniżej). Prosta metoda interwału percentyla jest rzadko stosowana w praktyce, ponieważ istnieją lepsze metody, takie jak skorygowany i przyspieszony bootstrap (BCa). Interwały BCa dostosowują się zarówno do odchylenia, jak i skosu w rozkładzie ładowania początkowego.

Odchylenie jest po prostu oszacowane jako różnica między średnią przechowywanych próbek bootstrap i oszacowanie pierwotnego (ów).$n$

Powtórzmy przykład ze strony internetowej, ale wykorzystując własną pętlę zawierającą pomysły, które nakreśliłem powyżej (ciągłe rysowanie z zamiennikiem):

#-----------------------------------------------------------------------------
# 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")

A oto nasza tabela podsumowań:

#-----------------------------------------------------------------------------
# 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

Kilka wyjaśnień

• Różnica między średnią szacunków bootstrap a pierwotnymi szacunkami jest tak zwana „stronniczość” w danych wyjściowych boot
• Wyjściowe bootwywołania „standardowego błędu” to standardowe odchylenie szacunków bootstrapped

Porównaj to z danymi wyjściowymi z 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

Porównaj kolumny „odchylenie” i „błąd standardowy” z kolumną „sd” naszej własnej tabeli podsumowań. Nasze przedziały ufności 95% są bardzo podobne do przedziałów ufności obliczanych przy boot.ciużyciu metody percentyla (choć nie wszystkie: spójrz na dolną granicę parametru o indeksie 9).

Powinieneś skupić się na funkcji, która jest przekazywana bootjako parametr „statystyczny” i zauważyć, jak jest zbudowana.

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]))
}

Argument „dane” otrzyma całą ramkę danych, ale argument „i” otrzyma próbkę indeksów wierszy wygenerowanych przez „boot” i pobranych z 1: NROW (dane). Jak widać z tego kodu, „i” jest następnie używane do utworzenia neopróbki, która jest przekazywana do, zeroinla następnie zwracane są tylko wybrane części jego wyników.

Wyobraźmy sobie, że „i” to {1,2,3,3,3,3,7,7,7,10}. Funkcja „[” zwróci tylko te wiersze z 3 kopiami wiersza 3 i 2 kopiami wiersza 7. To będzie podstawa do pojedynczego zeroinl()obliczenia, a następnie współczynniki zostaną zwrócone bootw wyniku tej repliki procesu. Liczba takich replikacji jest kontrolowana przez parametr „R”.

Ponieważ statisticw tym przypadku zwracane są tylko współczynniki regresji , bootfunkcja zwróci te skumulowane współczynniki jako wartość „t”. Dalsze porównania mogą być wykonywane przez inne funkcje pakietu rozruchowego.