如何强制bootMer忽略收敛警告，保留拟合模型用于bootstrap分析？

问题背景 #

我拟合了一个复杂的Beta-Binomial模型，响应变量的两个组成部分（Methylated_C和未甲基化变量）取值均很大，模型代码如下：

model <- glmmTMB(                            # Beta-Binomial模型
      cbind(Methylated_C, Total_C - Methylated_C) ~
      context * mattrat * offtrat +
      (1|Sample_Name) +  (1 | plantID) ,
      family = betabinomial(link = "logit"),
      data = globalml
)

运行后收到收敛警告：

In finalizeTMB(TMBStruc, obj, fit, h, data.tmb.old) :
Model convergence problem; false convergence (8). See vignette('troubleshooting'), help('diagnose')

经diagnose(model)诊断并参考Ben Bolker的建议后，确认这是参数估计精度过高导致的假收敛，模型摘要如下：

Family: betabinomial  ( logit )
 Formula:          cbind(Methylated_C, Total_C - Methylated_C) ~ context * mattrat *  
offtrat + (1 | Sample_Name) + (1 | plantID)
 Data: globalml

       AIC       BIC    logLik -2*log(L)  df.resid 
    5939.6    5989.8   -2954.8    5909.6       195 

 Random effects:

 Conditional model:
  Groups      Name        Variance Std.Dev.
  Sample_Name (Intercept) 0.1017   0.3188  
  plantID     (Intercept) 0.0355   0.1884  
 Number of obs: 210, groups:  Sample_Name, 70; plantID, 24
 
 Dispersion parameter for betabinomial family (): 4.61e+03 
 
 Conditional model:
                                       Estimate Std. Error z value Pr(>|z|)    
 (Intercept)                          -2.769080   0.108543  -25.51   <2e-16 ***
 contextCHG                           -0.197699   0.021492   -9.20   <2e-16 ***
 contextCHH                           -2.297768   0.045247  -50.78   <2e-16 ***
 mattratdrought                        0.353346   0.154621    2.29   0.0223 *  
 offtratdry                            0.009827   0.154812    0.06   0.9494    
 contextCHG:mattratdrought            -0.014047   0.028664   -0.49   0.6241    
 contextCHH:mattratdrought             0.006640   0.059616    0.11   0.9113    
 contextCHG:offtratdry                -0.014767   0.030834   -0.48   0.6320    
 contextCHH:offtratdry                -0.017955   0.065010   -0.28   0.7824    
 mattratdrought:offtratdry            -0.156159   0.218728   -0.71   0.4753    
 contextCHG:mattratdrought:offtratdry  0.023604   0.040963    0.58   0.5645    
 contextCHH:mattratdrought:offtratdry  0.022354   0.085255    0.26   0.7932    
 ---
 Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Bootstrap尝试遇到的问题 #

我尝试用Bootstrap方法估计p值，修改模拟数据后用update创建模型并调用bootMer：

simulML <- simulate(model)[[1]]
simulMLdf <- data.frame("Methylated_C"= simulML[,1], "Total_C"= simulML[,1]+simulML[,2])
newdata <- cbind(simulMLdf, model$frame[-1])
model_R <- update(model, data=newdata)
b1 <- lme4::bootMer(model_R, FUN=function(x) fixef(x)$cond, nsim=50, .progress="txt")

但所有50次模拟均因假收敛问题失败，提示：

Warning message:
In lme4::bootMer(model_R, FUN = function(x) fixef(x)$cond, nsim = 50,  :
some bootstrap runs failed (50/50)

解决方案 #

1. 自定义拟合函数忽略收敛警告（适配bootMer） #

bootMer没有直接强制使用带警告模型的参数，但可以通过自定义拟合函数，忽略收敛警告并强制返回结果：

# 自定义拟合函数，关闭收敛检查并忽略警告
my_fit_fun <- function(formula, data, start, control, ...) {
  # 设置glmmTMB控制参数，忽略收敛问题
  ctrl <- glmmTMBControl(
    ignoreConvergence = TRUE, # 直接返回拟合结果，不触发收敛错误
    optCtrl = list(maxfun = 1e5) # 增加迭代次数，减少真收敛失败概率
  )
  # 忽略拟合时的警告
  fit <- suppressWarnings(glmmTMB(formula, data = data, start = start, control = ctrl, ...))
  return(fit)
}

# 使用自定义函数调用bootMer
b1 <- lme4::bootMer(model_R, 
                    FUN = function(x) fixef(x)$cond, 
                    nsim = 50, 
                    .progress = "txt",
                    fitFun = my_fit_fun)

2. 调整glmmTMB拟合控制参数，减少假收敛 #

假收敛多因优化器迭代次数不足或阈值过严，可先调整原模型的控制参数，再做Bootstrap：

# 先拟合调整控制参数的模型
model_ctrl <- glmmTMB(
  cbind(Methylated_C, Total_C - Methylated_C) ~
    context * mattrat * offtrat +
    (1|Sample_Name) +  (1 | plantID) ,
  family = betabinomial(link = "logit"),
  data = globalml,
  control = glmmTMBControl(
    optCtrl = list(maxfun = 1e5, reltol = 1e-6), # 增加迭代次数、放宽收敛阈值
    optimizer = "bobyqa" # 换用更稳定的优化器
  )
)

# 基于调整后的模型做Bootstrap
simulML <- simulate(model_ctrl)[[1]]
simulMLdf <- data.frame("Methylated_C"= simulML[,1], "Total_C"= simulML[,1]+simulML[,2])
newdata <- cbind(simulMLdf, model_ctrl$frame[-1])
model_R_ctrl <- update(model_ctrl, data=newdata)

b1 <- lme4::bootMer(model_R_ctrl, 
                    FUN=function(x) fixef(x)$cond, 
                    nsim=50, 
                    .progress="txt")

3. 用boot包手动实现Bootstrap，完全控制拟合过程 #

如果bootMer限制过多，可使用boot包手动编写Bootstrap循环，自主处理收敛警告：

library(boot)

# 定义Bootstrap统计量计算函数
boot_stat <- function(data, indices) {
  # 抽取Bootstrap样本
  boot_data <- data[indices, ]
  # 拟合模型，忽略收敛警告并关闭收敛检查
  fit <- suppressWarnings(glmmTMB(
    cbind(Methylated_C, Total_C - Methylated_C) ~
      context * mattrat * offtrat +
      (1|Sample_Name) +  (1 | plantID) ,
    family = betabinomial(link = "logit"),
    data = boot_data,
    control = glmmTMBControl(ignoreConvergence = TRUE)
  ))
  # 返回固定效应参数
  return(fixef(fit)$cond)
}

# 执行Bootstrap
b1 <- boot(data = globalml, statistic = boot_stat, R = 50)

注意事项 #

即使强制使用带警告的模型，也要检查Bootstrap样本的参数估计是否稳定：比如对比原模型参数与Bootstrap样本的参数分布，避免因真拟合失败导致结果偏差。

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