5 Nov
2013
5 Nov
'13
13:02
Muito bem Walmes, o sinal das verossimilhanças!!!! O argumento
2*abs(diff(c(logLik(compl.mod), logLik(null.mod)))) Resolveu!
Obrigado,
Abraço,
--
======================================================================
Alexandre dos Santos
Proteção Florestal
IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso
Campus Cáceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
Cáceres - MT CEP: 78.200-000
Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO)
e-mails:alexandresantosbr@yahoo.com.br
alexandre.santos@cas.ifmt.edu.br
Lattes: http://lattes.cnpq.br/1360403201088680
======================================================================
Em 05/11/2013 09:40, walmes . escreveu:
> 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
> <http://www.leg.ufpr.br/%7Ewalmes>
> linux user number: 531218
> ==========================================================================
>
>
>
> _______________________________________________
> R-br mailing list
> R-br@listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
--
======================================================================
Alexandre dos Santos
Proteção Florestal
IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso
Campus Cáceres
Caixa Postal 244
Avenida dos Ramires, s/n
Bairro: Distrito Industrial
Cáceres - MT CEP: 78.200-000
Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO)
e-mails:alexandresantosbr@yahoo.com.br
alexandre.santos@cas.ifmt.edu.br
Lattes: http://lattes.cnpq.br/1360403201088680
======================================================================