如何在R中按ID分组基于标准差阈值计数幸福感峰值

这是《Count Peaks with R》问题的后续内容，目标是基于每个受试者（ID）的幸福感数据标准差推导阈值，完成峰值计数与超阈值时长统计。

模拟可复现数据 #

set.seed(9494)
Happiness <- round(runif(100, -100, 100))
ID <- rep(c("ID1", "ID2", "ID3", "ID4", "ID5"), 20)
Stimuli <- rep(1:4, 25)  # 修正：原代码仅生成4个值，调整为匹配100行数据
DF <- data.frame(ID, Stimuli, Happiness)

数据说明：DF记录了5名受试者观看刺激物时的幸福感数值，每行对应1秒的观测数据。

步骤1：按ID计算专属阈值 #

先按ID拆分数据集，再为每个ID计算基于均值±1/2/3倍标准差的阈值：

# 按ID拆分数据集
DF.id <- split(DF, DF$ID)

# 定义阈值计算函数
TR_SD <- function(Y){
  mean_val <- mean(Y)
  sd_val <- sd(Y)
  SD1_thresh <- mean_val + sd_val
  SD2_thresh <- mean_val + 2*sd_val
  SD3_thresh <- mean_val + 3*sd_val
  SD1_neg_thresh <- mean_val - sd_val
  SD2_neg_thresh <- mean_val - 2*sd_val
  SD3_neg_thresh <- mean_val - 3*sd_val
  return(data.frame(SD1_thresh, SD2_thresh, SD3_thresh, 
                    SD1_neg_thresh, SD2_neg_thresh, SD3_neg_thresh))
}

# 计算每个ID的专属阈值
SD.Thresh <- lapply(DF.id, function(x) TR_SD(x$Happiness))

步骤2：绑定ID阈值与数据，完成超阈值判断 #

核心修改：将每个ID的阈值与对应数据一起传入判断函数，确保阈值与个体匹配：

# 定义带阈值输入的超阈值判断函数
Thresh <- function(happiness_vec, thresholds){
  H_peaks_1a <- ifelse(happiness_vec >= thresholds$SD1_thresh, 1, 0)
  H_peaks_2a <- ifelse(happiness_vec >= thresholds$SD2_thresh, 1, 0)
  H_peaks_3a <- ifelse(happiness_vec >= thresholds$SD3_thresh, 1, 0)
  H_neg_peaks_1a <- ifelse(happiness_vec <= thresholds$SD1_neg_thresh, 1, 0)
  H_neg_peaks_2a <- ifelse(happiness_vec <= thresholds$SD2_neg_thresh, 1, 0)
  H_neg_peaks_3a <- ifelse(happiness_vec <= thresholds$SD3_neg_thresh, 1, 0)
  return(cbind(H_peaks_1a, H_peaks_2a, H_peaks_3a, 
               H_neg_peaks_1a, H_neg_peaks_2a, H_neg_peaks_3a))
}

# 为每个ID应用专属阈值
H_peaks.ID <- mapply(Thresh, 
                     happiness_vec = lapply(DF.id, function(x) x$Happiness),
                     thresholds = SD.Thresh,
                     SIMPLIFY = FALSE)

步骤3：峰值计数 #

统计每个ID在各阈值下的峰值数量（峰值定义为从超阈值状态回到非超阈值状态的次数）：

# 计算峰值数量
peaks <- t(sapply(H_peaks.ID, function(x) {
  apply(x, 2, function(y) sum(diff(c(y, 0)) < 0))
}))
peaks <- as.data.frame(peaks)
colnames(peaks) <- c("正峰值(1σ)", "正峰值(2σ)", "正峰值(3σ)", 
                     "负峰值(1σ)", "负峰值(2σ)", "负峰值(3σ)")

步骤4：超阈值时长统计 #

统计每个ID在各阈值下的总超阈值时长（即1的总和，对应秒数）：

# 计算超阈值总时长
time <- t(sapply(H_peaks.ID, function(x) apply(x, 2, sum)))
time <- as.data.frame(time)
colnames(time) <- c("正超阈值时长(1σ)", "正超阈值时长(2σ)", "正超阈值时长(3σ)", 
                    "负超阈值时长(1σ)", "负超阈值时长(2σ)", "负超阈值时长(3σ)")

合并汇总结果 #

将峰值计数与时长统计合并为最终结果：

# 合并结果，保留ID作为行名
final_result <- cbind(peaks, time)
final_result

内容的提问来源于stack exchange，提问作者Smuts94