lme模型收敛错误代码1：循环跳过方法及单独运行成功原因咨询

我对模型收敛问题的理解还不全面，但知道代码里有一些影响收敛的关键设置。由于数据量太大，无法提供完整数据库，现提出两个技术问题：

1. 如何在循环中跳过出现**收敛错误（错误代码1）**的变量？想要实现类似如下逻辑：

# 示例逻辑
if(error code = 1) next

2. 为什么单独对引发问题的变量rxra运行lme模型能得到正常结果，但在循环中运行却报错：nlminb problem, convergence error code = 1 message = false convergence (8)？

循环运行的代码及报错信息 #

genes_wide_s18 <- genes_s18 %>% pivot_wider(names_from = gen, values_from = dCt_s18)

genes_wide_s18$time [genes_wide_s18$time == "1"] <- 1 
genes_wide_s18$time [genes_wide_s18$time == "3"] <- 3 

genes_wide_s18$time<-as.numeric(genes_wide_s18$time)


# 变量名称
genes.names <- colnames(genes_wide_s18)[25:79]
gene.list$s18 <- NULL
no.genes <- length(genes.names)

# 创建用于存储拟合模型的命名列表
gene.list <- as.list(1:no.genes)
names(gene.list) <- genes.names

# 遍历基因名称
for(i in genes.names){ 

  # 打印运行状态
  print(paste("Running entity:", i, "which is", which(genes.names==i), "out of", no.genes))
  
  # 跳过观测值过少的变量
 if(i == "s18")next
     print(i)

  
  # 创建临时数据矩阵和模型公式
  # 第一种：调整多个变量（后续被覆盖，无实际作用）
  tmp <- genes_wide_s18[, c(i,"edad0","time", "grup_int", "id", "sexo", "peso1")]
  fml <- as.formula( paste( i, "~", paste(c("sexo:peso1", "edad0","time", "grup_int", "time:grup_int"), collapse="+")))

   
  # 第二种：仅按组调整（实际使用的模型）
  tmp <- genes_wide_s18[, c(i,"grup_int","time", "id")]
  fml <- as.formula( paste( i, "~", paste(c("time:grup_int"), collapse="+")))
  
  # 第一次拟合（后续被覆盖，无实际作用）
  gene.list[[i]] <- lme(fml, random= ~ time|id, method="REML", data=tmp, na.action = na.omit)
  
  # 第二次拟合（最终生效的模型）
  gene.list[[i]] <- lme(fml, random= ~ time|id, control = lmeControl(opt = "optim"), method="REML", data=tmp, na.action = na.omit)

 
}

运行到rxra时的报错：

[1] "Running entity: rxra which is 6 out of 55"
[1] "rxra"
Error in lme.formula(fml, random = ~time | id, method = "REML", data = tmp,  : 
  nlminb problem, convergence error code = 1
  message = false convergence (8)

单独运行rxra的代码及结果 #

lme(as.formula( paste( "rxra", "~", paste(c("time:grup_int"), collapse="+"))), random= ~ time|id, control = lmeControl(opt = "optim"), method="REML", data=tmp, na.action = na.omit)

输出结果：

Linear mixed-effects model fit by REML
  Data: tmp 
  Log-restricted-likelihood: -392.218
  Fixed: as.formula(paste("rxra", "~", paste(c("time:grup_int"), collapse = "+"))) 
   (Intercept) time:grup_int1 time:grup_int2 time:grup_int3 
     3.6723896     -0.2699351     -0.2630362     -0.2104078 

Random effects:
 Formula: ~time | id
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 1.6296951 (Intr)
time        0.7984285 -0.982
Residual    0.5870576       

Number of Observations: 279
Number of Groups: 140 

数据结构信息 #

str(tmp)

输出：

tibble [15,968 × 4] (S3: tbl_df/tbl/data.frame)
 $ rxra    : num [1:15968] NA NA NA NA NA ...
 $ grup_int: Factor w/ 3 levels "1","2","3": 3 3 3 3 3 3 3 3 3 3 ...
  ..- attr(*, "label")= chr "GENERAL: Grupo de intervención"
 $ time    : num [1:15968] 1 1 1 1 1 1 1 1 1 1 ...

问题1：循环中跳过收敛错误的变量 #

用tryCatch()捕获错误，匹配收敛错误代码后直接跳过当前迭代，修改后的循环核心代码如下：

for(i in genes.names){ 
  print(paste("Running entity:", i, "which is", which(genes.names==i), "out of", no.genes))
  
  if(i == "s18")next
  print(i)

  # 创建当前基因的专属临时数据
  tmp <- genes_wide_s18[, c(i,"grup_int","time", "id")]
  fml <- as.formula( paste( i, "~", paste(c("time:grup_int"), collapse="+")))
  
  # 捕获并处理收敛错误
  fit_result <- tryCatch(
    expr = lme(fml, random= ~ time|id, control = lmeControl(opt = "optim"), method="REML", data=tmp, na.action = na.omit),
    error = function(e) {
      # 匹配收敛错误代码1的信息
      if(grepl("convergence error code = 1", e$message)){
        message(paste("跳过变量", i, "：收敛错误代码1"))
        return(NULL)
      } else {
        # 其他错误正常抛出，避免遗漏问题
        stop(e)
      }
    }
  )
  
  # 仅当拟合成功时存入列表，否则跳过
  if(!is.null(fit_result)){
    gene.list[[i]] <- fit_result
  } else {
    next
  }
}

问题2：单独运行正常但循环中报错的原因及解决 #

核心原因 #

1. 临时变量污染：循环中先创建了包含多变量的tmp，之后重新赋值为4列的tmp，但R的变量环境可能存在残留，导致rxra的tmp数据和单独运行时不一致。
2. 优化器初始值复用：循环连续拟合模型时，lme可能复用前一个模型的优化器初始值，这些值对rxra模型不适用，导致收敛失败；单独运行时则使用默认初始值，恰好能收敛。
3. 冗余代码干扰：循环中连续两次调用lme，第一次的拟合过程可能影响了第二次的初始条件。

解决方法 #

1. 清理临时变量：每次循环开始时清除旧的tmp，避免残留：

for(i in genes.names){ 
  # 清除上一次循环的临时变量
  if(exists("tmp")) rm(tmp)
  
  # 后续代码不变...
}

2. 强制重置优化器初始值：在lmeControl中指定初始值为NULL，强制从默认值开始优化：

gene.list[[i]] <- lme(fml, random= ~ time|id, 
                      control = lmeControl(opt = "optim", optArgs = list(start = NULL)), 
                      method="REML", data=tmp, na.action = na.omit)

3. 删除冗余代码：直接删掉循环中第一次的lme调用，减少不必要的计算和干扰。

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