2

我想从数据框中选择列,以便生成的连续列序列尽可能长,而具有 NA 的行数尽可能少,因为之后必须删除它们。

(我想这样做的原因是,我想运行TraMineR::seqsubm()以自动获取转换成本矩阵(通过转换概率),然后cluster::agnes()在其上运行。TraMineR::seqsubm()不喜欢NA状态并且矩阵中cluster::agnes()NA状态不一定会使很有道理。)

为此,我已经编写了一个工作函数,它按原则计算所有可能的幂子集并检查它们是否为NAs。它适用于这个d代表 10x5 矩阵的玩具数据:

> 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 上进行了审查,但仍然无法应用于真实数据。

我知道使用这种方法会产生巨大的计算量。(对于非超级计算机来说可能太大了?!)

有谁知道更合适的方法?

我向您展示了来自 Code Review的 @minem 已经增强的功能:

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, X4X1, 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)
4

3 回答 3

1

转换为矩阵并计算每列的 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))
于 2018-12-04T14:31:59.667 回答
1

这需要不到一秒钟的海量数据

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])
于 2018-12-04T15:21:16.003 回答
1

澄清一下,seqsubm函数 fromTraMineR对 NA 完全没有问题,对于不同长度的序列也没有问题。但是,该函数需要一个状态序列对象(用 来创建seqdef)作为输入。

该功能seqsubm用于通过不同的方法计算状态之间的替代成本(即差异)。您可能会参考'TRATE'从观察到的转移概率得出成本的方法 ( ),即 2- p(i|j) - p(j|i),其中p(i|j)是处于状态的概率i in t当我们在t -1中处于状态j时。所以,我们所需要的只是转移概率,它可以很容易地从一组不同长度的序列或其中有间隙的序列中估计出来。

我在下面使用ex1随附的数据进行说明TraMineR。(由于您的玩具示例中有大量不同的状态,因此生成的替代成本矩阵对于此插图来说太大(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返回的替换成本矩阵),然后在聚类算法中输入结果距离矩阵。 seqsubmsm

于 2018-12-08T09:09:37.300 回答