如何按秒/分/时/日/周粒度对RawCH4与ConvCH4做回归及相关性分析？

针对你的大数据量时间序列数据集，我们可以用R的tidyverse、lubridate和broom包实现自动化的多粒度分组分析，无需手动子集化数据。以下是具体步骤：

1. 加载必要工具包 #

library(tidyverse)
library(lubridate)
library(broom)
# 可选：用并行加速大数据处理
# library(furrr)
# plan(multisession)

2. 清洗并转换时间列 #

首先将字符串格式的MeasTime转换为可操作的时间类型，适配示例中的特殊时间格式：

df <- df %>%
  mutate(
    # 解析"月/日/年 小时:分钟秒"格式（如0:0000对应00:00:00）
    datetime = parse_date_time(MeasTime, orders = "mdy HMOS")
  )

3. 创建多粒度分组变量 #

生成秒、分钟、小时、日期、周五个级别的分组标识：

df <- df %>%
  mutate(
    second_grp = datetime,          # 秒级（原始时间戳）
    minute_grp = floor_date(datetime, "minute"),  # 分钟级（截断到分钟）
    hour_grp = floor_date(datetime, "hour"),      # 小时级（截断到小时）
    date_grp = floor_date(datetime, "day"),       # 日期级（截断到天）
    # 周级：按周一为周起始，若需周日起始去掉week_start=1
    week_grp = floor_date(datetime, "week", week_start = 1)
  )

4. 定义分析函数 #

编写一个函数，对每个分组自动执行回归与相关性分析，并整理结构化结果：

analyze_group <- function(data) {
  # 跳过数据量不足2行的分组，避免报错
  if (nrow(data) < 2) {
    return(tibble(
      intercept_est = NA, intercept_se = NA, intercept_p = NA,
      convch4_est = NA, convch4_se = NA, convch4_p = NA,
      r_squared = NA, cor_coef = NA, cor_p = NA
    ))
  }
  
  # 线性回归分析
  lm_fit <- lm(RawCH4 ~ ConvCH4, data = data)
  lm_tidy <- tidy(lm_fit)
  
  # 相关性分析
  cor_test <- cor.test(data$RawCH4, data$ConvCH4)
  
  # 整理输出结果
  tibble(
    intercept_est = lm_tidy$estimate[lm_tidy$term == "(Intercept)"],
    intercept_se = lm_tidy$std.error[lm_tidy$term == "(Intercept)"],
    intercept_p = lm_tidy$p.value[lm_tidy$term == "(Intercept)"],
    convch4_est = lm_tidy$estimate[lm_tidy$term == "ConvCH4"],
    convch4_se = lm_tidy$std.error[lm_tidy$term == "ConvCH4"],
    convch4_p = lm_tidy$p.value[lm_tidy$term == "ConvCH4"],
    r_squared = glance(lm_fit)$r.squared,
    cor_coef = cor_test$estimate,
    cor_p = cor_test$p.value
  )
}

5. 批量执行多粒度分析 #

对每个时间粒度分组自动运行分析，最终合并所有结果：

# 秒级分析
second_results <- df %>%
  group_by(second_grp) %>%
  nest() %>%
  # 若用并行加速，将map替换为future_map
  mutate(analysis = map(data, analyze_group)) %>%
  unnest(analysis) %>%
  select(-data) %>%
  mutate(granularity = "second")

# 分钟级分析
minute_results <- df %>%
  group_by(minute_grp) %>%
  nest() %>%
  mutate(analysis = map(data, analyze_group)) %>%
  unnest(analysis) %>%
  select(-data) %>%
  mutate(granularity = "minute")

# 小时级分析
hour_results <- df %>%
  group_by(hour_grp) %>%
  nest() %>%
  mutate(analysis = map(data, analyze_group)) %>%
  unnest(analysis) %>%
  select(-data) %>%
  mutate(granularity = "hour")

# 日期级分析
date_results <- df %>%
  group_by(date_grp) %>%
  nest() %>%
  mutate(analysis = map(data, analyze_group)) %>%
  unnest(analysis) %>%
  select(-data) %>%
  mutate(granularity = "day")

# 周级分析
week_results <- df %>%
  group_by(week_grp) %>%
  nest() %>%
  mutate(analysis = map(data, analyze_group)) %>%
  unnest(analysis) %>%
  select(-data) %>%
  mutate(granularity = "week")

# 合并所有粒度结果到一个数据框
all_granularity_results <- bind_rows(
  second_results, minute_results, hour_results, date_results, week_results
)

结果说明 #

最终的all_granularity_results包含核心字段：

• 各粒度的时间分组标识（如minute_grp、date_grp）
• granularity：标记当前行所属的时间粒度
• 回归结果：截距/斜率的估计值、标准误、p值，以及R²
• 相关性结果：皮尔逊相关系数、p值

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