对于仅具有因子列的给定数据框,我想列出不出现在数据中的最多m属性的因子的所有组合。下面是一个简单的示例:

d <- expand.grid(w=factor(1:2), x=factor(1:2), y=factor(1:2),
                 z=factor(1:2))

# These combinations are removed by tail():
rmcomb <- 5; head(d, rmcomb)

##   w x y z
## 1 1 1 1 1
## 2 2 1 1 1
## 3 1 2 1 1
## 4 2 2 1 1
## 5 1 1 2 1


d <- tail(d, -rmcomb)
ftable(d, row.vars=c("w", "x"))

##     y 1   2
##     z 1 2 1 2
## w x
## 1 1   0 1 0 1
##   2   0 1 1 1
## 2 1   0 1 1 1
##   2   0 1 1 1

对于m == 3,我们考虑d中最多3个属性的所有4 + 6 + 4 = 14个组合:
m <- 3
library(plyr)
llply(
  1:m,
  function(i) combn(ncol(d), i, simplify=F)
) -> cc
unlist(cc, recursive=F) -> cc
length(cc)

## [1] 14

现在,我们可以使用tableuse which 将数据的选定列制成表格,以查找带有零的条目:
llply(
  cc,
  function(cols) {
    which(table(d[, cols]) == 0, arr.ind=T) -> z
    colnames(z) <- names(d)[cols]
    if (nrow(z) > 0) list(z) else NULL
  }
) -> zz
unlist(zz, recursive=F)

## [[1]]
##   y z
## 1 1 1
##
## [[2]]
##   w x z
## 1 1 1 1
##
## [[3]]
##   w y z
## 1 1 1 1
## 2 2 1 1
##
## [[4]]
##   x y z
## 1 1 1 1
## 2 2 1 1

但是,上面结果中的项目[[3]][[4]]是多余的,因为它们已被项目[[1]]覆盖(=没有对y == 1的观察,z == 1)。因此,解决方案应为(y,z) == (1,1); (w,x,z) == (1,1,1)

R中是否有一个内置工具可以用更少的编码来解决问题,也许包括删除冗余(=被覆盖的)元组?如果不是,您如何删除上面代码的这些多余项?

最佳答案

这是您如何继续您的算法来挑选那些序列的方法。首先,让我们将列表转换为矩阵,并填写NA。我发现这更容易处理,但是我敢肯定,您可以通过一些努力使它也可以与列表一起使用:

m = as.matrix(rbind.fill(lapply(zz, as.data.frame)))
#      y z  w  x
#[1,]  1 1 NA NA
#[2,] NA 1  1  1
#[3,]  1 1  1 NA
#[4,]  1 1  2 NA
#[5,]  1 1 NA  1
#[6,]  1 1 NA  2

现在让我们介绍一个函数,该函数将告诉我们subseq给定的矩阵的每一行是否是seq的“子序列”,这意味着按照OP的定义,它已经被seq覆盖了:
is.subsequence = function(seq, subseq) {
  comp = seq == t(subseq)

  rowSums(t(is.na(comp) == is.na(seq) &
            matrix(!(comp %in% FALSE), nrow = length(seq)))) == length(seq)
}

剩下的就是迭代矩阵并抛出覆盖的序列。由于OP中zz的自动排列,我们可以从上到下进行此操作。
i = 1
while(i < nrow(m)) {
  m = rbind(m[1:i,], tail(m, -i)[!is.subsequence(m[i,], tail(m, -i)),])

  i = i+1
}

m
#      y z  w  x
#[1,]  1 1 NA NA
#[2,] NA 1  1  1

如果愿意,您可以返回列表:
apply(m, 1, na.omit)

关于r - 列出数据帧中没有观察值的所有因素(相互作用)组合,直到给定维度为止,删除冗余,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/18278578/

10-12 17:20
查看更多