使用glmmTMB截断计数分布预测值偏低问题求助

DevOps 小满

2026-8-5

解决glmmTMB截断计数模型预测未考虑截断效应的问题

问题根源

glmmTMB中使用truncated_poisson/truncated_compois/truncated_nbinom1/truncated_nbinom2拟合模型时，默认predict(..., type="response")返回的是原始未截断分布的均值，而非截断后的条件均值。这就是你看到预测值低于观测均值的核心原因——观测数据已经剔除了截断范围内的样本，但默认预测没有对应调整。

解决方案：手动计算截断后预测值

对于左截断（例如截断值为0，仅保留≥1的观测），截断后的条件均值公式为：
$$\text{mu}{\text{truncated}} = \frac{\text{mu}{\text{original}}}{1 - P(Y \leq \text{trunc_val})}$$
其中$\text{mu}_{\text{original}}$是模型预测的原始分布均值，$P(Y \leq \text{trunc_val})$是原始分布在截断点及以下的累积概率。

以下是完整代码示例，覆盖所有四种截断分布，并包含95%置信区间的计算：

1. 准备环境与模拟数据

library(glmmTMB)
library(compoisson) # 用于compois分布的密度/累积概率计算
library(MASS)       # 用于参数自助法

# 模拟左截断（截断值=0，观测值≥1）的计数数据
set.seed(123)
n <- 500
x <- rnorm(n)
mu_original <- exp(1 + 0.5*x) # 原始泊松均值
y <- rpois(n, mu_original)
y[y == 0] <- NA # 剔除截断值（0）的样本
y <- na.omit(y)
x <- x[!is.na(y)]

2. 拟合四种截断模型

# 注意：truncation参数要与数据的截断规则一致（这里是左截断0）
mod_pois <- glmmTMB(y ~ x, family = truncated_poisson(link = "log", truncation = 0))
mod_compois <- glmmTMB(y ~ x, family = truncated_compois(link = "log", truncation = 0))
mod_nb1 <- glmmTMB(y ~ x, family = truncated_nbinom1(link = "log", truncation = 0))
mod_nb2 <- glmmTMB(y ~ x, family = truncated_nbinom2(link = "log", truncation = 0))

3. 定义函数计算截断后预测值与置信区间

该函数支持delta法近似置信区间，后续也可扩展参数自助法：

get_truncated_preds <- function(model, newdata = NULL, trunc_val = 0) {
  if (is.null(newdata)) newdata <- model.frame(model)
  
  # 获取线性预测值（log尺度）及标准误
  linpred_out <- predict(model, newdata = newdata, type = "link", se.fit = TRUE)
  fit_link <- linpred_out$fit
  se_link <- linpred_out$se.fit
  
  # 原始分布均值（response尺度）
  mu_original <- exp(fit_link)
  
  # 提取模型参数与分布类型
  family <- model$modelInfo$family$family
  theta <- plogis(model$fit$par["theta"]) # compois的z参数；nbinom1/2的size参数
  
  # 计算原始分布在截断点的累积概率P(Y ≤ trunc_val)
  p_trunc <- switch(family,
                    "truncated_poisson" = ppois(trunc_val, lambda = mu_original),
                    "truncated_compois" = pcompois(trunc_val, lambda = mu_original, z = theta),
                    "truncated_nbinom1" = pnbinom(trunc_val, size = theta, mu = mu_original),
                    "truncated_nbinom2" = pnbinom(trunc_val, size = theta, mu = mu_original)
  )
  
  # 计算截断后的条件均值
  mu_truncated <- mu_original / (1 - p_trunc)
  
  # delta法计算95%置信区间（近似）
  d_p_dmu <- switch(family,
                    "truncated_poisson" = dpois(trunc_val, lambda = mu_original),
                    "truncated_compois" = dcompois(trunc_val, lambda = mu_original, z = theta) * 
                      (1/trunc_val - 1/mu_original) * mu_original^2 / (mu_original + trunc_val*(theta-1)),
                    "truncated_nbinom1" = dnbinom(trunc_val, size = theta, mu = mu_original) * 
                      (trunc_val/mu_original^2) * (theta + mu_original)/(theta + trunc_val),
                    "truncated_nbinom2" = dnbinom(trunc_val, size = theta, mu = mu_original) * 
                      trunc_val/(mu_original^2)
  )
  d_log_trunc_d_link <- 1 - (d_p_dmu / (1 - p_trunc)) * mu_original
  se_log_trunc <- se_link * abs(d_log_trunc_d_link)
  
  ci_lower <- exp(log(mu_truncated) - 1.96 * se_log_trunc)
  ci_upper <- exp(log(mu_truncated) + 1.96 * se_log_trunc)
  
  return(data.frame(
    x = newdata$x,
    mu_original = mu_original,
    mu_truncated = mu_truncated,
    ci_lower_delta = ci_lower,
    ci_upper_delta = ci_upper
  ))
}

4. 生成预测结果

# 生成用于预测的新数据序列
new_x <- seq(min(x), max(x), length.out = 100)
newdata <- data.frame(x = new_x)

# 获取四个模型的截断后预测值
preds_pois <- get_truncated_preds(mod_pois, newdata)
preds_compois <- get_truncated_preds(mod_compois, newdata)
preds_nb1 <- get_truncated_preds(mod_nb1, newdata)
preds_nb2 <- get_truncated_preds(mod_nb2, newdata)

# 验证：观测均值与截断预测均值的匹配度
obs_mean <- mean(y)
cat("观测数据均值：", round(obs_mean, 2), "\n")
cat("泊松模型中间x的截断预测均值：", round(preds_pois$mu_truncated[50], 2), "\n")

5. 更准确的置信区间：参数自助法

如果需要更可靠的置信区间，可使用参数自助法（以泊松模型为例）：

# 参数自助法抽样预测值
boot_reps <- 1000
boot_preds_pois <- replicate(boot_reps, {
  # 从模型参数的多元正态分布中抽样
  params <- MASS::mvrnorm(1, coef(mod_pois)$cond, vcov(mod_pois)$cond)
  # 计算线性预测与原始均值
  linpred <- params[1] + params[2]*newdata$x
  mu_original <- exp(linpred)
  # 计算截断后均值
  p_trunc <- ppois(0, lambda = mu_original)
  mu_truncated <- mu_original/(1 - p_trunc)
  return(mu_truncated)
})

# 计算95%自助置信区间
preds_pois$ci_lower_boot <- apply(boot_preds_pois, 1, quantile, 0.025)
preds_pois$ci_upper_boot <- apply(boot_preds_pois, 1, quantile, 0.975)

关键注意事项

截断值一致性：拟合模型时的truncation参数必须与数据的截断规则完全一致（右截断需调整累积概率的计算逻辑）。
compois依赖：使用truncated_compois时需提前安装compoisson包，否则无法计算累积概率与密度。
置信区间选择：delta法计算速度快但为近似；参数自助法更准确但计算成本较高，适合样本量较小或分布复杂的场景。

内容的提问来源于stack exchange，提问作者user2602640