R语言嵌套循环高效替代：Wilcoxon检验加速实现方案问询

问题说明 #

• 需执行Wilcoxon检验，现有两个数据框列表：
  • Datalist：存储两年内不同观测的数量数据
  • Varlist：包含不同场景下的Case（病例）日与Control（对照）日标记
• 目标：检验各观测数量与不同病例/对照日场景之间的差异，使用wilcox.test()函数
• 期望输出：包含Case和Control的均值、p值，以及对应列表名、列名的数据框，用于结果匹配
• 当前问题：已有四层嵌套循环的可行代码，但运行速度极慢（预计至少10天），急需更高效的实现方式

示例数据 #

set.seed(42) 
n <- 365
df1 = data.frame(Date=seq.Date(as.Date("2017-01-01"), as.Date("2018-12-31"), "day"),
                  D1 = sample(18:30, n, replace=TRUE),
                  D2 = sample(0:7, n, replace=TRUE),
                  D3 = sample(0:10, n, replace=TRUE),
                  D4 = sample(0:4, n, replace=TRUE),
                  D5 = sample(0:23, n, replace=TRUE))
set.seed(7) 
n <- 365
df2 = data.frame(Date=seq.Date(as.Date("2017-01-01"), as.Date("2018-12-31"), "day"),
                 D1 = sample(18:30, n, replace=TRUE),
                 D2 = sample(0:7, n, replace=TRUE),
                 D3 = sample(0:10, n, replace=TRUE),
                 D4 = sample(0:4, n, replace=TRUE),
                 D5 = sample(0:23, n, replace=TRUE))

set.seed(9) 
n <- 365
df3 = data.frame(Date=seq.Date(as.Date("2017-01-01"), as.Date("2018-12-31"), "day"),
                 D1 = sample(18:30, n, replace=TRUE),
                 D2 = sample(0:7, n, replace=TRUE),
                 D3 = sample(0:10, n, replace=TRUE),
                 D4 = sample(0:4, n, replace=TRUE),
                 D5 = sample(0:23, n, replace=TRUE))

Datalist = list(df1, df2, df3)

set.seed(2) 
n <- 365
Var1 = data.frame(Date=seq.Date(as.Date("2017-01-01"), as.Date("2018-12-31"), "day"),
                 V1 = sample(c("Case", "Control", NA), n, replace=TRUE),
                 V2 = sample(c(NA, "Case", "Control"), n, replace=TRUE),
                 V3 = sample(c("Control", "Case", NA), n, replace=TRUE))

set.seed(6) 
n <- 365
Var2 = data.frame(Date=seq.Date(as.Date("2017-01-01"), as.Date("2018-12-31"), "day"),
                  V1 = sample(c("Case", "Control", NA), n, replace=TRUE),
                  V2 = sample(c(NA, "Case", "Control"), n, replace=TRUE),
                  V3 = sample(c("Control", "Case", NA), n, replace=TRUE))

set.seed(23) 
n <- 365
Var3 = data.frame(Date=seq.Date(as.Date("2017-01-01"), as.Date("2018-12-31"), "day"),
                  V1 = sample(c("Case", "Control", NA), n, replace=TRUE),
                  V2 = sample(c(NA, "Case", "Control"), n, replace=TRUE),
                  V3 = sample(c("Control", "Case", NA), n, replace=TRUE))

Varlist = list(Var1, Var2, Var3) 

原低效代码 #

Results = data.frame(matrix(ncol = 7, nrow = 0))
colnames(Results) = c("Code","ICD", "Cond", "Case", "Control", "pValue", "Ver")

for (a in 1:length(Datalist)) {
  print(names(Datalist)[a])
  for (b in 2:length(Datalist[[a]])) {
    for (c in 1:length(Varlist)) {
      for (d in 2:ncol(Varlist[[c]])){ 
        Ill = Datalist[[a]][,b]
        cutpoint = nrow(Datalist[[a]])
        Group = Varlist[[c]][,d]
        Group = Group[1:cutpoint]
        casecontrol = na.omit(data.frame(Ill, Group)) 
        wiltest = wilcox.test(casecontrol$Ill ~ casecontrol$Group)  
        stats = tapply(casecontrol$Ill,casecontrol$Group,mean) 
        
        Code = names(Datalist)[a]
        ICD = colnames(Datalist[[a]])[b]
        Cond = colnames(Varlist[[c]])[d]
        Case = round(stats[1],2)
        Control  = round(stats[2],2)
        pValue  = round(wiltest$p.value, 2)
        Ver  = names(Varlist)[c]
        
        addrow = c(Code, ICD, Case, Control, pValue, Ver)
        
        Results= rbind(Results,addrow)}}}}

代码慢因分析 #

1. 四层嵌套循环：R基础循环本身开销大，多次迭代累积耗时显著
2. 动态扩展数据框：每次用rbind添加行，会重复复制整个数据框，内存开销极大
3. 重复操作：每次循环都要做NA过滤、数据拼接，没有批量处理

优化方案 #

方案1：用tidyverse实现批量处理 #

利用长格式数据和分组操作替代循环，避免重复开销：

library(tidyverse)

# 给列表命名，方便标记来源
names(Datalist) <- paste0("df", 1:length(Datalist))
names(Varlist) <- paste0("Var", 1:length(Varlist))

# 将Datalist转换为长格式：Code(列表名)、Date、ICD(观测列名)、Value(观测值)
data_long <- imap_dfr(Datalist, ~ .x %>%
                        pivot_longer(cols = -Date, names_to = "ICD", values_to = "Value") %>%
                        mutate(Code = .y))

# 将Varlist转换为长格式：Ver(列表名)、Date、Cond(场景列名)、Group(Case/Control)
var_long <- imap_dfr(Varlist, ~ .x %>%
                       pivot_longer(cols = -Date, names_to = "Cond", values_to = "Group") %>%
                       mutate(Ver = .y))

# 按Date合并数据，过滤无效值
combined <- inner_join(data_long, var_long, by = "Date") %>%
  filter(!is.na(Group), Group %in% c("Case", "Control"))

# 分组执行检验和计算均值
results <- combined %>%
  group_by(Code, ICD, Cond, Ver) %>%
  summarise(
    Case = round(mean(Value[Group == "Case"]), 2),
    Control = round(mean(Value[Group == "Control"]), 2),
    pValue = round(wilcox.test(Value ~ Group)$p.value, 2),
    .groups = "drop"
  )

# 调整列顺序为目标格式
results <- results %>%
  select(Code, ICD, Cond, Case, Control, pValue, Ver)

方案2：用data.table实现（超大数据集推荐） #

data.table的分组操作基于C实现，内存管理更高效，适合大数据量：

library(data.table)

# 给列表命名
names(Datalist) <- paste0("df", 1:length(Datalist))
names(Varlist) <- paste0("Var", 1:length(Varlist))

# 转换Datalist为data.table长格式
data_dt <- rbindlist(Datalist, idcol = "Code") %>%
  melt(id.vars = c("Code", "Date"), variable.name = "ICD", value.name = "Value")

# 转换Varlist为data.table长格式，过滤无效值
var_dt <- rbindlist(Varlist, idcol = "Ver") %>%
  melt(id.vars = c("Ver", "Date"), variable.name = "Cond", value.name = "Group") %>%
  filter(!is.na(Group), Group %in% c("Case", "Control"))

# 合并数据
combined_dt <- merge(data_dt, var_dt, by = "Date")

# 分组计算统计量和p值
results_dt <- combined_dt[, .(
  Case = round(mean(Value[Group == "Case"]), 2),
  Control = round(mean(Value[Group == "Control"]), 2),
  pValue = round(wilcox.test(Value ~ Group)$p.value, 2)
), by = .(Code, ICD, Cond, Ver)]

# 调整列顺序
setcolorder(results_dt, c("Code", "ICD", "Cond", "Case", "Control", "pValue", "Ver"))

方案3：并行计算加速（多核CPU利用） #

如果机器有多个CPU核心，可使用furrr实现并行处理，进一步缩短时间：

library(tidyverse)
library(furrr)

# 开启多线程并行
plan(multisession)

# 重复方案1的前几步，得到combined数据框
names(Datalist) <- paste0("df", 1:length(Datalist))
names(Varlist) <- paste0("Var", 1:length(Varlist))

data_long <- imap_dfr(Datalist, ~ .x %>%
                        pivot_longer(cols = -Date, names_to = "ICD", values_to = "Value") %>%
                        mutate(Code = .y))

var_long <- imap_dfr(Varlist, ~ .x %>%
                       pivot_longer(cols = -Date, names_to = "Cond", values_to = "Group") %>%
                       mutate(Ver = .y))

combined <- inner_join(data_long, var_long, by = "Date") %>%
  filter(!is.na(Group), Group %in% c("Case", "Control"))

# 分组并行计算
results_parallel <- combined %>%
  group_nest(Code, ICD, Cond, Ver) %>%
  mutate(
    stats = future_map(data, ~ {
      list(
        Case = round(mean(.x$Value[.x$Group == "Case"]), 2),
        Control = round(mean(.x$Value[.x$Group == "Control"]), 2),
        pValue = round(wilcox.test(.x$Value ~ .x$Group)$p.value, 2)
      )
    })
  ) %>%
  unnest_wider(stats) %>%
  select(-data) %>%
  select(Code, ICD, Cond, Case, Control, pValue, Ver)

额外提速技巧 #

• 提前过滤无效分组：如果某个分组的Case或Control样本量过小（比如<5），可直接跳过检验，避免无效计算
• 减少精度：如果不需要保留两位小数，可适当减少，减少计算开销
• 内存优化：确保数据类型正确（比如把字符型Group转为因子型），减少内存占用

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