R中GLMMadaptive拟合二项重复数据报Hessian非正定如何解决

云上老周

2026-9-1

问题背景

在R环境中使用GLMMadaptive包分析二项分布重复测量数据，结局分别在基线、12周、24周3个时间点采集，数据存在随机缺失（MAR），需计算结局边际概率，初始建模设定为：固定效应纳入Time*Treatment + Duration，随机效应仅纳入个体水平随机截距，初始运行代码如下：

fm <- mixed_model(fixed = Outcome ~ Time * Treatment + Duration,
                  random = ~ 1 | ID,
                  data = data2,
                  family = binomial(),
                  nAGQ = 15,
                  iter_EM=0,
                  max_coef_value = 50,
                  initial_values = list(betas = rep(0, 5)),
                  na.action=na.exclude)

模型运行时触发警告：

Hessian matrix at convergence is not positive value: unstable solution.
后续调用summary(fm)查看结果时，标准误（SE）与P值均显示为NAs，以下为可落地的排查方向与解决方案。

排查方向与解决方案

优先修正初始值与迭代参数设置
当前设置iter_EM=0直接跳过了EM算法初始值寻优步骤，且手动指定的初始值betas = rep(0, 5)对二项混合模型来说极易让优化器落入参数空间的平坦区域，直接导致海森矩阵奇异。对应调整方式：
- 将iter_EM恢复为默认值30，EM算法的作用就是为后续牛顿迭代找到稳定的初始值，手动关闭会大幅提升收敛失败概率
- 暂时删除手动设置的initial_values参数，让包自动计算初始值；如果自动初值仍不收敛，可以先拟合不纳入随机效应的普通glm(Outcome ~ Time * Treatment + Duration, family = binomial(), data = data2)，将glm输出的系数作为betas初始值传入，不要使用全0初始值
- 初始调试阶段可将nAGQ降到7或11，15个节点的自适应高斯求积虽然精度更高，但初值不稳定时优化曲面复杂度也更高，更容易出现收敛问题，待模型稳定收敛后再逐步提升nAGQ验证结果一致性即可
排查数据层面的固有问题
- 检查二项结局分离问题：用table(data2$Outcome, data2$Time, data2$Treatment)逐个查看分组下的结局分布，如果存在某处理组某时间点结局全为0/全为1的完全/准完全分离情况，固定效应系数会趋向无穷大，直接导致海森矩阵非正定。出现该问题时可考虑对时间点做合理合并，或采用带弱信息先验的贝叶斯估计方式约束系数范围
- 检查共线性：计算固定效应项的方差膨胀因子，Time如果采用0/12/24的数值编码，和交互项的共线性会很高，建议先对Duration做个体均值中心化，Time采用以基线为参照的哑变量/0,1,2编码降低共线性，若VIF大于10需先处理共线性再建模
- 检查分组样本量：统计每个ID对应的观测条数，如果存在大量ID仅完成1次随访，随机截距的方差会因信息不足无法稳定估计，容易卡在方差为0的参数边界，导致海森矩阵奇异，可先过滤仅含1条观测的ID后重新拟合
调整优化器与数值计算参数
- 当前设置的max_coef_value = 50阈值过高，会允许优化器将系数迭代到完全不合理的极值区域，建议改回默认值10或设为15，及时拦截异常迭代路径
- 可在mixed_model()中加入参数optimizer = "optim"，换用L-BFGS-B优化器替代默认优化方法，不同优化器的收敛路径存在差异，多数情况下换用优化器即可解决非正定问题
- 如果确认模型实际已收敛（可通过fm$converged查看返回值是否为TRUE），SE为NA可能是海森矩阵数值计算误差导致，可加入参数hessian_args = list(eps = 1e-4, type = "central")，改用中心差分法、调整差分步长重新计算海森矩阵，即可得到正常的标准误结果
结合MAR缺失特征做验证
先统计整体缺失比例，如果缺失比例超过30%，直接用原始数据建模会因信息不足出现估计不稳定，建议先通过mice包做多重插补，对插补后的多个数据集分别拟合模型后再合并结果，符合MAR假设下的估计要求。

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