如何加速R语言lme4qtl包中线性混合模型的拟合过程？

问题背景 #

以亲缘矩阵作为随机效应拟合线性混合模型（LMM）在R环境中十分耗时，尤其是处理大规模标记数据集时。当前使用lme4qtl包的relmatLmer函数，处理500个样本、30K个因子型变异位点的数据集：

基准模型拟合 #

Model1 = relmatLmer(PHE ~ PC1+PC2+PC3+PC4+PC5+(1|ID), data, 
                    relmat = list(ID = Kin),
                    REML=F,
                    control = lmerControl(optimizer = "nloptwrap",optCtrl=list(xtol_abs=1e-6, ftol_abs=1e-6),calc.derivs = F))

批量检验变异位点 #

通过循环调用update加入单个变异位点G作为随机效应，再用anova比较模型获取显著性：

Model2 = update(Model1, .~.+(1|G))

即便用doSNOW开启40线程并行，整个过程仍耗时约20小时，现有优化手段效果有限，需更高效的加速方案。

可行加速方案 #

1. 换用专为GWAS/LMM优化的R包 #

放弃lme4qtl这类通用混合模型包，改用针对大规模标记优化的专用工具，这类包通常实现了更高效的REML算法（如平均信息AI-REML），避免重复计算冗余矩阵：

• sommer包：支持亲缘矩阵（加性遗传效应）和多随机效应，批量处理标记效率极高。示例代码：

  library(sommer)
  # 拟合基准模型（亲缘矩阵作为随机效应）
  base_mod = mmer2(PHE ~ PC1+PC2+PC3+PC4+PC5, 
                   random = ~ vsr(ID, Gu=Kin), 
                   data = data, REML=FALSE)
  # 批量检验所有变异位点（作为随机效应）
  # 可直接构建标记矩阵，用mmer2的批量模式或循环调用，底层已优化
  

• rrBLUP包：适合加性效应的GWAS分析，用mixed.solve快速拟合LMM，检验标记效应时无需重复拟合全模型：

  library(rrBLUP)
  # 先拟合亲缘关系对应的遗传方差
  fit = mixed.solve(y=data$PHE, X=model.matrix(~PC1+PC2+PC3+PC4+PC5, data=data), 
                    K=Kin, REML=FALSE)
  # 循环检验每个标记：利用残差计算得分统计量，大幅减少计算量
  

2. 改用得分检验替代全模型拟合 #

无需每次拟合加入G的完整模型，利用基准模型的残差和信息矩阵计算得分检验统计量，这是检验单个随机效应显著性的高效方法，计算量仅为全模型拟合的1/10甚至更低：

library(lmerTest)
# 对单个标记G，直接计算得分检验
score_res = scoreTest(Model1, ~ (1|G))
# 提取p值即可，无需拟合Model2和anova

此方法可将30K次模型拟合简化为30K次矩阵运算，直接缩短数小时耗时。

3. 优化模型拟合参数 #

调整lmerControl的优化器和收敛阈值，在保证结果稳定的前提下减少迭代次数：

• 换用bobyqa优化器：比nloptwrap更适合混合模型的参数优化，收敛速度更快

  control = lmerControl(optimizer = "bobyqa", optCtrl=list(maxfun=1e5), calc.derivs=F)
  

• 放宽收敛阈值：将xtol_abs和ftol_abs调整为1e-5，避免过度迭代（需验证结果稳定性）

4. 调整变异位点的效应类型 #

若研究目的是检验变异位点的加性效应显著性（而非方差组分），可将因子型标记转换为数值型（如0/1/2编码）作为固定效应纳入模型，而非随机效应：

# 转换标记为数值型
data$G_num = as.numeric(factor(data$G)) - 1
# 基准模型不变，检验时加入固定效应
Model2 = update(Model1, .~.+G_num)
# 用anova或summary提取p值

固定效应检验的计算量远低于随机效应，能大幅提升循环效率。

5. 优化并行计算逻辑 #

• 改用future+furrr的并行框架，比doSNOW的线程开销更小：

  library(furrr)
  plan(multisession, workers=40)
  # 批量处理标记
  results = future_map(markers_list, function(g) {
    # 这里用得分检验或简化模型
    scoreTest(Model1, ~ (1|get(g)))
  })
  

• 避免并行任务中重复复制基准模型：将Model1设为全局变量，或用foreach的.export参数仅传递必要对象，减少内存开销。

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