如何解决R中glmer的PIRLS循环NaN值错误及获取风险比？

一、解决PIRLS loop resulted in NaN value报错 #

这个错误通常源于log链接的二项式模型拟合时，预测概率趋近于0或1导致数值不稳定，或是随机效应分组存在极端结果（全0/全1）。以下是具体修正方案：

• 检查并处理极端分组：先查看每个随机效应组的结果分布：

  table(df$`random effect variable`, df$outcome)
  

  若存在某个组的outcome全为0或全为1，这类组会导致随机效应估计失效，可选择删除该组，或在业务逻辑允许的前提下与相似组合并。

• 更换优化器并调整迭代次数：默认优化器在log链接下稳定性不足，换用bobyqa优化器并增加最大迭代次数：

  library(lme4)
  model <- glmer(outcome ~ predictor + (1|`random effect variable`), 
                 data = df, 
                 family = binomial(link = "log"),
                 control = glmerControl(optimizer = "bobyqa", 
                                        optCtrl = list(maxfun = 20000)))
  summary(model, corr = FALSE)
  

• 改用广义估计方程（GEE）：若随机效应的个体间相关性不是核心研究目标，GEE对极端数据更稳健，可使用geepack包实现：

  library(geepack)
  model <- geeglm(outcome ~ predictor, 
                  data = df, 
                  family = binomial(link = "log"),
                  id = `random effect variable`)
  summary(model)
  

• 贝叶斯方法拟合：用brms包的贝叶斯模型，通过先验约束避免数值问题，同时能直接得到风险比的后验分布：

  library(brms)
  model <- brm(outcome ~ predictor + (1|`random effect variable`), 
               data = df, 
               family = binomial(link = "log"))
  summary(model)
  

二、关于风险比（RR）的计算 #

是的，指数化log链接二项模型的回归系数即可得到风险比。模型的核心形式为：

log(p) = β₀ + β₁×predictor

指数化后转化为：

p = exp(β₀) × exp(β₁)^predictor

其中exp(β₁)就是predictor每增加1单位时，事件发生概率的倍数，即风险比。

计算RR及95%置信区间的代码示例：

# 提取系数并指数化得到RR
rr <- exp(coef(model)$`random effect variable`[1, ])
# 计算95%置信区间
rr_ci <- exp(confint(model))

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