<div dir="ltr"><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">Você tá fazendo o teste de razão de verossimilhanças errado. Confusão com os sinais das verossimilhanças possivelmente (também me atrapalho com isso). Além do mais, na sua simulação suspeito que tenha esquecido de declarar 'trat' como fator,l mas isso é o de menos. Por outro lado, os graus de liberdade do teste de razão de verossimilhanças estão errados.<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>
> xtabs(~trat, da)<br>trat<br> 1 2 3 <br>100 100 100 <br>> head(da)<br> tempo trat<br>1 1 1<br>2 2 1<br>3 3 1<br>4 4 1<br>5 5 1<br>6 6 1<br>> <br>> ## Simulação de 3 distribuições inflacionadas de<br>
> ## Poisson usando pacote VGAM<br>> set.seed(54321)<br>> lambda = 10; phi = 0.1<br>> y1 <- rzipois(100, lambda, phi)<br>> lambda = 4; phi = 0.3<br>> y2 <- rzipois(100, lambda, phi)<br>> lambda = 8; phi = 0.5<br>
> y3 <- rzipois(100, lambda, phi)<br>> da$y <- c(y1,y2,y3)<br>> <br>> str(da)<br>'data.frame': 300 obs. of 3 variables:<br> $ tempo: int 1 2 3 4 5 6 7 8 9 10 ...<br> $ trat : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...<br>
$ y : num 9 7 13 7 12 5 14 11 17 10 ...<br> - attr(*, "out.attrs")=List of 2<br> ..$ dim : Named int 10 30<br> .. ..- attr(*, "names")= chr "tempo" "trat"<br> ..$ dimnames:List of 2<br>
.. ..$ tempo: chr "tempo= 1" "tempo= 2" "tempo= 3" "tempo= 4" ...<br> .. ..$ trat : chr "trat=1" "trat=1" "trat=1" "trat=1" ...<br>> xyplot(y~tempo|trat, data=da, jitter.x=TRUE)<br>
> <br>> is.factor(da$trat)<br>[1] TRUE<br>> <br>> #------------------------------------------------------------------<br>> # Modelo completo.<br>> compl.mod <- zeroinfl(y~trat+tempo|trat, data=da)<br>
> summary(compl.mod)<br><br>Call:<br>zeroinfl(formula = y ~ trat + tempo | trat, data = da)<br><br>Pearson residuals:<br> Min 1Q Median 3Q Max <br>-2.59962 -0.97565 -0.03836 0.68546 2.59079 <br>
<br>Count model coefficients (poisson with log link):<br> Estimate Std. Error z value Pr(>|z|) <br>(Intercept) 2.242577 0.056657 39.581 < 2e-16 ***<br>trat2 -1.051578 0.072962 -14.413 < 2e-16 ***<br>
trat3 -0.251806 0.057472 -4.381 1.18e-05 ***<br>tempo 0.015857 0.008327 1.904 0.0569 . <br><br>Zero-inflation model coefficients (binomial with logit link):<br> Estimate Std. Error z value Pr(>|z|) <br>
(Intercept) -2.9452 0.4592 -6.414 1.41e-10 ***<br>trat2 1.9999 0.5159 3.877 0.000106 ***<br>trat3 2.7437 0.5013 5.474 4.41e-08 ***<br>---<br>Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 <br>
<br>Number of iterations in BFGS optimization: 12 <br>Log-likelihood: -676.3 on 7 Df<br>> length(coef(compl.mod))<br>[1] 7<br>> <br>> ## Modelo nulo<br>> null.mod <- update(compl.mod, . ~ 1)<br>> summary(null.mod)<br>
<br>Call:<br>zeroinfl(formula = y ~ 1, data = da)<br><br>Pearson residuals:<br> Min 1Q Median 3Q Max <br>-1.3584 -1.3584 -0.1426 0.8299 3.0182 <br><br>Count model coefficients (poisson with log link):<br>
Estimate Std. Error z value Pr(>|z|) <br>(Intercept) 2.03004 0.02447 82.95 <2e-16 ***<br><br>Zero-inflation model coefficients (binomial with logit link):<br> Estimate Std. Error z value Pr(>|z|) <br>
(Intercept) -1.0135 0.1307 -7.752 9.05e-15 ***<br>---<br>Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 <br><br>Number of iterations in BFGS optimization: 9 <br>
Log-likelihood: -831.7 on 2 Df<br>> length(coef(null.mod))<br>[1] 2<br>> <br>> ## diferença em número de parâmetros<br>> ## (em dimensão dos espaços dos modelos)<br>> df <- length(coef(compl.mod))-length(coef(null.mod))<br>
> <br>> ## isso da número negativo para Deviance!!, montado errado<br>> D <- -2*(logLik(compl.mod)-logLik(null.mod))<br>> unclass(D)<br>[1] -310.8615<br>attr(,"df")<br>[1] 7<br>> pchisq(D, df=df, lower.tail=FALSE)<br>
'log Lik.' 1 (df=7)<br>> <br>> ## assim você evita preocupação com sinais<br>> D <- 2*abs(diff(c(logLik(compl.mod), logLik(null.mod))))<br>> pchisq(D, df=df, lower.tail=FALSE)<br>[1] 4.625179e-65<br>
> <br>> #-----------------------------------------------------------------------------<br>> </span><br></div><div class="gmail_extra"><br><div class="gmail_default" style="font-family:trebuchet ms,sans-serif">À disposição.<br>
Walmes.</div><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>