在R中基于队列研究数据计算疾病年度/月度发病率（含人月/人年）

针对你的队列数据，我们可以按以下步骤计算2014年1月1日至2021年8月31日的发病率：

1. 加载工具包并处理日期格式 #

首先把字符型日期转成R可识别的日期类型，用lubridate处理日期，dplyr做数据操作：

library(dplyr)
library(lubridate)

# 读取示例数据
df <- data.frame(patid=c("1","2","3","4","5","6","7"), 
                 start_date=c("01/03/1993","24/03/2000", 
                              "01/03/2020","24/03/2016", 
                              "24/03/2001","01/03/1999", 
                              "01/03/2016"), 
                 end_date=c("31/08/2021","31/08/2021", 
                            "23/08/2021","01/08/2019", 
                            "17/08/2020","04/08/2014", 
                            "31/08/2018"), 
                 disease=c("yes","no","yes","no", 
                           "no","yes","yes"), 
                 disease_date=c("15/11/2017",NA, 
                                "15/08/2020",NA,NA, 
                                "01/01/2014","18/03/2017") )

# 转换日期格式（原始为dd/mm/yyyy）
df <- df %>%
  mutate(
    start_date = dmy(start_date),
    end_date = dmy(end_date),
    disease_date = dmy(disease_date)
  )

2. 定义研究时段边界 #

明确分析的时间范围：

study_start <- ymd("2014-01-01")
study_end <- ymd("2021-08-31")

3. 确定每位患者的实际观察区间 #

队列研究中，患者的观察时间需考虑三个限制：

• 不能早于研究起始日
• 不能晚于研究终止日
• 若患者在研究期间发病，观察截止到发病当日（发病后不再计入暴露人群）

df <- df %>%
  mutate(
    # 实际观察起始日：取患者入组日和研究起始日的较晚者
    actual_start = pmax(start_date, study_start),
    # 实际观察终止日：分发病/未发病情况处理
    actual_end = case_when(
      disease == "yes" ~ pmin(disease_date, end_date, study_end),
      TRUE ~ pmin(end_date, study_end)
    ),
    # 排除实际起始日晚于终止日的无效记录
    valid = actual_start <= actual_end
  ) %>%
  filter(valid)

4. 拆分观察区间到年度/月度，计算人月/人年数 #

单个患者的观察时间可能跨多个时段，需拆分区间后计算对应时段的人月/人年数：

拆分到年度并计算人年数 #

# 生成年度序列
year_seq <- seq(study_start, study_end, by = "year") %>% floor_date("year")

# 拆分患者观察区间到年度并计算人年数
annual_data <- df %>%
  rowwise() %>%
  mutate(
    years = list(seq(floor_date(actual_start, "year"), floor_date(actual_end, "year"), by = "year")),
    annual_starts = list(pmax(actual_start, years)),
    annual_ends = list(pmin(actual_end, years %m+% years(1) - days(1)))
  ) %>%
  unnest(c(years, annual_starts, annual_ends)) %>%
  ungroup() %>%
  mutate(
    year = year(years),
    person_years = time_length(interval(annual_starts, annual_ends), "year") %>% round(2),
    person_months = time_length(interval(annual_starts, annual_ends), "month") %>% round(2)
  ) %>%
  group_by(year) %>%
  summarise(
    total_person_years = sum(person_years),
    total_person_months = sum(person_months)
  )

拆分到月度并计算人月数 #

# 生成月度序列
month_seq <- seq(study_start, study_end, by = "month") %>% floor_date("month")

# 拆分患者观察区间到月度并计算人月数
monthly_data <- df %>%
  rowwise() %>%
  mutate(
    months = list(seq(floor_date(actual_start, "month"), floor_date(actual_end, "month"), by = "month")),
    monthly_starts = list(pmax(actual_start, months)),
    monthly_ends = list(pmin(actual_end, months %m+% months(1) - days(1)))
  ) %>%
  unnest(c(months, monthly_starts, monthly_ends)) %>%
  ungroup() %>%
  mutate(
    year_month = format(months, "%Y-%m"),
    person_months = time_length(interval(monthly_starts, monthly_ends), "month") %>% round(2)
  ) %>%
  group_by(year_month) %>%
  summarise(total_person_months = sum(person_months))

5. 统计各时段的新发病例数 #

筛选研究期间发病的患者，按年度/月度分组统计：

# 年度发病数
annual_cases <- df %>%
  filter(disease == "yes", disease_date >= study_start, disease_date <= study_end) %>%
  mutate(year = year(disease_date)) %>%
  group_by(year) %>%
  summarise(new_cases = n())

# 月度发病数
monthly_cases <- df %>%
  filter(disease == "yes", disease_date >= study_start, disease_date <= study_end) %>%
  mutate(year_month = format(disease_date, "%Y-%m")) %>%
  group_by(year_month) %>%
  summarise(new_cases = n())

6. 计算发病率 #

合并发病数与总人年/人月数，计算以每1000人年/人月为单位的发病率：

年度发病率 #

annual_incidence <- annual_data %>%
  left_join(annual_cases, by = "year") %>%
  mutate(
    new_cases = replace_na(new_cases, 0),
    incidence_per_1000_py = (new_cases / total_person_years) * 1000 %>% round(2)
  )

月度发病率 #

monthly_incidence <- monthly_data %>%
  left_join(monthly_cases, by = "year_month") %>%
  mutate(
    new_cases = replace_na(new_cases, 0),
    incidence_per_1000_pm = (new_cases / total_person_months) * 1000 %>% round(2)
  )

结果说明 #

• annual_incidence包含各年度总人年数、新发病例数、每1000人年发病率
• monthly_incidence包含各月度总人月数、新发病例数、每1000人月发病率
• 人月数=人年数×12，可按需选择统计单位

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