R语言自定义ICC函数循环调用报错问题求助

自定义组内相关系数（ICC）计算函数icc9单独调用时正常（如icc9(y="Anhydro_1.5_D_glucitol", subject = "RegNum", dat=df)），但循环遍历代谢物列时出现报错。

循环尝试与报错 #

第一种循环代码 #

lapply(metabolite.names, function(response) {
 icc9(as.formula(paste0("y=",response)),subject="RegNum",dat=df)
 }) %>% do.call("rbind",.)

报错与警告：

Error in class(ff) <- "formula" : cannot set attribute on a symbol
In addition: Warning message:
invalid formula "y=Anhydro_1.5_D_glucitol": assignment is deprecated 

第二种循环代码 #

lapply(metabolite.names, function(response) {  
icc9(y=response,subject="RegNum",dat=df)  }) %>% do.call("rbind",.)

报错：

boundary (singular) fit: see help('isSingular') Error in 
asMethod(object) : not a positive definite matrix In addition: 
Warning message: In sqrt(m) : NaNs produced

错误原因分析 #

1. 第一种循环错误：icc9的y参数要求传入列名字符串，而非公式。传递as.formula(paste0("y=",response))不符合参数定义，导致公式解析失败。
2. 第二种循环错误：原函数存在变量名拼写错误（将iestimates误写为estimates），导致方差参数提取失败，进而引发奇异拟合、矩阵非正定等连锁错误；同时函数内data与dat变量混用，数据逻辑混乱。

解决方案 #

步骤1：修复icc9函数bug #

修正变量名错误与数据逻辑问题：

icc9 <- function(y, subject, fixed=NULL, dat, method="ML"){
  # 检查并安装依赖包
  p <- rownames(installed.packages())
  if (!"devtools" %in% p) install.packages("devtools")
  if (!"bootmlm" %in% p) devtools::install_github("marklhc/bootmlm")
  
  # 统一清理缺失值
  dat_clean <- dat[!is.na(dat[[subject]]) & !is.na(dat[[y]]), ]
  
  # 定义分析变量
  SUBJECT <- dat_clean[[subject]]
  Y <- dat_clean[[y]]
  if (!is.null(fixed)) {
    FIXED <- as.matrix(dat_clean[, fixed])
  } else {
    FIXED <- NULL
  }
  
  # 计算受试者频次统计
  tabx.f <- table(SUBJECT)
  tabx.p <- prop.table(tabx.f)
  tabx <- data.frame(subjid = names(tabx.f), 
                     COUNT = as.numeric(tabx.f), 
                     PERCENT = as.numeric(tabx.p)*100,
                     stringsAsFactors = FALSE)
  
  tab3 <- data.frame(numw2 = nrow(tabx),
                     obsw2 = sum(tabx$COUNT, na.rm = TRUE))
  tab3$avgrep <- tab3$obsw2 / tab3$numw2
  
  # 设置模型拟合方法
  modelmethod <- ifelse(method == "ML", FALSE, TRUE)
  
  # 构建模型公式
  if (!is.null(FIXED)) {
    f <- formula(paste0("Y ~ ", paste(c("(1|SUBJECT)", "FIXED"), collapse = "+")))
  } else {
    f <- formula("Y ~ (1|SUBJECT)")
  }
  fit <- lme4::lmer(f, REML = modelmethod, data = dat_clean)
  z <- summary(fit)
  
  # 提取渐近协方差矩阵
  asycov <- data.frame(bootmlm::vcov_vc(fit, FALSE))
  names(asycov) <- c("CovP1", "CovP2")
  rownames(asycov) <- c("Intercept", "Residual")
  
  # 提取方差参数（修正变量名错误）
  iestimates <- data.frame(z$varcor)
  estimates <- iestimates[, c("grp", "vcov", "sdcor")]
  names(estimates) <- c("Parm", "Estimate", "SD")
  
  # 提取固定效应
  fixed_effects <- data.frame(z$coefficients)[1, ]
  names(fixed_effects) <- c("Estimate", "SE", "T")
  
  # 提取计算所需参数
  subject_estimate <- estimates$Estimate[estimates$Parm == "SUBJECT"]
  residual_estimate <- estimates$Estimate[estimates$Parm == "Residual"]
  residual_covp2 <- asycov["Residual", "CovP2"]
  intercept_covp1 <- asycov["Intercept", "CovP1"]
  intercept_covp2 <- asycov["Intercept", "CovP2"]
  fixed_estimate <- fixed_effects$Estimate
  variance_fixed <- fixed_effects$SE^2
  
  # 格式化结果函数
  formatEstimate <- function(est, lcl, ucl){
    a <- format(round(est, 2), nsmall = 2)
    b <- format(round(lcl, 2), nsmall = 2)
    c <- format(round(ucl, 2), nsmall = 2)
    return(paste0(a, " (", b, ",", c, ")"))
  }
  
  # 计算ICC及置信区间
  ICC_Variance <- (subject_estimate^2 * residual_covp2 + residual_estimate^2 * intercept_covp1 - 
                     2*subject_estimate*residual_estimate*intercept_covp2) / subject_estimate^4
  ICC <- subject_estimate / sum(estimates$Estimate)
  logICC <- log(ICC / (1 - ICC))
  SE_ICC <- sqrt(((ICC^2) * ICC_Variance) / ((1 - ICC)^2))
  logLL <- logICC - 1.96*SE_ICC
  logUL <- logICC + 1.96*SE_ICC
  ICCLL <- exp(logLL) / (1 + exp(logLL))
  ICCUL <- exp(logUL) / (1 + exp(logUL))
  finalICC <- formatEstimate(ICC, ICCLL, ICCUL)
  
  # 计算CV及置信区间
  Variance_CV <- residual_covp2 / (4*fixed_estimate^2 * residual_estimate) +
    (residual_estimate / fixed_estimate^4) * variance_fixed
  CV <- sqrt(residual_estimate) / fixed_estimate
  logCV <- log(CV)
  varCV <- (1/CV^2) * Variance_CV
  logLL_CV <- logCV - 1.96*sqrt(varCV)
  logUL_CV <- logCV + 1.96*sqrt(varCV)
  CVLL <- exp(logLL_CV)
  CVUL <- exp(logUL_CV)
  finalCV <- formatEstimate(CV, CVLL, CVUL)
  
  # 构建结果数据框
  final <- data.frame(Response = y, tab3, ICC = finalICC, CV = finalCV, stringsAsFactors = FALSE)
  names(final) <- c("Response_Variable", "N_subjects", "N_Observations", 
                    "Average_Measurements", "ICC", "CV")
  return(final)
}

步骤2：正确的循环调用方式 #

使用修复后的函数，添加错误处理避免个别拟合失败中断循环：

# 定义代谢物列名范围
metab.start <- which(colnames(df) == "Anhydro_1.5_D_glucitol")
metab.stop <- which(colnames(df) == "Ethanolamine")
metabolite.names <- colnames(df)[metab.start:metab.stop]

# 循环调用并捕获错误
library(dplyr)
results <- lapply(metabolite.names, function(response) {
  tryCatch({
    icc9(y = response, subject = "RegNum", dat = df)
  }, error = function(e) {
    # 拟合失败时返回含错误信息的行
    data.frame(Response_Variable = response,
               N_subjects = NA,
               N_Observations = NA,
               Average_Measurements = NA,
               ICC = paste0("Error: ", e$message),
               CV = NA,
               stringsAsFactors = FALSE)
  })
}) %>% do.call("rbind", .)

# 查看结果
print(results)

关键修复点说明 #

• 修正了iestimates与estimates的变量名混淆问题，解决了核心拟合错误。
• 统一使用dat_clean处理数据，避免原函数中data与dat的逻辑混乱。
• 添加tryCatch错误处理，确保循环不会因个别代谢物拟合失败中断。
• 简化了频次统计代码，去掉不必要的变量清理操作。

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