[R-br] Contrastes em modelos de regressão de poisson inflacionados por zeros
ASANTOS
alexandresantosbr em yahoo.com.br
Sexta Novembro 8 17:23:21 BRST 2013
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 em 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 em yahoo.com.br
alexandre.santos em cas.ifmt.edu.br
Lattes: http://lattes.cnpq.br/1360403201088680
======================================================================
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20131108/37d02f9a/attachment.html>
Mais detalhes sobre a lista de discussão R-br