Zastosowanie standardowego błędu dystrybucji bootstrap

19

(w razie potrzeby zignoruj ​​kod R, ponieważ moje główne pytanie jest niezależne od języka)

Jeśli chcę spojrzeć na zmienność prostej statystyki (np. Średnia), wiem, że mogę to zrobić za pomocą teorii:

x = rnorm(50)

# Estimate standard error from theory
summary(lm(x~1))
# same as...
sd(x) / sqrt(length(x))

lub z bootstrap jak:

library(boot)

# Estimate standard error from bootstrap
(x.bs = boot(x, function(x, inds) mean(x[inds]), 1000))
# which is simply the standard *deviation* of the bootstrap distribution...
sd(x.bs$t)

Zastanawiam się jednak, czy użyteczne / prawidłowe (?) Może być sprawdzenie standardowego błędu dystrybucji bootstrap w określonych sytuacjach? Sytuacja, z którą mam do czynienia, to stosunkowo głośna funkcja nieliniowa, taka jak:

# Simulate dataset
set.seed(12345)
n   = 100
x   = runif(n, 0, 20)
y   = SSasymp(x, 5, 1, -1) + rnorm(n, sd=2)
dat = data.frame(x, y)

W tym przypadku model nawet nie jest zbieżny przy użyciu oryginalnego zestawu danych,

> (fit = nls(y ~ SSasymp(x, Asym, R0, lrc), dat))
Error in numericDeriv(form[[3L]], names(ind), env) : 
  Missing value or an infinity produced when evaluating the model

więc statystyki, którymi jestem zainteresowany, są bardziej ustabilizowanymi oszacowaniami tych parametrów nls - być może ich średnimi w wielu replikacjach bootstrap.

# Obtain mean bootstrap nls parameter estimates
fit.bs = boot(dat, function(dat, inds)
              tryCatch(coef(nls(y ~ SSasymp(x, Asym, R0, lrc), dat[inds, ])),
                       error=function(e) c(NA, NA, NA)), 100)
pars = colMeans(fit.bs$t, na.rm=T)

Oto one znajdują się w parku, w którym wykorzystałem symulację oryginalnych danych:

> pars
[1]  5.606190  1.859591 -1.390816

Wersja drukowana wygląda następująco:

# Plot
with(dat, plot(x, y))

newx = seq(min(x), max(x), len=100)
lines(newx, SSasymp(newx, pars[1], pars[2], pars[3]))

lines(newx, SSasymp(newx, 5, 1, -1), col='red')
legend('bottomright', c('Actual', 'Predicted'), bty='n', lty=1, col=2:1)

wprowadź opis zdjęcia tutaj

Teraz, jeśli chcę zmienność tych ustabilizowanych oszacowań parametrów, myślę, że mogę, zakładając normalność tego rozkładu ładowania początkowego, po prostu obliczyć ich standardowe błędy:

> apply(fit.bs$t, 2, function(x) sd(x, na.rm=T) / sqrt(length(na.omit(x))))
[1] 0.08369921 0.17230957 0.08386824

Czy to rozsądne podejście? Czy istnieje lepsze ogólne podejście do wnioskowania na temat parametrów takich niestabilnych modeli nieliniowych? (Przypuszczam, że mógłbym zamiast tego zrobić drugą warstwę ponownego próbkowania tutaj, zamiast polegać na teorii w ostatnim bicie, ale może to zająć dużo czasu w zależności od modelu. Mimo to nie jestem pewien, czy te standardowe błędy przydadzą się do wszystkiego, ponieważ zbliżyłyby się do 0, jeśli tylko zwiększę liczbę replik bootstrap.)

Wielkie dzięki, a przy okazji, jestem inżynierem, więc proszę wybacz mi, że jestem tutaj względnym nowicjuszem.

John Colby
źródło

Odpowiedzi:

13

W tym pytaniu jest kilka problemów. Po pierwsze, pojawia się pytanie, czy uśrednione wartości bootstrapów będą sensownymi estymatorami, nawet jeśli niektórych indywidualnych estymatorów bootstrapów nie da się obliczyć (brak zbieżności, brak rozwiązań). Po drugie, biorąc pod uwagę, że estymatory ładowania początkowego są sensowne, pojawia się pytanie, jak uzyskać przedziały ufności, a może po prostu standardowe błędy dla tych szacunków.

--

Celem tego pytania jest jednak uzyskanie oszacowań nawet w przypadkach, w których algorytm obliczania oszacowań może czasami zawieść lub gdy oszacowanie jest czasami niezdefiniowane. Ogólnym podejściem jest problem:

  • Uśrednianie oszacowań bootstrapowania podczas ślepego wyrzucania próbek bootstrapowania, dla których oszacowań nie da się obliczyć, ogólnie daje tendencyjne wyniki.

To, jak poważny jest ogólny problem, zależy od kilku rzeczy. Na przykład, jak często oszacowanie nie jest obliczalne i czy rozkład warunkowy próbki, biorąc pod uwagę, że oszacowanie nie jest obliczalne, różni się od rozkładu warunkowego próbki, biorąc pod uwagę, że oszacowanie jest obliczalne. Nie zalecałbym użycia tej metody.

Xθ^θ^(X)Y

θ~(X)=mi(θ^(Y)X,ZA(X))
ZA(X)Xθ^(Y)NA-XZA(X)θ~(X)

θ^(Y)XZA(X)θ~(X)

ZA(X)θ~(X)

Edytuj :

Bardzo ładny artykuł Oszacowanie i dokładność po wybraniu modelu przez Efron podaje ogólną metodę szacowania błędu standardowego w workowanym estymatorze bez użycia drugiej warstwy ładowania początkowego. Artykuł nie dotyczy wyraźnie estymatorów, które czasami nie są obliczalne.

NRH
źródło
Dzięki za wspaniałą odpowiedź. Szczególnie ważna jest kwestia stronniczości . Możesz sobie wyobrazić ekstremalny przypadek, w którym chmura punktów jest całkowicie jednolita, z wyjątkiem jednego zestawu odległych punktów, które bardzo dobrze pasują do modelu. Zdecydowana większość nlsnapadów może się nie powieść, ale z tych, które się zbiegają, odchylenie będzie ogromne, a przewidywane standardowe błędy / współczynniki CI fałszywie małe. nlsBootstosuje ad hoc wymóg 50% udanych dopasowań, ale zgadzam się z tobą, że (nie) podobieństwo rozkładów warunkowych jest równie ważne.
John Colby,
(Postaram się dać ci premię jutro, jeśli ta strona pozwoli mi na SO)
John Colby