嘿，你已经用glm.nb()搞定了基础的负二项回归模型，接下来要加交互项其实超简单，我给你拆解清楚操作方法和后续的分析要点：

负二项模型添加交互项的语法和普通线性模型完全一致，有两种常用写法：

• 自动包含主效应+交互项：用*符号，它会自动帮你包含两个变量的主效应，以及它们的交互项，代码如下：

# 拟合包含age、treatment主效应和交互项的负二项模型
ngbinmodel_interact <- glm.nb(seizure.rate ~ age * treatment, data = epilepsy_reduced)
summary(ngbinmodel_interact)

这个写法等价于seizure.rate ~ age + treatment + age:treatment，适合你想同时保留主效应和交互项的场景。

• 手动指定主效应+交互项：用:符号单独表示交互项，如果你已经明确要保留主效应，只想额外添加交互项，可以这么写：

# 手动声明主效应和交互项
ngbinmodel_interact <- glm.nb(seizure.rate ~ age + treatment + age:treatment, data = epilepsy_reduced)
summary(ngbinmodel_interact)

两种写法的模型结果完全一致，只是后者更直观展示所有变量。

添加完交互项后，重点要做这些分析：

2.1 先看交互项的统计显著性 #

查看summary(ngbinmodel_interact)输出中age:treatment对应的Pr(>|z|)值：

• 如果p值小于你的显著性水平（通常是0.05），说明交互项统计显著，意味着年龄对癫痫发作率的影响在不同治疗组之间存在差异；
• 如果p值不显著，说明年龄的效应在各治疗组中没有明显区别，可能不需要保留交互项。

2.2 解读交互项的实际含义（基于log链接） #

因为你用的是负二项模型的默认log链接，交互项的原始系数是对数尺度的，需要指数化后才能解读实际意义：

# 指数化交互项系数，得到倍率值
exp(coef(ngbinmodel_interact)["age:treatment"])

举个例子，如果指数化后的值是1.03，说明：在treatment组B（假设是二分类变量）中，年龄每增加1岁，发作率的变化倍率比treatment组A高3%。

2.3 可视化交互效应（更直观） #

用可视化能更清晰展示不同治疗组中年龄与发作率的关系差异，推荐用ggplot2绘制：

library(ggplot2)

# 生成覆盖全范围的预测数据
pred_data <- expand.grid(
  age = seq(min(epilepsy_reduced$age), max(epilepsy_reduced$age), length.out = 100),
  treatment = unique(epilepsy_reduced$treatment)
)

# 预测发作率（type="response"会自动反log转换到原始尺度）
pred_data$pred_seizure <- predict(ngbinmodel_interact, newdata = pred_data, type = "response")

# 绘制交互趋势图
ggplot(pred_data, aes(x = age, y = pred_seizure, color = treatment)) +
  geom_line(linewidth = 1) +
  labs(x = "年龄", y = "预测发作率", color = "治疗组") +
  theme_minimal()

从图里你能直接看到两条线的斜率差异，这就是交互效应的直观体现。

2.4 模型比较：确认交互项是否提升拟合效果 #

用似然比检验对比基础模型和带交互项的模型，判断交互项是否真的让模型更优：

# 似然比检验，比较两个嵌套模型
anova(ngbinmodel, ngbinmodel_interact, test = "Chisq")

如果检验的p值显著，说明带交互项的模型拟合效果显著优于基础模型，值得保留交互项；反之则可以考虑移除交互项，回到基础模型。

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