Comment obtenir les erreurs standard de la régression des données de comptage gonflées à zéro? [fermé]

9

Le code suivant

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

produit un data.frame--PredictNew à 3 colonnes , les valeurs ajustées, les erreurs standard et un terme d'échelle résiduelle.

Parfait ... Cependant en utilisant un modèle équipé de zeroinfl {pscl}:

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

ou

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

produire un vecteur de colonne unique de valeurs ajustées uniquement. Je serais cependant très désireux d'avoir des erreurs standard. Tout ce que j'ai lu dit qu'ils devraient être produits ..

(Le code a été quelque peu simplifié, j'ai en fait quatre variables et un décalage - pas de sondes avec predict.glmet se.fit = TRUEproduisant des SE.)

KalahariKev
la source
5
Jetez un œil à ce fil sur R-Help: stat.ethz.ch/pipermail/r-help/2008-December/thread.html#182806 (en particulier le message d'Achim Zeileis qui fournit du code pour faire ce que je pense que vous êtes essayer de faire). Il ne semble pas que des erreurs standard aient été implémentées dans la predict()fonction pour le zeroinfl()moment.
smillig
Merci, ce code a semblé produire des résultats assez raisonnables. D'autres devraient noter que le paramètre Predict () dans la nouvelle fonction zeroinfl.predict pour se.fit = TRUE a été changé en se = TRUE, afin d'extraire les intervalles prévus et se
KalahariKev

Réponses:

4

À ma connaissance, le predict méthode des résultats zeroinflne comprend pas les erreurs standard. Si votre objectif est de construire des intervalles de confiance, une alternative intéressante est d'utiliser le bootstrap. Je dis attrayant parce que le bootstrap a le potentiel d'être plus robuste (avec une perte d'efficacité si toutes les hypothèses pour les SE sont remplies).

Voici un code approximatif pour faire ce que vous voulez. Cela ne fonctionnera pas exactement, mais j'espère que vous pourrez apporter les corrections nécessaires.

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

J'ai dessiné ce code à partir de deux pages que j'ai écrites, un paramètres d'amorçage à partir d'une régression de poisson gonflé à zeroinfl zéro avec un poisson gonflé à zéro et un montrant comment obtenir des intervalles de confiance amorcés pour les valeurs prédites à partir d'un modèle binomial négatif à zéro tronqué Binôme négatif à zéro tronqué . Combiné, j'espère que cela vous fournit suffisamment d'exemples pour le faire fonctionner avec les valeurs prédites d'un poisson gonflé à zéro. Vous pouvez également obtenir des idées graphiques :)

Joshua
la source
J'ai essayé d'adapter votre code pour un modèle binomial négatif tronqué zéro dans le package VGAM, mais j'ai reçu une erreur. Dois-je créer une nouvelle question ici sur CV et un lien ici? J'apprécierais vraiment votre aide. Tel est précisément l'erreur que je reçois: Error in X.vlm.save %*% coefstart : non-conformable arguments.
Raphael K