[R-br] Estimativa do intervalo de confianca em um modelo ZINB

[Humberto Hazin (hotmail)]

Olá  a todos,

 

Tenho 2 problemas se alguém puder me ajudar (Um exemplo dos dados
encontra-se em anexo):

 

1-      como faço para estimar os intervalos de confiança dos coeficientes
de um modelo ZINB. Estou usando o package (pscl). Já tentei por bootstrap
mais ele não estima os intervalos. Segue abaixo o script que usei

 

data(bum2)

F1a<-formula(BUM ~ Year + Quarter + area + Strategy+offset(log(effort1)))

 

options(contrasts=c("contr.treatment","contr.poly"))

ZINB<-zeroinfl(F1a,dist="negbin",link = "logit",data=bum2)

 

boot.zinb <- function(data, indices){

                  require(pscl)

                  data <- data[indices,]      

                  try(mod <- zeroinfl(F1a,

                            data = data,      

                            link = "logit",

                            dist = "negbin",

                            trace = FALSE

                            EM = FALSE))

                  if (exists("mod")) {

                    coef(mod)

                    } else {

                    NA

                    }                                 }

 

system.time(zinb.boot.out <- boot(data = bum2, statistic = boot.zinb, R =
99,sim="permutation"))

 

zinb.boot.out    

boot.ci(zinb.boot.out,type="bca")

 

2-      como faço para estimar os predicts e os intervalos de confiança
apenas da variável  YEAR  de um modelo ZINB. Abaixo o código usado:

 

x<-seq(1980,2008)

newdata=cbind.data.frame(Year=as.factor(x),Quarter=as.factor(rep(1,length(x)
)),area=as.factor(rep(1,length(x))),Strategy=as.factor(rep(4,length(x))),eff
ort1=rep(1000,length(x)))

newdata

fm_zip <- zeroinfl(F1a,dist="negbin",link = "logit",data=bum2)

fit <- predict(fm_zip)

Pearson <- resid(fm_zip, type = "pearson")

VarComp <- resid(fm_zip, type = "response") / Pearson

year <- bum2$Year

bootstrap <- replicate(999, {

   yStar <- pmax(round(fit + sample(Pearson) * VarComp, 0), 0)

 

   predict(zeroinfl(yStar ~ year | 1), newdata = newdata)

})

newdata0 <- newdata

newdata0$fit <- predict(fm_zip, newdata = newdata, type = "response")

newdata0<- t(apply(bootstrap, 1, quantile, c(0.025, 0.975)))

 

 

Desde Ja agradeço a vossa atenção 

 

Humberto

-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20110324/1f968ff2/attachment-0001.html>
-------------- Próxima Parte ----------
Um anexo não-texto foi limpo...
Nome: bum2.csv
Tipo: application/vnd.ms-excel
Tamanho: 17781 bytes
Descrição: não disponível
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20110324/1f968ff2/attachment-0001.xlb>