陳文松 劉 曼 劉玉秀，4△ 許敏怡 熊 殷 鞏浩雯 李維勤

【提 要】 目的 探討觀察性研究中用于混雜偏倚控制的傾向性評分匹配、馬氏距離匹配和遺傳匹配三種方法的性能。方法 針對連續(xù)型結(jié)局變量，設(shè)定混雜變量與處理分組變量之間具有不同復(fù)雜度的回歸模型結(jié)構(gòu)，采用Monte-Carlo模擬方法比較三種匹配方法在處理組間效應(yīng)估計(jì)和匹配前后自變量均衡的區(qū)別，進(jìn)而對三種方法性能進(jìn)行評估。結(jié)果 在給定的模擬情形下，相比于傾向性評分匹配和馬氏距離匹配，遺傳匹配法得出的效應(yīng)估計(jì)偏差最小，匹配后兩組自變量均衡性最好。結(jié)論 遺傳匹配在三種匹配方法中表現(xiàn)出較好的統(tǒng)計(jì)性能，可考慮作為觀察性研究中控制混雜偏倚優(yōu)先推薦的匹配方法。

在非隨機(jī)觀察性研究中，當(dāng)需要考慮某目標(biāo)處理變量的效應(yīng)時，由于混雜偏倚的客觀存在，常常導(dǎo)致估計(jì)出的組間差別不能反映真實(shí)情況，此時混雜偏倚的控制就尤為關(guān)鍵。除已知的分層分析、多因素分析等傳統(tǒng)方法，近年來，利用距離度量(metric)對研究個體進(jìn)行匹配的傾向性評分匹配法(propensity score matching，PSM)以及馬氏距離匹配法(Mahalanobis distance matching，MDM)、遺傳匹配法(genetic matching，GM)在醫(yī)學(xué)研究中的應(yīng)用越來越廣泛[1-3]，然而有關(guān)這幾種匹配方法在統(tǒng)計(jì)推斷時的性能比較研究尚未見詳細(xì)文獻(xiàn)報道。

本文將在介紹基于進(jìn)化算法的遺傳匹配方法的基礎(chǔ)上，采用Monte-Carlo模擬方法對三種匹配方法的統(tǒng)計(jì)性能進(jìn)行比較研究，并通過一個實(shí)例進(jìn)行分析說明。

匹配方法

理論上，將同質(zhì)性相近的個體在不同的組間進(jìn)行匹配，組間的自變量的分布將趨于均衡，從而減少或抵消兩組混雜因素對治療效果的影響，這就是匹配的基本原理[4-5]。在醫(yī)學(xué)研究中，除規(guī)定匹配個體所有變量取值均相同的精確匹配外，以馬氏距離為基礎(chǔ)的多元變量匹配和將多維數(shù)據(jù)轉(zhuǎn)化為一維的傾向性評分匹配最為常見。有關(guān)傾向評分匹配已有較多文獻(xiàn)介紹，此不贅述。以下重點(diǎn)介紹基于馬氏距離發(fā)展的廣義馬氏距離，以及以廣義馬氏距離為度量的遺傳匹配。

1936年，P.C.Mahalanobis提出了馬氏距離(Mahalanobis distance，MD)，用以表示數(shù)據(jù)的協(xié)方差距離[6]：

(1)

其中Xi和Xj分別為第i個個體和第j個個體的自變量組成的向量，S-1為樣本協(xié)方差矩陣的逆。個體間馬氏距離越小，表示個體間越相似。

2013年，Diamond及其同事提出了基于進(jìn)化算法的遺傳匹配(genetic matching)來控制組間混雜因素的影響。遺傳匹配在馬氏距離的基礎(chǔ)上創(chuàng)立了更一般化的距離度量，即廣義馬氏距離(generalized Mahalanobis distance，GMD)[7]：

(2)

個體的自變量集合X除自變量以外，還可以包含傾向性評分。此時方程(1)中的X可以用Z代替，其中Z是一個由傾向性評分π(X)和自變量X組成的矩陣。如果將PS以外的其他變量的權(quán)重設(shè)為0，此時GMD等同于傾向性評分；如果將PS的權(quán)重設(shè)為0，其他自變量的各個權(quán)重設(shè)為1，此時GMD等價于馬氏距離。PSM和MDM都是遺傳匹配的極限情況。

常見的計(jì)算傾向性評分的方法有l(wèi)ogistic回歸、Probit回歸，以及機(jī)器學(xué)習(xí)估計(jì)中的神經(jīng)網(wǎng)絡(luò)、支持向量機(jī)、分類與回歸樹、隨機(jī)森林等方法[8-9]。其中利用logistic回歸計(jì)算傾向性評分最為常見，即利用組別為因變量，建立logistic模型，本文就采用這種方法計(jì)算傾向性評分。

遺傳匹配根據(jù)每個變量的相對重要性對其加權(quán)來獲取最佳的自變量平衡。遺傳匹配利用遺傳算法(genetic algorithms)搜索一系列距離度量，通過最小化損失函數(shù)，找到使得自變量均衡性最高的加權(quán)矩陣W。該算法允許使用者自定義損失函數(shù)，默認(rèn)損失函數(shù)定義為匹配后兩組自變量的Kolmogorov-Smirnov檢驗(yàn)的P值及配對t檢驗(yàn)的P值最小值，用以描述兩組自變量的最大差異。

遺傳算法是一種通過模擬達(dá)爾文生物進(jìn)化論的自然選擇和遺傳學(xué)機(jī)理的生物進(jìn)化過程搜索最優(yōu)解的方法[10-11]。圖1總結(jié)了遺傳算法的流程。該算法產(chǎn)生一批初代權(quán)重Ws，每個W對應(yīng)一個不同的距離度量。在每一代中，計(jì)算每個匹配樣本對應(yīng)的損失函數(shù)。選擇使得損失函數(shù)最小的權(quán)重作為父代，產(chǎn)生子代權(quán)重并不斷迭代這一流程，直至漸進(jìn)收斂于最優(yōu)解。每一代產(chǎn)生的權(quán)重?cái)?shù)量稱為種群規(guī)模，在迭代過程中種群規(guī)模保持不變。增大種群規(guī)模通常會改善遺傳匹配的效果，但會增加算法尋找最優(yōu)解所花費(fèi)的時間。

圖1 遺傳算法流程圖

當(dāng)損失函數(shù)滿足迭代標(biāo)準(zhǔn)時停止迭代，此時的權(quán)重即為最優(yōu)權(quán)重。利用最優(yōu)權(quán)重計(jì)算廣義馬氏距離，根據(jù)此距離度量進(jìn)行匹配，得到匹配后樣本。與傾向性評分一樣，此時可以使用不同的匹配策略，例如有無替換、有無卡鉗的最近鄰匹配或最優(yōu)匹配等。

值得注意的是，遺傳匹配與其他匹配方法的使用前提一樣，即要事先確定匹配的具體變量、選定測量自變量平衡性的指標(biāo)，以及指定合適的匹配方式。

模擬實(shí)驗(yàn)

本模擬實(shí)驗(yàn)基于模擬生成的數(shù)據(jù)集，采用Monte-Carlo方法比較遺傳匹配與經(jīng)典的馬氏距離匹配、傾向評分匹配在傾向性評分模型不同復(fù)雜程度下的效應(yīng)估計(jì)精度和匹配前后自變量均衡程度，評價其統(tǒng)計(jì)性能，為方法學(xué)選擇提供理論依據(jù)。

1.模擬數(shù)據(jù)

Setoguchi及其同事在探究機(jī)器學(xué)習(xí)技術(shù)構(gòu)建傾向性評分模型的研究中提出了一種模擬數(shù)據(jù)生成策略[12]，該策略產(chǎn)生的模擬數(shù)據(jù)包含了實(shí)際中可能出現(xiàn)的多數(shù)情景，具有較好的代表性。本研究在此策略的基礎(chǔ)上進(jìn)行了略微的修改。利用SAS軟件總共生成1000個模擬數(shù)據(jù)集，其中每個數(shù)據(jù)集的樣本量為1000。在每個數(shù)據(jù)集中，模擬生成10個自變量，其中4個混雜變量，3個暴露預(yù)測變量和3個結(jié)果預(yù)測變量。4個自變量(X2，X4，X7，X10)服從標(biāo)準(zhǔn)正態(tài)分布，其余6個二元變量(X1，X3，X5，X6，X8，X9)則由標(biāo)準(zhǔn)正態(tài)分布轉(zhuǎn)換得到，具體的模擬數(shù)據(jù)變量結(jié)構(gòu)及自變量間相關(guān)系數(shù)見圖2。設(shè)定group為組別變量，并使其在自變量均值上約等于0.5以保證兩組人數(shù)相同。Y作為連續(xù)型結(jié)局變量，由group和其他自變量的線性組合生成。在文后附錄中提供了計(jì)算傾向性評分、結(jié)局生成及治療分配的詳細(xì)過程。

圖2 模擬數(shù)據(jù)各變量關(guān)系的基本結(jié)構(gòu)

基于所設(shè)定的模擬數(shù)據(jù)結(jié)構(gòu)，本研究指定了7種復(fù)雜程度各不相同的情形來指代真正的傾向性模型，以便評估三種匹配方法在傾向性模型不同復(fù)雜程度下的表現(xiàn)情況：

A：線性可加(僅主效應(yīng))

B：輕微非線性(包含1個二次項(xiàng))

C：中等非線性(包含3個二次項(xiàng))

D：輕微非可加性(包含3個兩因素交互作用項(xiàng))

E：輕微非可加性和非線性(包含3個兩因素交互作用項(xiàng)和1個二次項(xiàng))

F：中等非可加性(包含10個兩因素交互作用項(xiàng))

G：中等非可加性和非線性(包含10個兩因素交互作用項(xiàng)和3個二次項(xiàng))

為評價模擬生成的數(shù)據(jù)集的統(tǒng)計(jì)特征，本研究比較了兩組間自變量的均衡性。圖3展示了傾向性評分模型為情形A時兩組間傾向性評分分布與自變量均衡情況，其他情形下的模擬數(shù)據(jù)表現(xiàn)與其類似，這里不予羅列。從圖3可知，試驗(yàn)組中傾向性評分大于0.5，對照組中傾向性評分小于0.5，兩組傾向性評分分布互有重疊；兩組間混雜變量與暴露預(yù)測變量的標(biāo)準(zhǔn)化差值大于0.1，混雜變量與暴露預(yù)測變量在兩組間分布不均衡。

圖3 情形A時模擬數(shù)據(jù)集兩組自變量均衡性

2.匹配方案

本研究使用MDM、PSM和GM來估計(jì)每個情形中的group的效應(yīng)。匹配方案均選擇有放回的最近鄰匹配(nearest neighbor matching)，匹配比例為1：1，同時PSM的卡鉗值設(shè)為0.02，GM中每代的種群規(guī)模(pop.size)設(shè)置為500。由于實(shí)際數(shù)據(jù)的復(fù)雜性，本研究很難確定哪些變量應(yīng)納入模型以及模型中變量間的具體關(guān)系，所以在PSM和GM計(jì)算傾向性評分時，采用兩種策略建立傾向性評分模型：在第一種策略中，默認(rèn)所有自變量均以線性可加的形式包含在傾向性評分模型中；而在第二種策略中，仍舊將所有自變量納入傾向性評分模型，但對于連續(xù)型變量，使用限制性立方樣條(restricted cubic spline，RCS)來考慮連續(xù)型變量的非線性關(guān)系。MDM使用所有自變量計(jì)算馬氏距離。三種匹配方法均應(yīng)用R軟件中的Matching包實(shí)現(xiàn)。此外，本研究還在遺傳匹配中去除傾向性評分，只使用所有自變量，探究其對結(jié)果的影響。其中，使用“PSM”表示使用線性可加模型進(jìn)行的傾向性評分匹配，“PSM+”表示使用非線性可加模型進(jìn)行的傾向性評分匹配，“MDM”表示馬氏距離匹配，“GM”表示使用線性可加模型得到的ps進(jìn)行的遺傳匹配，“GM-”表示不使用ps進(jìn)行的遺傳匹配，“GM+”表示使用非線性可加模型得到的ps進(jìn)行的遺傳匹配。

3.評價指標(biāo)

考慮采用2個度量值作為效應(yīng)估計(jì)性能的評價指標(biāo)，分別為：①相對偏倚(average percent relative bias，APRB)：估計(jì)值與真實(shí)值的相對差異百分比；其中采用配對t檢驗(yàn)對匹配后樣本計(jì)算估計(jì)值。②均方根誤差(root mean squared error，RMSE)：估計(jì)值與真實(shí)值偏差的平方與估計(jì)次數(shù)比值的平方根。另外，采用2個指標(biāo)來評價匹配前后自變量均衡性，分別為：①標(biāo)準(zhǔn)化差值(standardized mean difference，SMD)：用于定量評價均衡性。一般認(rèn)為標(biāo)準(zhǔn)化差值小于0.1時具有良好的均衡性。②配對t檢驗(yàn)或配對卡方檢驗(yàn)的P值。對于連續(xù)性變量使用配對t檢驗(yàn)，對于分類變量使用配對卡方檢驗(yàn)。考慮到篇幅所限，本文中自變量均衡性結(jié)果僅提供標(biāo)準(zhǔn)化差值。

4.模擬結(jié)果

在不同的傾向性評分模型下，使用不同方法對模擬數(shù)據(jù)集分析，比較不同方法的效應(yīng)估計(jì)性能評價指標(biāo)和匹配前后自變量均衡性。

(1)不同匹配方法后的處理效應(yīng)估計(jì)的比較

從表1可見，各方法在各情形下的總體偏倚從小到大依次為：GM≈GM+

表1 不同匹配方法效應(yīng)估計(jì)相對偏倚的比較

(2)不同匹配方法后組間自變量的均衡性比較

考慮結(jié)果的呈現(xiàn)性，本研究在此挑選了7種情形中的具有代表性的4種情形，并選擇了各類匹配方法中表現(xiàn)最優(yōu)的3種方法。從圖4可見，經(jīng)遺傳匹配處理后，兩組自變量的標(biāo)準(zhǔn)化差值最小，在各種情形中的自變量均衡情況表現(xiàn)均最優(yōu)。在指定了正確的傾向性評分模型后(即情形A)，PSM+與GM在整體上表現(xiàn)相近，但其SMD的分布略高于GM。隨著真實(shí)傾向性評分模型復(fù)雜度的增加，PSM+的兩組自變量均衡程度開始下降。MDM在控制分類變量的均衡性上表現(xiàn)較優(yōu)，但在連續(xù)型變量上表現(xiàn)最差，其中在匹配前最不均衡的變量X7上控制表現(xiàn)最差。

圖4 三種匹配方法匹配各自變量均衡性情況：SMD

(3)傾向性評分匹配中不同建模策略對統(tǒng)計(jì)性能的影響

圖5給出了在不同的真實(shí)傾向性評分模型下，使用兩種不同建模策略進(jìn)行傾向性評分匹配的自變量均衡情況。從圖5可見，使用非線性可加模型進(jìn)行傾向性評分匹配的組間自變量均衡程度更優(yōu)，在任何情形下各變量SMD的上四分位數(shù)均小于0.1。當(dāng)真實(shí)傾向性模型僅為簡單的線性可加模型時(情形A)，PSM與PSM+的表現(xiàn)一致；當(dāng)真實(shí)傾向性模型包含非線性項(xiàng)時(情形C、情形E、情形G)，在二元變量上，PSM與PSM+的表現(xiàn)一致。但在連續(xù)型變量上，PSM+的表現(xiàn)優(yōu)于PSM。從表1可見，在任何情形下，PSM+的療效相對偏倚小于PSM。

(4)遺傳匹配中使用傾向性評分與否對統(tǒng)計(jì)性能的影響

圖6給出了在不同的真實(shí)傾向性評分模型下，使用簡單線性相加模型得到的PS、使用非線性相加模型得到的PS與不使用PS進(jìn)行遺傳匹配的自變量均衡情況，可見，不使用PS作為自變量進(jìn)行遺傳匹配的組間自變量均衡程度更差。當(dāng)計(jì)算PS的模型與真實(shí)傾向性評分模型更相近時，使用PS的遺傳匹配效果更優(yōu)。GM與GM+的表現(xiàn)大體相似，在復(fù)雜情形下GM+的表現(xiàn)略優(yōu)于GM。從表1可見，使用PS進(jìn)行遺傳匹配的偏倚更小，從小到大依次為GM≈GM+

圖5 不同傾向性評分建模策略匹配各自變量均衡情況：SMD

圖6 不同策略遺傳匹配后各自變量均衡情況：SMD

案例分析

Dehejia等人創(chuàng)建的數(shù)據(jù)集“l(fā)alonde”基于一個全國性的職業(yè)培訓(xùn)研究，探討接受職業(yè)培訓(xùn)是否會增加學(xué)員收入，該數(shù)據(jù)集內(nèi)置于“Match”包中。該研究的研究對象共445例(試驗(yàn)組185例，對照組260例)，基線變量包含年齡、受教育時間、人種等8個變量。結(jié)局變量為研究對象在1978年的收入。本研究使用MDM、PSM與GM對此進(jìn)行分析，匹配方案與先前模擬試驗(yàn)中一致。匹配前后自變量均衡情況見表2及圖7。匹配前，有4個自變量在兩組間的SMD大于0.1(有無畢業(yè)文憑、族裔、受教育時間、年齡)，其中“有無畢業(yè)文憑”在兩組間的P值小于0.001。使用三種方法匹配后自變量均衡情況均有所改善，但GM匹配后兩組自變量均衡性最高，匹配后兩組SMD最小。

表2 匹配前后自變量特征

圖7 匹配前后兩組標(biāo)準(zhǔn)化差異

討 論

本文在簡要介紹利用進(jìn)化算法自動迭代地構(gòu)建匹配度量的遺傳匹配的基礎(chǔ)上，借助Monte-Carlo模擬方法，通過效應(yīng)估計(jì)的準(zhǔn)確性及匹配后組間自變量的均衡性兩個方面評價了馬氏距離匹配、傾向性評分匹配和遺傳匹配在數(shù)據(jù)不同復(fù)雜度下的統(tǒng)計(jì)性能。模擬結(jié)果顯示，當(dāng)指定的傾向性評分模型與真實(shí)的傾向性評分模型相近時，傾向性評分匹配表現(xiàn)出了較好的統(tǒng)計(jì)學(xué)性能。但當(dāng)模型復(fù)雜度上升后，即指定的傾向性評分模型與真實(shí)模型相差較大時，傾向性評分匹配的偏倚開始增大。而遺傳匹配不會受真實(shí)模型的影響，在不同模型復(fù)雜度下得出的效應(yīng)估計(jì)偏差最小，匹配后兩組自變量均衡性最高。實(shí)例分析結(jié)果也驗(yàn)證了，遺傳匹配在保持較多匹配例數(shù)下仍能達(dá)到最好的組間均衡。

近年來，傾向性評分匹配在醫(yī)學(xué)研究中得到了廣泛的應(yīng)用，但在實(shí)際使用中仍存在著如何確定合理傾向性評分模型的問題。Granger等人統(tǒng)計(jì)了在2014年到2016年發(fā)表在高影響因子的醫(yī)學(xué)期刊上使用傾向性評分進(jìn)行分析的研究，其中20.9%(187/894)的研究并未闡述如何評估所建立的傾向性評分模型的合理性[13]。通常情況下，當(dāng)研究人員根據(jù)數(shù)據(jù)特征并結(jié)合自身專業(yè)知識盡可能地納入混雜因素，確定傾向性評分模型后，可以通過評價匹配后自變量平衡情況來判斷是否指定正確并進(jìn)行調(diào)整。Rosenbaum和Rubin曾建議應(yīng)根據(jù)匹配后自變量均衡情況迭代修改傾向性評分模型已達(dá)到最優(yōu)均衡，但在實(shí)際情況中許多研究者缺乏這種建立合理的傾向性評分模型的意識或技能，這無疑會導(dǎo)致傾向性評分匹配在醫(yī)學(xué)研究中的錯用，得到并不真實(shí)的研究結(jié)論。

遺傳匹配將損失函數(shù)定義為匹配后組間均衡性指標(biāo)，利用遺傳算法自動地找到最優(yōu)的匹配組合，從而實(shí)現(xiàn)自變量的均衡。與傾向性評分匹配相比，遺傳匹配不需要研究者嚴(yán)格指定傾向性評分模型，不需要研究者根據(jù)匹配后自變量平衡情況手動調(diào)整傾向性評分模型，這極大地降低了研究者的工作難度。目前尚無評價匹配后自變量均衡性的最優(yōu)指標(biāo)，除最常見的組間檢驗(yàn)P值及標(biāo)準(zhǔn)化差值外，還有方差比、Hosmer-Lemeshow檢驗(yàn)統(tǒng)計(jì)量、KS檢驗(yàn)統(tǒng)計(jì)量及C統(tǒng)計(jì)量等[13]。在本文使用的R軟件包“Matching”中則可自定義損失函數(shù)，使用任何均衡性指標(biāo)進(jìn)行遺傳匹配。遺傳匹配獲得最優(yōu)的匹配度量后，匹配的方案與其他類似。傳統(tǒng)的最近鄰匹配、卡鉗匹配、1：n匹配等在遺傳匹配中仍可使用。

當(dāng)然，遺傳匹配與其他匹配方法一樣，只能處理已知的混雜因素，對未知混雜無法控制，也不適用于時依自變量所致混雜因素的處理[8-9]；當(dāng)其使用傾向性評分作為自變量進(jìn)行匹配時，也面臨傾向性評分模型自變量的選取問題。

本研究評價了使用不同傾向性評分建模策略對傾向性評分匹配和遺傳匹配的影響程度，研究結(jié)果表明，在建立傾向性評分模型時，考慮自變量的非線性關(guān)系會提高匹配的效果。另外，與不納入PS的遺傳匹配相比，納入PS進(jìn)行遺傳匹配具有更好的統(tǒng)計(jì)性能，且當(dāng)使用的傾向性評分模型越正確，遺傳匹配的效果越好，這恰好支持了Diamond等人的建議。這一結(jié)論也符合雙重穩(wěn)健估計(jì)的理論[14]，即納入PS的遺傳匹配提升了效應(yīng)估計(jì)的穩(wěn)健程度。

在真實(shí)世界研究中，為去除混雜因素的影響，得到無偏的效應(yīng)估計(jì)，眾多學(xué)者提出了許多方法，例如傳統(tǒng)的分層分析、多因素回歸分析、傾向性得分匹配、逆概率加權(quán)匹配、重疊加權(quán)等方法。近年來一些更為復(fù)雜的方法也被提出，例如G估計(jì)、雙重穩(wěn)健估計(jì)等方法。本文因篇幅所限，僅比較了三種匹配方法的統(tǒng)計(jì)表現(xiàn)，遺傳匹配與其他方法的比較有待進(jìn)一步研究。

本研究只比較了在設(shè)定的模擬情形下三種匹配方法的結(jié)果，此假設(shè)下模擬數(shù)據(jù)盡管能滿足組間傾向性評分的重疊假定(overlap assumption)，但并未對方法的匹配度進(jìn)行探討比較，也未探討兩組傾向性評分的重疊程度對匹配效能的影響；另外，本研究的模擬試驗(yàn)只考慮了1∶1有放回的匹配策略，未進(jìn)行其他比例的匹配，這些問題均有待進(jìn)一步研究。