Plot poisson mixed models with ggplot2

I try to make a plot for standard purposes with zero inflated model and zero inflated mixed model using ggplot2 without success. For this, I try:

#Packages
library(pscl)
library(glmmTMB)
library(ggplot2)
library(gridExtra)


# 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 thal 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)
#-

Doesn't work of sure. Is possible to make this using ggplot2? Thanks in advance

Solution

I'm not sure, but it looks to me that you're looking for marginal effects. You can do this with the ggeffects-package. Here are two examples, using your simulated data, that create a ggplot-object, one with and one w/o raw data.

library(glmmTMB)
library(ggeffects)

mZIPmix<- glmmTMB(y~poly(time,2)+sex+(1|id), data=DF, ziformula=~1,family=poisson)

# compute marginal effects and create a plot.
# the tag "[all]" is useful for polynomial terms, to produce smoother plots
ggpredict(mZIPmix, c("time [all]", "sex")) %>% plot(rawdata = TRUE, jitter = .01)

ggpredict(mZIPmix, c("time [all]", "sex")) %>% plot(rawdata = FALSE)

^{Created on 2019-05-16 by the reprex package (v0.2.1)}

Note that sex only has an "additive" effect. Maybe you want to model an intercation between time and sex?

mZIPmix<- glmmTMB(y~poly(time,2)*sex+(1|id), data=DF, ziformula=~1,family=poisson)

ggpredict(mZIPmix, c("time [all]", "sex")) %>% plot(rawdata = TRUE, jitter = .01)

ggpredict(mZIPmix, c("time [all]", "sex")) %>% plot()

^{Created on 2019-05-16 by the reprex package (v0.2.1)}