基于多列中的直接和间接相似性对变量进行分组的快速方法

Posted 2023-03-05

【中文标题】基于多列中的直接和间接相似性对变量进行分组的快速方法

【英文标题】：Fast way to group variables based on direct and indirect similarities in multiple columns

【发布时间】：2019-11-06 12:30:57

【问题描述】：

我有一个相对较大的数据集（1,750,000 行，5 列），其中包含具有唯一 ID 值的记录（第一列），由四个标准（其他 4 列）描述。一个小例子是：

# example
library(data.table)
dt <- data.table(id=c("a1","b3","c7","d5","e3","f4","g2","h1","i9","j6"), 
                 s1=c("a","b","c","l","l","v","v","v",NA,NA), 
                 s2=c("d","d","e","k","k","o","o","o",NA,NA),
                 s3=c("f","g","f","n","n","s","r","u","w","z"),
                 s4=c("h","i","j","m","m","t","t","t",NA,NA))

看起来像这样：

   id   s1   s2 s3   s4
 1: a1    a    d  f    h
 2: b3    b    d  g    i
 3: c7    c    e  f    j
 4: d5    l    k  n    m
 5: e3    l    k  n    m
 6: f4    v    o  s    t
 7: g2    v    o  r    t
 8: h1    v    o  u    t
 9: i9 <NA> <NA>  w <NA>
10: j6 <NA> <NA>  z <NA>

我的最终目标是在 any 描述列（不考虑 NA）上找到所有具有相同字符的记录，并将它们分组到一个新的 ID 下，以便我可以轻松识别重复记录。这些 ID 是通过连接每一行的 ID 构成的。

事情变得更加混乱，因为我可以直接和间接地找到具有重复描述的那些记录。所以，我目前分两步做这个操作。

第 1 步 - 根据直接重复构造重复 ID

# grouping ids with duplicated info in any of the columns
#sorry, I could not find search for duplicates using multiple columns simultaneously...
dt[!is.na(dt$s1),ids1:= paste(id,collapse="|"), by = list(s1)]
dt[!is.na(dt$s1),ids2:= paste(id,collapse="|"), by = list(s2)]
dt[!is.na(dt$s1),ids3:= paste(id,collapse="|"), by = list(s3)]
dt[!is.na(dt$s1),ids4:= paste(id,collapse="|"), by = list(s4)]

# getting a unique duplicated ID for each row
dt$new.id <- apply(dt[,.(ids1,ids2,ids3,ids4)], 1, paste, collapse="|")
dt$new.id <- apply(dt[,"new.id",drop=FALSE], 1, function(x) paste(unique(strsplit(x,"\\|")[[1]]),collapse="|"))

此操作产生以下结果，唯一的重复 ID 定义为“new.id”：

   id   s1   s2 s3   s4     ids1     ids2  ids3     ids4   new.id
 1: a1    a    d  f    h       a1    a1|b3 a1|c7       a1 a1|b3|c7
 2: b3    b    d  g    i       b3    a1|b3    b3       b3    b3|a1
 3: c7    c    e  f    j       c7       c7 a1|c7       c7    c7|a1
 4: d5    l    k  n    m    d5|e3    d5|e3 d5|e3    d5|e3    d5|e3
 5: e3    l    k  n    m    d5|e3    d5|e3 d5|e3    d5|e3    d5|e3
 6: f4    v    o  s    t f4|g2|h1 f4|g2|h1    f4 f4|g2|h1 f4|g2|h1
 7: g2    v    o  r    t f4|g2|h1 f4|g2|h1    g2 f4|g2|h1 f4|g2|h1
 8: h1    v    o  u    t f4|g2|h1 f4|g2|h1    h1 f4|g2|h1 f4|g2|h1
 9: i9 <NA> <NA>  w <NA>     <NA>     <NA>  <NA>     <NA>       NA
10: j6 <NA> <NA>  z <NA>     <NA>     <NA>  <NA>     <NA>       NA

请注意，记录“b3”和“c7”是通过“a1”间接复制的（所有其他示例都是直接复制，应该保持不变）。这就是我们需要下一步的原因。

第 2 步 - 根据间接重复更新重复的 ID

#filtering the relevant columns for the indirect search
dt = dt[,.(id,new.id)]

#creating the patterns to be used by grepl() for the look-up for each row
dt[,patt:= .(paste(paste("^",id,"\\||",sep=""),paste("\\|",id,"\\||",sep=""),paste("\\|",id,"$",sep=""),collapse = "" ,sep="")), by = list(id)]

#Transforming the ID vector into factor and setting it as a 'key' to the data.table (speed up the processing)
dt$new.id = as.factor(dt$new.id)
setkeyv(dt, c("new.id"))

#Performing the loop using sapply
library(stringr)
for(i in 1:nrow(dt)) 
  pat = dt$patt[i] # retrieving the research pattern
  tmp = dt[new.id %like% pat] # searching the pattern using grepl()
  if(dim(tmp)[1]>1) 
    x = which.max(str_count(tmp$new.id, "\\|"))
    dt$new.id[i] = as.character(tmp$new.id[x])
  


#filtering the final columns 
dt = dt[,.(id,new.id)]

决赛桌长这样：

   id   new.id
 1: a1 a1|b3|c7
 2: b3 a1|b3|c7
 3: c7 a1|b3|c7
 4: d5    d5|e3
 5: e3    d5|e3
 6: f4 f4|g2|h1
 7: g2 f4|g2|h1
 8: h1 f4|g2|h1
 9: i9       NA
10: j6       NA

请注意，现在前三个记录（“a1”、“b3”、“c7”）被分组在一个更广泛的重复 ID 下，其中包含直接记录和间接记录。

一切正常，但我的代码非常慢。运行一半的数据集（~800,0000）花了整整 2 天的时间。我可以将循环并行化到不同的内核中，但这仍然需要几个小时。而且我几乎可以肯定我可以以更好的方式使用 data.table 功能，也许在循环中使用“set”。我今天花了几个小时尝试使用 data.table 实现相同的代码，但我对它的语法不熟悉，我在这里真的很难。关于如何优化此代码的任何建议？

注意：代码中最慢的部分是循环，而循环内部效率最低的步骤是 data.table 中模式的 grepl()。似乎为 data.table 设置一个“键”可以加快这个过程，但在我的例子中，我没有改变执行 grepl() 的时间。

【问题讨论】：

您是否考虑过矢量化操作？由于只有 26 个字母，您可以为每个字母/标准列创建一个变量。然后，您可以为每个在 s1 中具有 a 的观察创建一个 id。但是，我不清楚这些新 ID 会有什么帮助。

感谢您的评论@spazznolo。我正在使用此代码在数百个植物标本馆中搜索重复的植物标本。在真实数据中，id 和描述实际上是更大且可变的字符，并且有很多缺失的描述（这就是我寻找间接重复的原因）。我确实尝试过使用 'sapply' 来矢量化循环，但由于数据的大小，它仍然需要将近 2 天的计算时间。再次感谢！

描述是否标准化？您有正在评估的样本清单吗？

再次感谢@spazznolo！在真实数据中，我的 id 看起来像：“alcb_237”或“p_120356”。描述列是包含每个记录的连接描述的向量，它们看起来像：“Melastomataceae_jesus_34_sao sebastiao passe”或“Fabaceae_costa_1022_juazeiro”。所以，我不会说它们是标准化的。此外，在某些记录的某些描述列中存在 NA，而对于其他记录，它们是完整的。

为了正确理解，描述出现在哪一列有关系吗？例如，s2 列中的字母“k”是否与s4 列中出现的含义相同。例如，如果第一行 (id == "a1") 中 s4 列中的“h”变为“k”，那么预期结果是否为 c(rep('a1|b3|c7|d5|e3', 5), rep('f4|g2|h1', 3))？

【参考方案1】：

您可以将此视为网络问题。这里我使用了igraph 包中的函数。基本步骤：

melt数据转长格式。

使用graph_from_data_frame 创建图表，其中“id”和“value”列被视为边列表。

使用components 获取图形的连接组件，即通过其标准直接或间接连接哪些“id”。

选择membership元素获取“每个顶点所属的簇id”。

加入原始数据的成员。

连接按集群成员分组的“id”。

---

library(igraph)

# melt data to long format, remove NA values
d <- melt(dt, id.vars = "id", na.rm = TRUE)

# convert to graph
g <- graph_from_data_frame(d[ , .(id, value)])

# get components and their named membership id 
mem <- components(g)$membership

# add membership id to original data
dt[.(names(mem)), on = .(id), mem := mem] 

# for groups of length one, set 'mem' to NA
dt[dt[, .I[.N == 1], by = mem]$V1, mem := NA]

如果需要，将“id”与“mem”列连接（对于非NA“mem”）（恕我直言，这只会使进一步的数据操作更加困难；））。无论如何，我们开始吧：

dt[!is.na(mem), id2 := paste(id, collapse = "|"), by = mem]

#     id   s1   s2 s3   s4  mem      id2
#  1: a1    a    d  f    h    1 a1|b3|c7
#  2: b3    b    d  g    i    1 a1|b3|c7
#  3: c7    c    e  f    j    1 a1|b3|c7
#  4: d5    l    k  l    m    2    d5|e3
#  5: e3    l    k  l    m    2    d5|e3
#  6: f4    o    o  s    o    3 f4|g2|h1
#  7: g2    o    o  r    o    3 f4|g2|h1
#  8: h1    o    o  u    o    3 f4|g2|h1
#  9: i9 <NA> <NA>  w <NA>   NA     <NA>
# 10: j6 <NA> <NA>  z <NA>   NA     <NA>

---

这个小例子中的基本图，只是为了说明连接的组件：

plot(g, edge.arrow.size = 0.5, edge.arrow.width = 0.8, vertex.label.cex = 2, edge.curved = FALSE)

【讨论】：

亲爱的@Henrik，感谢您的回答。但与上一个答案类似，我省略了数据集在 sp1:4 中有很多 NA，不应视为相等。我已经编辑了我的问题以包含带有 NA 的示例。我正在尝试在这里调整您的想法。如果我能做到，我会告诉你的。再次感谢。

你认为简单的 'na.omit(d,cols = "value")' 会是一个快速的解决方案吗？

哦，我喜欢这个，应该比重复连接更快。

感激不尽，@Henrik 和 Alexis。我刚刚用我的整个数据集（1,758,253 行）测试了这两种解决方案。两者都工作得很好，但是 Henrik 解决方案大约需要 30 秒，而 Alexis 解决方案大约需要 3 分钟！从 2 天到不到 3 分钟的计算时间简直太神奇了！ Henrik 感谢额外的连接代码，因为它们将用于集合之间的存储和检索目的。它终于是代码中最慢的部分，大约占总时间的 20 秒，但无论如何它都非常快。再次感谢大家！

也感谢@R.Lima 发布了一个很好的第一个问题，其中包含玩具数据和您尝试过的代码。干杯

【参考方案2】：

我认为这种递归方法可以满足您的需求。 基本上，它对每一列执行自连接， 一次一个， 如果匹配多于一行 （即考虑的行以外的行）， 它保存匹配中的所有唯一 ID。 它通过利用secondary indices 避免使用带有NA 的行。 诀窍是我们进行两次递归， 一次使用ids，再次使用新创建的new_ids。

dt[, new_id := .(list(character()))]

get_ids <- function(matched_ids, new_id) 
  if (length(matched_ids) > 1L) 
    list(unique(
      c(new_id[[1L]], unlist(matched_ids))
    ))
   else 
    new_id
  


find_recursively <- function(dt, cols, pass) 
  if (length(cols) == 0L) return(invisible())

  current <- cols[1L]
  next_cols <- cols[-1L]

  next_dt <- switch(
    pass,

    first = dt[!list(NA_character_),
               new_id := dt[.SD, .(get_ids(x.id, i.new_id)), on = current, by = .EACHI]$V1,
               on = current],

    second = dt[!list(NA_character_),
                new_id := dt[.SD, .(get_ids(x.new_id, i.new_id)), on = current, by = .EACHI]$V1,
                on = current]
  )

  find_recursively(next_dt, next_cols, pass)


find_recursively(dt, paste0("s", 1:4), "first")
find_recursively(dt, paste0("s", 1:4), "second")

dt[, new_id := sapply(new_id, function(nid) 
  ids <- unlist(nid)
  if (length(ids) == 0L) 
    NA_character_
   else 
    paste(ids, collapse = "|")
  
)]

print(dt)
    id   s1   s2 s3   s4   new_id
 1: a1    a    d  f    h a1|b3|c7
 2: b3    b    d  g    i a1|b3|c7
 3: c7    c    e  f    j a1|c7|b3
 4: d5    l    k  l    m    d5|e3
 5: e3    l    k  l    m    d5|e3
 6: f4    o    o  s    o f4|g2|h1
 7: g2    o    o  r    o f4|g2|h1
 8: h1    o    o  u    o f4|g2|h1
 9: i9 <NA> <NA>  w <NA>     <NA>
10: j6 <NA> <NA>  z <NA>     <NA>

连接使用this idiom。

【讨论】：

我觉得这可以改进，不确定它是否足够快。我会考虑的。

感谢一百万 @Alexis。我进行了一个测试，看看它在原始数据样本上的表现如何。

亲爱的@Alexis，我正要测试你提供的代码的性能，但后来我意识到我遗漏了一个非常重要的事情：数据集在 sp1:4 中有很多 NA，这不应该被视为平等。我已经编辑了我的问题以包含带有 NA 的示例，但我无法调整您的代码来做同样的事情。我应该在您的代码中声明类似于 !is.na(sp1:4) 的任何想法？再次感谢。

@R.Lima 我已经修改了代码，无论如何我使用的成语是错误的。我仍然不确定这是否是最好的方法，但一定要试一试并告诉我。

@R.Lima 出于好奇，你为什么不能使用 Henrik 的解决方案？