如何预测装有 glmmTMB 函数的模型的概率

Question

我使用了下面的代码，但我得到了 po.gp 的错误，因为它是一个装有 glmmTMB 包的广义泊松模型。我得到错误：

#MODEL COMPARISONS
library(pscl)
po.p <- predprob(model.p) %>% colMeans
po.nb <- predprob(model.nb) %>% colMeans
po.gp <- predprob(model.gp) %>% colMeans

df <- data.frame(x = 0:max(na.omit(bo$Grav)), Poisson = po.p, 
                 NegBin = po.nb, Gen_Poisson = po.gp)

obs <- table(bo$Grav) %>% prop.table() %>% data.frame #Observed
names(obs) <- c("x", 'Observed')

p1 <- predict(linear) %>% round() %>% table %>% prop.table %>% data.frame #for OLS
names(p1) <- c('x', 'OLS')

tmp <- merge(p1, obs, by = 'x', all = T)
tmp$x <- as.numeric(as.character(tmp$x))

comb <- merge(tmp, df, by = 'x', all = T)
comb[is.na(comb)] <- 0

comb2 <- comb[1:11, ] #just for the first 11 results, including zero

mm <- melt(comb2, id.vars = 'x', value.name = 'prob', variable.name = 'Model')
mm <- filter(mm, Model != "OLS") #can include the linear model too if you want
#the SAS note does not, so I am not including it

ggplot(mm, aes(x = x, y = prob, group = Model, col = Model)) +
  geom_line(aes(lty = Model), lwd = 1) +
  theme_bw() +
  labs(x = "Number of pregnancies", y = 'Probability',
       title = "Models for number of pregnancies") +
  scale_color_manual(values = c('black', 'blue', 'red', 'green')) +
  scale_linetype_manual(values = c('solid', 'dotted', 'dotted', 'dotted')) +
  theme(legend.position=c(.2, .80), axis.title.y = element_text(angle = 0))

我试图从装有 glmmTMB 包的广义泊松模型中预测和绘制概率，以便将其与其他计数模型进行比较，但函数 predprob 不适用于 glmmTMB 包，但我正在寻找合适的函数来使用。

Answer 1

这不是

predprob

的完全替代品，但您可以使用

VGAM::dgenpois1()

生成相应的概率。 （要检查

glmmTMB

和

VGAM

的参数化之间的对应关系，请参阅

?glmmTMB::family_glmmTMB

，其中指出方差-均值关系为

V = mu*exp(eta)

；这对应于

VGAM

的

dgenpois1()

函数。

library(VGAM)
library(glmmTMB)
set.seed(101)
z <- rgenpois1(1000, meanpar = exp(2), dispind = 2)
m <- glmmTMB(z ~ 1, data = data.frame(z_nb), family = genpois)


xvec <- 0:23
probs <- dgenpois1(xvec, meanpar = exp(predict(m)[1]),
          dispind = sigma(m), log = FALSE)


par(las = 1)
plot(prop.table(table(z)))
points(xvec, probs, col = 2, pch = 16)