Jak uzyskać standardowe błędy z regresji danych zliczanych z zerową liczbą R? [Zamknięte]

9

Poniższy kod

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

tworzy 3-kolumnowy - data.framePredictNew, dopasowane wartości, błędy standardowe i rezydualny element skali.

Idealne ... Jednak przy użyciu modelu wyposażonego w zeroinfl {pscl}:

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

lub

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

tworzy tylko jeden wektor kolumn z dopasowanymi wartościami. Byłbym jednak bardzo zainteresowany standardowymi błędami. Wszystko, co przeczytałem, mówi, że powinny być produkowane.

(Kod został nieco uproszczony, w rzeczywistości mam cztery zmienne i przesunięcie - bez problemów z SE predict.glmi se.fit = TRUEprodukującymi SE).

KalahariKev
źródło
5
Spójrz na ten wątek w R-Help: stat.ethz.ch/pipermail/r-help/2008-December/thread.html#182806 (szczególnie wiadomość od Achima Zeileisa, który zapewnia kod do robienia tego, co myślę, że jesteś próbuje zrobić). Nie wygląda na to, predict()aby zeroinfl()w tej chwili zaimplementowano w tej funkcji standardowe błędy .
smillig
Dzięki, ten kod wydawał się całkiem rozsądny. Inni powinni zauważyć, że parametr przewidywana () w nowej funkcji zeroinfl.predict dla se.fit = TRUE został zmieniony na se = TRUE, aby wyodrębnić przewidywane przedziały i se
KalahariKev

Odpowiedzi:

4

Według mojej wiedzy predictmetoda uzyskiwania wyników zeroinflnie obejmuje standardowych błędów. Jeśli Twoim celem jest zbudowanie przedziałów ufności, jedną z atrakcyjnych alternatyw jest użycie ładowania początkowego. Mówię atrakcyjnie, ponieważ ładowanie początkowe może być bardziej niezawodne (przy utracie wydajności, jeśli wszystkie założenia dotyczące SE zostaną spełnione).

Oto trochę z grubsza kod do robienia tego, co chcesz. Nie będzie działać dokładnie, ale mam nadzieję, że możesz wprowadzić niezbędne poprawki.

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

Narysowałem ten kod z dwóch stron, które napisałem: jednego parametru ładowania początkowego z regresji poissona z zeroinfl zerowym napełnieniem zerowym poissonem z napompowaniem zerowym i jednego demonstrującego, jak uzyskać przedziały ufności bootstrapowania dla przewidywanych wartości z ujemnego dwumianowego modelu zerowego obciętego Zero obcięty ujemny dwumianowy . W połączeniu, mam nadzieję, że daje to wystarczające przykłady, aby uruchomić go z przewidywanymi wartościami z zerowo napompowanego poissona. Możesz również uzyskać pomysły na grafikę :)

Jozuego
źródło
Próbowałem dostosować kod do modelu dwumianowego o obciętym zeru w pakiecie VGAM, ale otrzymałem błąd. Czy powinienem utworzyć nowe pytanie tutaj w CV i link tutaj? Byłbym bardzo wdzięczny za twoją pomoc. W szczególności, jest to błąd pojawia się: Error in X.vlm.save %*% coefstart : non-conformable arguments.
Raphael K,