在R中利用起止日期计算活跃暴发疫情数的技术求助

问题背景 #

我有一份机构呼吸道疾病暴发数据集，每条记录包含暴发通知日期、结束日期，以及COVID-19、流感、RSV的存在状态（多病原体共存为混合暴发）。暴发从通知日到结束日视为活跃状态，目标是绘制从最早通知日到当日，按病原体分类的每日活跃暴发数图表，但统计每日总活跃暴发数时遇到问题。

现有代码问题 #

第一版代码问题 #

第一版代码尝试为每条暴发生成活跃日期序列，但统计病原体分类时逻辑错误——用is.na判断其他病原体不存在，而实际数据中不存在是用0表示，导致仅计数一次而非每日计数：

ari_test <- ari_data %>%
  select(record_id, notification_date, declaration_date, c_cov_present, c_flu_present, c_rsv_present) |> 
  mutate(notification_date = as.Date(notification_date),
         declaration_date = as.Date(declaration_date)) %>%
  filter(!is.na(notification_date) & !is.na(declaration_date)) %>%
  # Generate a sequence of dates from notification_date to declaration_date for each facility
  rowwise() %>%
  mutate(date = list(seq(notification_date, declaration_date, by = "day"))) %>%
  unnest(date) %>%
  select(-notification_date, -declaration_date) %>%
  # Count the number of active outbreaks per day for each pathogen
  group_by(date) %>%
  summarise(active_covid = sum(c_cov_present == 1 & is.na(c_flu_present) & is.na(c_rsv_present)),
            active_influenza = sum(is.na(c_cov_present) & c_flu_present == 1 & is.na(c_rsv_present)),
            active_rsv = sum(is.na(c_cov_present) & is.na(c_flu_present) & c_rsv_present == 1),
            active_mixed = sum(rowSums(cbind(c_cov_present, c_flu_present, c_rsv_present), na.rm = TRUE) >= 2))

第二版代码问题 #

第二版代码用complete扩展日期时，因分组后min(notification_date)和max(declaration_date)是每组的范围，而非全局范围，导致record_id找不到的错误：

ari_test <- ari_data %>%
  mutate(notification_date = as.Date(notification_date),
         declaration_date = as.Date(declaration_date)) %>%
  filter(!is.na(notification_date) & !is.na(declaration_date)) %>%
  mutate(across(starts_with("c_"), ~if_else(is.na(.), 0, 1))) %>%  # Convert NA to 0 for presence/absence
  group_by(record_id) %>%
  mutate(active_covid = +(any(c_cov_present == 1 & is.na(c_flu_present) & is.na(c_rsv_present))),
         active_influenza = +(any(is.na(c_cov_present) & c_flu_present == 1 & is.na(c_rsv_present))),
         active_rsv = +(any(is.na(c_cov_present) & is.na(c_flu_present) & c_rsv_present == 1)),
         active_mixed = +(any(rowSums(select(., starts_with("c_"))) >= 2))) %>%
  complete(record_id, date = seq.Date(min(notification_date), max(declaration_date), by = "day"), fill = list(active_covid = 0, active_influenza = 0, active_rsv = 0, active_mixed = 0)) %>%
  ungroup()

修正后的解决方案 #

核心思路：

• 先统一将病原体状态的NA转为0，明确0表示不存在，1表示存在
• 为每条暴发生成完整的活跃日期序列
• 对每条暴发先判定其属于哪类（单一COVID、单一流感、单一RSV、混合）
• 按日期分组统计各类活跃暴发数

# 加载所需包
library(tidyverse)

# 处理数据
ari_test <- ari_data %>%
  # 转换日期格式并过滤有效记录
  mutate(notification_date = as.Date(notification_date),
         declaration_date = as.Date(declaration_date)) %>%
  filter(!is.na(notification_date) & !is.na(declaration_date)) %>%
  # 将NA转为0，统一状态标识
  mutate(across(starts_with("c_"), ~replace_na(., 0))) %>%
  # 为每条暴发生成活跃日期序列
  rowwise() %>%
  mutate(date = list(seq(notification_date, declaration_date, by = "day"))) %>%
  unnest(date) %>%
  # 判定每条暴发的类型
  mutate(
    total_pathogens = c_cov_present + c_flu_present + c_rsv_present,
    outbreak_type = case_when(
      total_pathogens >= 2 ~ "mixed",
      c_cov_present == 1 ~ "covid",
      c_flu_present == 1 ~ "influenza",
      c_rsv_present == 1 ~ "rsv",
      TRUE ~ "none" # 无病原体的情况，可根据实际需求处理
    )
  ) %>%
  # 过滤无病原体的记录（如果有的话）
  filter(outbreak_type != "none") %>%
  # 按日期和暴发类型分组计数，再转为宽格式
  group_by(date, outbreak_type) %>%
  summarise(count = n(), .groups = "drop") %>%
  pivot_wider(names_from = outbreak_type, values_from = count, values_fill = 0) %>%
  # 重命名列以匹配需求
  rename(
    active_covid = covid,
    active_influenza = influenza,
    active_rsv = rsv,
    active_mixed = mixed
  )

验证示例数据 #

使用提供的示例数据测试：

# 示例数据
ari_data <- structure(list(record_id = c(1, 2, 5, 6, 7, 8, 10, 11, 12, 13
), notification_date = structure(c(19523, 19524, 19535, 19535, 
19535, 19535, 19536, 19536, 19542, 19542), class = "Date"), declaration_date = structure(c(19544, 
19537, 19548, 19559, 19542, 19555, 19548, 19549, 19550, 19569
), class = "Date"), c_cov_present = c(1, 1, 1, 1, 0, 1, 1, 1, 
1, 0), c_flu_present = c(1, 0, 0, 0, 
0, 0, 0, 0, 0, 1), 
    c_rsv_present = c(0, 0, 0, 0, 1, 0, 0, 1, 0, 0)), row.names = c(NA, 
-10L), class = c("tbl_df", "tbl", "data.frame"))

# 运行修正后的代码
ari_test <- ari_data %>%
  mutate(notification_date = as.Date(notification_date),
         declaration_date = as.Date(declaration_date)) %>%
  filter(!is.na(notification_date) & !is.na(declaration_date)) %>%
  mutate(across(starts_with("c_"), ~replace_na(., 0))) %>%
  rowwise() %>%
  mutate(date = list(seq(notification_date, declaration_date, by = "day"))) %>%
  unnest(date) %>%
  mutate(
    total_pathogens = c_cov_present + c_flu_present + c_rsv_present,
    outbreak_type = case_when(
      total_pathogens >= 2 ~ "mixed",
      c_cov_present == 1 ~ "covid",
      c_flu_present == 1 ~ "influenza",
      c_rsv_present == 1 ~ "rsv",
      TRUE ~ "none"
    )
  ) %>%
  filter(outbreak_type != "none") %>%
  group_by(date, outbreak_type) %>%
  summarise(count = n(), .groups = "drop") %>%
  pivot_wider(names_from = outbreak_type, values_from = count, values_fill = 0) %>%
  rename(
    active_covid = covid,
    active_influenza = influenza,
    active_rsv = rsv,
    active_mixed = mixed
  )

# 查看结果
head(ari_test)

运行后会得到每日各类型活跃暴发数，可直接用于绘图。

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