如何获取一行在 R 中高于或低于临界值的次数
How to obtain counts of the number of times a row changes from being above or below a critical value in R
我有一个正在使用的数据框，它是 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
我的问题是，虽然这会生成我正在寻找的答案（每行总和的向量偏离 above/below 临界值），但这在较大的数据集上花费的时间太长。我有几组数据框，它们都是大约 6,000 行乘以 46,000 列。我知道 R 对于 for 循环效率低下，但我对 R 相当缺乏经验，对 bash 的经验略多一些，一般来说是编码的新手。任何可以优化此过程的帮助将不胜感激。抱歉，如果这个问题的格式不符合标准，或者如果它在其他地方被问过，这是我的第一个 post 并且我无法在之前的问题中找到解决方案。
更新
小假设数据框和预期输出
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} 或 {df[x,y]<=crit_value & df[x,y+1]>crit_value}，我需要一个计数，这样我就可以得到相对于给定 crit_value.
的符号变化总和
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))
你可以试试：
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)
我有一个正在使用的数据框，它是 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
我的问题是，虽然这会生成我正在寻找的答案（每行总和的向量偏离 above/below 临界值），但这在较大的数据集上花费的时间太长。我有几组数据框，它们都是大约 6,000 行乘以 46,000 列。我知道 R 对于 for 循环效率低下，但我对 R 相当缺乏经验，对 bash 的经验略多一些，一般来说是编码的新手。任何可以优化此过程的帮助将不胜感激。抱歉，如果这个问题的格式不符合标准，或者如果它在其他地方被问过，这是我的第一个 post 并且我无法在之前的问题中找到解决方案。
更新 小假设数据框和预期输出
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} 或 {df[x,y]<=crit_value & df[x,y+1]>crit_value}，我需要一个计数，这样我就可以得到相对于给定 crit_value.
的符号变化总和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))
你可以试试：
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)