【问题标题】:Subsetting within a subset子集中的子集
【发布时间】:2013-10-08 22:36:22
【问题描述】:

我想知道我是否可以使用data.table 有效地做到这一点。我有一个由不同样本组成的数据集,用于不同的时期(日期)和不同的组(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)))

我想应用以下功能。

    #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)))
        }
    }
}

}

这个函数比较不同的样本;鉴于我们有三个样本 A、B、C,它允许例如计算样本 A 的抽取大于样本 B 的抽取的概率,假设样本 A 的抽取大于样本 C 的抽取。我想使用 data.tables 以某种方式应用此函数。下面的例子应该能说明我想要做什么:

    #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

我想通过对所有 id 和所有日期使用 data.table 来自动化上述过程。换句话说:我想计算从 id=1 的变量 A 的平局大于从变量 B 的平局的概率,因为从 id=1 的变量 A 的平局大于从 id 的变量的平局!=1 (expand.grid 的使用意味着使用蛮力方法查看所有可能的组合,上面的 prob() 函数使用更优雅的秩和方法)。

这意味着我需要子集中的某种子集。直觉上我玩过类似的东西:

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

但是,这种方法会导致错误消息。谁能帮我解决这个问题?任何评论都非常感谢!

【问题讨论】:

  • 能否用一个更小的函数来说明您的尝试?
  • 我很困惑你为什么使用 expand.grid...它将 ex 与 dt 中的数据行分开。我认为您的条件概率可能没有明确定义。
  • 您的dt[.... 行中存在一些问题,但最重要的是,您希望prob 的第6 个参数是什么? (即,.SD[!.BY,A] 对您来说代表什么?)
  • @RicardoSaporta 我使用 by=date 按日期(月)进行子集化。然后我想再次按 = id 拆分日期子集。但是,在这里我必须将一个 id 的样本与所有其他 id 的样本进行比较。这意味着通过使用 .SD[!.BY,A] 我想在相应的日期子样本 .SD 中排除所有不属于相应 id 的 A 的所有观察值(此时 alogirthm 看起来)
  • @Frank prob 函数工作得非常好。它是产生问题的子集的子集。我正在尝试找到正确的 data.table 语法。

标签: r data.table


【解决方案1】:

重要的是:在上面的示例中,请注意您正在回收 A 值以匹配 B 值的长度。目前尚不清楚这是否是您的实际意图,答案是否错误,或者答案是否正确,但更多的是由于对称性而不是实际方法。您可能需要仔细检查您的示例。 同时,这可以有效地完成上述操作


## 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
  ]

基准测试:

你是对的,prob 函数更快(顺便说一句,不是很多)。

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

【讨论】:

  • 感谢您的回复。不幸的是,您的建议中的机器人提供了不正确的结果。这就是为什么我添加了上面的示例以提供可以比较的参考。我想我想做的事情一定是有可能的。
  • @chameau13: library(fortunes); fortune(109)
  • 我只是想运行一个我认为昨天还在工作的示例。 prob 函数是比较样本的秩和检验。它节省了计算时间,因为不必每次都计算所有组合。我希望我能尽快更新我的问题。
  • 我想用.N 替换sum(ow &gt; ot) 并设置i=ow &gt; ot 应该会快一点。
  • @Frank,这是个好建议
猜你喜欢
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
  • 2013-01-08
  • 2017-12-10
  • 1970-01-01
  • 1970-01-01
  • 1970-01-01
相关资源
最近更新 更多