如何正确将SAS负二项回归模型转换为R代码复现结果

你用glmer.nb拟合的是带受试者随机截距的负二项广义线性混合模型（GLMM），但SAS中PROC GENMOD搭配repeated SUBJECT = ID语句拟合的是负二项分布广义估计方程（GEE）模型，两类模型存在本质区别：

• GLMM属于条件模型，估计的是控制个体随机效应后的受试者特异性处理效应，基于随机效应分布假设做极大似然估计
• SAS的GENMOD GEE属于边际模型，估计的是总体平均水平的处理效应，采用拟似然估计+三明治稳健标准误计算
  二者估计目标、计算逻辑都不相同，出现结果差异是模型选型不匹配导致的，不是代码语法错误。

使用geepack包的GEE拟合函数即可完全匹配SAS输出，注意对齐SAS的默认设置：负二项为NB2分布、默认独立工作相关结构、默认稳健标准误，同时要保证分类变量参考组和SAS一致。

# 加载依赖包
library(geepack)
library(MASS)

# 提前处理TRT变量，保证参考组和SAS一致，示例以TRT="A"为参考组
你的数据集名$TRT <- relevel(factor(你的数据集名$TRT), ref = "A")

# 预拟合常规负二项模型，估计离散参数theta，匹配SAS的NB2分布设定
nb_pre <- glm.nb(EVENTS ~ TRT + offset(lnTIME), data = 你的数据集名)
sas_theta <- nb_pre$theta

# 拟合和SAS完全一致的GEE负二项模型
sas_matched_model <- geeglm(
  formula = EVENTS ~ TRT + offset(lnTIME),
  family = negative.binomial(theta = sas_theta),
  data = 你的数据集名,
  id = ID, # 指定聚类/受试者ID
  corstr = "independence", # 匹配SAS PROC GENMOD repeated默认的独立工作相关结构
  std.err = "san.se" # 匹配SAS默认的三明治稳健标准误
)

# 输出结果
summary(sas_matched_model)

• 如果你在SAS的repeated语句中手动指定了相关结构，比如type=exch（交换相关）、type=ar1（一阶自相关），只需要把上述代码中corstr参数改为对应值即可：交换结构对应"exchangeable"，一阶自相关对应"ar1"。
• 不要尝试用GLMM类函数（包括glmer.nb、glmmTMB随机效应模型）匹配SAS的GENMOD GEE结果，二者的参数解释本身就存在差异，即使数值接近也不是同一估计目标的结果。
• 如果结果仍有极小幅数值差，可检查离散参数的估计值是否和SAS输出完全对齐，也可以手动将SAS输出的离散参数赋值给theta参数，消除迭代估计的精度差。

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