Penggunaan kesalahan standar distribusi bootstrap

19

(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)

masukkan deskripsi gambar di sini

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.

John Colby
sumber

Jawaban:

13

Ada beberapa masalah dalam pertanyaan ini. Pertama, ada pertanyaan apakah rata-rata bootstrap akan menjadi penaksir yang masuk akal bahkan ketika beberapa penaksir bootstrap individual tidak dapat dihitung (kurangnya konvergensi, tidak adanya solusi). Kedua, mengingat bahwa estimator yang bootstrap masuk akal, ada pertanyaan tentang bagaimana mendapatkan interval kepercayaan atau mungkin hanya kesalahan standar untuk estimasi ini.

--

Namun, tujuan dari pertanyaan ini adalah untuk menghasilkan estimasi bahkan dalam kasus di mana algoritma untuk menghitung estimasi mungkin gagal sesekali atau di mana estimator kadang-kadang tidak terdefinisi. Sebagai pendekatan umum ada masalah:

  • Rata-rata perkiraan bootstrap sementara secara membabi buta membuang sampel bootstrap yang estimasi tidak dapat dihitung secara umum akan memberikan hasil yang bias.

Seberapa parah masalah umumnya tergantung pada beberapa hal. Misalnya, seberapa sering estimasi tersebut tidak dapat dihitung dan apakah distribusi kondisional dari sampel mengingat bahwa estimasi tersebut tidak dapat dihitung berbeda dari distribusi kondisional sampel yang diberikan bahwa estimasi tersebut dapat dihitung. Saya tidak akan merekomendasikan untuk menggunakan metode ini.

Xθ^θ^(X)Y

θ~(X)=E(θ^(Y)X,SEBUAH(X))
SEBUAH(X)Xθ^(Y)NA-XSEBUAH(X)θ~(X)

θ^(Y)XSEBUAH(X)θ~(X)

SEBUAH(X)θ~(X)

Edit :

Estimasi kertas yang sangat bagus dan Akurasi Setelah Pemilihan Model oleh Efron memberikan metode umum untuk memperkirakan kesalahan standar estimator kantong tanpa menggunakan lapisan kedua bootstrap. Makalah ini tidak berurusan secara eksplisit dengan penduga yang terkadang tidak dapat dihitung.

NRH
sumber
Terima kasih atas jawaban yang luar biasa. Poin bias biasanya diambil dengan baik. Anda dapat membayangkan kasus ekstrim di mana cloud titik benar-benar seragam, simpan untuk satu set titik jauh yang cocok dengan model dengan sangat baik. Sebagian besar nlskecocokan mungkin gagal, tetapi, dari pencocokan itu, bias akan sangat besar dan kesalahan / CI standar yang diprediksi sangat kecil. nlsBootmenggunakan persyaratan ad hoc 50% berhasil, tetapi saya setuju dengan Anda bahwa (dis) kesamaan distribusi kondisional sama-sama memprihatinkan.
John Colby
(Saya akan mencoba memberi Anda bonus besok jika situs ini mengizinkan saya seperti SO)
John Colby