[R-br] Cálculo do deviance das variáveis explicativas em modelo de Poisson inflacionado de zeros

walmes . walmeszeviani em gmail.com
Segunda Dezembro 15 17:15:16 BRST 2014


Seguem sugestões no CMR.

#------------------------------------------------------------------
# Dados artificiais.

da <- expand.grid(trat=gl(4,1), tempo=1:10)
X <- model.matrix(~trat+tempo, da); ncol(X)
## betas <- c(0.1,0.9,0.6,0.3,0.7)
betas <- c(0.1,0.9,0.6,0.3,0) ## Simulando sem efeito de tempo.
eta <- X%*%betas
y1 <- rpois(da$trat, lambda=exp(eta))
y2 <- rbinom(y1, size=1, prob=0.7)
da$y <- y1*y2
str(da)

#------------------------------------------------------------------
# Ajuste do modelo.

m0 <- zeroinfl(y~trat+tempo|trat, data=da)
summary(m0)

length(coef(m0))==dim(vcov(m0))
names(coef(m0))==colnames(vcov(m0))

##-----------------------------------------------------------------------------
## Teste por meio de comparações entre modelos aninhados.

## Testar o efeito de tempo na parte Poisson.
## ~trat+tempo|trat vs ~trat|trat.

m1 <- zeroinfl(y~trat|trat, data=da)

class(m1)
anova.zeroinf <- function(m0, m1){
    ll0 <- c(logLik(m0))
    ll1 <- c(logLik(m1))
    np0 <- length(coef(m0))
    np1 <- length(coef(m1))
    dnp <- abs(np0-np1)
    dll <- 2*abs(ll0-ll1)
    pchi <- pchisq(dll, df=dnp, lower.tail=FALSE)
    list("twice difference em log-likelihood"=dll,
         "difference in number of parameters"=dnp,
         "Pr(>chi)"=pchi)
}

## Por razão de log-verossimilhanças.
anova.zeroinf(m0, m1)

## Usando Wald para testar a mesma hipótese.
names(coef(m0))
L <- rbind(c(0,0,0,0,1,0,0,0,0))

require(multcomp)

## Pela aproximação quadrática da log-verossimilhança.
summary(glht(m0, linfct=L), test=Chisqtest())

## Em modelos lineares de efeito fixo, as duas aboradagens dão mesmo
## resultado porque a aproximação quadrática de uma log-verossimilhança.

​
​À disposição.
Walmes.​
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20141215/2952bc6e/attachment.html>


Mais detalhes sobre a lista de discussão R-br