CODEWITHSUNDEEP
rdata.table
chameau13
chameau13

Reputation: 644

Subsetting within a subset

I am wondering if I could do this efficiently with a data.table. I have got a data set which consists of different samples, for different periods (date) and different groups (id). 

    #the data
    require(data.table)
    dt <- data.table(id=c(rep(1,50),rep(2,50),rep(1,50),rep(2,50)),date=c(rep("2004-01-01",100),rep("2004-02-01",100)),A=c(rnorm(50,1,3),rnorm(50,2,3),rnorm(50,1,4),rnorm(50,1.5,3)),
             B=c(rnorm(50,1.3,2.9),rnorm(50,1.8,3.1),rnorm(50,1.6,4),rnorm(50,1.7,2.4)))

I want to apply the following function.

    #the function which should be applied
    function(a, ie1, b, a1, ie2, b2, ...) {
    ipf <- function(a, b, ...) {
    m <- length(a)
    n <- length(b)
    if (m < n) {
        r <- rank(c(a, b), ...)[1:m] - 1:m
    } else {
        r <- rank(c(a, b), ...)[(m + 1):(m + n)] - 1:n
    }
    s <- ifelse((n + m)^2 > 2^31, sum(as.double(r)), sum(r))/(as.double(m) * n)
    return(ifelse(m < n, s, 1 - s))
}

expand.grid.alt <- function(seq1, seq2) {
    cbind(rep.int(seq1, length(seq2)), c(t(matrix(rep.int(seq2, length(seq1)), nrow = length(seq2)))))
}

if (missing(a1) | missing(b2) | missing(ie2)) {
    if (ie1 == ">") {
        return(ipf(a, b))
    } else {
        return(ipf(b, a))
    }
} else {
    if (ie1 == ">") {
        if (ie2 == ">") {
            return(ipf(a, apply(expand.grid.alt(b, b2), 1, max))/ipf(a1, b2))
        } else {
            return(1 - ipf(apply(expand.grid.alt(b, b2), 1, min), a)/(1 - ipf(a1, b2)))
        }
    } else {
        if (ie2 == ">") {
            return(1 - ipf(a, apply(expand.grid.alt(b, b2), 1, max))/ipf(a1, b2))
        } else {
            return(ipf(apply(expand.grid.alt(b, b2), 1, min), a)/(1 - ipf(a1, b2)))
        }
    }
}

}

This function compares different samples; Given we have three samples A, B, C it allows e.g. to compute the probability that a draw from sample A is greater than a draw from sample B given that the draw from sample A is greater than a draw from sample C. I want to apply this function in a certain manner using data.tables. The following example should illustrate you what I want to do:

    #example - what I want to do
    dt1 <-  dt[date=="2004-01-01"]
    ow <-   dt1[id==1,A]
    ot <-   dt1[id!=1,A]
    cs  <-  dt1[,B]
    ex <- expand.grid(unique(ow),unique(ot),unique(cs))
    names(ex) <- c("ow","ot","cs")
    sum(ex$ow > ex$ot & ex$ow > ex$cs)/sum(ex$ow > ex$ot)

    #check if the result is correct
    all.equal(prob(ow,">",cs,ow,">",ot),sum(ex$ow > ex$ot & ex$ow > ex$cs)/sum(ex$ow > ex$ot))
    [1] TRUE

I want to automatize the procedure above with the use of data.table for all ids and all dates. In words: I want to compute the probability that a draw from variable A of id=1 is greater than a draw from variable B given that a draw from variable A of id=1 is greater than a draw from variable of id!=1 (the use of expand.grid implies the brute force method which looks at all possible combinations, the prob() function above use a more elegant rank-sum approach).

This means I need some kind of subset within a subset. Intuitively I have played around with something like that:

    dt[,.SD[,prob(A,">",B,A,">",.SD[!.BY,A]),key=id],key=date]

This approach however leads to an error messages. Who can help me with this problem? Any comment is highly appreciated!

Upvotes: 0

Views: 270

Answers (1)

Ricardo Saporta
Ricardo Saporta

Reputation: 55420

Importantly: In your example above, note that you are recycling your A values to match the length of the B values. It's not clear if this is what you actually intend, if the answer is wrong, or if the answer is correct, but moreso due to a symmetry than to the actual method. You might want to double check your example. Meanwhile, this does what you have above, in an efficient manner

## USING CJ
setkey(dt, id)
dt[, {
      .SD1 <- .SD;
      .SD1[, {.B <- unlist(.BY);
              CJ( ow=.SD1[.(.B)][["A"]], 
                  ot=.SD1[!.(.B)][["A"]], 
                  cs=.SD1[["B"]]
                )[
                  , sum(ow>ot & ow>cs) / sum(ow > ot)] 
             }
    , by=id ]
    }
  , by=date
  ]

## USING PROB
setkey(dt, id)
dt[, {
      .SD1 <- .SD;
      .SD1[, {.B <- unlist(.BY);
              ow <- .SD1[.(.B)][["A"]] 
              ot <- .SD1[!.(.B)][["A"]]
              cs <- .SD1[["B"]]
              prob(ow,">",cs,ow,">",ot)
             }
    , by=id ]
    }
  , by=date
  ]

Benchmarks:

You are right, the prob function is faster (incidentally, not by much). 

usingProb <- quote(dt[, {.SD1 <- .SD;.SD1[, {.B <- unlist(.BY);ow <- .SD1[.(.B)][["A"]] ;ot <- .SD1[!.(.B)][["A"]];cs <- .SD1[["B"]];prob(ow,">",cs,ow,">",ot)}, by=id ]}, by=date  ])
usingCJ <- quote(dt[, {.SD1 <- .SD;.SD1[, {.B <- unlist(.BY);CJ( ow=.SD1[.(.B)][["A"]], ot=.SD1[!.(.B)][["A"]], cs=.SD1[["B"]])[, sum(ow>ot & ow>cs) / sum(ow > ot)] }, by=id ]}, by=date])

eval(usingProb)
eval(usingCJ)
all.equal(eval(usingProb), eval(usingCJ))

library(microbenchmark)
microbenchmark(PROB=eval(usingProb), CJ=eval(usingCJ), times=20L)

Unit: milliseconds
 expr      min       lq   median       uq      max neval
 PROB 50.59504 53.62986 62.78143 80.64911 106.2133    20
   CJ 67.63520 69.59654 74.56110 79.45636 136.6357    20

Upvotes: 1

Related Questions