(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)
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.
źródło
nls
napadó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.nlsBoot
stosuje 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.