Contrastes em modelos de regressão de poisson inflacionados por zeros
Caros membros, Gostaria de fazer a junção de termos qualitativos não significativos, através de análises de contraste de modelos, para verificar a semelhança entre tratamentos, no entanto, usando modelos de regressão de poisson inflacionados por zeros, aproveitando um modelo do Walmes fiz: #------------------------------------------------------------------ # Definições da sessão. rm(list=ls()) require(pscl) require(multcomp) require(lattice) require(latticeExtra) #------------------------------------------------------------------ # Dados artificiaiscom 3 tratamentos e 10 tempos da <- expand.grid(trat=gl(3,1), tempo=1:10) X <- model.matrix(~trat+tempo, da); ncol(X) betas <- c(0.1,0.1,0.3,0.5) eta <- X%*%betas y1 <- rpois(da$trat, lambda=exp(eta)) y2 <- rbinom(y1, size=1, prob=0.7) da$y <- y1*y2 str(da) xyplot(y~tempo|trat, data=da, jitter.x=TRUE) #------------------------------------------------------------------ # Ajuste do modelo. compl.mod <- zeroinfl(y~trat+tempo|trat, data=da) summary(compl.mod) #------------------------------------------------------------------ # Contrate de modelos de poisson inflacionado com zeros sort(tapply(da$y,da$trat,mean))#Ordena as médias trat2<-da$trat levels(trat2) levels(trat2)[1]<-"T1eT3" levels(trat2)[3]<-"T1eT3" levels(trat2) reduc.mod<-zeroinfl(y~trat2+tempo|trat, data=da) # Comparando o modelo completo e o ajunte com junção dos tratamentos 1 e 3 pchisq(-2*(logLik(compl.mod)-logLik(reduc.mod)),df=3,lower.tail=FALSE) # Na minha simulação: #[1] 1 #attr(,"df") #[1] 7 #attr(,"class") #[1] "logLik" Então pergunto: T1 e T3 são iguais e podem ser representados por um único modelo, isto esta correto? Ou existe alguma outra abordagem? Obrigado, -- ====================================================================== 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 ======================================================================
Eu não corri o CMR mas a princípio está correto. Só o grau de liberdade do teste de razão de verossimilhanças que deveria ser 1 e não 3, certo? Porque a diferença entre os espaços dos modelos é de 1 parâmetro. À 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 ==========================================================================
Caros membros, Tentei criar dois tratamentos diferentes simulados através de rzipois() do pacote VGAM, mas estranhamente o modelo nulo não é diferente do modelo completo, o que pode estar acontecendo? Segue CRM: #------------------------------------------------------------------ # Definições da sessão. rm(list=ls()) require(pscl) require(VGAM) require(multcomp) require(lattice) require(latticeExtra) #------------------------------------------------------------------ # Dados artificiais. trat <- gl(3,100) ## 3 tratamentos tempo<- rep(sort(rep(1:10,10)),3) ## 10 tempos #Simulação de 3 distribuições inflacionadas de Poisson usando pacote VGAM 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) y <- c(y1,y2,y3) da<-as.data.frame(cbind(trat,tempo,y)) str(da) xyplot(y~tempo|trat, data=da, jitter.x=TRUE) #------------------------------------------------------------------ # Ajuste do modelo. compl.mod <- zeroinfl(y~trat+tempo|trat, data=da) summary(compl.mod) ## Testando nulo e o completo compl.mod <- zeroinfl(y~trat+tempo|trat, data=da) null.mod <- update(compl.mod, . ~ 1) pchisq(-2*(logLik(compl.mod)-logLik(null.mod)),df=2,lower.tail=FALSE) ## A não significância indica que os tratamentos são iguais, por que? Obrigado, -- ====================================================================== 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 03/11/2013 23:46, walmes . escreveu:
Eu não corri o CMR mas a princípio está correto. Só o grau de liberdade do teste de razão de verossimilhanças que deveria ser 1 e não 3, certo? Porque a diferença entre os espaços dos modelos é de 1 parâmetro.
À 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.
Caros membros, Tentei criar três tratamentos diferentes simulados através de rzipois() do pacote VGAM, mas estranhamente o modelo nulo não é diferente do modelo completo, o que pode estar acontecendo? Segue CRM: #------------------------------------------------------------------ # Definições da sessão. rm(list=ls()) require(pscl) require(VGAM) require(multcomp) require(lattice) require(latticeExtra) #------------------------------------------------------------------ # Dados artificiais. trat <- gl(3,100) ## 3 tratamentos tempo<- rep(sort(rep(1:10,10)),3) ## 10 tempos #Simulação de 3 distribuições inflacionadas de Poisson usando pacote VGAM 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) y <- c(y1,y2,y3) da<-as.data.frame(cbind(trat,tempo,y)) str(da) xyplot(y~tempo|trat, data=da, jitter.x=TRUE) #------------------------------------------------------------------ # Ajuste do modelo. compl.mod <- zeroinfl(y~trat+tempo|trat, data=da) summary(compl.mod) ## Testando nulo e o completo compl.mod <- zeroinfl(y~trat+tempo|trat, data=da) null.mod <- update(compl.mod, . ~ 1) pchisq(-2*(logLik(compl.mod)-logLik(null.mod)),df=2,lower.tail=FALSE) ## A não significância indica que os tratamentos são iguais, por que? Obrigado, -- ====================================================================== 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 04/11/2013 16:15, ASANTOS escreveu:
Em 03/11/2013 23:46, walmes . escreveu:
Eu não corri o CMR mas a princípio está correto. Só o grau de liberdade do teste de razão de verossimilhanças que deveria ser 1 e não 3, certo? Porque a diferença entre os espaços dos modelos é de 1 parâmetro.
À 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.
_______________________________________________ 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 ======================================================================
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
==========================================================================
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
======================================================================
Walmes,
Achei que o gráfico do lattice não permite visualizar bem a
parte binomial do modelo, então fiz um gráfico separado, porém percebi
que a parte binomial do modelo forma uma reta, mesmo fazendo o
curve(exp(a+bx)/(1+exp(a+bx))) com os coeficientes não ficou correto e
pergunto não deveria ser a parte binomial um modelo parecido com a forma
de um S com os valores médios observados em torno, conforme dados
artificiais do CRM abaixo?
Obrigado,
Segue CRM:
#------------------------------------------------------------------
# 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çã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 = 10; phi = 0.1
y3 <- rzipois(100, lambda, phi)
da$y <- c(y1,y2,y3)
str(da)
xyplot(y~tempo|trat, data=da, jitter.x=TRUE)
is.factor(da$trat)
#------------------------------------------------------------------
# Modelo completo.
compl.mod <- zeroinfl(y~trat+tempo|trat, data=da)
## Modelo nulo
null.mod <- update(compl.mod, . ~ 1)
## 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))
D <- 2*abs(diff(c(logLik(compl.mod), logLik(null.mod))))
pchisq(D, df=df, lower.tail=FALSE)
#------------------------------------------------------------------
# Contraste de modelos
### Observando as médias
sort(tapply(da$y,da$trat,mean, na.rm = T))
# Juntando T1 e T3
trat1<-da$trat
levels(trat1)
levels(trat1)[1]<-"T1_T3"
levels(trat1)[3]<-"T1_T3"
levels(trat1)
reduc.mod1<-zeroinfl(y~trat1+tempo|trat1, data=da)
# Comparando ao modelo completo
df <- length(coef(compl.mod))-length(coef(reduc.mod1))
D <- 2*abs(diff(c(logLik(compl.mod), logLik(reduc.mod1))))
pchisq(D, df=df, lower.tail=FALSE)
## São iguais, tenho por tanto uma curva para T2 e outra para T1 e T3
#------------------------------------------------------------------
# Predição do modelo considerando as duas porções.
X <- model.matrix(~trat1+tempo, da)
i <- grep("^count\\_", names(coef(reduc.mod1)))
eta <- X%*%coef(reduc.mod1)[i]
da$y.pois <- exp(eta)
X <- model.matrix(~trat1, da)
i <- grep("^zero\\_", names(coef(reduc.mod1)))
eta <- X%*%coef(reduc.mod1)[i]
da$y.zero <- exp(eta)/(1+exp(eta))
xyplot(y~tempo|trat1, data=da, jitter.x=TRUE)+
as.layer(xyplot(y.pois~tempo|trat1, data=da, type="l"))+
as.layer(xyplot(y.zero~tempo|trat1, data=da,
type="l", lty=2, lwd=2))+
layer(panel.abline(h=1, lty=2))
# contínua: média da contagem ~ Poisson.
# tracejada: probabilidade de um zero não Poisson.
# abline: linha no 1, referência.
#------------------------------------------------------------------
##Gráfico 2
x<-c(1,2,3,4,5,6,7,8,9,10)
dados2<- with(da, tapply(y.pois, list(trat1, tempo), mean, na.rm = T))
### Medias estimadas de contagem
dados3<- with(da, tapply(y.zero, list(trat1, tempo), mean, na.rm = T))
### Medias estimadas binomial
dados4a<-da[da[,3]!=0,]### Somente dados de contagem
trat.p<-dados4a$trat
levels(trat.p)
levels(trat.p)[1]<-"T1eT3"
levels(trat.p)[3]<-"T1eT3"
levels(trat.p)
dados4<- with(dados4a, tapply(y, list(trat.p, tempo), mean, na.rm = T))
binom<-rep(1,length(dados4a[,1])) ### Transformando os dados observados
em 1 e 0
dados4b<-da[da[,3]==0,]
da$binom<-ifelse(da[,3]>0,1,0)
dados5<- with(da, tapply(binom, list(trat1, tempo), mean, na.rm = T))
#Gráfico para binomial "
plot(x, dados3[1,], ylim=c(0,1), xlim=c(0,10), ylab="Probabilidade de
ocorrência", xlab="tempo",lty=1,type="l")
lines(x, dados3[2,],lty=2)
points(x, dados5[1,],pch=18)
points(x, dados5[2,],pch=22)
#
#-------------------------------------------------------------------------------
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
======================================================================
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 ==========================================================================
Obrigado Walmes, Esqueci de y~trat+tempo|trat+tempo, desapercebidamente não reparei que faltava o coeficiente do tempo na parte binomial, Redobrados agradecimentos, -- ====================================================================== 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 08/11/2013 17:32, walmes . escreveu:
xyplot(p.zero~tempo|trat, data=pred, type="l", lty=2, lwd=2)
Alexandre, A menos que a versão que você está usando do pacote pscl mudou a interface de usuário, se você quer os mesmos coeficientes para parte binomial e a parte inflacionada com zeros, não carece de repetir os mesmos depois da barra vertical. 2013/11/8 ASANTOS <alexandresantosbr@yahoo.com.br>
Obrigado Walmes,
Esqueci de y~trat+tempo|trat+tempo, desapercebidamente não reparei que faltava o coeficiente do tempo na parte binomial,
Redobrados agradecimentos,
-- ====================================================================== 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 08/11/2013 17:32, walmes . escreveu:
xyplot(p.zero~tempo|trat, data=pred, type="l", lty=2, lwd=2)
_______________________________________________ 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.
participantes (3)
-
ASANTOS -
Cesar Rabak -
walmes .