在R中按组匹配而不替换

Question

我有一个包含约 20,000 个病例和约 50,000 个对照的数据集。每个案例都使用 SQL 中的连接过程与 3 个对照进行了初步匹配，该过程要求对某些协变量（例如，该人是否在 2019 年参加研究的二元变量）进行精确匹配，并允许在一个范围内匹配其他协变量（例如，在 +/- 25 天内登记的总天数）。因此，一个控件可以被列为多个案例的可能控件。

在 R 中，我想从每个案例的 3 个可能匹配项中为每个案例选择控件，同时确保两件事：1）控件不重复，2）对于没有完全匹配的案例，控件在协变量上尽可能地匹配。

我认为我需要的方法是按 case_id 对数据进行分组，并从每个组中选择 control_id，其中协变量之间的差异最小化。问题是我需要确保一旦为一组中的一个案例选择了一个控件，它就不能用作第二个案例的控件，无论协变量匹配有多好。

这是一些示例数据：

case_id <- c(1,1,1,2,2,2,3,3,3,4,4,4)
control_id <- c(5,14,7,8,9,10,11,12,13,14,15,5)
case_enrolled_2019 <- c(1,1,1,0,0,0,0,0,0,1,1,1)
control_enrolled_2019 <- c(1,1,1,0,0,0,0,0,0,1,1,1)
case_enrollment_count <- c(1,1,1,3,3,3,1,1,1,1,1,1)
control_enrollment_count <- c(2,1,1,3,3,3,1,1,1,1,1,2)
case_enrolled_days <- c(300,300,300,200,200,200,600,600,600,300,300,300)
control_enrolled_days <- c(300,300,280,200,200,210,601,610,598,300,301,300)

cbind(case_id, control_id, case_enrolled_2019, control_enrolled_2019, case_enrollment_count, control_enrollment_count, case_enrolled_days, control_enrolled_days)

我需要这样的输出：

case_id <- c(1,2,3,4)
control_id <- c(14,8,11,15)
case_enrolled_2019 <- c(1,0,0,1)
control_enrolled_2019 <- c(1,0,0,1)
case_enrollment_count <- c(1,3,1,1)
control_enrollment_count <- c(1,3,1,1)
case_enrolled_days <- c(300,200,600,300)
control_enrolled_days <- c(300,200,601,301)

cbind(case_id, control_id, case_enrolled_2019, control_enrolled_2019, case_enrollment_count, control_enrollment_count, case_enrolled_days, control_enrolled_days)

Answer 1

这里有一个解决方案。它很慢，因为它必须遍历每个案例，计算它与块中尚未匹配的每个控制单元之间的距离，然后选择最近的。使用包含 72000 行的此数据集的复制版本，我的计算机需要 183 秒。您可以使用并行化来加快速度。

第一步是创建一个数据集，其中每个唯一行（包括最初与每次处理的多个匹配的控制单元）都有一个条目，就好像它们是它们自己的单元一样。可能有一种更优雅、更通用的方法来做到这一点，但老实说，它在这里工作得很好。

 df2 <- setNames(as.data.frame(rbind(cbind(1, unique(df[,c(1,1,3,5,7)])), 
                                     cbind(0, df[,c(1,2,4,6,8)]))),
                 c("cc", "block", "id", "enrolled_2019", "enrollment_count", "enrolled_days"))

head(df2)
#>   cc block id enrolled_2019 enrollment_count enrolled_days
#> 1  1     1  1             1                1           300
#> 2  1     2  2             0                3           200
#> 3  1     3  3             0                1           600
#> 4  1     4  4             1                1           300
#> 5  0     1  5             1                2           300
#> 6  0     1 14             1                1           300

这给我们留下了一个有 6 列的数据集：第一个 cc 是单位是案例 (1) 还是控制 (0)。第二个block 是每个单元所属的初始匹配块。第三个id 是单位的ID。其他的是变量值。

接下来，我们初始化要匹配的案例和控件的 ID，并在df2 中添加一列，其中包含每个单元是否匹配。这对于确保对应于同一控件 ID 的两行不匹配非常重要。

cases <- unique(df2$id[df2$cc == 1])
controls <- rep(NA, length(cases))
df2$matched <- FALSE

接下来我们遍历每个案例并找到它的匹配项。在循环中，我们首先确保该案例有任何可用的控制单元。如果不是，我们跳过这个案例，它仍然是不匹配的。如果是这样，我们使用optmatch 包中的match_on 计算它与块中剩余控制单元之间的距离。这只是为每个案例创建一个距离矩阵，为每个控件创建一个列。在这种情况下，我们有一个案例，最多三个控件。我们取距离最小的值，找到它对应的控件ID，并将其存储在controls向量中。我们还确保将 df2 中具有相同 ID 的任何行标记为“匹配”，这样控制单元就不会被重复使用。

for (i in seq_along(cases)) {
    if (any(df2$block==cases[i] & !df2$matched)) {
        dist <- optmatch::match_on(cc ~ enrolled_2019 + enrollment_count + enrolled_days, 
                                   data = df2[df2$block==cases[i] & !df2$matched,])
        controls[i] <- df2$id[as.numeric(colnames(dist)[which.min(dist)])]
        df2$matched[df2$id == controls[i]] <- TRUE
    }
}

最后，我们从匹配的控件中提取对。我们使用na.omit，以防有不匹配的情况被丢弃。

(pairs <- na.omit(data.frame(cases, controls)))
#>   cases controls
#> 1     1       14
#> 2     2        8
#> 3     3       11
#> 4     4       15

matched_df <- df[interaction(df[,1], df[,2]) %in% interaction(pairs[,1], pairs[,2]),]

^{由reprex package (v0.3.0) 于 2020 年 5 月 28 日创建}