如何获取行从高于或低于R中的临界值变化的次数的计数

Posted 2023-03-11

【中文标题】如何获取行从高于或低于R中的临界值变化的次数的计数

【英文标题】：How to obtain counts of the number of times a row changes from being above or below a critical value in R

【发布时间】：2020-02-12 14:01:01

【问题描述】：

我有一个正在使用的数据框，它是 HMM 输出的一系列概率。我想知道概率从高于任意临界值切换到低于该值的次数，反之亦然。我对 R 很陌生，虽然我开发了一个产生输出的代码，但它相当耗时。

> Haplo                         #Subset of original dataframe
chr2L_502618 chr2L_502999 chr2L_504449 chr2L_504509 chr2L_504686 chr2L_504688 chr2L_504690 chr2L_504706 chr2L_505918 chr2L_506002
3       0.04865      0.04864       0.0486       0.0486       0.0486       0.0486       0.0486       0.0486      0.04857      0.04856
4       0.04769      0.04767      0.04764      0.04764      0.04764      0.04764      0.04764      0.04764      0.04761       0.0476
5       0.04817      0.04817      0.04813      0.04813      0.04813      0.04813      0.04813      0.04813      0.04808      0.04807
6        0.0612      0.06118      0.06114      0.06114      0.06114      0.06114      0.06113      0.06113      0.06112      0.06112
7       0.41175      0.41178      0.41193      0.41194      0.41194      0.41194      0.41194      0.41194      0.41206       0.4121
8       0.04754      0.04752      0.04749      0.04749      0.04749      0.04749      0.04749      0.04749      0.04746      0.04745
9       0.27742      0.27742      0.27751      0.27751      0.27751      0.27751      0.27751      0.27751      0.27756      0.27759
10      0.05761       0.0576      0.05757      0.05757      0.05756      0.05756      0.05756      0.05756      0.05753      0.05753
11      0.00067      0.00065      0.00059      0.00059      0.00059      0.00059      0.00059      0.00059      0.00055      0.00053
12      0.00075      0.00073      0.00067      0.00067      0.00067      0.00067      0.00067      0.00067      0.00063      0.00061
> probs <- array(0,dim=dim(Haplo))
> for (i in 1:ncol(probs)) probs[,i] <- as.character(Haplo[,i])
> crits <- matrix(as.numeric(probs>0.27751),nrow=nrow(probs),ncol=ncol(probs))
> crits              
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    0    0    0    0    0    0    0    0    0     0
 [2,]    0    0    0    0    0    0    0    0    0     0
 [3,]    0    0    0    0    0    0    0    0    0     0
 [4,]    0    0    0    0    0    0    0    0    0     0
 [5,]    1    1    1    1    1    1    1    1    1     1
 [6,]    0    0    0    0    0    0    0    0    0     0
 [7,]    0    0    0    0    0    0    0    0    1     1
 [8,]    0    0    0    0    0    0    0    0    0     0
 [9,]    0    0    0    0    0    0    0    0    0     0
[10,]    0    0    0    0    0    0    0    0    0     0

这给了我一个数据框，其中任何高于临界值的值为 1，低于临界值的任何值为 0，然后我可以将其输入嵌套的 for 循环以判断行何时从 0 更改为 1，反之亦然

> shifts <- c()
> for (g in 1:nrow(crits))
+     for (i in 1:(ncol(crits)-1))
+         shifts <- c(shifts, sapply(crits[g,i], identical, y=crits[g,i+1]))
+      
+  
> shifts2 <- matrix(as.numeric(!shifts), nrow=nrow(crits), ncol=(ncol(crits)-1), byrow=TRUE)
> shifts2                   #Times a column isn't identical to previous by row
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
 [1,]    0    0    0    0    0    0    0    0    0
 [2,]    0    0    0    0    0    0    0    0    0
 [3,]    0    0    0    0    0    0    0    0    0
 [4,]    0    0    0    0    0    0    0    0    0
 [5,]    0    0    0    0    0    0    0    0    0
 [6,]    0    0    0    0    0    0    0    0    0
 [7,]    0    0    0    0    0    0    0    1    0
 [8,]    0    0    0    0    0    0    0    0    0
 [9,]    0    0    0    0    0    0    0    0    0
[10,]    0    0    0    0    0    0    0    0    0
> sums <- c()
> for (i in 1:nrow(shifts2))
+      sums <- c(sums, sum(shifts2[i,]))
+      
> sums
 [1] 0 0 0 0 0 0 1 0 0 0

我的问题是，虽然这会生成我正在寻找的答案（每行总和的向量从高于/低于临界值移动），但在较大的数据集上花费的时间太长。我有几组数据框，它们都是大约 6,000 行 x 46,000 列。我知道 R 在 for 循环方面效率低下，但我对 R 相当缺乏经验，而且对 bash 的经验也稍多一些，通常是编码新手。任何可以优化此过程的帮助将不胜感激。抱歉，如果此问题未按标准格式化或已在其他地方提出过，这是我的第一篇文章，我无法在之前的提问中找到解决方案。

更新 小假设数据框和预期输出

          X1         X2         X3        X4         X5
1  0.9650217 0.07409232 0.22213328 0.3121305 0.31466359
2  0.1475712 0.06802015 0.63699272 0.2434809 0.17147398
3  0.2951922 0.65086116 0.09405872 0.2389092 0.10440221
4  0.6780534 0.73516696 0.62324000 0.9203979 0.89965700
5  0.4788420 0.16794910 0.13661247 0.5266925 0.52919389
6  0.6738885 0.68843836 0.17165125 0.2478758 0.94910386
7  0.8461378 0.74790781 0.16186888 0.8145674 0.13336087
8  0.3557357 0.65646290 0.21965522 0.6859082 0.55574490
9  0.5262744 0.74453676 0.18037489 0.2106494 0.01274704
10 0.9694096 0.41149759 0.03084501 0.8243646 0.42332927
critical_value=0.3
#expected output: 2, 2, 2, 0, 2, 2, 3, 2, 1, 2

澄清一下，随时 df[x,y]>crit_value & df[x,y+1] crit_value，我需要一个计数，这样我才能得到相对于给定 crit_value 的符号变化的总和。

【问题讨论】：

只有值 [5,3] 高于，[7,3] 完全相等，其余低于。基于此为该子集选择临界值，但它似乎适用于测试中的任何给定值。

@HowlArgwen 添加一个小数据集（例如 10 行 x 5 列）并显示您的预期输出。

【参考方案1】：

你可以试试：

colSums(diff(t(as.matrix(df) > .3)) != 0)

 1  2  3  4  5  6  7  8  9 10 
 2  2  2  0  2  2  3  2  1  2    

数据：

df <- df <- read.table(text = "          X1         X2         X3        X4         X5
1  0.9650217 0.07409232 0.22213328 0.3121305 0.31466359
2  0.1475712 0.06802015 0.63699272 0.2434809 0.17147398
3  0.2951922 0.65086116 0.09405872 0.2389092 0.10440221
4  0.6780534 0.73516696 0.62324000 0.9203979 0.89965700
5  0.4788420 0.16794910 0.13661247 0.5266925 0.52919389
6  0.6738885 0.68843836 0.17165125 0.2478758 0.94910386
7  0.8461378 0.74790781 0.16186888 0.8145674 0.13336087
8  0.3557357 0.65646290 0.21965522 0.6859082 0.55574490
9  0.5262744 0.74453676 0.18037489 0.2106494 0.01274704
10 0.9694096 0.41149759 0.03084501 0.8243646 0.42332927", header = TRUE)

【参考方案2】：

R 中的经验法则是，如果您想编写快速代码，则必须使用向量化的 R 函数，而不是循环。根据我对您问题的理解，我编写了一个函数来满足您的要求：

find_switch <- function(test_ds, crit_val)
 m <- sapply(test_ds, function(x) as.integer(x > crit_val))
 tm <- t(m)
 nrtm <- nrow(tm)
 colSums(tm - rbind(tm[1,], tm[1:(nrtm-1),]) != 0)

注意我在矩阵上使用向量化操作。

我将你的代码包装成一个函数：

find_switch2 <- function(test_ds, crit_val)
  crits <- matrix(as.numeric(test_ds > crit_val),nrow=nrow(test_ds),ncol=ncol(test_ds))
  shifts <- c()
  for (g in 1:nrow(crits))
    for (i in 1:(ncol(crits)-1))
      shifts <- c(shifts, sapply(crits[g,i], identical, y=crits[g,i+1]))
      
  

  shifts2 <- matrix(as.numeric(!shifts), nrow=nrow(crits), ncol=(ncol(crits)-1), byrow=TRUE)

  sums <- c()
  for (i in 1:nrow(shifts2))
    sums <- c(sums, sum(shifts2[i,]))
    
  sums

并提出了一些模拟数据集来对以下两个函数进行基准测试：

set.seed(123)
n_row <- 5e2

crit_val <- 0.3

test_ds <- data.frame(p1 = runif(n_row),
                      p2 = runif(n_row),
                      p3 = runif(n_row),
                      p4 = runif(n_row))

临界值设置为0.3。

然后我对两个实现都进行了计时：

microbenchmark::microbenchmark(find_switch(test_ds, crit_val), find_switch2(test_ds, crit_val))

 #Unit: microseconds expr       min         lq       mean    median         uq       max neval
 #find_switch(test_ds, crit_val)    96.265   121.8295   177.7687   176.132   206.4575   352.265   100
 #find_switch2(test_ds, crit_val) 27499.848 31556.8755 36564.2898 34315.394 40223.6580 93957.460   100

速度差异为 250 倍。所以，这就是为什么使用矢量化函数很重要的原因。

最后，让我们确保这两个函数产生相同的输出：

identical(find_switch(test_ds, 0.3), find_switch2(test_ds, 0.3))

【讨论】：

我尝试在我的原始数据集上运行您的函数。您的函数似乎给了我该行超过临界值的次数。我试图告诉它从超过变为低于或反之亦然的次数。口头上，如果一半的数字连续高于或低于临界值，我需要知道它们是否被分组为两个块（即 100 个值，首先高于 50，然后低于 50），或者切换每个值（高于，低于，上面，下面... x25)，或介于两者之间。我也许可以适应 ifelse 功能。不过，这有一些有用的想法，所以谢谢。

@HowlArgwen 包括一个小的测试数据集（5 x 10）和它的预期输出，我会更容易理解你想要什么

@HowlArgwen 请查看上面的测试数据集，并在需要时更正预期输出。

对于您的测试数据集，预期输出应为：1、0、0、0、1。我将尝试创建一个数据集，更清楚地说明我需要什么并发布。

@HowlArgwen 你能解释一下为什么第 1 行的输出为 1 吗？两个元素（col 1 和 col 2）超过 0.3 的临界值。