如何在sjPlot的plot_model()中处理offset变量并调整Y轴为比率

问题背景 #

我正在构建含offset变量的Poisson分布广义线性模型，响应变量为狍表现警戒行为的时长，offset为在样地停留的总时长。参考Zuur《Mixed Effect Models for Ecology》修改的可复现代码如下：

library(glmmTMB)
Owls$NCalls <- Owls$SiblingNegotiation

library(lme4)
Formula <- formula(NCalls ~ offset(logBroodSize) + SexParent * FoodTreatment + SexParent * ArrivalTime)
fit1 <- glm(Formula, data = Owls, family = poisson(link="log"))

（注：因习惯使用lme4，故用glmmTMB的数据集在lme4中建模）

尝试用sjPlot包的plot_model()函数绘制结果：

library(sjPlot)
plot_model(fit1, type="pred")

但绘图结果未纳入logBroodSize变量，而是单独展示offset。我希望每个回归图的Y轴为NCalls/BroodSize（警戒时长/总时长的比率），这对研究至关重要。

曾尝试将offset变量直接加入响应变量，但出现“值为非整数”的警告：

Formula.2 <- formula(NCalls/logBroodSize ~ SexParent * FoodTreatment + SexParent * ArrivalTime)
fit2 <- glm(Formula.2, data = Owls, family = poisson(link="log"))

考虑过转换为高斯分布，但数据含大量0值，效果不理想。

解决方案 #

方法1：基于原Poisson模型，手动转换预测值为比率 #

你的原模型fit1的本质是拟合：
$$\log(E[NCalls]) = \text{协变量效应} + \log(\text{BroodSize})$$
整理后等价于：
$$\log\left(E\left[\frac{NCalls}{\text{BroodSize}}\right]\right) = \text{协变量效应}$$
因此，只需将模型的预测值（对数尺度）转换为原始尺度后，除以BroodSize即可得到比率的期望。

可以通过以下步骤实现sjPlot绘图的调整：

1. 提取plot_model()的预测数据，手动计算比率：

# 获取预测数据
pred_data <- plot_model(fit1, type="pred", return_data = TRUE)
# 计算NCalls/BroodSize的期望：exp(predicted) / BroodSize
pred_data$ratio <- exp(pred_data$predicted) / pred_data$logBroodSize

# 用ggplot重新绘制（sjPlot基于ggplot2，可直接复用格式）
library(ggplot2)
ggplot(pred_data, aes(x = x, y = ratio, color = group)) +
  geom_line() +
  geom_ribbon(aes(ymin = exp(conf.low)/logBroodSize, ymax = exp(conf.high)/logBroodSize), alpha = 0.2) +
  facet_wrap(~group) +
  labs(y = "NCalls/BroodSize", x = "") +
  theme_sjplot()

2. 或者直接在plot_model()中通过transform参数调整（部分版本支持）：

plot_model(fit1, type="pred", transform = function(x) exp(x)/Owls$BroodSize)

注：若transform参数无法直接关联分组数据，优先用第一种手动提取数据的方法。

方法2：改用负二项模型（若存在过度离散） #

如果你的实际数据（狍的警戒时长）存在过度离散（Poisson模型的残差方差远大于均值），可以改用负二项模型，同样保留offset，后续绘图调整逻辑与Poisson模型一致：

library(MASS)
fit_nb <- glm.nb(NCalls ~ offset(logBroodSize) + SexParent * FoodTreatment + SexParent * ArrivalTime, data = Owls)
# 绘图调整同方法1

方法3：零膨胀Poisson模型（针对大量0值） #

若数据中0值占比极高，普通Poisson模型拟合效果差，可使用零膨胀Poisson（ZIP）模型，同样保留offset：

library(pscl)
fit_zip <- zeroinfl(NCalls ~ offset(logBroodSize) + SexParent * FoodTreatment + SexParent * ArrivalTime | 1, data = Owls, dist = "poisson")
# 提取预测值时需同时考虑计数部分和零膨胀概率，计算比率的期望：
pred_zip <- predict(fit_zip, type = "response")
ratio_zip <- pred_zip / Owls$BroodSize
# 再结合协变量绘制分组图

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