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
[email protected]
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:[email protected]
[email protected]
Lattes: http://lattes.cnpq.br/1360403201088680
======================================================================
_______________________________________________
R-br mailing list
[email protected]
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.
