隨機對照試驗是評估干預措施效果的金標準，主要用于估計研究人群整體的平均干預效果。但同一種干預對不同特征人群的效果可能存在顯著差異，即存在處理效應異質性。傳統的亞組分析和交互作用分析對處理效應異質性的檢驗功效通常較低，難以有效發現異質性來源。近年來，隨著機器學習技術的發展，因果森林成為評估處理效應異質性的新方法，可有效解決傳統分析的局限性。然而，該方法在醫學研究領域的應用尚處于起步階段。為促進該方法的合理使用，本文在醫學干預效果評價的情境下，介紹因果森林方法的用途、基本原理和實現步驟，并結合實例解讀及代碼實現樣例，討論使用因果森林方法的注意事項。
引用本文： 周文岳, 易飛, 李冰麗, 孫鳳, 楊智榮. 因果森林在醫學處理效應異質性評價中的基本原理與應用. 中國循證醫學雜志, 2023, 23(4): 485-491. doi: 10.7507/1672-2531.202212074 復制
1 前言
隨機對照試驗是評估干預效果的金標準，主要用于估計干預在研究人群整體中的平均效果，但試驗平均效果不一定適用于每一個體，同種干預對于不同特征人群的效果可能存在顯著差異，即存在處理效應異質性。如何針對處理效應異質性進行評價是推動精準醫學計劃、促進個體化醫療的關鍵。
為確定處理效應異質性的來源，研究者通常根據不同人群特征分組（一般根據預定的假設）進行亞組分析，或在回歸分析的傳統模型中加入人群特征與干預措施的交互項并檢驗其統計學顯著性。然而，每一次亞組分析通常只根據一種人群特征進行分組，若同時有多個亞組分析，可能會增加假陽性的機會[1, 2]。但當多種人群特征同時疊加才會影響干預效果時，則按單個特征分組的亞組分析難以奏效[3]。而在對交互項的分析中，當需要同時疊加多種人群特征時，交互項個數會呈指數增加，研究者手動設置這些交互項是非常費時費力的[4]。此外，相比平均效應的估計，不管是亞組分析還是交互項檢驗都需要更大樣本量的數據支持，否則難以保證足夠的檢驗效能[5]。
為克服上述傳統方法的局限性，因果森林（causal forest）作為一種新興技術開始應用于評估處理效應異質性。該方法不需要預先設定異質性來源，而是通過機器學習，幫助研究者找出傳統因果分析中被掩蓋的異質性。該算法由Atheya（2021年諾貝爾經濟學獎得主）等[6]于2016年在決策樹的基礎上提出，首次將分類回歸樹引入到因果識別框架，形成因果樹。2018年Wager等[7]進一步探討了如何使用生成的因果樹執行隨機森林的算法，即因果森林。該方法在2019年首次引進國內[8]。
本文將進一步介紹因果樹和因果森林的算法、實例解讀和軟件實現，探討因果森林在醫學干預措施的效應異質性評價中的應用。
2 基本原理
2.1 因果樹
因果樹是在決策樹的基礎上發展而來的機器學習方法。決策樹以信息增益、信息增益比、基尼系數或誤差平方和等指標作為分裂準則，將樣本層層分組[9]。為避免單一一棵決策樹過擬合的情況，可通過隨機選取變量和樣本的方式來生成許多不同的決策樹，形成隨機森林，然后由每一棵樹進行“投票”來決定最后的分組[10]。
與決策樹不同，因果樹按照某特征分組后將組間處理效應值之差最大作為分裂準則，分別估算出各組的處理效應。此時組內個體的平均處理效應就是該組個體處理效應的預測值。根據因果推斷的“反事實”原理，按特征X分組后的平均處理效應值
可以用如下等式來表示：
 |
（公式1：個體平均處理效應值）
公式1假設個體i能夠接受兩種干預，其中Xi=x代表某個人群特征x，Yi(1)表示個體i接受干預的結局，Yi(0)表示個體i不接受干預（或接受其他干預）的結局。即按照某個特征分組后，分別計算各組
的差值，再找出分組后各組差值最大的特征x作為分裂的節點。
由于實際情況下只能獲得到個體i的觀測結局Yi(1)或Yi(0)，而無法同時獲得兩者。因此在隨機對照試驗的背景下，我們一般用人群的平均處理效應（即E[Y(1)]和E[Y(0)]）來代替個體的處理效應，進而衡量兩組人群結局間差異：
 |
（公式2：人群分組平均處理效應值）
根據公式2，在構建因果樹時算法自動找出當前人群按某個特征值分成兩組后的組間處理效應值之差最大（分裂準則），并按此特征值分組；然后在各組中重復上述步驟，直到生成符合預設要求的因果樹。
2.2 因果森林
與隨機森林相比，因果森林在節點分叉、模型擬合和處理效應估計三個方面都有自己的特點[6]。對于決策樹-隨機森林模型，我們可以用驗證集對其進行驗證，使用如均方誤差之類的指標來度量預測的準確性。但在因果樹-因果森林模型中，由于無法觀察到真正的處理效應，因而無法實現如決策樹-隨機森林模型那樣的驗證。為此，需要通過檢查偏差、標準誤差和相關估計的可信區間來評估因果森林的準確性。為實現這種統計推斷，Atheya和Wager兩人引入了誠實（honesty）算法[7]。
誠實算法要求樹的分裂和
的估計分開進行，即每次從訓練集中隨機抽取一定比例的樣本用于構建單個樹模型，將抽樣后的樣本分為兩部分，一部分用于估計樹模型參數（分裂子樣本），另一部分用于計算節點內部的處理效應
（估計子樣本），這樣做可以減少對
的估計誤差。估計子樣本填充在分裂子樣本擬合的樹中，根據末端葉子節點中估計子樣本的E[Y(1)]和E[Y(0)]之間的差異來對處理效應做預測。在構建完森林模型后，根據特征變量在森林的每個層次上進行分裂的次數，計算每個特征變量的重要性并根據重要性排序。最后，基于因果森林中的最優樹模型，在驗證集中采用傳統因果分析方法對不同亞組的處理效應進行估計。
因果森林實現步驟如下：
① 隨機抽樣：將原始數據集隨機分成訓練集和驗證集，然后無放回地再從訓練集里隨機抽取樣本（一般比例為50%，可根據原始樣本量進行調整）作為分裂子樣本，將剩下部分作為估計子樣本。由于訓練集（包括分裂子樣本和估計子樣本）和驗證集是隨機抽樣形成的，理論上均可代表研究的源人群。
② 生成因果樹和因果森林：先用分裂子樣本的數據，以特征變量X基于遞歸分區的方式構造因果樹模型，即從根節點開始自頂向下對樣本進行劃分。利用以下公式計算估計子樣本每個葉子節點上所有個體的平均處理效應值，即：
 |
（公式3：平均處理效應估計值）
其中公式3的前半部分表示在葉子節點L中干預變量W為1時所有個體的結局變量Y的均值；后半部分表示葉子節點L中干預變量W為0時所有的個體結局變量Y的均值。然后基于處理效應值之差最大化的分裂準則選擇特征值x后，將父節點分裂為Xi<=x和Xi>x左右兩個子節點。例如，在考慮按哪個年齡閾值來分組時，≤30歲的人群處理效應值為10，>30歲的人群處理效應值為5，效應值之差為5；而≤40歲的人群處理效應值為6，>40歲的人群處理效應值為3，效應值之差為3；那么算法會優先考慮以30歲作為分類閾值，以此類推。將子節點按照相同的準則繼續分割，直到最終不再生成新的節點為止。同時我們需要確定模型的處理效應估計值的準確性，即各個組內的個體處理效應值和組內平均處理效應值的均方誤差最小化，使得組內效應估計方差最小化，組間異質性最大化。
重復上述步驟B次，最終形成有B棵樹的因果森林，此時第i個個體的處理效應估計值為綜合B棵樹的均值進行計算，公式為：
 |
（公式4：個體平均處理效應值的估計值）
③ 生成最佳模型：根據因果森林的特征重要性和處理效應異質性，算法會自動選擇出最能代表整個森林的處理效應異質性的一棵最優樹。
④ 估計各組處理效應：根據最優樹模型的每個終端節點把驗證集的人群劃分成相應的亞組，利用傳統回歸模型來估計每個人群亞組的平均處理效應并作組間比較，以此驗證因果森林的劃分結果。若組間比較結果具有統計學差異，提示處理效應可能存在異質性。
上述原理適用于隨機對照試驗的數據，這類數據在干預組和對照組的人群特征上具有可比性，即使按因果森林進行分組后，每一個亞組內部仍然保留原有的隨機分組。若是采用觀察性研究的數據，還需對潛在混雜因素進行校正[6, 11]。
3 實例分析
下面以Look AHEAD（Action for Health in Diabetes）研究為例進行解讀[12-13]。這是一項在美國16個臨床中心開展的隨機對照試驗，將超重或肥胖的2型糖尿病患者隨機分配到干預組或對照組，干預組通過促進健康飲食和增加身體活動來減輕體重，對照組接受糖尿病常規支持及教育，主要結局是心血管復合事件，包括心血管死亡、非致命性心肌梗死、非致命性卒中和心絞痛住院等。該研究發現，與糖尿病常規支持和教育相比，在總體人群中進行生活方式強化干預并未降低2型糖尿病患者的心血管疾病和死亡風險[12]。考慮到傳統的亞組分析不一定能有效發現處理效應異質性，可能會掩蓋某些人群亞組的獲益或風險，因此研究者基于該試驗的數據，利用因果森林進行了二次分析[13]。
該二次分析先將總體樣本分為訓練集和驗證集，然后利用訓練集建立因果森林模型，再按照模型的分類將驗證集分為四類人群，最后用傳統回歸分析驗證干預措施在不同的人群亞組中的效果是否相同（圖1）。
圖1
基于Look AHEAD實驗的數據進行因果森林分析原理*
*：亞組3和4均未表現出明顯異質性，原文未提供具體ARR值；HbA1C：糖化血紅蛋白；SF-36：SF-36健康調查簡表。
在因果森林分析中，研究者一共納入84個基線因素作為特征變量，涵蓋了社會人口學特征、疾病史、實驗室檢查、行為特征四個方面，具體包括糖化血紅蛋白（glycated hemoglobin，HbA1c）水平、SF-36健康調查簡表（36-item short form survey，SF-36）、年齡、種族、吸煙史、體重等。4 901名患者被隨機分為訓練集（2 450名）和驗證集（2 451名），分別用于訓練模型與驗證模型。由于數據來自隨機對照試驗，訓練集和驗證集內部以及每一個亞組內部仍符合干預措施的隨機分配原則。
具體步驟包括：
① 隨機抽樣：無放回地隨機選擇50%的數據作為訓練集，將另外50%作為驗證集。
② 生成因果樹和因果森林：在訓練集中，根據特征變量對處理效應的影響大小（處理效應的均方誤差最小）分裂結點，最大限度地增大不同組之間的干預與發生主要結局風險的關聯性差異，以此得到1 000棵不同的因果樹。
③ 生成最佳模型：從1 000棵因果樹中根據特征重要性和處理效應異質性選出擬合最佳的因果樹，以此作為分組依據。
④ 估計各組處理效應：按照因果森林的結果分組，在驗證集里采用傳統Cox回歸分析得出：整體ARR值為1.89［95%CI（?1.08，4.85），P=0.212］，HR值為0.87［95%CI（0.71，1.05），P=0.15］。可以看出對于被隨機分到驗證集的2 451名患者而言，干預組與對照組之間的主要結局發生風險在整體上并沒有統計學差異。而研究者通過因果森林訓練出的模型將驗證集分為4個亞組后發現：亞組1（基線HbA1c<6.8%且基線SF-36得分<48）的ARR值為?7.41［95%CI（?14.22，?0.60），P=0.033］，HR值為1.99［95%CI（1.06，3.75），P=0.03］，亞組1以外其余人群的ARR值為3.46［95%CI（0.21，6.73），P=0.038］，HR值為0.78［95%CI（0.63，0.96），P=0.02］，干預措施和亞組特征存在交互作用（PInteraction=0.006）；亞組2（基線HbA1c<6.8%且基線SF-36得分≥48）的ARR值為8.22［95%CI（2.17，14.27），P=0.007］，HR值為0.55［95%CI（0.35，0.86），P=0.01］，亞組2以外其余人群的ARR值為0.05［95%CI（?3.35，3.45），P=0.977］，HR值為0.97［95%CI（0.78，1.20），P=0.78］，干預措施和亞組特征之間同樣存在交互作用（PInteraction=0.025）。
由此可見，對于整體人群來說，與糖尿病常規支持和教育相比，減重干預并未降低超重或肥胖2型糖尿病患者的心血管疾病和死亡風險。但對于HbA1c<6.8%且SF-36得分<48的患者來說，減重干預反而使心血管疾病和死亡的風險增加；而對于HbA1c<6.8%且SF-36得分≥48的患者來說，減重干預能使心血管疾病和死亡的風險降低。然而，這一重要結果并沒有在Look AHEAD試驗的傳統亞組分析中發現[12]。通過因果森林，研究者不再需要手動調節交互項或進行多次亞組分析就可以找出這些具有處理效應異質性的人群，從而針對不同人群采取不同的干預措施。
除了Look AHEAD試驗以外，因果森林在其他臨床干預領域也有成功應用的研究案例。比如，在SPRINT試驗結果顯示，強化降壓能預防心血管事件，傳統亞組分析并未發現異質性亞組的存在[14]，但因果森林結果發現，對于基線時收縮壓>144 mmHg且吸煙的患者來說強化降壓反而使其心血管事件風險上升[15]。近年來，因果森林也開始應用于觀察性研究。比如，在MESA前瞻性隊列研究中，研究者在傾向評分匹配的基礎上構建因果森林，發現冠狀動脈鈣化與心血管疾病風險之間的關聯存在異質性，男性、西班牙裔、具有心血管危險因素的冠狀動脈鈣患者是心血管疾病的高危人群[16]。
4 代碼實現
本案例采用表1的隨機對照試驗模擬數據，利用R軟件進行因果森林建模分析，代碼及解釋如下。
表1
隨機對照試驗模擬數據庫的結構（只截取前六名患者數據）
隨機對照試驗模擬數據覆蓋3 582例腦卒中患者，結局是癡呆（dement），干預是他汀用藥（statin），協變量包括連續性指標的年齡（age）、體質指數（bmi），二分類指標的性別（gender）、慢性阻塞性肺病（copd）和冠心病（mi）。
setwd('路徑') # 設置工作路徑
# install.packages("R包的名字") # 首先安裝要用到的R包
#導入包
library(scorecard) # 用于訓練和驗證集分割
library(tidyverse) # 數據處理
library(grf) # 廣義隨機森林，需要grf0.10.2版本或以上
# 讀取數據，本案例所使用的數據文件為'statin.csv'，格式如表1所示
data <- read.csv("statin.csv")
data <- na.omit(data) # 刪除有缺失值的個體數據
data <- data.frame(lapply(data, function(x) as.numeric(as.character(x))))
# 將數據集按1∶1劃分訓練集和驗證集
data_split <- split_df(data, ratio = 0.5)
train_data <- data_split$train # 訓練集（50%，用于因果森林建模）
test_data <- data_split$test # 驗證集（50%，用于邏輯回歸驗證）
### 因果森林建模 ###
W <- train_data$statin # 干預W（如表1，本例中干預變量名稱為statin）
Y <- as.numeric(train_data$dement) # 結局Y（如表1，本例中結局變量名稱為dement）
covars <- c("age", "bmi", "copd", "gender", "mi") # 協變量
X <- select(train_data, covars) # 從訓練集中選擇協變量
# 生成因果森林
set.seed(42)
cf <- causal_forest(X, Y, W,
num.trees = 2000,
sample.fraction=0.5, # 構建每棵樹時從訓練集中抽樣的比例
min.node.size = 30,
honesty = TRUE,
honesty.fraction = 0.5, # 分裂子樣本的比例
tune.parameters = c("mtry", "alpha", "imbalance.penalty",
"honesty.prune.leaves"))
cf$tunable.params # 查看因果森林參數
# 查看變量重要性排序
var_imp <- c(variable_importance(cf))
names(var_imp) <- covars
sort(var_imp, decreasing=TRUE)
# 計算平均處理效應值
average_treatment_effect(cf, target.sample = "treated")
# ---------- 補充grf函數: find_best_tree ---------- #
# 由于版本更新，請從以下網址復制源碼至此處并運行
# https://gist.github.com/msegar/c15af2bfc067e5319eecfd67d785fde5
# ------------------------------------------------- #
# 在森林中找到最佳樹
best_trees <- find_best_tree(cf, "causal")
plot(tree <- get_tree(cf, best_trees$best_tree))
# 估計CATE的最佳線性擬合，可通過系數評估模型的性能
# 若mean.forest.prediction=1，則森林產生的平均預測是正確的；
# 若differential.forest.prediction = 1，則森林的異質性已被校準。
test_calibration(cf)
# 使用驗證集進行預測
X_test <- select(test_data, covars)
test.hat <- predict(cf, newdata = X_test, estimate.variance=TRUE)$predictions
### 邏輯回歸驗證 ###
# 根據最佳樹的分裂條件將驗證集分成4個亞組
group1 <- test_data[test_data$bmi <= 28.43 & test_data$bmi <= 24.79]
group2 <- test_data[test_data$bmi <= 28.43 & test_data$bmi > 24.79]
group3 <- test_data[test_data$bmi > 28.43 & test_data$age <= 61]
group4 <- test_data[test_data$bmi > 28.43 & test_data$age > 61]
grps <- list(group1, group2, group3, group4)
logistic_res <- matrix(1, 4, 3)
# 擬合邏輯回歸模型
for (i in 1:4){
glm_fit <- glm(dement ~ statin, data =grps[[i]], family ='binomial')
# print(glm_fit)
logistic_res[i, ] <- exp(cbind(OR = coef(glm_fit),confint(glm_fit)))[2, 1:3]
}
logistic_res # 查看擬合結果
通過運行上述代碼，可繪制出因果森林中最優的分組結果（圖2）。
圖2
模擬數據的因果森林最佳樹分組結果
圖中訓練樣本被切割為4個亞組，分別為
① BMI≤28.43且BMI≤24.79；
② BMI≤28.43且BMI>24.79；
③ BMI>28.43且年齡≤61歲；
④ BMI>28.43且年齡>61歲；
然后對四組進行Logistic回歸，在驗證集中計算出四組OR值（他汀組與非他汀組腦卒中后癡呆的比值比）分別為0.80［95%CI（0.48，1.31）］、0.66［95%CI（0.44，0.99）］、1.30［95%CI（0.83，2.05）］、0.33［95%CI（0.19，0.56）］。從結果可以看出，③組的處理效應相較其他亞組存在明顯異質性，即提示對于BMI>28.43且年齡≤61歲的患者來說，他汀治療可能無效，甚至存在潛在風險。
5 總結與展望
處理效應異質性是醫學實踐中關注的熱點問題，也是與精準醫學密切相關的重要問題。在異質性人群的發現和個體治療效果的預測方面，因果森林具有重要意義和應用價值[17]。
作為機器學習與因果推斷相結合的方法，因果森林可彌補傳統亞組分析和交互項檢驗的局限性。它通過數據驅動的方式識別各種特征組合的亞組，而不需要事先預設分組特征，在盡可能保證各類亞組人群都能被覆蓋的同時也不會增加假陽性的概率。因果森林不僅適用于隨機對照試驗，也適用于觀察性研究，但觀察性研究在構建和驗證因果森林時都需校正可能的混雜因素，比如采用傾向評分匹配，以滿足無混雜的假設[6, 11]。因果森林一般作為探索性分析，其結果需要后續研究的進一步驗證。
除了本文所介紹的因果森林以外，還有其他更復雜的機器學習算法也能用于處理效應異質性的評價，比如，貝葉斯加性回歸樹、神經網絡等[18]。這些新型方法比傳統方法的實施難度大，而且存在報告不規范的問題，在使用這些機器學習方法開展相關研究時，需要多學科交叉團隊的緊密合作。未來隨著統計方法學的日趨完善，機器學習與因果推斷相結合的方法會越來越多地應用于醫學處理效應異質性的評價[19]，推動精準醫學的發展。
1 前言
隨機對照試驗是評估干預效果的金標準，主要用于估計干預在研究人群整體中的平均效果，但試驗平均效果不一定適用于每一個體，同種干預對于不同特征人群的效果可能存在顯著差異，即存在處理效應異質性。如何針對處理效應異質性進行評價是推動精準醫學計劃、促進個體化醫療的關鍵。
為確定處理效應異質性的來源，研究者通常根據不同人群特征分組（一般根據預定的假設）進行亞組分析，或在回歸分析的傳統模型中加入人群特征與干預措施的交互項并檢驗其統計學顯著性。然而，每一次亞組分析通常只根據一種人群特征進行分組，若同時有多個亞組分析，可能會增加假陽性的機會[1, 2]。但當多種人群特征同時疊加才會影響干預效果時，則按單個特征分組的亞組分析難以奏效[3]。而在對交互項的分析中，當需要同時疊加多種人群特征時，交互項個數會呈指數增加，研究者手動設置這些交互項是非常費時費力的[4]。此外，相比平均效應的估計，不管是亞組分析還是交互項檢驗都需要更大樣本量的數據支持，否則難以保證足夠的檢驗效能[5]。
為克服上述傳統方法的局限性，因果森林（causal forest）作為一種新興技術開始應用于評估處理效應異質性。該方法不需要預先設定異質性來源，而是通過機器學習，幫助研究者找出傳統因果分析中被掩蓋的異質性。該算法由Atheya（2021年諾貝爾經濟學獎得主）等[6]于2016年在決策樹的基礎上提出，首次將分類回歸樹引入到因果識別框架，形成因果樹。2018年Wager等[7]進一步探討了如何使用生成的因果樹執行隨機森林的算法，即因果森林。該方法在2019年首次引進國內[8]。
本文將進一步介紹因果樹和因果森林的算法、實例解讀和軟件實現，探討因果森林在醫學干預措施的效應異質性評價中的應用。
2 基本原理
2.1 因果樹
因果樹是在決策樹的基礎上發展而來的機器學習方法。決策樹以信息增益、信息增益比、基尼系數或誤差平方和等指標作為分裂準則，將樣本層層分組[9]。為避免單一一棵決策樹過擬合的情況，可通過隨機選取變量和樣本的方式來生成許多不同的決策樹，形成隨機森林，然后由每一棵樹進行“投票”來決定最后的分組[10]。
與決策樹不同，因果樹按照某特征分組后將組間處理效應值之差最大作為分裂準則，分別估算出各組的處理效應。此時組內個體的平均處理效應就是該組個體處理效應的預測值。根據因果推斷的“反事實”原理，按特征X分組后的平均處理效應值可以用如下等式來表示：
|
（公式1：個體平均處理效應值）
公式1假設個體i能夠接受兩種干預，其中Xi=x代表某個人群特征x，Yi(1)表示個體i接受干預的結局，Yi(0)表示個體i不接受干預（或接受其他干預）的結局。即按照某個特征分組后，分別計算各組的差值，再找出分組后各組差值最大的特征x作為分裂的節點。
由于實際情況下只能獲得到個體i的觀測結局Yi(1)或Yi(0)，而無法同時獲得兩者。因此在隨機對照試驗的背景下，我們一般用人群的平均處理效應（即E[Y(1)]和E[Y(0)]）來代替個體的處理效應，進而衡量兩組人群結局間差異：
|
（公式2：人群分組平均處理效應值）
根據公式2，在構建因果樹時算法自動找出當前人群按某個特征值分成兩組后的組間處理效應值之差最大（分裂準則），并按此特征值分組；然后在各組中重復上述步驟，直到生成符合預設要求的因果樹。
2.2 因果森林
與隨機森林相比，因果森林在節點分叉、模型擬合和處理效應估計三個方面都有自己的特點[6]。對于決策樹-隨機森林模型，我們可以用驗證集對其進行驗證，使用如均方誤差之類的指標來度量預測的準確性。但在因果樹-因果森林模型中，由于無法觀察到真正的處理效應，因而無法實現如決策樹-隨機森林模型那樣的驗證。為此，需要通過檢查偏差、標準誤差和相關估計的可信區間來評估因果森林的準確性。為實現這種統計推斷，Atheya和Wager兩人引入了誠實（honesty）算法[7]。
誠實算法要求樹的分裂和的估計分開進行，即每次從訓練集中隨機抽取一定比例的樣本用于構建單個樹模型，將抽樣后的樣本分為兩部分，一部分用于估計樹模型參數（分裂子樣本），另一部分用于計算節點內部的處理效應（估計子樣本），這樣做可以減少對的估計誤差。估計子樣本填充在分裂子樣本擬合的樹中，根據末端葉子節點中估計子樣本的E[Y(1)]和E[Y(0)]之間的差異來對處理效應做預測。在構建完森林模型后，根據特征變量在森林的每個層次上進行分裂的次數，計算每個特征變量的重要性并根據重要性排序。最后，基于因果森林中的最優樹模型，在驗證集中采用傳統因果分析方法對不同亞組的處理效應進行估計。
因果森林實現步驟如下：
① 隨機抽樣：將原始數據集隨機分成訓練集和驗證集，然后無放回地再從訓練集里隨機抽取樣本（一般比例為50%，可根據原始樣本量進行調整）作為分裂子樣本，將剩下部分作為估計子樣本。由于訓練集（包括分裂子樣本和估計子樣本）和驗證集是隨機抽樣形成的，理論上均可代表研究的源人群。
② 生成因果樹和因果森林：先用分裂子樣本的數據，以特征變量X基于遞歸分區的方式構造因果樹模型，即從根節點開始自頂向下對樣本進行劃分。利用以下公式計算估計子樣本每個葉子節點上所有個體的平均處理效應值，即：
|
（公式3：平均處理效應估計值）
其中公式3的前半部分表示在葉子節點L中干預變量W為1時所有個體的結局變量Y的均值；后半部分表示葉子節點L中干預變量W為0時所有的個體結局變量Y的均值。然后基于處理效應值之差最大化的分裂準則選擇特征值x后，將父節點分裂為Xi<=x和Xi>x左右兩個子節點。例如，在考慮按哪個年齡閾值來分組時，≤30歲的人群處理效應值為10，>30歲的人群處理效應值為5，效應值之差為5；而≤40歲的人群處理效應值為6，>40歲的人群處理效應值為3，效應值之差為3；那么算法會優先考慮以30歲作為分類閾值，以此類推。將子節點按照相同的準則繼續分割，直到最終不再生成新的節點為止。同時我們需要確定模型的處理效應估計值的準確性，即各個組內的個體處理效應值和組內平均處理效應值的均方誤差最小化，使得組內效應估計方差最小化，組間異質性最大化。
重復上述步驟B次，最終形成有B棵樹的因果森林，此時第i個個體的處理效應估計值為綜合B棵樹的均值進行計算，公式為：
|
（公式4：個體平均處理效應值的估計值）
③ 生成最佳模型：根據因果森林的特征重要性和處理效應異質性，算法會自動選擇出最能代表整個森林的處理效應異質性的一棵最優樹。
④ 估計各組處理效應：根據最優樹模型的每個終端節點把驗證集的人群劃分成相應的亞組，利用傳統回歸模型來估計每個人群亞組的平均處理效應并作組間比較，以此驗證因果森林的劃分結果。若組間比較結果具有統計學差異，提示處理效應可能存在異質性。
上述原理適用于隨機對照試驗的數據，這類數據在干預組和對照組的人群特征上具有可比性，即使按因果森林進行分組后，每一個亞組內部仍然保留原有的隨機分組。若是采用觀察性研究的數據，還需對潛在混雜因素進行校正[6, 11]。
3 實例分析
下面以Look AHEAD（Action for Health in Diabetes）研究為例進行解讀[12-13]。這是一項在美國16個臨床中心開展的隨機對照試驗，將超重或肥胖的2型糖尿病患者隨機分配到干預組或對照組，干預組通過促進健康飲食和增加身體活動來減輕體重，對照組接受糖尿病常規支持及教育，主要結局是心血管復合事件，包括心血管死亡、非致命性心肌梗死、非致命性卒中和心絞痛住院等。該研究發現，與糖尿病常規支持和教育相比，在總體人群中進行生活方式強化干預并未降低2型糖尿病患者的心血管疾病和死亡風險[12]。考慮到傳統的亞組分析不一定能有效發現處理效應異質性，可能會掩蓋某些人群亞組的獲益或風險，因此研究者基于該試驗的數據，利用因果森林進行了二次分析[13]。
該二次分析先將總體樣本分為訓練集和驗證集，然后利用訓練集建立因果森林模型，再按照模型的分類將驗證集分為四類人群，最后用傳統回歸分析驗證干預措施在不同的人群亞組中的效果是否相同（圖1）。
圖1
基于Look AHEAD實驗的數據進行因果森林分析原理*
*：亞組3和4均未表現出明顯異質性，原文未提供具體ARR值；HbA1C：糖化血紅蛋白；SF-36：SF-36健康調查簡表。
在因果森林分析中，研究者一共納入84個基線因素作為特征變量，涵蓋了社會人口學特征、疾病史、實驗室檢查、行為特征四個方面，具體包括糖化血紅蛋白（glycated hemoglobin，HbA1c）水平、SF-36健康調查簡表（36-item short form survey，SF-36）、年齡、種族、吸煙史、體重等。4 901名患者被隨機分為訓練集（2 450名）和驗證集（2 451名），分別用于訓練模型與驗證模型。由于數據來自隨機對照試驗，訓練集和驗證集內部以及每一個亞組內部仍符合干預措施的隨機分配原則。
具體步驟包括：
① 隨機抽樣：無放回地隨機選擇50%的數據作為訓練集，將另外50%作為驗證集。
② 生成因果樹和因果森林：在訓練集中，根據特征變量對處理效應的影響大小（處理效應的均方誤差最小）分裂結點，最大限度地增大不同組之間的干預與發生主要結局風險的關聯性差異，以此得到1 000棵不同的因果樹。
③ 生成最佳模型：從1 000棵因果樹中根據特征重要性和處理效應異質性選出擬合最佳的因果樹，以此作為分組依據。
④ 估計各組處理效應：按照因果森林的結果分組，在驗證集里采用傳統Cox回歸分析得出：整體ARR值為1.89［95%CI（?1.08，4.85），P=0.212］，HR值為0.87［95%CI（0.71，1.05），P=0.15］。可以看出對于被隨機分到驗證集的2 451名患者而言，干預組與對照組之間的主要結局發生風險在整體上并沒有統計學差異。而研究者通過因果森林訓練出的模型將驗證集分為4個亞組后發現：亞組1（基線HbA1c<6.8%且基線SF-36得分<48）的ARR值為?7.41［95%CI（?14.22，?0.60），P=0.033］，HR值為1.99［95%CI（1.06，3.75），P=0.03］，亞組1以外其余人群的ARR值為3.46［95%CI（0.21，6.73），P=0.038］，HR值為0.78［95%CI（0.63，0.96），P=0.02］，干預措施和亞組特征存在交互作用（PInteraction=0.006）；亞組2（基線HbA1c<6.8%且基線SF-36得分≥48）的ARR值為8.22［95%CI（2.17，14.27），P=0.007］，HR值為0.55［95%CI（0.35，0.86），P=0.01］，亞組2以外其余人群的ARR值為0.05［95%CI（?3.35，3.45），P=0.977］，HR值為0.97［95%CI（0.78，1.20），P=0.78］，干預措施和亞組特征之間同樣存在交互作用（PInteraction=0.025）。
由此可見，對于整體人群來說，與糖尿病常規支持和教育相比，減重干預并未降低超重或肥胖2型糖尿病患者的心血管疾病和死亡風險。但對于HbA1c<6.8%且SF-36得分<48的患者來說，減重干預反而使心血管疾病和死亡的風險增加；而對于HbA1c<6.8%且SF-36得分≥48的患者來說，減重干預能使心血管疾病和死亡的風險降低。然而，這一重要結果并沒有在Look AHEAD試驗的傳統亞組分析中發現[12]。通過因果森林，研究者不再需要手動調節交互項或進行多次亞組分析就可以找出這些具有處理效應異質性的人群，從而針對不同人群采取不同的干預措施。
除了Look AHEAD試驗以外，因果森林在其他臨床干預領域也有成功應用的研究案例。比如，在SPRINT試驗結果顯示，強化降壓能預防心血管事件，傳統亞組分析并未發現異質性亞組的存在[14]，但因果森林結果發現，對于基線時收縮壓>144 mmHg且吸煙的患者來說強化降壓反而使其心血管事件風險上升[15]。近年來，因果森林也開始應用于觀察性研究。比如，在MESA前瞻性隊列研究中，研究者在傾向評分匹配的基礎上構建因果森林，發現冠狀動脈鈣化與心血管疾病風險之間的關聯存在異質性，男性、西班牙裔、具有心血管危險因素的冠狀動脈鈣患者是心血管疾病的高危人群[16]。
4 代碼實現
本案例采用表1的隨機對照試驗模擬數據，利用R軟件進行因果森林建模分析，代碼及解釋如下。
表1
隨機對照試驗模擬數據庫的結構（只截取前六名患者數據）
隨機對照試驗模擬數據覆蓋3 582例腦卒中患者，結局是癡呆（dement），干預是他汀用藥（statin），協變量包括連續性指標的年齡（age）、體質指數（bmi），二分類指標的性別（gender）、慢性阻塞性肺病（copd）和冠心病（mi）。
setwd('路徑') # 設置工作路徑
# install.packages("R包的名字") # 首先安裝要用到的R包
#導入包
library(scorecard) # 用于訓練和驗證集分割
library(tidyverse) # 數據處理
library(grf) # 廣義隨機森林，需要grf0.10.2版本或以上
# 讀取數據，本案例所使用的數據文件為'statin.csv'，格式如表1所示
data <- read.csv("statin.csv")
data <- na.omit(data) # 刪除有缺失值的個體數據
data <- data.frame(lapply(data, function(x) as.numeric(as.character(x))))
# 將數據集按1∶1劃分訓練集和驗證集
data_split <- split_df(data, ratio = 0.5)
train_data <- data_split$train # 訓練集（50%，用于因果森林建模）
test_data <- data_split$test # 驗證集（50%，用于邏輯回歸驗證）
### 因果森林建模 ###
W <- train_data$statin # 干預W（如表1，本例中干預變量名稱為statin）
Y <- as.numeric(train_data$dement) # 結局Y（如表1，本例中結局變量名稱為dement）
covars <- c("age", "bmi", "copd", "gender", "mi") # 協變量
X <- select(train_data, covars) # 從訓練集中選擇協變量
# 生成因果森林
set.seed(42)
cf <- causal_forest(X, Y, W,
num.trees = 2000,
sample.fraction=0.5, # 構建每棵樹時從訓練集中抽樣的比例
min.node.size = 30,
honesty = TRUE,
honesty.fraction = 0.5, # 分裂子樣本的比例
tune.parameters = c("mtry", "alpha", "imbalance.penalty",
"honesty.prune.leaves"))
cf$tunable.params # 查看因果森林參數
# 查看變量重要性排序
var_imp <- c(variable_importance(cf))
names(var_imp) <- covars
sort(var_imp, decreasing=TRUE)
# 計算平均處理效應值
average_treatment_effect(cf, target.sample = "treated")
# ---------- 補充grf函數: find_best_tree ---------- #
# 由于版本更新，請從以下網址復制源碼至此處并運行
# https://gist.github.com/msegar/c15af2bfc067e5319eecfd67d785fde5
# ------------------------------------------------- #
# 在森林中找到最佳樹
best_trees <- find_best_tree(cf, "causal")
plot(tree <- get_tree(cf, best_trees$best_tree))
# 估計CATE的最佳線性擬合，可通過系數評估模型的性能
# 若mean.forest.prediction=1，則森林產生的平均預測是正確的；
# 若differential.forest.prediction = 1，則森林的異質性已被校準。
test_calibration(cf)
# 使用驗證集進行預測
X_test <- select(test_data, covars)
test.hat <- predict(cf, newdata = X_test, estimate.variance=TRUE)$predictions
### 邏輯回歸驗證 ###
# 根據最佳樹的分裂條件將驗證集分成4個亞組
group1 <- test_data[test_data$bmi <= 28.43 & test_data$bmi <= 24.79]
group2 <- test_data[test_data$bmi <= 28.43 & test_data$bmi > 24.79]
group3 <- test_data[test_data$bmi > 28.43 & test_data$age <= 61]
group4 <- test_data[test_data$bmi > 28.43 & test_data$age > 61]
grps <- list(group1, group2, group3, group4)
logistic_res <- matrix(1, 4, 3)
# 擬合邏輯回歸模型
for (i in 1:4){
glm_fit <- glm(dement ~ statin, data =grps[[i]], family ='binomial')
# print(glm_fit)
logistic_res[i, ] <- exp(cbind(OR = coef(glm_fit),confint(glm_fit)))[2, 1:3]
}
logistic_res # 查看擬合結果
通過運行上述代碼，可繪制出因果森林中最優的分組結果（圖2）。
圖2
模擬數據的因果森林最佳樹分組結果
圖中訓練樣本被切割為4個亞組，分別為
① BMI≤28.43且BMI≤24.79；
② BMI≤28.43且BMI>24.79；
③ BMI>28.43且年齡≤61歲；
④ BMI>28.43且年齡>61歲；
然后對四組進行Logistic回歸，在驗證集中計算出四組OR值（他汀組與非他汀組腦卒中后癡呆的比值比）分別為0.80［95%CI（0.48，1.31）］、0.66［95%CI（0.44，0.99）］、1.30［95%CI（0.83，2.05）］、0.33［95%CI（0.19，0.56）］。從結果可以看出，③組的處理效應相較其他亞組存在明顯異質性，即提示對于BMI>28.43且年齡≤61歲的患者來說，他汀治療可能無效，甚至存在潛在風險。
5 總結與展望
處理效應異質性是醫學實踐中關注的熱點問題，也是與精準醫學密切相關的重要問題。在異質性人群的發現和個體治療效果的預測方面，因果森林具有重要意義和應用價值[17]。
作為機器學習與因果推斷相結合的方法，因果森林可彌補傳統亞組分析和交互項檢驗的局限性。它通過數據驅動的方式識別各種特征組合的亞組，而不需要事先預設分組特征，在盡可能保證各類亞組人群都能被覆蓋的同時也不會增加假陽性的概率。因果森林不僅適用于隨機對照試驗，也適用于觀察性研究，但觀察性研究在構建和驗證因果森林時都需校正可能的混雜因素，比如采用傾向評分匹配，以滿足無混雜的假設[6, 11]。因果森林一般作為探索性分析，其結果需要后續研究的進一步驗證。
除了本文所介紹的因果森林以外，還有其他更復雜的機器學習算法也能用于處理效應異質性的評價，比如，貝葉斯加性回歸樹、神經網絡等[18]。這些新型方法比傳統方法的實施難度大，而且存在報告不規范的問題，在使用這些機器學習方法開展相關研究時，需要多學科交叉團隊的緊密合作。未來隨著統計方法學的日趨完善，機器學習與因果推斷相結合的方法會越來越多地應用于醫學處理效應異質性的評價[19]，推動精準醫學的發展。
