<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body text="#000000" bgcolor="#FFFFFF">
Walmes,<br>
<br>
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?<br>
<br>
Obrigado,<br>
<br>
Segue CRM:<br>
#------------------------------------------------------------------<br>
# Definições da sessão.<br>
<br>
rm(list=ls())<br>
require(pscl)<br>
require(VGAM)<br>
require(multcomp)<br>
require(lattice)<br>
require(latticeExtra)<br>
<br>
#------------------------------------------------------------------<br>
# Dados artificiais.<br>
<br>
da <- expand.grid(tempo=rep(1:10), trat=gl(3,10))<br>
<br>
## Simulação de 3 distribuições inflacionadas de<br>
## Poisson usando pacote VGAM<br>
set.seed(54321)<br>
lambda = 10; phi = 0.1<br>
y1 <- rzipois(100, lambda, phi)<br>
lambda = 4; phi = 0.3<br>
y2 <- rzipois(100, lambda, phi)<br>
lambda = 10; phi = 0.1<br>
y3 <- rzipois(100, lambda, phi)<br>
da$y <- c(y1,y2,y3)<br>
str(da)<br>
xyplot(y~tempo|trat, data=da, jitter.x=TRUE)<br>
is.factor(da$trat)<br>
<br>
#------------------------------------------------------------------<br>
# Modelo completo.<br>
compl.mod <- zeroinfl(y~trat+tempo|trat, data=da)<br>
<br>
## Modelo nulo<br>
null.mod <- update(compl.mod, . ~ 1)<br>
<br>
## diferença em número de parâmetros<br>
## (em dimensão dos espaços dos modelos)<br>
df <- length(coef(compl.mod))-length(coef(null.mod))<br>
<br>
D <- 2*abs(diff(c(logLik(compl.mod), logLik(null.mod))))<br>
pchisq(D, df=df, lower.tail=FALSE)<br>
<br>
#------------------------------------------------------------------<br>
# Contraste de modelos<br>
<br>
### Observando as médias<br>
sort(tapply(da$y,da$trat,mean, na.rm = T))<br>
<br>
# Juntando T1 e T3<br>
trat1<-da$trat<br>
levels(trat1)<br>
levels(trat1)[1]<-"T1_T3"<br>
levels(trat1)[3]<-"T1_T3"<br>
levels(trat1)<br>
reduc.mod1<-zeroinfl(y~trat1+tempo|trat1, data=da)<br>
<br>
# Comparando ao modelo completo<br>
<br>
df <- length(coef(compl.mod))-length(coef(reduc.mod1))<br>
D <- 2*abs(diff(c(logLik(compl.mod), logLik(reduc.mod1))))<br>
pchisq(D, df=df, lower.tail=FALSE)<br>
<br>
## São iguais, tenho por tanto uma curva para T2 e outra para T1 e
T3<br>
<br>
<br>
#------------------------------------------------------------------<br>
# Predição do modelo considerando as duas porções.<br>
<br>
X <- model.matrix(~trat1+tempo, da)<br>
i <- grep("^count\\_", names(coef(reduc.mod1)))<br>
eta <- X%*%coef(reduc.mod1)[i]<br>
da$y.pois <- exp(eta)<br>
<br>
X <- model.matrix(~trat1, da)<br>
i <- grep("^zero\\_", names(coef(reduc.mod1)))<br>
eta <- X%*%coef(reduc.mod1)[i]<br>
da$y.zero <- exp(eta)/(1+exp(eta))<br>
<br>
xyplot(y~tempo|trat1, data=da, jitter.x=TRUE)+<br>
as.layer(xyplot(y.pois~tempo|trat1, data=da, type="l"))+<br>
as.layer(xyplot(y.zero~tempo|trat1, data=da,<br>
type="l", lty=2, lwd=2))+<br>
layer(panel.abline(h=1, lty=2))<br>
<br>
# contínua: média da contagem ~ Poisson.<br>
# tracejada: probabilidade de um zero não Poisson.<br>
# abline: linha no 1, referência.<br>
<br>
#------------------------------------------------------------------<br>
<br>
##Gráfico 2 <br>
<br>
x<-c(1,2,3,4,5,6,7,8,9,10)<br>
<br>
dados2<- with(da, tapply(y.pois, list(trat1, tempo), mean, na.rm
= T)) ### Medias estimadas de contagem<br>
dados3<- with(da, tapply(y.zero, list(trat1, tempo), mean, na.rm
= T)) ### Medias estimadas binomial<br>
dados4a<-da[da[,3]!=0,]### Somente dados de contagem<br>
trat.p<-dados4a$trat<br>
levels(trat.p)<br>
levels(trat.p)[1]<-"T1eT3"<br>
levels(trat.p)[3]<-"T1eT3"<br>
levels(trat.p)<br>
dados4<- with(dados4a, tapply(y, list(trat.p, tempo), mean, na.rm
= T))<br>
binom<-rep(1,length(dados4a[,1])) ### Transformando os dados
observados em 1 e 0<br>
dados4b<-da[da[,3]==0,]<br>
da$binom<-ifelse(da[,3]>0,1,0)<br>
dados5<- with(da, tapply(binom, list(trat1, tempo), mean, na.rm =
T))<br>
<br>
#Gráfico para
binomial
"<br>
plot(x, dados3[1,], ylim=c(0,1), xlim=c(0,10), ylab="Probabilidade
de ocorrência", xlab="tempo",lty=1,type="l")<br>
lines(x, dados3[2,],lty=2)<br>
points(x, dados5[1,],pch=18)<br>
points(x, dados5[2,],pch=22) <br>
#<br>
#-------------------------------------------------------------------------------<br>
<br>
<br>
<br>
<div class="moz-cite-prefix">Em 05/11/2013 09:40, walmes . escreveu:<br>
</div>
<blockquote
cite="mid:CAFU=EkbZw6ZRrjGTJU0q1VpYb_LZ5fvdp15vcmsTpSs3Z-P+tg@mail.gmail.com"
type="cite">
<div dir="ltr">
<div class="gmail_default" style="font-family:trebuchet
ms,sans-serif">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.<br>
<br>
<span style="font-family:courier new,monospace">>
#------------------------------------------------------------------<br>
> # Definições da sessão.<br>
> <br>
> rm(list=ls())<br>
> require(pscl)<br>
> require(VGAM)<br>
> require(multcomp)<br>
> require(lattice)<br>
> require(latticeExtra)<br>
> <br>
>
#------------------------------------------------------------------<br>
> # Dados artificiais.<br>
> <br>
> da <- expand.grid(tempo=rep(1:10), trat=gl(3,10))<br>
> xtabs(~trat, da)<br>
trat<br>
1 2 3 <br>
100 100 100 <br>
> head(da)<br>
tempo trat<br>
1 1 1<br>
2 2 1<br>
3 3 1<br>
4 4 1<br>
5 5 1<br>
6 6 1<br>
> <br>
> ## Simulação de 3 distribuições inflacionadas de<br>
> ## Poisson usando pacote VGAM<br>
> set.seed(54321)<br>
> lambda = 10; phi = 0.1<br>
> y1 <- rzipois(100, lambda, phi)<br>
> lambda = 4; phi = 0.3<br>
> y2 <- rzipois(100, lambda, phi)<br>
> lambda = 8; phi = 0.5<br>
> y3 <- rzipois(100, lambda, phi)<br>
> da$y <- c(y1,y2,y3)<br>
> <br>
> str(da)<br>
'data.frame': 300 obs. of 3 variables:<br>
$ tempo: int 1 2 3 4 5 6 7 8 9 10 ...<br>
$ trat : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1
1 ...<br>
$ y : num 9 7 13 7 12 5 14 11 17 10 ...<br>
- attr(*, "out.attrs")=List of 2<br>
..$ dim : Named int 10 30<br>
.. ..- attr(*, "names")= chr "tempo" "trat"<br>
..$ dimnames:List of 2<br>
.. ..$ tempo: chr "tempo= 1" "tempo= 2" "tempo= 3"
"tempo= 4" ...<br>
.. ..$ trat : chr "trat=1" "trat=1" "trat=1" "trat=1" ...<br>
> xyplot(y~tempo|trat, data=da, jitter.x=TRUE)<br>
> <br>
> is.factor(da$trat)<br>
[1] TRUE<br>
> <br>
>
#------------------------------------------------------------------<br>
> # Modelo completo.<br>
> compl.mod <- zeroinfl(y~trat+tempo|trat, data=da)<br>
> summary(compl.mod)<br>
<br>
Call:<br>
zeroinfl(formula = y ~ trat + tempo | trat, data = da)<br>
<br>
Pearson residuals:<br>
Min 1Q Median 3Q Max <br>
-2.59962 -0.97565 -0.03836 0.68546 2.59079 <br>
<br>
Count model coefficients (poisson with log link):<br>
Estimate Std. Error z value Pr(>|z|) <br>
(Intercept) 2.242577 0.056657 39.581 < 2e-16 ***<br>
trat2 -1.051578 0.072962 -14.413 < 2e-16 ***<br>
trat3 -0.251806 0.057472 -4.381 1.18e-05 ***<br>
tempo 0.015857 0.008327 1.904 0.0569 . <br>
<br>
Zero-inflation model coefficients (binomial with logit
link):<br>
Estimate Std. Error z value Pr(>|z|) <br>
(Intercept) -2.9452 0.4592 -6.414 1.41e-10 ***<br>
trat2 1.9999 0.5159 3.877 0.000106 ***<br>
trat3 2.7437 0.5013 5.474 4.41e-08 ***<br>
---<br>
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' '
1 <br>
<br>
Number of iterations in BFGS optimization: 12 <br>
Log-likelihood: -676.3 on 7 Df<br>
> length(coef(compl.mod))<br>
[1] 7<br>
> <br>
> ## Modelo nulo<br>
> null.mod <- update(compl.mod, . ~ 1)<br>
> summary(null.mod)<br>
<br>
Call:<br>
zeroinfl(formula = y ~ 1, data = da)<br>
<br>
Pearson residuals:<br>
Min 1Q Median 3Q Max <br>
-1.3584 -1.3584 -0.1426 0.8299 3.0182 <br>
<br>
Count model coefficients (poisson with log link):<br>
Estimate Std. Error z value Pr(>|z|) <br>
(Intercept) 2.03004 0.02447 82.95 <2e-16 ***<br>
<br>
Zero-inflation model coefficients (binomial with logit
link):<br>
Estimate Std. Error z value Pr(>|z|) <br>
(Intercept) -1.0135 0.1307 -7.752 9.05e-15 ***<br>
---<br>
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' '
1 <br>
<br>
Number of iterations in BFGS optimization: 9 <br>
Log-likelihood: -831.7 on 2 Df<br>
> length(coef(null.mod))<br>
[1] 2<br>
> <br>
> ## diferença em número de parâmetros<br>
> ## (em dimensão dos espaços dos modelos)<br>
> df <- length(coef(compl.mod))-length(coef(null.mod))<br>
> <br>
> ## isso da número negativo para Deviance!!, montado
errado<br>
> D <- -2*(logLik(compl.mod)-logLik(null.mod))<br>
> unclass(D)<br>
[1] -310.8615<br>
attr(,"df")<br>
[1] 7<br>
> pchisq(D, df=df, lower.tail=FALSE)<br>
'log Lik.' 1 (df=7)<br>
> <br>
> ## assim você evita preocupação com sinais<br>
> D <- 2*abs(diff(c(logLik(compl.mod),
logLik(null.mod))))<br>
> pchisq(D, df=df, lower.tail=FALSE)<br>
[1] 4.625179e-65<br>
> <br>
>
#-----------------------------------------------------------------------------<br>
> </span><br>
</div>
<div class="gmail_extra"><br>
<div class="gmail_default" style="font-family:trebuchet
ms,sans-serif">À disposição.<br>
Walmes.</div>
<br clear="all">
<div>
<div dir="ltr"><span style="font-family:trebuchet
ms,sans-serif">==========================================================================</span><br
style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">Walmes
Marques Zeviani</span><br style="font-family:trebuchet
ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">LEG
(Laboratório de Estatística e Geoinformação, 25.450418
S, 49.231759 W)</span><br style="font-family:trebuchet
ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">Departamento
de Estatística - Universidade Federal do Paraná</span><br
style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">fone:
(+55) 41 3361 3573</span><br
style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">skype:
walmeszeviani<br style="font-family:trebuchet
ms,sans-serif">
</span><span style="font-family:trebuchet ms,sans-serif">homepage:
<a moz-do-not-send="true"
href="http://www.leg.ufpr.br/%7Ewalmes"
target="_blank">http://www.leg.ufpr.br/~walmes</a></span><br
style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">linux
user number: 531218</span><br
style="font-family:trebuchet ms,sans-serif">
<span style="font-family:trebuchet ms,sans-serif">==========================================================================</span><br>
</div>
</div>
<br>
</div>
</div>
<br>
<fieldset class="mimeAttachmentHeader"></fieldset>
<br>
<pre wrap="">_______________________________________________
R-br mailing list
<a class="moz-txt-link-abbreviated" href="mailto:R-br@listas.c3sl.ufpr.br">R-br@listas.c3sl.ufpr.br</a>
<a class="moz-txt-link-freetext" href="https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br">https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br</a>
Leia o guia de postagem (<a class="moz-txt-link-freetext" href="http://www.leg.ufpr.br/r-br-guia">http://www.leg.ufpr.br/r-br-guia</a>) e forneça código mínimo reproduzível.</pre>
</blockquote>
<br>
<pre class="moz-signature" cols="72">--
======================================================================
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)
<a class="moz-txt-link-abbreviated" href="mailto:e-mails:alexandresantosbr@yahoo.com.br">e-mails:alexandresantosbr@yahoo.com.br</a>
<a class="moz-txt-link-abbreviated" href="mailto:alexandre.santos@cas.ifmt.edu.br">alexandre.santos@cas.ifmt.edu.br</a>
Lattes: <a class="moz-txt-link-freetext" href="http://lattes.cnpq.br/1360403201088680">http://lattes.cnpq.br/1360403201088680</a>
======================================================================
</pre>
</body>
</html>