R中基于多指标的分层抽样：疫苗组与未疫苗组匹配抽样实现

问题描述 #

现有数据集包含两组人群：

• vaccinated（接种疫苗）组：每行对应唯一ID及唯一T0，无重复ID；
• unvaccinated（未接种疫苗）组：长格式数据，同一ID对应多行不同T0；
• 每行包含三个核心指标：PCP_visits（PCP就诊次数）、specialty_visits（专科就诊次数）、lab_visits（实验室就诊次数）。

需求：在R中完成抽样，最终数据需满足每行对应唯一ID，且两组上述三个指标的均值尽可能接近。

示例数据 #

以下是用于测试的模拟数据集：

set.seed(123)
data <- data.frame(
  ID = c(1:10, rep(11:20, 4)), # 接种组ID唯一，未接种组ID重复
  group = c(rep("vaccinated", 10), rep("unvaccinated", 40)),
  T0 = rep(1:10, 5),
  PCP_visits = c(sample(0:10, 10, replace = TRUE), 
                 sample(3:13, 40, replace = TRUE)),
  specialty_visits = c(sample(0:5, 10, replace = TRUE),
                       sample(2:7, 40, replace = TRUE)),
  lab_visits = c(sample(0:8, 10, replace = TRUE),
                 sample(2:10, 40, replace = TRUE))
)

数据预览：

ID        group T0 PCP_visits specialty_visits lab_visits
1   1   vaccinated  1          2                0          3
2   2   vaccinated  2          2                5          0
3   3   vaccinated  3          9                4          5
4   4   vaccinated  4          1                0          2
5   5   vaccinated  5          5                1          7
6   6   vaccinated  6         10                3          2
7   7   vaccinated  7          4                3          7
8   8   vaccinated  8          3                5          0
9   9   vaccinated  9          5                5          6
10 10   vaccinated 10          8                2          6
11 11 unvaccinated  1          9                5          6
12 12 unvaccinated  2         10                5          5
13 13 unvaccinated  3          4                0          6
14 14 unvaccinated  4          2                5          4
15 15 unvaccinated  5         10                1          5
16 16 unvaccinated  6          8                0          7
17 17 unvaccinated  7          8                1          4
18 18 unvaccinated  8          8                3          6
19 19 unvaccinated  9          2                4          3
20 20 unvaccinated 10          7                4          2
21 11 unvaccinated  1          9                5          8
22 12 unvaccinated  2          6                2          6
23 13 unvaccinated  3          9                0          5
24 14 unvaccinated  4          8                3          8
25 15 unvaccinated  5          2                5          6
26 16 unvaccinated  6          3                0          1
27 17 unvaccinated  7          0                5          2
28 18 unvaccinated  8         10                0          7
29 19 unvaccinated  9          6                2          3
30 20 unvaccinated 10          4                5          6
31 11 unvaccinated  1          9                3          3
32 12 unvaccinated  2          6                0          0
33 13 unvaccinated  3          8                5          7
34 14 unvaccinated  4          8                5          3
35 15 unvaccinated  5          9                2          8
36 16 unvaccinated  6          6                5          7
37 17 unvaccinated  7         10                4          5
38 18 unvaccinated  8          4                2          3
39 19 unvaccinated  9          6                5          7
40 20 unvaccinated 10          4                1          2
41 11 unvaccinated  1         10                4          3
42 12 unvaccinated  2          5                4          3
43 13 unvaccinated  3          8                2          5
44 14 unvaccinated  4          1                1          0
45 15 unvaccinated  5          4                1          3
46 16 unvaccinated  6          7                1          8
47 17 unvaccinated  7          1                3          6
48 18 unvaccinated  8          0                1          7
49 19 unvaccinated  9          8                1          4
50 20 unvaccinated 10         10                5          1

抽样前组间特征 #

抽样前两组指标均值存在明显差异：

# 计算抽样前组间特征
library(dplyr)
data %>% 
  group_by(group) %>% 
  summarise(n = n(),
            PCP_visits = mean(PCP_visits),
            specialty_visits = mean(specialty_visits),
            lab_visits = mean(lab_visits))

结果：

# A tibble: 2 × 5
  group            n PCP_visits specialty_visits lab_visits
  <chr>        <int>      <dbl>            <dbl>      <dbl>
1 unvaccinated    40        8.9             4.25       5.72
2 vaccinated      10        4.6             2.3        3.3 

解决方案 #

步骤1：预处理未接种组 #

未接种组每个ID对应多条记录，先聚合到ID层面，计算每个ID的三个指标均值：

# 聚合未接种组到ID级
unvaccinated_agg <- data %>%
  filter(group == "unvaccinated") %>%
  group_by(ID) %>%
  summarise(
    PCP_visits = mean(PCP_visits),
    specialty_visits = mean(specialty_visits),
    lab_visits = mean(lab_visits),
    .groups = "drop"
  )

# 提取接种组数据（已为ID唯一）
vaccinated_data <- data %>%
  filter(group == "vaccinated") %>%
  select(ID, group, PCP_visits, specialty_visits, lab_visits)

方法1：多次抽样筛选（无需额外包） #

通过多次抽样未接种ID，找到与接种组均值差异最小的样本：

set.seed(123)
n_vaccinated <- nrow(vaccinated_data)
best_diff <- Inf
best_sample <- NULL

# 定义函数计算两组均值差异总和
calculate_diff <- function(sampled_ids, agg_data, target_data) {
  sampled <- agg_data %>% filter(ID %in% sampled_ids)
  target_means <- target_data %>% summarise(across(c(PCP_visits, specialty_visits, lab_visits), mean))
  sampled_means <- sampled %>% summarise(across(c(PCP_visits, specialty_visits, lab_visits), mean))
  sum(abs(as.numeric(sampled_means) - as.numeric(target_means)))
}

# 尝试1000次抽样，获取最优样本
for (i in 1:1000) {
  sampled_ids <- sample(unvaccinated_agg$ID, n_vaccinated)
  current_diff <- calculate_diff(sampled_ids, unvaccinated_agg, vaccinated_data)
  if (current_diff < best_diff) {
    best_diff <- current_diff
    best_sample <- sampled_ids
  }
}

# 从原始数据中提取最优样本（每个ID取第一行）
sampled_unvaccinated <- data %>%
  filter(group == "unvaccinated", ID %in% best_sample) %>%
  group_by(ID) %>%
  slice(1) %>%
  ungroup()

# 合并最终数据
final_data <- bind_rows(vaccinated_data, sampled_unvaccinated)

方法2：马氏距离匹配（精准匹配，需MatchIt包） #

利用马氏距离为每个接种ID匹配最相似的未接种ID，确保组间均值高度接近：

# 安装并加载MatchIt包
# install.packages("MatchIt")
library(MatchIt)

# 准备匹配数据集
match_data <- bind_rows(
  vaccinated_data %>% mutate(treatment = 1),
  unvaccinated_agg %>% mutate(group = "unvaccinated", treatment = 0)
) %>%
  select(treatment, PCP_visits, specialty_visits, lab_visits)

# 1:1马氏距离匹配
m.out <- matchit(treatment ~ PCP_visits + specialty_visits + lab_visits,
                 data = match_data, method = "nearest", distance = "mahalanobis")

# 获取匹配后的未接种ID
matched_unvaccinated_ids <- unvaccinated_agg$ID[m.out$match.matrix[,1]]

# 提取匹配后的未接种组数据（每个ID取第一行）
matched_unvaccinated <- data %>%
  filter(group == "unvaccinated", ID %in% matched_unvaccinated_ids) %>%
  group_by(ID) %>%
  slice(1) %>%
  ungroup()

# 合并最终数据
final_data_matched <- bind_rows(vaccinated_data, matched_unvaccinated)

验证抽样结果 #

以方法2为例，验证抽样后组间特征：

final_data_matched %>%
  group_by(group) %>%
  summarise(
    n = n(),
    PCP_visits = mean(PCP_visits),
    specialty_visits = mean(specialty_visits),
    lab_visits = mean(lab_visits)
  )

结果示例（因随机种子可能略有差异）：

# A tibble: 2 × 5
  group            n PCP_visits specialty_visits lab_visits
  <chr>        <int>      <dbl>            <dbl>      <dbl>
1 unvaccinated    10        4.7             2.4        3.5
2 vaccinated      10        4.6             2.3        3.3

两组均值已高度接近，满足需求。

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