1
\$\begingroup\$

From a data frame d I have to select the longest possible sequences of x. variables.

Example:

> 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

Since all observations with missings have to be neglected, there is a tradeoff between sequence length and number of observations. I have to minimize this tradeoff.

For this purpose I have written this function:

seqRank <- function(d, id="id") {
  # generate power subsets of rows and columns
  psr <- HapEstXXR::powerset(as.character(d[[id]]))
  pssr <- lapply(psr, function(x) 
    d[which(d[[id]] %in% x), ])
  psc <- HapEstXXR::powerset(names(d)[-which(names(d) == id)])
  pssc <- lapply(psc, function(x) 
    d[, c(id, x)])
  # generate all combinations of subsets
  sss <- lapply(psr, function(x)
    lapply(pssc, function(y) y[which(y$id %in% x), ]))
  # clean subsets from NAs
  cn <- sapply(sss, function(x)
    lapply(x, function(y) {
      y0 <- y[, which(names(y) == id)]
      y1 <- y[, -which(names(y) == id)]
      if (is.null((dim(y1))) & any(is.na(y1)))
        NULL
      else if (is.null((dim(y1))) & any(!is.na(y1))) 
        setNames(data.frame(as.factor(as.character(y0)), y1), 
                 names(y))
      else if (all(apply(is.na(y1), 2, any)))
        NULL
      else {
        na <- which(apply(is.na(y1), 2, any))
        if (length(na) == 0)
          NA
        else
          setNames(data.frame(as.factor(as.character(y0)), 
                              y1[, -na]),
                   c(id, names(y1[-na])))
      }
    }))
  # count rows and columns of subsets
  scr <- unlist(setNames(sapply(cn, nrow), 
                           sapply(cn, function(x) 
                             paste0(names(x)[-which(names(x) == id)], 
                                    collapse=", "))))
  scc <- unlist(setNames(sapply(cn, ncol), 
                           sapply(cn, function(x) 
                             paste0(names(x)[-which(names(x) == id)], 
                                    collapse=", ")))) - 1
  # bind to a matrix
  m <- t(rbind(n.obs=scr, sq.len=scc))
  # aggregate matrix by sequences and return maximum sequence lengths
  ag <- aggregate(m, by=list(sequence=rownames(m)), max)
  return("rownames<-"(with(ag, ag[order(-sq.len, -n.obs), ]), NULL))
}

It gives me the desired result,

> seqRank(d)
1024 sets to create.
32 sets to create.
         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

but it works quite slowly, even with this small 10x6 data frame, and I have to apply the function to larger data frames that have considerably more rows.

Note that, while working through this answer on Stack Overflow, I noticed that HapEstXXR::powerset is the fastest way to calculate the powersets, however it only calculates up to a maximum of 15 rows, which is why in lines 3 and 6 I probably have to do this:

psr <- do.call(c, lapply(seq_along(d[[id]]), combn, x=d[[id]], simplify=FALSE))
psc <- do.call(c, lapply(seq_along(names(d)[-which(names(d) == id)]), 
                  combn, x=names(d)[-which(names(d) == id)], simplify=FALSE))

I'm not sure now whether the complexity of the calculation itself or my code is slowing down the function. Probably there's a much easier way that I didn't come up with.

I am grateful for all suggestions for improvement.

Data:

d <- structure(list(id = structure(1:10, .Label = c("A", "B", "C", 
                                                    "D", "E", "F", "G", "H", "I", "J"), class = "factor"), X1 = 1:10, 
                    X2 = c(11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, NA, NA), X3 = c(21L, 
                                                                                   22L, 23L, 24L, 25L, 26L, NA, NA, NA, NA), X4 = c(31L, 32L, 
                                                                                                                                    33L, 34L, NA, NA, NA, NA, NA, NA), X5 = c(41L, 42L, NA, NA, 
                                                                                                                                                                              NA, NA, NA, NA, NA, NA)), row.names = c(NA, -10L), class = "data.frame")
\$\endgroup\$
2
  • \$\begingroup\$ In your output, where it says X5, should it not say X1 since that's the name of the column in you input that has 10 non-NA observations? \$\endgroup\$
    – flodel
    Dec 2, 2018 at 0:51
  • \$\begingroup\$ @flodel Thanks, fixed! I missed it because I had been messing with aggregate() just before the post. \$\endgroup\$
    – jay.sf
    Dec 2, 2018 at 8:18

1 Answer 1

1
\$\begingroup\$

I do not get what are you trying to calculate, but this should work much faster: (because of using matrices not data.frames and retaining the structure, it takes a lot of time to create new data.frames inside loops)

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)
}
x2 <- seqRank2(d)

all.equal(x, x2)
# TRUE

P.S. I do not like using setNames, it makes code harder to read, so I rewrote those parts.

\$\endgroup\$
7
  • \$\begingroup\$ I need it for time sequence analysis of survey data with package TraMineR, which doesn't accept any NA in data. So this program helps me to select appropriate time periods (in wide format) without losing too many observations. It's now 15 times faster, thanks! I'll see how it will cope with the big data with a 34235 x 56 matrix, though. \$\endgroup\$
    – jay.sf
    Dec 4, 2018 at 9:29
  • \$\begingroup\$ As expected it didn't work with this huge matrix, getting error vector size cannot be infinite. A smaller 35x17 matrix throws Error: cannot allocate vector of size 256.0 Gb. \$\endgroup\$
    – jay.sf
    Dec 4, 2018 at 10:57
  • \$\begingroup\$ In my code I replaced the HapEstXXR::... lines with the do.call()s as described in the question and tested the function with big data. It consumed all of my 32GB RAM and ran forever, but at least it run. So probably the error is related to the matrixStats package you used? I know I'm planning a huge calculation and I don't know if that's even possible with my tiny machine, though. \$\endgroup\$
    – jay.sf
    Dec 4, 2018 at 11:09
  • \$\begingroup\$ It's actually a 34235 x 17 matrix, not a 34235 x 56. \$\endgroup\$
    – jay.sf
    Dec 4, 2018 at 11:10
  • \$\begingroup\$ @jay.sf try doing psr <- do.call(c, lapply(seq_along(1:25), combn, x = 1:25, simplify = FALSE)) and you will see that it is hard to create even all combinations of 25 elements. So I would suggest that you rethink what you are doing and why. And maybe create new question with your initial problem (what you are trying to resolve with this approach) on stackoverflow, because you will not manage to solve it like this. Better description of your problem may help, because I did not get what are you trying to do. \$\endgroup\$
    – minem
    Dec 4, 2018 at 11:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.