在dplyr中按ID计算阈值后，按ID和Stimuli统计峰值

这是《Counting peaks in r per group》的跟进问题。

可复现数据 #

set.seed(949494)
Happiness <- round(runif(100, -100, 100))
ID <- rep(c("ID1", "ID2", "ID3", "ID4", "ID5"), 20)
Stimuli <- rep(1:4, 1)
DF <- data.frame(ID, Stimuli, Happiness)

现有函数 #

计算ID级1SD阈值的函数 #

# 1SD
f.SD1 <- function(y) {
  SD1_thresh <- mean(y) + (1*sd(y))
  return(SD1_thresh)
}

判断是否为峰值的函数 #

# SD1
f.Peaks_SD1 <- function(X, thresh) {
  H_peaks_1 <- ifelse(X >= thresh ,1,0)
  return(H_peaks_1)
}

需求 #

需先按**ID（跨所有Stimuli）**计算统一阈值，再按ID和Stimuli分组统计每个刺激下的峰值数量、总观测数（对应总时长）等指标，同时保留ID+Stimuli的分组输出结构，避免同一ID因子分组计算导致阈值不一致。

dplyr简洁实现方案 #

可以通过"先全局计算ID阈值，再分组统计"的链式操作实现，以下提供两种方式：

方式一：复用现有函数 #

library(dplyr)

result <- DF %>%
  # 为每个ID计算统一阈值并添加到所有行
  group_by(ID) %>%
  mutate(SD1_thresh = f.SD1(Happiness)) %>%
  ungroup() %>%
  # 标记每个观测是否为峰值
  mutate(peak = f.Peaks_SD1(Happiness, SD1_thresh)) %>%
  # 按ID+Stimuli分组统计指标
  group_by(ID, Stimuli) %>%
  summarise(
    total_duration = n(),  # 该分组总时长（观测数）
    peak_count = sum(peak),  # 峰值数量
    peak_avg_ratio = mean(peak)  # 峰值占总时长的比例
  ) %>%
  ungroup()

print(result)

方式二：整合逻辑简化代码（无需额外函数） #

如果不需要保留原有函数，可直接将计算逻辑嵌入管道，代码更紧凑：

library(dplyr)

result <- DF %>%
  group_by(ID) %>%
  mutate(SD1_thresh = mean(Happiness) + sd(Happiness)) %>%
  ungroup() %>%
  mutate(peak = as.integer(Happiness >= SD1_thresh)) %>%
  group_by(ID, Stimuli) %>%
  summarise(
    total_duration = n(),
    peak_count = sum(peak),
    peak_avg_ratio = mean(peak)
  ) %>%
  ungroup()

print(result)

核心逻辑说明 #

1. 全局ID阈值计算：通过group_by(ID) + mutate，为每个ID的所有观测行赋予相同的阈值，确保同一ID下阈值统一。
2. 峰值标记：逐行判断Happiness是否超过对应ID的阈值，生成二进制标记列。
3. 分组统计：最后按ID+Stimuli聚合，输出你需要的统计指标和分组结构。

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