从成对矩阵中，找到等于某个值的最大个体群体

Question

我有一个 39x39 的成对相关性矩阵，其中包含 39 个个体的所有成对组合的相关性值。我想找到完全不相关的最大个体组，即该组中所有成对相关性值都等于 0。

在 R 中是否有一种简单的方法可以做到这一点？

一个更简单的例子：

set.seed(420)

#Create the matrix
relatedness.matrix <- matrix(data = sample(x = c(0.5, 1, 0,0), size = 25, replace = TRUE), nrow = 5, ncol = 5)

# Matrix has the same upper and lower triangles
relatedness.matrix[upper.tri(relatedness.matrix)] <- relatedness.matrix[lower.tri(relatedness.matrix)]

# Add names for simplicity of reference
colnames(relatedness.matrix) <- letters[1:5]
rownames(relatedness.matrix) <- letters[1:5]

# Relatedness between the same individual does not count
diag(relatedness.matrix) <- NA

在这种情况下，存在三种可能的解决方案：仅包含 和 的 2x2 矩阵e、仅包含和 的b2x2 矩阵以及仅包含和 的2x2 矩阵。将任何其他个体添加到任何这些矩阵中都会添加一个相关个体。cdae

编辑：补充说上三角形和下三角形是相同的，并且此示例中实际上有多个 2x2 解决方案。

Answer 1

对称邻接矩阵中不相关个体的集合描述了一个独立的集合。我们可以用来igraph::largest_ivs找到最大的这样的集合。

我将使用一个实际上对称的更大矩阵。

set.seed(420)

m <- matrix(sample(0:2, 26^2, 1), 26, 26, 0, rep(list(letters), 2))
m[lower.tri(m)] <- t(m)[lower.tri(m)]
diag(m) <- NA

检查矩阵是否对称

isSymmetric(m)
#> [1] TRUE

m
#>    a  b  c  d  e  f  g  h  i  j  k  l  m  n  o  p  q  r  s  t  u  v  w  x  y  z
#> a NA  0  1  2  1  2  0  2  2  0  1  2  0  1  2  2  0  0  0  0  0  2  2  2  1  1
#> b  0 NA  0  1  2  2  0  2  0  2  0  2  2  2  1  2  2  1  1  0  1  2  1  2  0  1
#> c  1  0 NA  0  1  0  2  1  2  1  0  1  0  1  2  2  2  2  1  2  2  0  2  0  1  0
#> d  2  1  0 NA  2  2  2  2  2  2  1  1  0  1  2  1  2  2  1  2  1  0  1  0  2  1
#> e  1  2  1  2 NA  2  1  0  1  0  1  0  0  0  1  2  0  2  0  2  2  1  2  1  2  2
#> f  2  2  0  2  2 NA  2  2  2  1  1  2  1  2  0  2  0  2  2  0  1  1  0  2  2  2
#> g  0  0  2  2  1  2 NA  0  2  1  2  2  2  2  0  1  2  0  2  1  0  0  1  1  2  1
#> h  2  2  1  2  0  2  0 NA  2  2  1  0  2  2  1  0  1  1  1  1  2  1  1  1  1  2
#> i  2  0  2  2  1  2  2  2 NA  1  2  1  0  2  2  0  2  2  2  0  2  0  0  0  0  2
#> j  0  2  1  2  0  1  1  2  1 NA  1  1  2  2  0  0  1  1  2  2  2  1  0  0  2  2
#> k  1  0  0  1  1  1  2  1  2  1 NA  2  2  1  0  0  2  0  2  0  0  1  1  1  1  2
#> l  2  2  1  1  0  2  2  0  1  1  2 NA  1  1  2  0  2  2  1  2  1  0  0  2  1  1
#> m  0  2  0  0  0  1  2  2  0  2  2  1 NA  0  2  2  0  2  1  1  1  1  0  2  1  1
#> n  1  2  1  1  0  2  2  2  2  2  1  1  0 NA  1  0  1  2  1  2  0  1  0  1  1  2
#> o  2  1  2  2  1  0  0  1  2  0  0  2  2  1 NA  2  2  0  1  2  1  2  2  1  1  0
#> p  2  2  2  1  2  2  1  0  0  0  0  0  2  0  2 NA  2  2  2  1  0  2  0  0  1  2
#> q  0  2  2  2  0  0  2  1  2  1  2  2  0  1  2  2 NA  1  0  1  2  2  1  0  1  1
#> r  0  1  2  2  2  2  0  1  2  1  0  2  2  2  0  2  1 NA  1  1  2  1  2  2  2  1
#> s  0  1  1  1  0  2  2  1  2  2  2  1  1  1  1  2  0  1 NA  2  1  1  2  1  1  1
#> t  0  0  2  2  2  0  1  1  0  2  0  2  1  2  2  1  1  1  2 NA  0  0  1  2  2  0
#> u  0  1  2  1  2  1  0  2  2  2  0  1  1  0  1  0  2  2  1  0 NA  2  2  0  2  0
#> v  2  2  0  0  1  1  0  1  0  1  1  0  1  1  2  2  2  1  1  0  2 NA  2  0  1  1
#> w  2  1  2  1  2  0  1  1  0  0  1  0  0  0  2  0  1  2  2  1  2  2 NA  0  2  0
#> x  2  2  0  0  1  2  1  1  0  0  1  2  2  1  1  0  0  2  1  2  0  0  0 NA  1  2
#> y  1  0  1  2  2  2  2  1  0  2  1  1  1  1  1  1  1  2  1  2  2  1  2  1 NA  0
#> z  1  1  0  1  2  2  1  2  2  2  2  1  1  2  0  2  1  1  1  0  0  1  0  2  0 NA

获取最大的独立集：

library(igraph)

is <- largest_ivs(graph_from_adjacency_matrix(m, "undirected"))
is
#> [[1]]
#> + 4/26 vertices, named, from 272900e:
#> [1] i p w x
#> 
#> [[2]]
#> + 4/26 vertices, named, from 272900e:
#> [1] c d v x
#> 
#> [[3]]
#> + 4/26 vertices, named, from 272900e:
#> [1] j p w x

验证独立集之间不存在边。

lapply(is, \(i) m[i, i])
#> [[1]]
#>    i  p  w  x
#> i NA  0  0  0
#> p  0 NA  0  0
#> w  0  0 NA  0
#> x  0  0  0 NA
#> 
#> [[2]]
#>    c  d  v  x
#> c NA  0  0  0
#> d  0 NA  0  0
#> v  0  0 NA  0
#> x  0  0  0 NA
#> 
#> [[3]]
#>    j  p  w  x
#> j NA  0  0  0
#> p  0 NA  0  0
#> w  0  0 NA  0
#> x  0  0  0 NA

附带说明一下，由邻接矩阵 形成的图中的独立集m将与由 形成的图中的派系相同!m。有趣的是，对于我的小例子来说，比找到所需的答案largest_cliques更高效：largest_ivs

microbenchmark::microbenchmark(
  cliques = largest_cliques(graph_from_adjacency_matrix(!m, "undirected")),
  ivs = largest_ivs(graph_from_adjacency_matrix(m, "undirected"))
)
#> Unit: microseconds
#>     expr   min    lq    mean median     uq    max neval
#>  cliques 319.7 348.6 372.581 368.90 388.55  555.0   100
#>      ivs 560.8 589.6 629.992 616.55 654.35 1187.6   100

随着矩阵变大，性能差异也随之增大：

m <- matrix(sample(0:2, 1e4, 1), 100, 100, 0)
m[lower.tri(m)] <- t(m)[lower.tri(m)]
diag(m) <- NA

microbenchmark::microbenchmark(
  cliques = largest_cliques(graph_from_adjacency_matrix(!m, "undirected")),
  ivs = largest_ivs(graph_from_adjacency_matrix(m, "undirected"))
)
#> Unit: milliseconds
#>     expr      min       lq       mean   median       uq      max neval
#>  cliques   2.5735   2.7651   3.275977   2.9013   3.3138   7.9742   100
#>      ivs 161.9572 182.3812 191.595736 191.2344 202.1377 243.5654   100

m <- matrix(sample(0:2, 4e4, 1), 200, 200, 0)
m[lower.tri(m)] <- t(m)[lower.tri(m)]
diag(m) <- NA

system.time(cl <- largest_cliques(graph_from_adjacency_matrix(!m, "undirected")))
#>    user  system elapsed 
#>    0.05    0.00    0.05
system.time(is <- largest_ivs(graph_from_adjacency_matrix(m, "undirected")))
#>    user  system elapsed 
#>   10.14    0.00   10.15

检查答案是否相同。

library(data.table)

identical(
  setorder(as.data.table(t(mapply(sort, cl)))),
  setorder(as.data.table(t(mapply(sort, is))))
)
#> [1] TRUE