相对风险多变量回归报错求助：logbin模型MLE边界问题

问题背景 #

尝试通过log-binomial模型进行相对风险多变量回归时，出现拟合警告，无法计算参数的p值和置信区间，具体表现为NaN结果及参数空间边界相关警告。

所用代码 #

fit.glm <- glm(RemDAS ~ Sex + Age, data = db, family=binomial())

fit.logbin <- logbin(formula(fit.glm), data = db)
coef(summary(fit.logbin))[,'Pr(>|z|)']
exp(coef(fit.logbin))
exp(confint(fit.logbin))

数据概览 #

> db$RemDAS
 [1] 1 1 1 1 1 1 0 1 1 1 0 1 1 1 0 0 1 0 1 1 1 1 0 1 0 0 1 1 0 1 0 1 1 0 0 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 0 1 1 0
[61] 1 0 1 1 0 1 1 0 0 1 1 0 1 0 1 1 1 0 1 0 1 0 1 0 1 1 1 0 1
Levels: 0 1
> db$Sex
 [1] 0 0 0 0 1 0 0 1 1 0 1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 1 1 1 0 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 1 0 1 0 0 0 0
[61] 1 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 1 0 1 0 0 1 0 0 1 0 0 1 1
Levels: 0 1
> db$Age
 [1] 68 35 40 41 63 46 65 65 49 54 85 65 39 61 55 62 21 46 58 49 71 46 37 31 59 49 62 52 71 34 62 59 65 58 74 58 48 51 76 53
[41] 53 28 77 55 50 42 55 48 32 36 40 41 51 64 54 58 20 67 70 52 53 58 49 42 58 48 43 60 77 69 49 60 38 73 45 72 46 60 39 47
[81] 59 58 51 47 67 71 37 74 46

报错信息 #

> fit.logbin <- logbin(formula(fit.glm), data = db)
Warning message:
nplbin: fitted probabilities numerically 1 occurred 
> 
> coef(summary(fit.logbin))[,'Pr(>|z|)']
(Intercept)        Sex1         Age 
        NaN         NaN         NaN 
Warning message:
MLE on boundary of parameter space, cannot use asymptotic covariance matrix 
> 
> exp(coef(fit.logbin))
(Intercept)        Sex1         Age 
  1.0640070   1.1193059   0.9913005 
> exp(confint(fit.logbin))
            2.5 % 97.5 %
(Intercept)   NaN    NaN
Sex1          NaN    NaN
Age           NaN    NaN
Warning message:
MLE on boundary of parameter space, cannot use asymptotic covariance matrix 

解决方案 #

1. 核心问题分析 #

出现警告的原因是结局变量RemDAS中1的占比过高，导致log-binomial模型拟合时预测概率接近1，最大似然估计（MLE）落在参数空间边界，无法计算渐近协方差矩阵，进而导致p值和置信区间为NaN。

2. 替代模型：泊松回归+稳健标准误 #

当log-binomial模型拟合失败时，可改用泊松回归配合稳健标准误，同样能得到相对风险估计，且不受高事件率影响：

# 拟合泊松模型
fit.pois <- glm(RemDAS ~ Sex + Age, data = db, family = poisson())

# 加载计算稳健标准误的包
library(sandwich)
library(lmtest)

# 输出带稳健标准误的系数检验结果
coeftest(fit.pois, vcov = vcovHC(fit.pois, type = "HC3"))

# 计算相对风险（RR）
exp(coef(fit.pois))

# 计算稳健置信区间
exp(confint.default(fit.pois))

3. 调整log-binomial模型拟合参数 #

尝试调整logbin函数的拟合控制参数，增加迭代次数或调整收敛阈值，可能解决边界估计问题：

fit.logbin_adjusted <- logbin(
  formula(fit.glm), 
  data = db, 
  control = list(maxit = 1000, epsilon = 1e-8)
)

# 重新查看结果
coef(summary(fit.logbin_adjusted))
exp(confint(fit.logbin_adjusted))

4. 改用拟二项log链接模型 #

使用quasibinomial族配合log链接，可避免严格的二项分布假设，有时能解决高事件率下的拟合问题：

fit.qlog <- glm(RemDAS ~ Sex + Age, data = db, family = quasibinomial(link = "log"))

# 查看结果
coef(summary(fit.qlog))
exp(coef(fit.qlog))
exp(confint(fit.qlog))

5. 检查数据分离情况 #

检查自变量与结局变量的交叉表，确认是否存在完全分离（某自变量组合下结局全为1/0）：

# 查看Sex与RemDAS的交叉表
table(db$Sex, db$RemDAS)

# 将Age分组后查看与RemDAS的关系
db$Age_group <- cut(db$Age, breaks = 3)
table(db$Age_group, db$RemDAS)

若存在完全分离，需考虑合并分组或移除相关变量。

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