(abaikan kode R jika perlu, karena pertanyaan utama saya adalah bahasa-independen)
Jika saya ingin melihat variabilitas statistik sederhana (mis: rata-rata), saya tahu saya bisa melakukannya melalui teori seperti:
x = rnorm(50)
# Estimate standard error from theory
summary(lm(x~1))
# same as...
sd(x) / sqrt(length(x))
atau dengan bootstrap seperti:
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)
Namun, yang saya pikirkan adalah, bisakah berguna / valid (?) Untuk melihat kesalahan standar dari distribusi bootstrap dalam situasi tertentu? Situasi yang saya hadapi adalah fungsi nonlinear yang relatif bising, seperti:
# 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)
Di sini model bahkan tidak konvergen menggunakan set data asli,
> (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
jadi statistik yang saya minati adalah perkiraan yang lebih stabil dari parameter nls ini - mungkin artinya di sejumlah replikasi 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)
Inilah ini, memang, di taman bola dari apa yang saya gunakan untuk mensimulasikan data asli:
> pars
[1] 5.606190 1.859591 -1.390816
Versi yang diplot terlihat seperti:
# 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)
Sekarang, jika saya ingin variabilitas estimasi parameter yang distabilkan ini , saya rasa saya bisa, dengan asumsi normalitas distribusi bootstrap ini, cukup hitung kesalahan standar mereka:
> apply(fit.bs$t, 2, function(x) sd(x, na.rm=T) / sqrt(length(na.omit(x))))
[1] 0.08369921 0.17230957 0.08386824
Apakah ini pendekatan yang masuk akal? Apakah ada pendekatan umum yang lebih baik untuk inferensi pada parameter model nonlinier tidak stabil seperti ini? (Saya kira saya bisa melakukan layer kedua dari resampling di sini, daripada mengandalkan teori untuk bit terakhir, tapi itu mungkin membutuhkan banyak waktu tergantung pada model. Bahkan, saya tidak yakin apakah kesalahan standar ini akan berguna untuk apa pun, karena mereka akan mendekati 0 jika saya hanya menambah jumlah replikasi bootstrap.)
Terima kasih banyak, dan, ngomong-ngomong, saya seorang insinyur jadi tolong maafkan saya karena saya seorang pemula di sekitar sini.
sumber
nls
kecocokan mungkin gagal, tetapi, dari pencocokan itu, bias akan sangat besar dan kesalahan / CI standar yang diprediksi sangat kecil.nlsBoot
menggunakan persyaratan ad hoc 50% berhasil, tetapi saya setuju dengan Anda bahwa (dis) kesamaan distribusi kondisional sama-sama memprihatinkan.