Блог

Главная - Вопросы по программированию - Построить смешанные модели Пуассона с помощью ggplot2

Построить смешанные модели Пуассона с помощью ggplot2

#r #ggplot2 #glm #mixed-models #poisson

#r #ggplot2 #glm #смешанные модели #poisson

Вопрос:

Я безуспешно пытаюсь построить график для стандартных целей с использованием модели с нулевым увеличением и смешанной модели с нулевым увеличением, используя ggplot2. Для этого я пытаюсь:

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

Используя модель Пуассона с нулевым увеличением с пакетом pscl

 time2<-(DF$time)^2
mZIP <- zeroinfl(y~time time2 sex|time sex, data=DF)
summary(mZIP)

Если я представляю, что все коэффициенты значимы

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

Теперь используется смешанная модель Пуассона с нулевым увеличением с пакетом glmmTMB

 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")

Построить модель Пуассона с нулевым увеличением и смешанную модель Пуассона 

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

предварительный

Конечно, не работает. Возможно ли это сделать с помощью ggplot2?
Заранее спасибо

Комментарии:

1. Есть ли проблема на обоих графиках или только на одном? Для GLMM вам понадобятся прогнозы численности населения вместо прогнозов по каждому из них id . Вы можете получить это, установив id в вашем наборе данных прогнозирования значение NA (согласно документации для predict.glmmTMB ).

2. проблема заключается только в построении модели GLMM. Мне нравится график, чтобы посмотреть на различия в скорректированных моделях.

Ответ №1:

Я не уверен, но мне кажется, что вы ищете маргинальные эффекты. Вы можете сделать это с помощью ggeffects-package. Вот два примера, использующие ваши смоделированные данные, которые создают ggplot-объект, один с исходными данными, а другой без них.

 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)

Создано 2019-05-16 пакетом reprex (версия 0.2.1)

Обратите внимание, что 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()

Создано 2019-05-16 пакетом reprex (версия 0.2.1)