<div dir="ltr"><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">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".<br>
<br><span style="font-family:courier new,monospace">#------------------------------------------------------------------<br># Definições da sessão.<br><br>rm(list=ls())<br>require(pscl)<br>require(VGAM)<br>require(multcomp)<br>
require(lattice)<br>require(latticeExtra)<br><br>#------------------------------------------------------------------<br># Dados artificiais.<br><br>da <- expand.grid(tempo=rep(1:10), trat=gl(3,10))<br><br>## simula porção Poisson<br>
X <- model.matrix(~trat+tempo, da)<br>colnames(X)<br>beta.pois <- c(-1,0.25,0.5,0.25)<br>eta.pois <- X%*%beta.pois<br>da$y.pois <- rpois(nrow(X), lambda=exp(eta.pois))<br>xyplot(y.pois~tempo|trat, da)<br><br>## simula porção binomial<br>
X <- model.matrix(~tempo, da)<br>beta.bern <- c(0,0.3)<br>eta.bern <- X%*%beta.bern<br>da$y.bern <- rbinom(nrow(X), size=1, prob=1/(1+exp(-eta.bern)))<br><br>xyplot(y.bern~tempo|trat, da)<br><br>## monta a contagem com inflação de zeros<br>
da$y <- with(da, ifelse(y.bern==0, 0, y.pois))<br>xyplot(y~tempo|trat, da)<br><br>#------------------------------------------------------------------<br># Modelo completo.<br><br>compl.mod <- zeroinfl(y~trat+tempo|trat+tempo, data=da)<br>
coef(compl.mod)<br><br>#------------------------------------------------------------------<br># Predição do modelo considerando as duas porções.<br><br>pred <- expand.grid(tempo=rep(1:10), trat=gl(3,1))<br>X <- model.matrix(~trat+tempo, pred)<br>
i <- grep("^count\\_", names(coef(compl.mod)))<br>eta <- X%*%coef(compl.mod)[i]<br>pred$m.pois <- exp(eta)<br>i <- grep("^zero\\_", names(coef(compl.mod)))<br>eta <- X%*%coef(compl.mod)[i]<br>
pred$p.zero <- exp(eta)/(1+exp(eta))<br><br>xyplot(y~tempo|trat, data=da, jitter.x=TRUE)+<br> as.layer(xyplot(m.pois~tempo|trat, data=pred, type="l"))+<br> as.layer(xyplot(p.zero~tempo|trat, data=pred,<br>
type="l", lty=2, lwd=2))+<br> layer(panel.abline(h=1, lty=2))<br><br>xyplot(p.zero~tempo|trat, data=pred,<br> type="l", lty=2, lwd=2)<br><br>#------------------------------------------------------------------<br>
</span><br><br>À disposição.<br>Walmes.<br></div><div class="gmail_extra"><br clear="all"><div><div dir="ltr"><span style="font-family:trebuchet ms,sans-serif">==========================================================================</span><br style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">Walmes Marques Zeviani</span><br style="font-family:trebuchet ms,sans-serif"><span style="font-family:trebuchet ms,sans-serif">LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)</span><br style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">Departamento de Estatística - Universidade Federal do Paraná</span><br style="font-family:trebuchet ms,sans-serif"><span style="font-family:trebuchet ms,sans-serif">fone: (+55) 41 3361 3573</span><br style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">skype: walmeszeviani<br style="font-family:trebuchet ms,sans-serif"></span><span style="font-family:trebuchet ms,sans-serif">homepage: <a href="http://www.leg.ufpr.br/%7Ewalmes" target="_blank">http://www.leg.ufpr.br/~walmes</a></span><br style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">linux user number: 531218</span><br style="font-family:trebuchet ms,sans-serif"><span style="font-family:trebuchet ms,sans-serif">==========================================================================</span><br>
</div></div><br></div></div>