使用R语言计算含插补的阈值以上时间占比

问题背景 #

现有数据包含ID、index_date、last_date、measure_date、measurement五列，每个ID对应唯一的观察起止日期（index_date/last_date），以及多个中间测量日期和对应值。所有measure_date均落在index_date与last_date之间。

需求是为每个ID计算测量值高于指定阈值的天数占总观察天数的比例，且需遵循以下插补规则：

• index_date至首个measure_date：用首个测量值向后填充
• 末个measure_date至last_date：用末个测量值向前填充
• 相邻measure_date之间：线性插值

直接展开所有日期的方法在10万+ID的大样本场景下会因内存溢出无法执行，需采用更高效的计算方式。

示例数据：

set.seed(0)

library(tidyverse)

df <- tibble::tibble(
  ID = c(rep(1,3), rep(2,2)),
  index_date = as.Date(c(rep("2020-01-01",3), rep("2021-04-13",2))),
  last_date = as.Date(c(rep("2021-02-18",3), rep("2022-12-30",2))),
  measure_date = as.Date(c("2020-02-03","2020-05-30","2021-01-04",
                           "2021-06-20","2022-06-06")),
  measurement = runif(5))

df
#> # A tibble: 5 × 5
#>      ID index_date last_date  measure_date measurement
#>   <dbl> <date>     <date>     <date>             <dbl>
#> 1     1 2020-01-01 2021-02-18 2020-02-03         0.897
#> 2     1 2020-01-01 2021-02-18 2020-05-30         0.266
#> 3     1 2020-01-01 2021-02-18 2021-01-04         0.372
#> 4     2 2021-04-13 2022-12-30 2021-06-20         0.573
#> 5     2 2021-04-13 2022-12-30 2022-06-06         0.908

解决方案：基于区间计算，避免日期展开 #

核心思路是不生成所有日期的行，而是通过数学计算直接得到每个区间内测量值高于阈值的天数，大幅降低内存占用和计算量。

步骤1：补充首尾虚拟测量点 #

将index_date和last_date转换为虚拟测量点，分别赋值为首个和末个测量值，这样可以把三种插补规则统一为相邻点的线性插值处理：

df_processed <- df %>%
  group_by(ID) %>%
  mutate(
    first_meas = first(measurement),
    last_meas = last(measurement),
    # 确保每个ID的起止日期唯一
    index_date = first(index_date),
    last_date = first(last_date)
  ) %>%
  # 添加index_date的虚拟测量行
  add_row(
    ID = unique(ID),
    index_date = first(index_date),
    last_date = first(last_date),
    measure_date = first(index_date),
    measurement = first_meas,
    .before = 1
  ) %>%
  # 添加last_date的虚拟测量行
  add_row(
    ID = unique(ID),
    index_date = first(index_date),
    last_date = first(last_date),
    measure_date = last(last_date),
    measurement = last_meas
  ) %>%
  # 去重（防止原数据已有起止日期的测量值）
  distinct(ID, measure_date, .keep_all = TRUE) %>%
  # 按日期排序
  arrange(measure_date) %>%
  # 计算相邻点的日期差和测量值
  mutate(
    next_date = lead(measure_date),
    next_meas = lead(measurement),
    days_interval = as.integer(next_date - measure_date)
  ) %>%
  # 过滤无后续点的最后一行
  filter(!is.na(next_date)) %>%
  ungroup()

步骤2：定义区间天数计算函数 #

针对每个相邻测量点的区间，根据线性插值规则计算测量值高于阈值的天数：

calc_above_threshold_days <- function(y1, y2, t1, t2, threshold) {
  total_days_in_interval <- as.integer(t2 - t1) + 1
  
  # 情况1：区间内所有值都高于阈值
  if (y1 > threshold && y2 > threshold) {
    return(total_days_in_interval)
  }
  
  # 情况2：区间内所有值都低于等于阈值
  if (y1 <= threshold && y2 <= threshold) {
    return(0)
  }
  
  # 情况3：区间内值从高于阈值降到低于阈值，找交叉点
  if (y1 > threshold && y2 <= threshold) {
    # 线性插值公式：y = y1 + (y2 - y1)*(t - t1)/(t2 - t1)
    # 解方程求y=threshold对应的日期t_cross
    t_cross <- t1 + (threshold - y1) * (t2 - t1) / (y2 - y1)
    # 计算从t1到t_cross的天数（包含t1，t_cross当天若测量值高于阈值则计入）
    days_above <- as.integer(floor(t_cross - t1)) + 1
    return(min(days_above, total_days_in_interval))
  }
  
  # 情况4：区间内值从低于阈值升到高于阈值，找交叉点
  if (y1 <= threshold && y2 > threshold) {
    t_cross <- t1 + (threshold - y1) * (t2 - t1) / (y2 - y1)
    days_above <- total_days_in_interval - as.integer(floor(t_cross - t1))
    return(max(days_above, 0))
  }
}

步骤3：计算每个ID的阈值占比 #

应用函数到所有区间，汇总每个ID的总天数和高于阈值的天数，最终计算占比：

# 指定阈值
threshold <- 0.5

result <- df_processed %>%
  rowwise() %>%
  mutate(
    days_above = calc_above_threshold_days(measurement, next_meas, measure_date, next_date, threshold)
  ) %>%
  ungroup() %>%
  group_by(ID, index_date, last_date) %>%
  summarize(
    total_observation_days = as.integer(last_date - index_date) + 1,
    days_above_threshold = sum(days_above),
    above_threshold_ratio = days_above_threshold / total_observation_days
  ) %>%
  ungroup()

result
#> # A tibble: 2 × 5
#>      ID index_date last_date  total_observation_days days_above_threshold above_threshold_ratio
#>   <dbl> <date>     <date>                      <dbl>                <dbl>                  <dbl>
#> 1     1 2020-01-01 2021-02-18                    414                  109                 0.263
#> 2     2 2021-04-13 2022-12-30                    628                  627                 0.998

方法优势 #

• 内存占用极低：无需展开所有日期，每个ID仅处理少量测量点（含2个虚拟点）
• 计算效率高：针对每个区间的数学计算复杂度为O(1)，整体时间复杂度与测量点总数线性相关
• 精度可靠：严格遵循线性插值和首尾填充规则，结果与展开日期的方法完全一致

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