R数据表：对于每个组成员，在子集向量上计算函数

Question

我有一个数据表格，类似于以下内容：

set.seed(1)

dt<-data.table(med=sample(letters,50,T),
    diag=sample(LETTERS[1:7],50,T),
    val=sample(1:100,50,F))

我希望计算任何val在相同diag中大于val的概率，并将其分配给表的新列，例如prob（我知道概率不一定是正常的。我可以接受这种情况）。

我可以使用for循环来完成：

for (i in 1:50){
    dg<-dt[i,diag]
    vl<-dt[i,val]
    dt$prob[i]<-pnorm(vl,
                      mean(dt[diag==dg,val]),
                      sd(dt[diag==dg,val]),
                      lower.tail = F)
}

但是我的数据相当大（dt大约有800k行，对角线上有一些2k级别的数据），因此我想矢量化而不是循环。

我尝试过

dt[,
   .(lapply(.SD,function(x) 
                pnorm(x[1],
                mean(x),
                sd(x),
                lower.tail = F))),
   by=diag,
   .SDcols="val"]

当然，按diag分组只会产生一个概率，因此没有什么用处。我也尝试过

dt[,
   .(lapply(.SD,function(x) 
                pnorm(x[1],
                mean(x),
                sd(x),
                lower.tail = F))),
   by=.EACHI,
   .SDcols="val"]

但是它会产生一个错误：

Error in `[.data.table`(dt, , .(lapply(.SD, function(x) pnorm(x[1], mean(x),  : 
  logicial error. i is not data.table, but mult='all' and 'by'=.EACHI

如何通过向量化生成所需结果的代码？

由于我正在逐渐熟悉 data.table，因此我更喜欢使用该包的解决方案，但我完全可以接受其他解决方案（plyr、dplyr等）。

谢谢。

- PavoDive

---

似乎dt[,prob2:= 1 - pnorm(val,mean(val),sd(val),lower.tail=FALSE),by=diag]或类似的语句可以匹配您的结果。不确定为什么需要执行1 - pnorm操作。 - thelatemail

结果不应该是长度为7吗？我想要一个长度为50的结果，对于每一行都计算prob，只是它是通过考虑所有共享相同diag值的所有值的分布来计算的。 - PavoDive

不是的，by= 默认会对每个组中的每一行进行计算。这就是整个操作的意义所在。 - thelatemail

我对你的示例有几个问题。50比你需要的要多，如果在答案中显示结果会占用很多空间。而且你的示例是随机生成的，但你没有设置种子，因此它是不可重复的。 - Frank

现在dt[,prob2:= pnorm(val,mean(val),sd(val),lower.tail=FALSE),by=diag]会得到相同的结果 - 只有在您使用了lower.tail=F而不是lower.tail=FALSE并且我有一个名为F的变量时，才需要1-。全写出FALSE和TRUE是一个好习惯。 - thelatemail

就性能而言，使用apply函数并不能胜过for循环。如果你查看apply的代码实现，会发现其本质就是一个for循环。 - Dean MacGregor

Answer 1

在 data.table 中：

dt[, prob2 := pnorm(val, mean(val), sd(val), lower.tail=FALSE), by=diag]

似乎符合您的要求：

head(dt)
#   med diag val       prob      prob2
#1:   p    E  91 0.04713131 0.04713131
#2:   f    E   3 0.92991675 0.92991675
#3:   o    B  26 0.83792988 0.83792988
#4:   t    C  38 0.70877125 0.70877125
#5:   g    E  71 0.16909178 0.16909178
#6:   i    E  25 0.75428819 0.75428819

- thelatemail

---

我甚至没有考虑尝试这个，这就是为什么我选择了.SD。我不完全理解在调用pnorm时如何区分第一个val是值，而第二个和第三个是向量。我认为如果我在val上调用pnorm，那么我会得到一个向量作为结果。如果您能帮助我理解为什么不是这样，我将不胜感激。谢谢。 - PavoDive

@PavoDive - pnorm是向量化的。给定pnorm(1:3,mean=mean(1:10),sd=sd(1:10))，它将循环遍历1:3，对于每个计算，'循环使用'均值(5.5)和标准差(3.02)的单个值。这本质上与pnorm(1,mean(1:10),sd(1:10)); pnorm(2,mean(1:10),sd(1:10)); pnorm(3,mean(1:10),sd(1:10))相同。 - thelatemail

Answer 2

Here's a dplyr solution:

dt %>% group_by(diag) %>% 
       mutate(prob = pnorm(val, mean(val), sd(val), lower.tail = FALSE))

- jeremycg