sommer包mmer函数非结构化方差模型奇异及负方差问题咨询

问题背景 #

• 分析场景：4个环境下水稻育种系多点试验分析，此前已在lme4、nlme包中实现非结构化方差模型拟合，本次计划使用sommer包完成同类分析
• 初始拟合代码：

m3 <- mmer(RDM ~ ENV ,
              random=~ vsr(usr(ENV),GEN),
              rcov=~ units,
              data=d)

• 初始报错信息：

System is singular (V). Stopping the job. Try a bigger number of tolparinv.

• 临时调整现象：将tolparinv参数调至1000时模型可运行，但输出多个不符合方差数值逻辑的负方差分量。

问题1：为什么模型会快速出现奇异问题？ #

核心原因是当前指定的无结构（US）方差协方差结构待估参数过多，和数据可提供的有效自由度不匹配：

• 4个环境下的US结构，仅基因型与环境互作项就需要估计4个环境的遗传方差、6个环境间的两两遗传协方差，共10个（n*(n+1)/2）方差协方差参数，再加上残差方差，待估参数总量较高。
• 如果数据存在部分基因型在个别环境无观测值、单环境内未设置生物学重复、部分基因型多环境表型共线性过高、单环境内表型变异极小的情况，会直接导致混合模型的V矩阵（边缘方差协方差矩阵）不满秩，迭代初始阶段就触发奇异检测，直接终止任务。
• 此前在lme4、nlme包中未遇到即时报错，是因为两个包默认的矩阵奇异容忍阈值更宽松，不会在迭代初期直接终止计算，不代表模型本身不存在参数识别问题。

问题2：将tolparinv调至1000的大取值会产生什么影响？ #

tolparinv是sommer包控制矩阵求逆奇异值容忍度的核心参数，默认取值为1e-6量级，作用是判定矩阵特征值是否为0（即矩阵是否奇异）：

• 参数取值越大，求逆过程中会将越小的特征值判定为非零值，可强行对不满秩的矩阵完成广义逆计算，跳过奇异终止判断。
• 取值为1000属于极端宽松的阈值，相比默认值放大了近9个数量级，相当于完全放弃迭代过程中对V矩阵正定性的校验，最终得到的参数估计标准误会被严重低估，固定效应、随机效应的显著性检验完全不可靠，方差协方差参数估计值也会大幅偏离真实值。

问题3：tolparinv设置过大是否是导致负方差结果的原因？ #

是直接诱因，但不是根本原因：

• 根本原因是当前模型复杂度超过了数据可支撑的参数识别上限，参数空间内不存在满足正定要求的方差协方差矩阵合法解。
• 大阈值的tolparinv相当于移除了sommer默认对方差分量的边界约束（方差≥0、协方差矩阵正定），允许迭代过程进入参数空间的非法区域，最终输出负方差分量。如果保持默认阈值，模型只会触发奇异报错终止，不会输出不符合统计逻辑的估计值。

问题4：该问题是否可修复？具体修复方法有哪些？ #

该问题完全可修复，按优先级从高到低推荐以下方案：

• 校正模型结构错误：当前代码中rcov=~units默认假设所有环境残差方差同质，和多点试验的实际数据特征不符，是诱发矩阵奇异的高频错误。拟合非结构化遗传协方差时，需要同步指定异质残差结构，将rcov项修改为rcov=~vsr(dsr(ENV), units)，允许每个环境独立估计残差方差，使用默认tolparinv参数重新运行，绝大多数场景下可直接解决奇异问题。
• 核查数据结构：确认每个供试基因型在4个环境中均有表型观测值，剔除仅在单个环境出现的基因型；确认每个环境内设置了足够的生物学重复，若单环境无重复，无法独立分离残差与G×E互作方差，不要强行拟合全非结构化协方差结构。
• 合理降低协方差结构复杂度：如果数据量确实不足以支撑10个参数的全US结构估计，可替换为更简约的协方差结构，比如cs(ENV)（复合对称结构）、ar1(ENV)（一阶自回归结构）或fa(ENV, k=1/2)（低阶因子分析结构），这类结构待估参数更少，对数据量要求更低，不易出现奇异问题。
• 参数调整兜底方案：如果确实需要拟合全US结构，不要直接将tolparinv拉到1000这类极端值，可尝试设置tolparinv=1e-2，同时添加参数constraint="P"强制方差协方差矩阵保持正定，避免出现负方差；该方法得到的结果需要和简约模型结果交叉比对，确认估计偏差在可接受范围内再使用。

内容的提问来源于stack exchange，提问作者Sergio Antonio Castro Pacheco