[R-sig-ME] Comparative plot between zero inflated poisson model (ZIP) and zero inflated poisson mixed model

ASANTOS @|ex@ndre@@nto@br @end|ng |rom y@hoo@com@br
Wed Mar 13 23:02:03 CET 2019


Dear Members,

I've like to plot zero inflated poisson model (ZIP) and zero inflated 
poisson mixed model for look the differences in the models adjusted. For 
this I try:

#Packages
library(pscl) #For zero inflated poisson
library(glmmTMB)# For mixed zero inflated poisson
library(ggplot2) #Plot the results
library(gridExtra) ## Put 2 ggplot together


# Artificial data set
set.seed(007)
n <- 100 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 5 # maximum follow-up time
DF <- data.frame(id = rep(seq_len(n), each = K),
                  time = c(replicate(n, c(0, sort(runif(K - 1, 0, 
t_max))))),
                  sex = rep(gl(2, n/2, labels = c("male", "female")), 
each = K))
DF$y <- rnbinom(n * K, size = 2, mu = exp(1.552966))
str(DF)

# Using zero inflated poisson model with pscl package
time2<-(DF$time)^2
mZIP <- zeroinfl(y~time+time2+sex|time+sex, data=DF)
summary(mZIP)

#If I imagine that all coefficients are significant

# Y estimated
pred.data1 = data.frame(
time<-DF$time,
time2<-(DF$time)^2,
sex<-DF$sex)
pred.data1$y = predict(mZIP, newdata=pred.data1, type="response")


# Now using zero inflated poisson mixed model with glmmTMB package
mZIPmix<- glmmTMB(y~time+time2+sex+(1|id),
data=DF, ziformula=~1,family=poisson)
summary(mZIPmix)
#

# new Y estimated
pred.data2 = data.frame(
time<-DF$time,
time2<-(DF$time)^2,
sex<-DF$sex,
id<-DF$id)
pred.data2$y = predict(mZIPmix, newdata=pred.data2, type="response")


# Plot zero inflated poisson model and mixed poisson model
par(mfrow=c(1,2))
plot1<-ggplot(DF, aes(time, y, colour=sex)) +
   labs(title="Zero inflated model") +
   geom_point() +
   geom_line(data=pred.data1) +
   stat_smooth(method="glm", family=poisson(link="log"), formula = 
y~poly(x,2),fullrange=TRUE)

plot2<-ggplot(DF, aes(time, y, colour=sex)) +
   labs(title="Zero inflated mixed model") +
   geom_point() +
   geom_line(data=pred.data2) +
   stat_smooth(method="glm", family=poisson(link="log"), formula = 
y~poly(x,2),fullrange=TRUE)## here a don't find any method to mixed glm
grid.arrange(plot1, plot2, ncol=2)
#-

But doesn't work of sure. Please, how I could make this two plots for 
direct comparing purposes?

Thanks in advanced,

Alexandre

-- 
======================================================================
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 99686-6970 (VIVO) (+55) 65 3221-2674 (FIXO)

         alexandre.santos using cas.ifmt.edu.br
Lattes: http://lattes.cnpq.br/1360403201088680
OrcID: orcid.org/0000-0001-8232-6722
Researchgate: www.researchgate.net/profile/Alexandre_Santos10
LinkedIn: br.linkedin.com/in/alexandre-dos-santos-87961635
Mendeley:www.mendeley.com/profiles/alexandre-dos-santos6/



More information about the R-sig-mixed-models mailing list