Reputation: 868
I'm running a glmmTMB model with various truncated count distributions (truncated_poisson
, truncated_compois
, truncated_nbinom1
, truncated_nbinom2
). When I predict from the model, the values seem to be lower than expected, as if the prediction is not accounting for the truncation. Where am I going wrong? A toy example is provided, showing that predicted values are lower than observed means.
Any advice would be appreciated. Extra points if the advice can extend to the other truncated count distributions (see above) and if it shows how to correctly get the 95% confidence band around the estimated values in these cases.
library(dplyr)
library( extraDistr)
library(glmmTMB)
set.seed(1)
df <- data.frame(Group = rep(c("a", "b"), each = 20), N = rtpois(40, 1, a = 0), ran = "a") %>%
mutate(N = ifelse(N == 0, 1, N))
m <- glmmTMB(N ~ Group + (1|ran), data = df, family = "truncated_poisson")
df %>% group_by(Group) %>% summarize(mean(N))
predict(m, newdata = data.frame(Group = c("a", "b"), ran = NA), type = "response")
Upvotes: 1
Views: 574
Reputation: 7832
Maybe to add another alternative to get the results shown in Ben's answer, you could use the ggeffects package to get predictions for those models. However, in order to correctly detect the zero-inflation, you have to specify ziformula = ~1
, then type = "zero_inflated"
will give you predictions on the response scale, taking the zero-inflation part into account.
library(ggeffects)
library(extraDistr)
library(glmmTMB)
set.seed(1)
df <- data.frame(Group = rep(c("a", "b"), each = 20), Np = rtpois(40, 1, a = 0))
## clunky trunc nbinom generator
tnb <- rep(0, 40)
z <- (tnb==0)
while (any(z)) {
tnb[z] <- rnbinom(sum(z), mu = 1, size = 1)
z <- (tnb == 0)
}
df$Nnb <- tnb
m1 <- glmmTMB(Np ~ Group, ziformula = ~1, data = df, family = "truncated_poisson")
m2 <- update(m1, Nnb ~ ., family = truncated_nbinom2)
ggpredict(m1, "Group", type = "zero_inflated")
#> # Predicted counts of Np
#>
#> Group | Predicted | 95% CI
#> --------------------------------
#> a | 1.75 | [0.77, 2.73]
#> b | 1.45 | [0.78, 2.12]
ggpredict(m2, "Group", type = "zero_inflated")
#> # Predicted counts of Nnb
#>
#> Group | Predicted | 95% CI
#> --------------------------------
#> a | 1.80 | [0.83, 2.77]
#> b | 2.35 | [0.91, 3.79]
Created on 2023-11-03 with reprex v2.0.2
Upvotes: 1
Reputation: 226861
I think the main issue is probably that you're using a slightly older version of glmmTMB
(< 1.1.5, where a bug was fixed, see e.g. e.g. https://github.com/glmmTMB/glmmTMB/issues/860).
streamlined slightly (we don't need to include a random effect for this example), and adding a truncated nbinom2.
library(dplyr)
library(extraDistr)
library(glmmTMB)
set.seed(1)
df <- data.frame(Group = rep(c("a", "b"), each = 20),
Np = rtpois(40, 1, a = 0))
## clunky trunc nbinom generator
tnb <- rep(0, 40)
z <- (tnb==0)
while(any(z)) {
tnb[z] <- rnbinom(sum(z), mu = 1, size = 1)
z <- (tnb==0)
}
df$Nnb <- tnb
## summarize
df %>% group_by(Group) %>% summarize(across(starts_with("N"), mean))
## Group Np Nnb
## 1 a 1.75 1.8
## 2 b 1.45 2.35
m1 <- glmmTMB(Np ~ Group, data = df, family = "truncated_poisson")
m2 <- update(m1, Nnb ~ ., family = truncated_nbinom2)
Predicting with se.fit = TRUE
will give you standard errors for the predictions, from which you can compute confidence intervals (assuming Normality/Wald intervals/blah blah blah ...) ...
pfun <- function(m, level = 0.95) {
pp <- predict(m, newdata = data.frame(Group = c("a", "b")),
type = "response",
se.fit = TRUE)
list(est = unname(pp$fit),
lwr = unname(pp$fit + qnorm((1-level)/2)*pp$se.fit),
upr = unname(pp$fit + qnorm((1+level)/2)*pp$se.fit))
}
pfun(m1)
pfun(m2)
Upvotes: 2