R语言计算二元变量存在/缺失对score的贡献值

测试数据集构建 #

使用如下代码生成测试用分析数据集：

set.seed(123)

var1 = sample(0:1, 10000, replace=T)
var2 = sample(0:1, 10000, replace=T)
var3 = sample(0:1, 10000, replace=T)
var4 = sample(0:1, 10000, replace=T)
score = rnorm(10000,10,5)

my_data = data.frame(var1,var2, var3,var4, score)

协变量组合计数统计 #

可通过data.table语法快速统计4个二元变量所有唯一取值组合的样本量，4个二元变量共存在16种取值组合：

dt = my_data[,c(1,2,3,4)]
setDT(dt)[,list(Count=.N) ,names(dt)]

运行后输出计数结果如下：

var1 var2 var3 var4 Count
 1:    0    0    0    0   667
 2:    0    1    0    0   601
 3:    1    1    1    1   651
 4:    0    1    1    1   608
 5:    1    0    1    1   613
 6:    1    1    0    1   588
 7:    0    1    1    0   607
 8:    0    0    1    1   607
 9:    1    0    1    0   625
10:    0    1    0    1   661
11:    1    1    1    0   635
12:    0    0    1    0   640
13:    1    1    0    0   608
14:    1    0    0    0   607
15:    0    0    0    1   626
16:    1    0    0    1   656

计算需求说明 #

核心计算逻辑：在其余协变量取值固定的前提下，计算目标二元变量取1（存在）和取0（缺失）时结果变量score的均值差值，即为该变量在对应场景下对score的贡献值。
计算示例：

• Var4在其余变量全为1场景下的贡献：(var1=1,var2=1,var3=1,var4=1)组score均值 减去 (var1=1,var2=1,var3=1,var4=0)组score均值
• Var2在其余变量全为1场景下的贡献：(var1=1,var2=1,var3=1,var4=1)组score均值 减去 (var1=1,var2=0,var3=1,var4=1)组score均值
• 所有场景逻辑通用，例如Var4在其余变量全为0场景下的贡献，可直接对比(1,0,0,0)和(1,0,0,1)两组的score均值差得到。

手动硬编码逐个子集筛选的实现方式重复度高、效率低，示例代码如下：

var1_present <- my_data[which(my_data$var1 == 1 & my_data$var2 == 1 & my_data$var3 == 1 & my_data$var4 == 1 ), ]
var1_present_score = mean(var1_present$score)

var1_absent <- my_data[which(my_data$var1 == 0 & my_data$var2 == 1 & my_data$var3 == 1 & my_data$var4 == 1 ), ]
var1_absent_score = mean(var1_absent$score)

var_1_contribution = var1_present_score - var1_absent_score
# 输出结果：[1] 0.1288283

通用实现方案 #

实现思路：预计算所有协变量组合对应的score均值生成映射表，再通过参数传入目标变量和固定协变量取值，自动匹配两组均值做差，无需重复编写子集筛选逻辑。

函数定义 #

# 加载data.table用于高效分组计算
library(data.table)

# 转换原始数据为data.table格式
setDT(my_data)

# 预计算所有协变量组合的score均值，生成查询映射表
combo_mean <- my_data[, .(mean_score = mean(score)), by = .(var1, var2, var3, var4)]

#' 计算指定场景下目标变量对score的贡献值
#' @param target_var 字符串，待计算贡献的目标变量名，例如"var1"
#' @param fixed_cov 命名向量，其余协变量的固定取值，例如c(var2=1,var3=1,var4=1)
#' @return 包含目标变量名、固定协变量场景、贡献值的数据框
calc_contribution <- function(target_var, fixed_cov) {
  # 构造存在组（目标变量=1）、缺失组（目标变量=0）的筛选条件
  present_cond <- c(setNames(1, target_var), fixed_cov)
  absent_cond <- c(setNames(0, target_var), fixed_cov)
  
  # 从预计算映射表中匹配两组均值
  present_mean <- combo_mean[as.list(present_cond), on = names(present_cond)]$mean_score
  absent_mean <- combo_mean[as.list(absent_cond), on = names(absent_cond)]$mean_score
  
  # 返回结构化结果
  return(data.frame(
    target_var = target_var,
    fixed_scenario = paste(names(fixed_cov), fixed_cov, sep = "=", collapse = ","),
    contribution = present_mean - absent_mean
  ))
}

单场景调用示例 #

计算其余变量全为1时var1的贡献，和手动硬编码结果完全一致：

calc_contribution(
  target_var = "var1",
  fixed_cov = c(var2=1, var3=1, var4=1)
)

输出：

target_var             fixed_scenario contribution
1       var1 var2=1,var3=1,var4=1    0.1288283

批量计算所有场景贡献 #

如果需要一次性输出所有固定协变量场景下4个变量的贡献值，可通过循环批量生成：

all_vars <- c("var1", "var2", "var3", "var4")
result_list <- list()
idx <- 1

for (target in all_vars) {
  # 提取目标变量之外的其余协变量
  other_vars <- setdiff(all_vars, target)
  # 生成其余协变量所有可能的0/1取值组合
  other_combos <- expand.grid(rep(list(0:1), length(other_vars)), KEEP.OUT.ATTRS = FALSE)
  names(other_combos) <- other_vars
  
  # 遍历每个组合计算贡献
  for (i in seq_len(nrow(other_combos))) {
    fixed_val <- as.numeric(other_combos[i, ])
    names(fixed_val) <- other_vars
    result_list[[idx]] <- calc_contribution(target, fixed_val)
    idx <- idx + 1
  }
}

# 合并为完整结果表，总计32条记录（4个变量，每个变量对应8种固定协变量场景）
all_contributions <- rbindlist(result_list)

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