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

[walmes .]

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.

> #------------------------------------------------------------------
> # 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))
> xtabs(~trat, da)
trat
  1   2   3
100 100 100
> head(da)
  tempo trat
1     1    1
2     2    1
3     3    1
4     4    1
5     5    1
6     6    1
>
> ## Simulação de 3 distribuições inflacionadas de
> ## Poisson usando pacote VGAM
> set.seed(54321)
> lambda = 10; phi = 0.1
> y1 <- rzipois(100, lambda, phi)
> lambda = 4; phi = 0.3
> y2 <- rzipois(100, lambda, phi)
> lambda = 8; phi = 0.5
> y3 <- rzipois(100, lambda, phi)
> da$y <- c(y1,y2,y3)
>
> str(da)
'data.frame':    300 obs. of  3 variables:
 $ tempo: int  1 2 3 4 5 6 7 8 9 10 ...
 $ trat : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
 $ y    : num  9 7 13 7 12 5 14 11 17 10 ...
 - attr(*, "out.attrs")=List of 2
  ..$ dim     : Named int  10 30
  .. ..- attr(*, "names")= chr  "tempo" "trat"
  ..$ dimnames:List of 2
  .. ..$ tempo: chr  "tempo= 1" "tempo= 2" "tempo= 3" "tempo= 4" ...
  .. ..$ trat : chr  "trat=1" "trat=1" "trat=1" "trat=1" ...
> xyplot(y~tempo|trat, data=da, jitter.x=TRUE)
>
> is.factor(da$trat)
[1] TRUE
>
> #------------------------------------------------------------------
> # Modelo completo.
> compl.mod <- zeroinfl(y~trat+tempo|trat, data=da)
> summary(compl.mod)

Call:
zeroinfl(formula = y ~ trat + tempo | trat, data = da)

Pearson residuals:
     Min       1Q   Median       3Q      Max
-2.59962 -0.97565 -0.03836  0.68546  2.59079

Count model coefficients (poisson with log link):
             Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.242577   0.056657  39.581  < 2e-16 ***
trat2       -1.051578   0.072962 -14.413  < 2e-16 ***
trat3       -0.251806   0.057472  -4.381 1.18e-05 ***
tempo        0.015857   0.008327   1.904   0.0569 .

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -2.9452     0.4592  -6.414 1.41e-10 ***
trat2         1.9999     0.5159   3.877 0.000106 ***
trat3         2.7437     0.5013   5.474 4.41e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 12
Log-likelihood: -676.3 on 7 Df
> length(coef(compl.mod))
[1] 7
>
> ## Modelo nulo
> null.mod <- update(compl.mod, . ~ 1)
> summary(null.mod)

Call:
zeroinfl(formula = y ~ 1, data = da)

Pearson residuals:
    Min      1Q  Median      3Q     Max
-1.3584 -1.3584 -0.1426  0.8299  3.0182

Count model coefficients (poisson with log link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  2.03004    0.02447   82.95   <2e-16 ***

Zero-inflation model coefficients (binomial with logit link):
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -1.0135     0.1307  -7.752 9.05e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Number of iterations in BFGS optimization: 9
Log-likelihood: -831.7 on 2 Df
> length(coef(null.mod))
[1] 2
>
> ## diferença em número de parâmetros
> ## (em dimensão dos espaços dos modelos)
> df <- length(coef(compl.mod))-length(coef(null.mod))
>
> ## isso da número negativo para Deviance!!, montado errado
> D <- -2*(logLik(compl.mod)-logLik(null.mod))
> unclass(D)
[1] -310.8615
attr(,"df")
[1] 7
> pchisq(D, df=df, lower.tail=FALSE)
'log Lik.' 1 (df=7)
>
> ## assim você evita preocupação com sinais
> D <- 2*abs(diff(c(logLik(compl.mod), logLik(null.mod))))
> pchisq(D, df=df, lower.tail=FALSE)
[1] 4.625179e-65
>
>
#-----------------------------------------------------------------------------
>

À 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/20131105/84476916/attachment.html>