如何从一个巨大的矩阵中以尽可能少的行 NA 获得尽可能大的列序列？
How to get the largest possible column sequence with the least possible row NAs from a huge matrix?
我想从数据框中提取 select 列，以便生成的 连续 列序列尽可能长，而带有 NA 的行数是尽可能小，因为之后必须删除它们。
（我想这样做的原因是，我想 运行 TraMineR::seqsubm() 自动获取转换成本矩阵（按转换概率），然后 运行 cluster::agnes() 在上面。TraMineR::seqsubm() 不喜欢 NA 状态并且 cluster::agnes() 矩阵中有 NA 状态不一定有多大意义。）
为此，我已经写了一个工作 function，它原则上计算所有可能的幂子集并检查它们是否有 NA。它适用于代表 10x5 矩阵的玩具数据 d：
> d
id X1 X2 X3 X4 X5
1 A 1 11 21 31 41
2 B 2 12 22 32 42
3 C 3 13 23 33 NA
4 D 4 14 24 34 NA
5 E 5 15 25 NA NA
6 F 6 16 26 NA NA
7 G 7 17 NA NA NA
8 H 8 18 NA NA NA
9 I 9 NA NA NA NA
10 J 10 NA NA NA NA
11 K NA NA NA NA NA
现在的问题是我实际上想将该算法应用于代表 34235 x 17 矩阵的调查数据！
我的代码已经在Code Review上审核过了，但是仍然不能应用于真实数据。
我知道使用这种方法会有大量的计算。 （对于非超级计算机来说可能太大了？！）
有谁知道更合适的方法吗？
我给你看已经 enhanced function by @minem 来自 Code Review：
seqRank2 <- function(d, id = "id") {
require(matrixStats)
# change structure, convert to matrix
ii <- as.character(d[, id])
dm <- d
dm[[id]] <- NULL
dm <- as.matrix(dm)
rownames(dm) <- ii
your.powerset = function(s){
l = vector(mode = "list", length = 2^length(s))
l[[1]] = numeric()
counter = 1L
for (x in 1L:length(s)) {
for (subset in 1L:counter) {
counter = counter + 1L
l[[counter]] = c(l[[subset]], s[x])
}
}
return(l[-1])
}
psr <- your.powerset(ii)
psc <- your.powerset(colnames(dm))
sss <- lapply(psr, function(x) {
i <- ii %in% x
lapply(psc, function(y) dm[i, y, drop = F])
})
cn <- sapply(sss, function(x)
lapply(x, function(y) {
if (ncol(y) == 1) {
if (any(is.na(y))) return(NULL)
return(y)
}
isna2 <- matrixStats::colAnyNAs(y)
if (all(isna2)) return(NULL)
if (sum(isna2) == 0) return(NA)
r <- y[, !isna2, drop = F]
return(r)
}))
scr <- sapply(cn, nrow)
scc <- sapply(cn, ncol)
namesCN <- sapply(cn, function(x) paste0(colnames(x), collapse = ", "))
names(scr) <- namesCN
scr <- unlist(scr)
names(scc) <- namesCN
scc <- unlist(scc)
m <- t(rbind(n.obs = scr, sq.len = scc))
ag <- aggregate(m, by = list(sequence = rownames(m)), max)
ag <- ag[order(-ag$sq.len, -ag$n.obs), ]
rownames(ag) <- NULL
return(ag)
}
产量：
> seqRank2(d)
sequence n.obs sq.len
1 X1, X2, X3, X4 4 4
2 X1, X2, X3 6 3
3 X1, X2, X4 4 3
4 X1, X3, X4 4 3
5 X2, X3, X4 4 3
6 X1, X2 8 2
7 X1, X3 6 2
8 X2, X3 6 2
9 X1, X4 4 2
10 X2, X4 4 2
11 X3, X4 4 2
12 X1 10 1
13 X2 8 1
14 X3 6 1
15 X4 4 1
16 X5 2 1
> system.time(x <- seqRank2(d))
user system elapsed
1.93 0.14 2.93
在这种情况下，我会选择 X1, X2, X3, X4、X1, X2, X3 或 X2, X3, X4，因为它们 连续 并产生适当数量的观察结果.
预期输出：
所以对于玩具数据 d 预期的输出类似于：
> seqRank2(d)
sequence n.obs sq.len
1 X1, X2, X3, X4 4 4
2 X1, X2, X3 6 3
3 X2, X3, X4 4 3
4 X1, X2 8 2
5 X2, X3 6 2
6 X3, X4 4 2
7 X1 10 1
8 X2 8 1
9 X3 6 1
10 X4 4 1
11 X5 2 1
最后函数应该 运行 在巨大的矩阵 d.huge 上正确地导致错误：
> seqRank2(d.huge)
Error in vector(mode = "list", length = 2^length(s)) :
vector size cannot be infinite
玩具资料d:
d <- structure(list(id = structure(1:11, .Label = c("A", "B", "C",
"D", "E", "F", "G", "H", "I", "J", "K"), class = "factor"), X1 = c(1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, NA), X2 = c(11L, 12L, 13L,
14L, 15L, 16L, 17L, 18L, NA, NA, NA), X3 = c(21L, 22L, 23L, 24L,
25L, 26L, NA, NA, NA, NA, NA), X4 = c(31L, 32L, 33L, 34L, NA,
NA, NA, NA, NA, NA, NA), X5 = c(41L, 42L, NA, NA, NA, NA, NA,
NA, NA, NA, NA)), row.names = c(NA, -11L), class = "data.frame")
玩具资料d.huge:
d.huge <- setNames(data.frame(matrix(1:15.3e5, 3e4, 51)),
c("id", paste0("X", 1:50)))
d.huge[, 41:51] <- lapply(d.huge[, 41:51], function(x){
x[which(x %in% sample(x, .05*length(x)))] <- NA
x
})
附录（见评论最新回复）：
d.huge <- read.csv("d.huge.csv")
d.huge.1 <- d.huge[sample(nrow(d.huge), 3/4*nrow(d.huge)), ]
d1 <- seqRank3(d.huge.1, 1.27e-1, 1.780e1)
d2 <- d1[complete.cases(d1), ]
dim(d2)
names(d2)
转换为矩阵并计算每列的 Na 计数：
dm <- is.na(d[, -1])
na_counts <- colSums(dm)
x <- data.frame(na_counts = na_counts, non_na_count = nrow(dm) - na_counts)
x <- as.matrix(x)
# create all combinations for column indexes:
nx <- 1:nrow(x)
psr <- do.call(c, lapply(seq_along(nx), combn, x = nx, simplify = FALSE))
# test if continuous:
good <- sapply(psr, function(y) !any(diff(sort.int(y)) != 1L))
psr <- psr[good == T] # remove non continuous
# for each combo count nas and non NA:
s <- sapply(psr, function(y) colSums(x[y, , drop = F]))
# put all together in table:
res <- data.frame(var_count = lengths(psr), t(s))
res$var_indexes <- sapply(psr, paste, collapse = ',')
res
# var_count na_counts non_na_count var_indexes
# 1 1 1 10 1
# 2 1 3 8 2
# 3 1 5 6 3
# 4 1 7 4 4
# 5 1 9 2 5
# 6 2 4 18 1,2
# 7 2 8 14 2,3
# 8 2 12 10 3,4
# 9 2 16 6 4,5
# 10 3 9 24 1,2,3
# 11 3 15 18 2,3,4
# 12 3 21 12 3,4,5
# 13 4 16 28 1,2,3,4
# 14 4 24 20 2,3,4,5
# 15 5 25 30 1,2,3,4,5
# choose
由于 var 索引是排序的，为了速度我们可以简单地使用：
good <- sapply(psr, function(y) !any(diff(y) != 1L))
海量数据用时不到一秒
l1 = combn(2:length(d), 2, function(x) d[x[1]:x[2]], simplify = FALSE)
# If you also need "combinations" of only single columns, then uncomment the next line
# l1 = c(d[-1], l1)
l2 = sapply(l1, function(x) sum(complete.cases(x)))
score = sapply(1:length(l1), function(i) NCOL(l1[[i]]) * l2[i])
best_score = which.max(score)
best = l1[[best_score]]
问题不清楚如何对各种组合进行排名。我们可以使用不同的评分公式来产生不同的偏好。例如，要分别对行数和列数进行加权，我们可以这样做
col_weight = 2
row_weight = 1
score = sapply(1:length(l1), function(i) col_weight*NCOL(l1[[i]]) + row_weight * l2[i])
澄清一下，TraMineR 中的 seqsubm 函数对于 NA 和不同长度的序列都没有任何问题。但是，该函数需要一个状态序列对象（使用 seqdef 创建）作为输入。
函数seqsubm用于通过不同的方法计算状态之间的替代成本（即差异）。您可能指的是从观察到的转移概率中得出成本的方法 ('TRATE')，即 2-p(i|j) - p (j|i)，其中 p(i|j) 是在 [=29] 中处于状态 i 的概率=]t 当我们在 t-1 中处于状态 j 时。所以，我们所需要的只是转移概率，它可以很容易地从一组不同长度的序列或其中有间隙的序列中估计出来。
我在下面使用 TraMineR 附带的 ex1 数据进行说明。 （由于您的玩具示例中有大量不同的状态，因此生成的替代成本矩阵对于此说明来说太大 (28 x 28)。）
library(TraMineR)
data(ex1)
sum(is.na(ex1))
# [1] 38
sq <- seqdef(ex1[1:13])
sq
# Sequence
# s1 *-*-*-A-A-A-A-A-A-A-A-A-A
# s2 D-D-D-B-B-B-B-B-B-B
# s3 *-D-D-D-D-D-D-D-D-D-D
# s4 A-A-*-*-B-B-B-B-D-D
# s5 A-*-A-A-A-A-*-A-A-A
# s6 *-*-*-C-C-C-C-C-C-C
# s7 *-*-*-*-*-*-*-*-*-*-*-*-*
sm <- seqsubm(sq, method='TRATE')
round(sm,digits=3)
# A-> B-> C-> D->
# A-> 0 2.000 2 2.000
# B-> 2 0.000 2 1.823
# C-> 2 2.000 0 2.000
# D-> 2 1.823 2 0.000
现在，我不清楚您想对状态差异做什么。在聚类算法中输入它们，您将对状态进行聚类。如果你想对序列进行聚类，那么你应该首先计算序列之间的差异（使用 seqdist 并可能将 seqsubm 返回的替代成本矩阵作为 sm 参数传递）然后输入聚类算法中的结果距离矩阵。
我想从数据框中提取 select 列，以便生成的 连续 列序列尽可能长，而带有 NA 的行数是尽可能小，因为之后必须删除它们。
（我想这样做的原因是，我想 运行 TraMineR::seqsubm() 自动获取转换成本矩阵（按转换概率），然后 运行 cluster::agnes() 在上面。TraMineR::seqsubm() 不喜欢 NA 状态并且 cluster::agnes() 矩阵中有 NA 状态不一定有多大意义。）
为此，我已经写了一个工作 function，它原则上计算所有可能的幂子集并检查它们是否有 NA。它适用于代表 10x5 矩阵的玩具数据 d：
> d
id X1 X2 X3 X4 X5
1 A 1 11 21 31 41
2 B 2 12 22 32 42
3 C 3 13 23 33 NA
4 D 4 14 24 34 NA
5 E 5 15 25 NA NA
6 F 6 16 26 NA NA
7 G 7 17 NA NA NA
8 H 8 18 NA NA NA
9 I 9 NA NA NA NA
10 J 10 NA NA NA NA
11 K NA NA NA NA NA
现在的问题是我实际上想将该算法应用于代表 34235 x 17 矩阵的调查数据！
我的代码已经在Code Review上审核过了，但是仍然不能应用于真实数据。
我知道使用这种方法会有大量的计算。 （对于非超级计算机来说可能太大了？！）
有谁知道更合适的方法吗？
我给你看已经 enhanced function by @minem 来自 Code Review：
seqRank2 <- function(d, id = "id") {
require(matrixStats)
# change structure, convert to matrix
ii <- as.character(d[, id])
dm <- d
dm[[id]] <- NULL
dm <- as.matrix(dm)
rownames(dm) <- ii
your.powerset = function(s){
l = vector(mode = "list", length = 2^length(s))
l[[1]] = numeric()
counter = 1L
for (x in 1L:length(s)) {
for (subset in 1L:counter) {
counter = counter + 1L
l[[counter]] = c(l[[subset]], s[x])
}
}
return(l[-1])
}
psr <- your.powerset(ii)
psc <- your.powerset(colnames(dm))
sss <- lapply(psr, function(x) {
i <- ii %in% x
lapply(psc, function(y) dm[i, y, drop = F])
})
cn <- sapply(sss, function(x)
lapply(x, function(y) {
if (ncol(y) == 1) {
if (any(is.na(y))) return(NULL)
return(y)
}
isna2 <- matrixStats::colAnyNAs(y)
if (all(isna2)) return(NULL)
if (sum(isna2) == 0) return(NA)
r <- y[, !isna2, drop = F]
return(r)
}))
scr <- sapply(cn, nrow)
scc <- sapply(cn, ncol)
namesCN <- sapply(cn, function(x) paste0(colnames(x), collapse = ", "))
names(scr) <- namesCN
scr <- unlist(scr)
names(scc) <- namesCN
scc <- unlist(scc)
m <- t(rbind(n.obs = scr, sq.len = scc))
ag <- aggregate(m, by = list(sequence = rownames(m)), max)
ag <- ag[order(-ag$sq.len, -ag$n.obs), ]
rownames(ag) <- NULL
return(ag)
}
产量：
> seqRank2(d)
sequence n.obs sq.len
1 X1, X2, X3, X4 4 4
2 X1, X2, X3 6 3
3 X1, X2, X4 4 3
4 X1, X3, X4 4 3
5 X2, X3, X4 4 3
6 X1, X2 8 2
7 X1, X3 6 2
8 X2, X3 6 2
9 X1, X4 4 2
10 X2, X4 4 2
11 X3, X4 4 2
12 X1 10 1
13 X2 8 1
14 X3 6 1
15 X4 4 1
16 X5 2 1
> system.time(x <- seqRank2(d))
user system elapsed
1.93 0.14 2.93
在这种情况下，我会选择 X1, X2, X3, X4、X1, X2, X3 或 X2, X3, X4，因为它们 连续 并产生适当数量的观察结果.
预期输出：
所以对于玩具数据 d 预期的输出类似于：
> seqRank2(d)
sequence n.obs sq.len
1 X1, X2, X3, X4 4 4
2 X1, X2, X3 6 3
3 X2, X3, X4 4 3
4 X1, X2 8 2
5 X2, X3 6 2
6 X3, X4 4 2
7 X1 10 1
8 X2 8 1
9 X3 6 1
10 X4 4 1
11 X5 2 1
最后函数应该 运行 在巨大的矩阵 d.huge 上正确地导致错误：
> seqRank2(d.huge)
Error in vector(mode = "list", length = 2^length(s)) :
vector size cannot be infinite
玩具资料d:
d <- structure(list(id = structure(1:11, .Label = c("A", "B", "C",
"D", "E", "F", "G", "H", "I", "J", "K"), class = "factor"), X1 = c(1L,
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, NA), X2 = c(11L, 12L, 13L,
14L, 15L, 16L, 17L, 18L, NA, NA, NA), X3 = c(21L, 22L, 23L, 24L,
25L, 26L, NA, NA, NA, NA, NA), X4 = c(31L, 32L, 33L, 34L, NA,
NA, NA, NA, NA, NA, NA), X5 = c(41L, 42L, NA, NA, NA, NA, NA,
NA, NA, NA, NA)), row.names = c(NA, -11L), class = "data.frame")
玩具资料d.huge:
d.huge <- setNames(data.frame(matrix(1:15.3e5, 3e4, 51)),
c("id", paste0("X", 1:50)))
d.huge[, 41:51] <- lapply(d.huge[, 41:51], function(x){
x[which(x %in% sample(x, .05*length(x)))] <- NA
x
})
附录（见评论最新回复）：
d.huge <- read.csv("d.huge.csv")
d.huge.1 <- d.huge[sample(nrow(d.huge), 3/4*nrow(d.huge)), ]
d1 <- seqRank3(d.huge.1, 1.27e-1, 1.780e1)
d2 <- d1[complete.cases(d1), ]
dim(d2)
names(d2)
转换为矩阵并计算每列的 Na 计数：
dm <- is.na(d[, -1])
na_counts <- colSums(dm)
x <- data.frame(na_counts = na_counts, non_na_count = nrow(dm) - na_counts)
x <- as.matrix(x)
# create all combinations for column indexes:
nx <- 1:nrow(x)
psr <- do.call(c, lapply(seq_along(nx), combn, x = nx, simplify = FALSE))
# test if continuous:
good <- sapply(psr, function(y) !any(diff(sort.int(y)) != 1L))
psr <- psr[good == T] # remove non continuous
# for each combo count nas and non NA:
s <- sapply(psr, function(y) colSums(x[y, , drop = F]))
# put all together in table:
res <- data.frame(var_count = lengths(psr), t(s))
res$var_indexes <- sapply(psr, paste, collapse = ',')
res
# var_count na_counts non_na_count var_indexes
# 1 1 1 10 1
# 2 1 3 8 2
# 3 1 5 6 3
# 4 1 7 4 4
# 5 1 9 2 5
# 6 2 4 18 1,2
# 7 2 8 14 2,3
# 8 2 12 10 3,4
# 9 2 16 6 4,5
# 10 3 9 24 1,2,3
# 11 3 15 18 2,3,4
# 12 3 21 12 3,4,5
# 13 4 16 28 1,2,3,4
# 14 4 24 20 2,3,4,5
# 15 5 25 30 1,2,3,4,5
# choose
由于 var 索引是排序的，为了速度我们可以简单地使用：
good <- sapply(psr, function(y) !any(diff(y) != 1L))
海量数据用时不到一秒
l1 = combn(2:length(d), 2, function(x) d[x[1]:x[2]], simplify = FALSE)
# If you also need "combinations" of only single columns, then uncomment the next line
# l1 = c(d[-1], l1)
l2 = sapply(l1, function(x) sum(complete.cases(x)))
score = sapply(1:length(l1), function(i) NCOL(l1[[i]]) * l2[i])
best_score = which.max(score)
best = l1[[best_score]]
问题不清楚如何对各种组合进行排名。我们可以使用不同的评分公式来产生不同的偏好。例如，要分别对行数和列数进行加权，我们可以这样做
col_weight = 2
row_weight = 1
score = sapply(1:length(l1), function(i) col_weight*NCOL(l1[[i]]) + row_weight * l2[i])
澄清一下，TraMineR 中的 seqsubm 函数对于 NA 和不同长度的序列都没有任何问题。但是，该函数需要一个状态序列对象（使用 seqdef 创建）作为输入。
函数seqsubm用于通过不同的方法计算状态之间的替代成本（即差异）。您可能指的是从观察到的转移概率中得出成本的方法 ('TRATE')，即 2-p(i|j) - p (j|i)，其中 p(i|j) 是在 [=29] 中处于状态 i 的概率=]t 当我们在 t-1 中处于状态 j 时。所以，我们所需要的只是转移概率，它可以很容易地从一组不同长度的序列或其中有间隙的序列中估计出来。
我在下面使用 TraMineR 附带的 ex1 数据进行说明。 （由于您的玩具示例中有大量不同的状态，因此生成的替代成本矩阵对于此说明来说太大 (28 x 28)。）
library(TraMineR)
data(ex1)
sum(is.na(ex1))
# [1] 38
sq <- seqdef(ex1[1:13])
sq
# Sequence
# s1 *-*-*-A-A-A-A-A-A-A-A-A-A
# s2 D-D-D-B-B-B-B-B-B-B
# s3 *-D-D-D-D-D-D-D-D-D-D
# s4 A-A-*-*-B-B-B-B-D-D
# s5 A-*-A-A-A-A-*-A-A-A
# s6 *-*-*-C-C-C-C-C-C-C
# s7 *-*-*-*-*-*-*-*-*-*-*-*-*
sm <- seqsubm(sq, method='TRATE')
round(sm,digits=3)
# A-> B-> C-> D->
# A-> 0 2.000 2 2.000
# B-> 2 0.000 2 1.823
# C-> 2 2.000 0 2.000
# D-> 2 1.823 2 0.000
现在，我不清楚您想对状态差异做什么。在聚类算法中输入它们，您将对状态进行聚类。如果你想对序列进行聚类，那么你应该首先计算序列之间的差异（使用 seqdist 并可能将 seqsubm 返回的替代成本矩阵作为 sm 参数传递）然后输入聚类算法中的结果距离矩阵。