使用R统计不同阈值下Happiness的峰值数量及持续时长

问题背景与需求 #

我是R语言初学者，之前一直用Excel手动分析数据，现在想转用R实现自动化分析。

模拟数据 #

set.seed(94)
Happiness <- round(runif(60, -100, 100))
ID <- rep(1:3, 20)
Stimuli <- rep(1:3, 1)
DF <- data.frame(ID, Stimuli, Happiness)

数据集说明：记录3名受试者观看3种不同图片时的Happiness情绪值（每行对应1秒的时间段）。

分析目标 #

1. 按ID（受试者）和Stimuli（刺激物）统计Happiness在不同正阈值（20/50/70）及负阈值（-20/-50/-70）下的峰值数量；
2. 统计Happiness处于各阈值区间的总时长（秒）；
   最终生成对应的汇总表。

当前疑问 #

1. 是否需要在按ID拆分后再按Stimuli拆分，再运行lapply()？我的目标是对比不同Stimuli下受试者的平均Happiness情况；
2. 代码中function(id)的含义是什么？
3. 如何统计矩阵中连续1的集群数以得到峰值数量？
4. 如何统计矩阵中1的总数以得到总时长，以及如何将矩阵转换为向量或数据框？

问题解答与实现代码 #

1. 关于数据拆分的方式 #

不需要先按ID拆分再按Stimuli拆分，直接按ID和Stimuli分组处理更高效。可以用dplyr的group_by(ID, Stimuli)或者split(DF, list(DF$ID, DF$Stimuli))一次性拆分到最细的分组，之后用lapply()遍历每个分组计算，这样能直接得到每个受试者-刺激物组合的统计结果，正好满足你对比不同Stimuli的需求。

2. 代码中function(id)的含义 #

这是自定义匿名函数的写法，id是函数的参数名（可换成任意合法变量名，比如group_data）。当你用lapply()遍历拆分后的分组时，每个分组的数据会被传入这个参数，函数内部就可以对该分组的数据进行计算。示例：

grouped_data <- split(DF, list(DF$ID, DF$Stimuli))
lapply(grouped_data, function(group) {
  # group就是每个ID+Stimuli组合的子数据集
  mean(group$Happiness)
})

3. 统计连续1的集群数（峰值数量） #

连续的1代表情绪值处于阈值之上/之下的连续时间段，一个连续集群就是一个峰值。用rle()（游程编码）函数即可实现：

• 先把情绪值转换成0/1逻辑向量（满足阈值为1，否则为0）；
• 用rle()提取连续值的长度和对应值；
• 筛选出值为1的游程，统计其数量就是峰值数。

示例函数：

count_peaks <- function(vec, threshold, direction = "above") {
  # direction可选"above"（超过正阈值）或"below"（低于负阈值）
  binary_vec <- if (direction == "above") {
    as.integer(vec >= threshold)
  } else {
    as.integer(vec <= threshold)
  }
  rle_result <- rle(binary_vec)
  sum(rle_result$values == 1)
}

4. 统计1的总数（总时长）与格式转换 #

• 统计1的总数直接用sum()函数，0/1向量的和就是1的个数，对应总秒数；
• 矩阵转向量用as.vector()，转数据框用as.data.frame()，或者在计算时直接整理成数据框格式。

完整实现代码（含汇总表生成） #

用dplyr和tidyr实现简洁的自动化分析：

library(dplyr)
library(tidyr)

# 修正模拟数据的Stimuli生成（原代码仅生成3个值，不符合每行1秒的场景）
set.seed(94)
Happiness <- round(runif(60, -100, 100))
ID <- rep(1:3, 20)
Stimuli <- rep(rep(1:3, each = 7), 2)[1:60] # 每个ID对应每种Stimuli有6-7条数据
DF <- data.frame(ID, Stimuli, Happiness)

# 定义阈值列表
thresholds <- list(
  positive = c(20, 50, 70),
  negative = c(-20, -50, -70)
)

# 编写分组分析函数
analyze_group <- function(data) {
  result <- list()
  
  # 处理正阈值
  for (thresh in thresholds$positive) {
    binary <- data$Happiness >= thresh
    peak_count <- sum(rle(binary)$values == TRUE)
    duration <- sum(binary)
    result[[paste0("above_", thresh)]] <- data.frame(peak_count, duration)
  }
  
  # 处理负阈值
  for (thresh in thresholds$negative) {
    binary <- data$Happiness <= thresh
    peak_count <- sum(rle(binary)$values == TRUE)
    duration <- sum(binary)
    result[[paste0("below_", thresh)]] <- data.frame(peak_count, duration)
  }
  
  # 合并结果并整理格式
  bind_cols(result) %>%
    mutate(ID = unique(data$ID), Stimuli = unique(data$Stimuli)) %>%
    relocate(ID, Stimuli)
}

# 分组计算并生成汇总表
summary_table <- DF %>%
  group_by(ID, Stimuli) %>%
  group_modify(~ analyze_group(.x)) %>%
  ungroup()

# 查看结果
print(summary_table)

代码说明 #

• 修正了原模拟数据中Stimuli的生成逻辑，确保每个ID对应每种刺激有多个时间点数据；
• group_modify()直接对group_by后的分组应用分析函数，比split()+lapply()更简洁；
• 最终的summary_table包含每个ID+Stimuli组合在所有阈值下的峰值数和时长，方便对比不同刺激物的差异。

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