因果森林（R包grf）-治疗异质性探索

一、概念介绍

传统方法的局限：一个“平均”的答案：

在做很多分析时，最常问的问题是：“这个新干预/措施到底有没有用？”

传统的因果分析方法，比如我们熟知的A/B测试，能很好地回答这个问题。它会给我们一个平均处理效应（Average Treatment Effect, ATE）。比如，在A干预相较于B干预死亡率降低2.2%。

随着个体化医疗和精准医学相关研究的进展，研究者越来越关心研究效应的异质性处理效应(treatment efectheterogeneity)。即由于个体异质性的存在, 采用相同药物或治疗方法的患者往往会在治疗效果上产生差异,这也使得将临床试验结果直接转化到个体患者变得困难，并非所有患者的治疗效果都与参加试验患者的平均效果相似, 因此对个体处理效应的估计变得尤为重要。

因此，只知道“平均效果”是远远不够的。我们更想知道：这个措施，到底对“谁”最有效？对“谁”没用，甚至有反作用？

因果森林：

因果森林是一种基于随机森林的机器学习算法，它的核心目标不是进行“预测”，而是估计“干预”（Treatment）对不同个体产生的“因果效应”，特别是找出这种效应在不同人群中的差异（即异质性处理效应）。

一些概念：

因果森林可以估计处理或干预对结果变量的平均因果效应（ATE）以及个体化的因果效应（ITE）。
# ITE: Individual Treatment Effect (个体处理效应) 干预（Treatment）对于某一个特定个体 i 的真实效果
# ATE: Average Treatment Effect (平均处理效应) 在整个目标群体中，处理效应的平均值。它不考虑任何个体或亚群体的特征。
# CATE: Conditional Average Treatment Effect (条件平均处理效应) 在给定一组协变量（特征）X=x 的条件下，处理效应的平均值。

二、文献理解（以下面几篇文献为例，仅列举可参考进行分析的内容）

文献1：利用因果森林估计异质性人群下个体的处理效应（https://html.rhhz.net/zhlxbx/20190620.htm）

数据情况：
- W=1代表处理组,即使用RHC,W=0代表对照组,为未使用RHC。
- X：协变量
- Y：180d生存结局

分析内容：比较干预对个体处理效应的影响/分布

个体处理效应估计r(x)可以解释为在给定协变量下使用RHC 和未使用RHC180d生存结局之间差异的条件期望值,当结局变量Y为二分类时,个体处理效应解释为个体接受处理与未接受处理时结局的率差。

个体处理效应预测为负值意味着患者使用RHC时死亡率比较小,因为RHC的使用降低了相同协变量情况下该患者的死亡可能性,效应值为正则为未使用RHC时180d时死亡可能性较小。

结果解释：上图横轴代表个体，纵坐标为个体处理效应预测值。可以看出：预测值为正的个体占大多数，即使用RHC时死亡可能性增高的个体占绝大多数，从另一个方面说，使用RHC会使该样本人群180 d死亡率增加。同时该结果显示，约22.9%的个体效应值<0，表明使用RHC对这部分个体起到了保护作用，使用RHC时死亡率最高可降低0.0655（在所有患者中，RHC治疗能带来的最佳效果就是使某个患者的死亡率降低 0.0655（即降低 6.55%）。）.

文献2：Targeting weight loss interventions to reduce cardiovascular complications of type 2 diabetes: a machine learning-based post-hoc analysis of heterogeneous treatment effects in the Look AHEAD trial

数据情况：
- W：W=1为健康饮食和增加体力活动干预组；W=0为对照组
- X：协变量
- Y： 复合终点（首次发生心血管原因死亡、非致命性心肌梗死、非致命性中风或因心绞痛住院等）

分析内容：因果森林模型指导下的亚型异质性分析-基于特征重要性分析筛选出的协变量

因果森林模型揭示，基线 HbA1c 和总体健康状况（通过 SF-36 健康调查自我报告）是区分从强化减肥干预中获益高与低的主要协变量。

进一步，基于因果树推断的协变量cutoff进行亚组划分。

（亚组分析，Cox分析）针对划分出的四个亚组，比如，在“HbA1c低且健康评分低”的这部分人和“剩余的”那部分人中，分别控制协变量后，比较干预组和对照组，事件的发生率（或绝对风险降低百分比），以及干预治疗是否显著增加了风险（HR值和p值）。

文献3：Estimating individual treatment effects on COPD exacerbations by causal machine learning on randomised controlled trials ( https://mp.weixin.qq.com/s/thHt6L97qcZk05D6J1e8eg)

数据情况：
- W: FF/VI吸入剂 vs 安慰剂
- X: 协变量
- Y: 急性加重

分析内容：因果森林模型指导下的亚型异质性分析-基于个体处理效应的分位数分组

文献2是基于协变量分组讨论治疗异质性，文献3是直接基于个体处理效应的分位数（Q1-Q5）分组，探索治疗异质性（同时对每个分组比较比较基线差异）

上图结果表明Q1亚群的患者对治疗是显著有效的。同时，作者对不同亚型（Q1-Q5）中治疗组(FFNI)和安慰剂组的基线特征进行分析，表明患者的基线特征(如年龄、肺功能等)没有显著差异，说明了各组内的疗效差异不是因为随机分配不均造成的，而是真实的。

临床解释：从Q1到Q5 患者的肺功能(FEV1)是越来越好的(呈上升趋势)，而患者既往的急性加重次数和药物使用情况是越来越轻的(呈下降趋势)。表明模型预测的“疗效最好”的Q1组，恰恰是那些基线时肺功能最差、病情最重的患者。

分析内容：特征重要性分析

三、代码示例（主要参考其他文章，特异性分析还需要进一步敲码）

set.seed(39)  # 设定随机数
rm(list = ls()) library(grf)          #  causal_forest()
library(ggplot2)
library(patchwork)    #  Combine ggplots
library(hstats)       #  Friedman's H, PDP
library(kernelshap)   #  General SHAP
library(shapviz)      #  SHAP plots
library(dplyr)
library(tidyverse)# 此处假设我们的数据为一个数据框df，包含下述列：
# X：特征变量（协变量）
# W：处理变量（通常是0或1）
# Y：结果变量，连续变量或分类变量；当结局变量为二分类时，个体处理效应解释为个体接受处理与未接受处理时结局的率差；W <- as.integer(df$treat) #
table(W)Y <- as.numeric(df$OS_time)
mean(Y)xvars <- c('性别','年龄','吸烟年支','糖尿病史')
X<-df[,xvars]
X <- as_tibble(X)
head(X)

Step1：模型拟合

# 用grf包来拟合因果森林，这是一个试图估计条件平均处理效果(CATE) E[Y(1) - Y(0) | X = X]的树集合。因此，它可以用于研究治疗效果的不均匀性。
fit <- causal_forest(X = X,Y = Y,W = W,num.trees = 1000,mtry = 4,sample.fraction = 0.7,seed = 1,ci.group.size = 1,
)

Step2：个体处理效应（ITE）估计

## 使用因果森林来估计个体处理效应（ITE）。
ite_estimates <- predict(fit)$predictions
df$ITE <- ite_estimates
head(df,2)# 绘制ITE的直方图
ggplot(df, aes(x = ITE)) +geom_histogram(binwidth = 0.05, fill = "grey", color = "black", alpha = 0.5) +labs(title = "Individual Treatment Effect Estimates",x = "Estimated Individual Treatment Effect (ITE)",y = "Frequency") +theme_minimal()

Step2：Average Treatment Effect (平均处理效应) 估计

# 计算ATE
ate_result <- average_treatment_effect(fit)
estimate <- ate_result[1]
std_err <- ate_result[2]# --- 计算95%置信区间 ---
confidence_level <- 0.95
alpha <- 1 - confidence_level
critical_value <- qnorm(1 - alpha / 2) # 计算误差范围
margin_of_error <- critical_value * std_err# 计算置信区间的上下限
ci_lower <- estimate - margin_of_error
ci_upper <- estimate + margin_of_errorcat("Estimate Average Treatment Effect (ATE):", estimate, "\n")
cat("Standard Error:", std_err, "\n")
cat(confidence_level * 100, "% Confidence Interval: [", ci_lower, ", ", ci_upper, "]\n", sep="")# 计算p-value (双边检验)
# 计算Z-score
z_score <- estimate / std_err# pnorm(-abs(z_score)) 计算的是单边的p值，乘以2得到双边的p值
p_value <- 2 * pnorm(-abs(z_score))
cat("P-value:", p_value, "\n")

Step3：模型解释-grf包（自带函数）

变量的重要性

因果森林的可变重要性可以通过每个特征被用来分割的相对计数来衡量。

imp <- sort(setNames(variable_importance(fit), xvars))
par(mai = c(0.7, 2, 0.2, 0.2))
barplot(imp, horiz = TRUE, las = 1, col = "orange")

主要效应

为了研究对CATE的主要影响，我们考虑了部分依赖图(PDP)。这样的图显示了平均预测如何依赖于一个特征的值，保持所有其他特征值不变(可能是不自然的)。

par(mai = c(0.7, 0.1, 0.8, 0.2)) # 下左上右pred_fun <- function(object, newdata, ...) {predict(object, newdata, ...)$predictions
}
pdps <- lapply(xvars, function(v) plot(partial_dep(fit, v, X = X, pred_fun = pred_fun)))
wrap_plots(pdps, guides = "collect", ncol = 3) &
#   ylim(c(-0.11, -0.06)) &ylab("Treatment effect")
#当1为治疗组，预测术后并发症发生，负值意味着更强的(积极的)治疗效果。   

Step3：模型解释-shapviz包（SHAP）

用瀑布图解释第一例患者特征值对应的CATE

# Explaining one CATE
kernelshap(fit, X = X[1, ], bg_X = X, pred_fun = pred_fun) |> shapviz() |> sv_waterfall() +xlab("Prediction")# Explaining all CATEs globally
system.time(  # 13 minks <- kernelshap(fit, X = X, pred_fun = pred_fun)  
)
shap_values <- shapviz(ks)

sv_importance(shap_values)
sv_importance(shap_values, kind = "bee")

sv_dependence(shap_values, v = xvars) +plot_layout(ncol = 2)  

Step4：亚型的平均ATE估计

library(grf)
library(dplyr)# 1. 计算特征重要性，并排序
var_imp <- variable_importance(fit)
importance_scores <- data.frame(Feature = xvars,Importance = var_imp
)
sorted_importance <- importance_scores[order(importance_scores$Importance, decreasing = TRUE), ]cat("--- 特征重要性排名 ---\n")
print(sorted_importance)if (nrow(sorted_importance) < 2) {stop("特征数量少于3个，无法选择排名第一和第三的特征。")
}feature1_name <- as.character(sorted_importance$Feature[1])
feature3_name <- as.character(sorted_importance$Feature[2])feature1_id <- which(xvars == feature1_name)
feature3_id <- which(xvars == feature3_name)cat(paste("\n将使用以下两个特征进行交叉分组:\n"))
cat(paste(" - 排名第一:", feature1_name))
cat(paste(" - 排名第二:", feature3_name))# 2: 分别提取这两个特征的分割点# 定义一个函数来提取指定特征的分割点，避免代码重复
get_split_points_for_feature <- function(fit, feature_id) {points <- c()for (i in 1:fit$`_num_trees`) {tree <- get_tree(fit, i)for (node in tree$nodes) {if (!is.null(node$split_variable) && node$split_variable == feature_id) {points <- c(points, node$split_value)}}}return(sort(unique(points)))
}# 提取特征的分割点
split_points1 <- get_split_points_for_feature(fit, feature1_id)
split_points3 <- get_split_points_for_feature(fit, feature3_id)cat(paste("\n为特征 '", feature1_name, "' 找到了", length(split_points1), "个独特的分割点。\n"))
cat(paste("为特征 '", feature3_name, "' 找到了", length(split_points3), "个独特的分割点。\n"))# 3: 根据两个特征的中位数分割点创建交叉分组# 4. 计算两个特征各自的中位数分割值
if (length(split_points1) == 0 || length(split_points3) == 0) {stop("其中一个或两个选定的特征从未被用作分割，无法进行交叉分组。")
}
median_split1 <- median(split_points1)
median_split3 <- median(split_points3)cat(paste("\n使用 '", feature1_name, "' 的中位数分割点:", round(median_split1, 4), "\n"))
cat(paste("使用 '", feature3_name, "' 的中位数分割点:", round(median_split3, 4), "\n"))# 5.在数据框 df 中创建两个临时的分组列
df$group1 <- ifelse(df[[feature1_name]] <= median_split1, "Low", "High")
df$group3 <- ifelse(df[[feature3_name]] <= median_split3, "Low", "High")# 6. 将两个临时分组结合成一个最终的交叉分组列
df$cross_group <- paste0(feature1_name, "_", df$group1, " & ",feature3_name, "_", df$group3
)
df$cross_group <- as.factor(df$cross_group)cat("\n--- 交叉分组后各组的样本数量 ---\n")
print(table(df$cross_group))# 7. 分析并比较各交叉分组的平均处理效应 (ATE)  #  为每个交叉分组计算 ATE
group_names <- levels(df$cross_group)
ate_results_list <- list()for (group_name in group_names) {# 获取当前组的行索引subset_indices <- which(df$cross_group == group_name)if (length(subset_indices) > 0) {# 使用 average_treatment_effect 计算该子集的 ATEate_est <- average_treatment_effect(fit, subset = subset_indices)ate_results_list[[group_name]] <- ate_est} else {cat("该组没有样本，跳过分析。\n")}
}if (length(ate_results_list) > 0) {results_df <- do.call(rbind, lapply(names(ate_results_list), function(name) {data.frame(Group = name,ATE_Estimate = ate_results_list[[name]][1],Std_Error = ate_results_list[[name]][2])}))cat("\n--- 各交叉分组的平均处理效应 (ATE) 及统计显著性检验结果 ---\n")# Z-score = (估计值 - 假设值) / 标准误。这里假设值 H0: ATE = 0。results_df$Z_Score <- results_df$ATE_Estimate / results_df$Std_Error# 计算双边检验的 P-value# P-value = 2 * P(Z > |z_score|)results_df$P_Value <- 2 * pnorm(-abs(results_df$Z_Score))# 添加一个显著性判断列 (使用常见的 alpha = 0.05)results_df$Significance <- ifelse(results_df$P_Value < 0.05, "显著 (p < 0.05)", "不显著 (p >= 0.05)")# 打印最终的、包含P值和显著性判断的表格results_df$Z_Score <- round(results_df$Z_Score, 3)results_df$P_Value <- format.pval(results_df$P_Value, eps = 0.001, digits = 3)    
}
results_df

四、参考

方法：

https://mp.weixin.qq.com/s/PI_QkxEe4igDH9P3Rf5PiA

https://mp.weixin.qq.com/s/NefDA3h-JJka9ysWpmcuvg

写作或文章参考：

https://mp.weixin.qq.com/s/thHt6L97qcZk05D6J1e8eg

https://mp.weixin.qq.com/s/6IBq6gHAE4A4uumvFKH_Eg

利用因果森林估计异质性人群下个体的处理效应

Targeting weight loss interventions to reduce cardiovascular complications of type 2 diabetes: a machine learning-based post-hoc analysis of heterogeneous treatment effects in the Look AHEAD trial