如何基于含出生/死亡日期的数据集绘制带病例数和人年的Lexis图？

问题背景 #

现有个体级数据集，每行对应一名个体，包含出生日期与死亡日期。已筛选出2005-2010年间去世、年龄0-18岁的个体，目标：

• 绘制以年龄为Y轴、日历时间为X轴的Lexis图
• 在每个Lexis三角区域标注病例数（该区域内的死亡人数）和人年数
• 不清楚如何将个体级数据适配Epi包的N2Y函数，且参考示例未覆盖病例数统计

示例数据集 #

library(dplyr)

set.seed(40)

df <- data.frame(
  Id = seq(1, 100, by = 1),
  DateOfBirth = sample(seq(as.Date("1990-01-01"), as.Date("2007-01-01"), by = "days"), 100, replace = TRUE),
  DateOfDeath = sample(seq(as.Date("2005-01-01"), as.Date("2010-12-31"), by = "days"), 100, replace = TRUE)
)

df_filtered <- df %>%
  mutate(AgeAtDeath = round(as.numeric(DateOfDeath - DateOfBirth) / 365.25, 1)) %>%
  filter(AgeAtDeath >= 0 & AgeAtDeath <= 18)

print(df_filtered)

解决方案 #

1. 加载依赖包 #

除dplyr外，还需Epi处理Lexis图和人年计算、lubridate简化日期操作：

library(Epi)
library(dplyr)
library(lubridate)

2. 将个体级数据转换为N2Y要求的格式 #

N2Y需要**按年龄组(A)、日历年份(P)分组的人口数(N)**数据集，因此先将个体数据聚合到对应单元格：

# 提取年份与年龄组信息
df_prep <- df_filtered %>%
  mutate(
    # 死亡年份（日历时间P）
    P = year(DateOfDeath),
    # 整数年龄组（与示例格式对齐）
    A = floor(AgeAtDeath)
  )

# 统计每个(A, P)单元格的个体数，补充缺失单元格避免报错
nx_df <- df_prep %>%
  count(A, P, name = "N") %>%
  complete(A = 0:18, P = 2005:2010, fill = list(N = 0))

3. 用N2Y计算Lexis三角区域的人年数 #

运行N2Y生成每个三角区域的人年数据：

# 计算三角区域人年，返回数据框格式
nt_df <- N2Y(data = nx_df, return.dfr = TRUE)

4. 统计每个三角区域的病例数 #

根据死亡时间在年份中的位置，将死亡个体分配到对应Lexis三角并计数：

# 为每个死亡个体匹配所属三角区域
df_triangles <- df_prep %>%
  mutate(
    # 死亡日期在当年的占比
    prop_year = yday(DateOfDeath) / ifelse(leap_year(P), 366, 365),
    # 判断归属三角：上半年归下三角，下半年归上三角
    triangle_type = ifelse(prop_year <= 0.5, "bottom", "top"),
    # 三角区域的坐标（用于后续标注）
    tri_A = ifelse(triangle_type == "bottom", A + 0.25, A + 0.75),
    tri_P = ifelse(triangle_type == "bottom", P + 0.25, P + 0.75)
  )

# 统计每个三角区域的病例数
case_counts <- df_triangles %>%
  count(tri_A, tri_P, name = "Cases")

5. 绘制Lexis图并标注数据 #

绘制基础Lexis图，添加人年数和病例数标注：

# 绘制Lexis网格图
Lexis.diagram(
  age = c(0, 19), 
  date = c(2005, 2011), 
  int = 1, 
  coh.grid = TRUE,
  main = "Lexis图：0-18岁死亡个体（2005-2010）"
)

# 标注人年数（红色）
with(nt_df, text(P, A, formatC(Y, format = "f", digits = 1), col = "red", cex = 0.8))

# 标注病例数（蓝色）
with(case_counts, text(tri_P, tri_A, paste("Cases:", Cases), col = "blue", cex = 0.7, pos = 4))

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