Bagaimana cara mendapatkan kesalahan standar dari regresi data hitung R nol? [Tutup]

9

Kode berikut

PredictNew <- predict (glm.fit, newdata = Predict, X1 =X1, Y1= Y1, 
                       type = "response", se.fit = TRUE)

menghasilkan 3 kolom data.frame PredictNew, nilai-nilai yang dipasang, kesalahan standar dan istilah skala residual.

Sempurna ... Namun menggunakan model yang dilengkapi dengan zeroinfl {pscl}:

PredictNew <- predict (zeroinfl.fit, newdata = Predict, X1 =X1, Y1= Y1, 
                       type = "response", se.fit = TRUE)

atau

PredictNew <- predict (zeroinfl.fit, newdata = Predict, X1 =X1, Y1= Y1, 
                       type = "response", se.fit = TRUE, MC = 2500, conf = .95))

menghasilkan vektor kolom tunggal dengan nilai yang dipasang saja. Namun, saya sangat ingin memiliki kesalahan standar. Semua yang saya baca mengatakan bahwa itu harus diproduksi ..

(Kode agak disederhanakan, saya benar-benar memiliki empat variabel dan offset - tidak ada masalah dengan predict.glmdan se.fit = TRUEmenghasilkan UK).

KalahariKev
sumber
5
Lihat utas ini di Bantuan-R: stat.ethz.ch/pipermail/r-help/2008-December/thread.html#182806 (khususnya pesan dari Achim Zeileis yang menyediakan kode untuk melakukan apa yang saya pikir Anda lakukan mencoba melakukan). Sepertinya tidak ada kesalahan standar yang diimplementasikan ke dalam predict()fungsi untuk zeroinfl()saat ini.
smillig
Terima kasih, kode itu sepertinya memberikan hasil yang cukup masuk akal. Yang lain harus mencatat bahwa parameter prediksikan () di fungsi zeroinfl.predict baru untuk se.fit = TRUE diubah menjadi se = TRUE, untuk mengekstraksi interval yang diprediksi dan se
KalahariKev

Jawaban:

4

Sepengetahuan saya, predictmetode untuk hasil dari zeroinfltidak termasuk kesalahan standar. Jika tujuan Anda adalah membangun interval kepercayaan, salah satu alternatif menarik adalah menggunakan bootstrap. Saya katakan menarik karena bootstrap berpotensi menjadi lebih kuat (dengan kehilangan efisiensi jika semua asumsi untuk SE terpenuhi).

Berikut ini beberapa kode kasar untuk melakukan apa yang Anda inginkan. Ini tidak akan bekerja dengan tepat, tetapi mudah-mudahan Anda dapat melakukan koreksi yang diperlukan.

## load boot package
require(boot)
## output coefficients from your original model
## these can be used as starting values for your bootstrap model
## to help speed up convergence and the bootstrap
dput(round(coef(zeroinfl.fit, "count"), 3))
dput(round(coef(zeroinfl.fit, "zero"), 3))

## function to pass to the boot function to fit your model
## needs to take data, an index (as the second argument!) and your new data
f <- function(data, i, newdata) {
  require(pscl)
  m <- zeroinfl(count ~ child + camper | persons, data = data[i, ], start = list(count = c(1.598, -1.0428, 0.834), zero = c(1.297, -0.564)))
  mparams <- as.vector(t(do.call(rbind, coef(summary(m)))[, 1:2]))
  yhat <- predict(m, newdata, type = "response")
  return(c(mparams, yhat))    
}

## set the seed and do the bootstrap, make sure to set your number of cpus
## note this requires a fairly recent version of R
set.seed(10)
res <- boot(dat, f, R = 1200, newdata = Predict, parallel = "snow", ncpus = 4)

## get the bootstrapped percentile CIs
## the 10 here is because in my initial example, there were 10 parameters before predicted values
yhat <- t(sapply(10 + (1:nrow(Predict)), function(i) {
  out <- boot.ci(res, index = i, type = c("perc"))
  with(out, c(Est = t0, pLL = percent[4], pUL = percent[5]))
}))

## merge CIs with predicted values
Predict<- cbind(Predict, yhat)

Aku menarik kode ini dari dua halaman saya menulis, salah satu parameter bootstrap dari poisson regresi zero-meningkat dengan zeroinfl poisson Zero-meningkat dan satu berdemonstrasi bagaimana untuk dinyalakan interval kepercayaan untuk nilai-nilai diprediksi dari model binomial negatif nol-terpotong Zero-terpotong binomial negatif . Gabungan, mudah-mudahan itu memberi Anda contoh yang cukup untuk membuatnya bekerja dengan nilai yang diprediksi dari poisson nol-meningkat. Anda juga dapat memperoleh beberapa ide grafik :)

Joshua
sumber
Saya mencoba mengadaptasi kode Anda untuk model binomial negatif terpotong nol dalam paket VGAM, tetapi menerima kesalahan. Haruskah saya membuat pertanyaan baru di sini di CV dan tautan di sini? Saya akan sangat menghargai bantuan Anda dengan ini. Secara khusus, ini adalah kesalahan saya: Error in X.vlm.save %*% coefstart : non-conformable arguments.
Raphael K