[R-br] Contrastes em modelos de regressão de poisson inflacionados por zeros

[walmes .]

Você não deve esperar ver uma curva exponencial para a média das contagens
em função do tempo e uma sigmoide para a proporção de zeros não Poisson
para os zeros se não simulou os dados com efeito de tempo e não declarou
efeito de tempo na porção binomial. Para que isso aconteça a simulação tem
que estar condizente, tem que ter os efeitos para que as curvas não sejam
as "curvas nulas".

#------------------------------------------------------------------
# Definições da sessão.

rm(list=ls())
require(pscl)
require(VGAM)
require(multcomp)
require(lattice)
require(latticeExtra)

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

da <- expand.grid(tempo=rep(1:10), trat=gl(3,10))

## simula porção Poisson
X <- model.matrix(~trat+tempo, da)
colnames(X)
beta.pois <- c(-1,0.25,0.5,0.25)
eta.pois <- X%*%beta.pois
da$y.pois <- rpois(nrow(X), lambda=exp(eta.pois))
xyplot(y.pois~tempo|trat, da)

## simula porção binomial
X <- model.matrix(~tempo, da)
beta.bern <- c(0,0.3)
eta.bern <- X%*%beta.bern
da$y.bern <- rbinom(nrow(X), size=1, prob=1/(1+exp(-eta.bern)))

xyplot(y.bern~tempo|trat, da)

## monta a contagem com inflação de zeros
da$y <- with(da, ifelse(y.bern==0, 0, y.pois))
xyplot(y~tempo|trat, da)

#------------------------------------------------------------------
# Modelo completo.

compl.mod <- zeroinfl(y~trat+tempo|trat+tempo, data=da)
coef(compl.mod)

#------------------------------------------------------------------
# Predição do modelo considerando as duas porções.

pred <- expand.grid(tempo=rep(1:10), trat=gl(3,1))
X <- model.matrix(~trat+tempo, pred)
i <- grep("^count\\_", names(coef(compl.mod)))
eta <- X%*%coef(compl.mod)[i]
pred$m.pois <- exp(eta)
i <- grep("^zero\\_", names(coef(compl.mod)))
eta <- X%*%coef(compl.mod)[i]
pred$p.zero <- exp(eta)/(1+exp(eta))

xyplot(y~tempo|trat, data=da, jitter.x=TRUE)+
    as.layer(xyplot(m.pois~tempo|trat, data=pred, type="l"))+
    as.layer(xyplot(p.zero~tempo|trat, data=pred,
                    type="l", lty=2, lwd=2))+
    layer(panel.abline(h=1, lty=2))

xyplot(p.zero~tempo|trat, data=pred,
                    type="l", lty=2, lwd=2)

#------------------------------------------------------------------


À disposição.
Walmes.

==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
skype: walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20131108/19a0813a/attachment.html>