合理進(jìn)行多重Logistic回歸分析
——結(jié)合平均處理效應(yīng)分析

胡純嚴(yán) ，胡良平 ，2*

（1.軍事科學(xué)院研究生院，北京 100850；2.世界中醫(yī)藥學(xué)會(huì)聯(lián)合會(huì)臨床科研統(tǒng)計(jì)學(xué)專業(yè)委員會(huì)，北京 100029*通信作者：胡良平，E-mail：lphu927@163.com）

在進(jìn)行多重Logistic回歸分析時(shí)，若自變量中包含“處理變量”，研究者不僅要關(guān)注處理變量對(duì)因變量的影響是否有統(tǒng)計(jì)學(xué)意義，還希望了解處理變量取“處理水平”與“對(duì)照水平”所產(chǎn)生的因果效應(yīng)，這兩個(gè)水平下的因果效應(yīng)分別被稱為潛在結(jié)果Y（1）與Y（0）。進(jìn)一步，還需求出兩種潛在結(jié)果期望值的差值，該差值被稱為平均處理效應(yīng)（ATE）。本文將介紹如何結(jié)合平均處理效應(yīng)分析，合理進(jìn)行多重Logistic回歸分析的基本思想和具體方法。

1 基本概念

1.1 潛在結(jié)果與潛在結(jié)果均值

潛在結(jié)果（Potential Outcomes，PO）被定義為與受試者相關(guān)的獨(dú)特和內(nèi)在價(jià)值，它描述了一個(gè)理想化的數(shù)據(jù)來源，研究者可以觀察受試者對(duì)所有可能處理分配的響應(yīng)［1］。例如，如果數(shù)據(jù)分析者正在估計(jì)一個(gè)二值處理變量T的因果效應(yīng)（設(shè)結(jié)果變量為Y），其中，T=0表示受試者接受的處理為對(duì)照水平，T=1表示受試者接受的處理為真實(shí)的處理水平，那么每個(gè)受試者都有兩個(gè)潛在結(jié)果Y（0）與Y（1）。受試者水平（也稱為單位水平）上的處理效應(yīng)定義為潛在結(jié)果Y（1）-Y（0）的差值。

潛在結(jié)果Y（1）與Y（0）各自的期望值被稱為潛在結(jié)果均值（Potential Outcome Mean，POM）［2-3］，POM的計(jì)算見式（1）。

式（1）中，“E［·］”代表求“·”的期望值（即總體均值）。

也可以說，潛在結(jié)果Y（1）與Y（0）分別是T=1與T=0條件下，Y的樣本均值；而潛在結(jié)果均值μ1=E[Y(1)]與μ0=E[Y(0)]分別是T=1與T=0條件下，Y的總體均值。

1.2 平均處理效應(yīng)

整個(gè)總體的平均處理效應(yīng)有時(shí)也稱為平均因果效應(yīng)（ACE）［1］，其計(jì)算公式見式（2）。

由式（2）可知，平均處理效應(yīng)在本質(zhì)上反映了兩種處理在總體上的因果效應(yīng)之差。此差值的絕對(duì)值越大，表明處理變量的一個(gè)水平相對(duì)于另一個(gè)水平的作用越大。

1.3 處理組中的平均處理效應(yīng)

處理組中的平均處理效應(yīng)（Average Treatment Effect for the Treated，ATET或ATT）是處理組中的個(gè)體之間的平均因果效應(yīng)［1］，其計(jì)算公式見式（3）。

式（3）中 ，μ1|T=1=E[Y(1)|T=1]與 μ0|T=1=E[Y(0)|T=1]是處理組中分別接受真實(shí)處理與未接受真實(shí)處理的潛在結(jié)果。

1.4 逆概率權(quán)重

在估計(jì)上述提及的統(tǒng)計(jì)量（POM、ATE和ATT）時(shí)，需要對(duì)每個(gè)個(gè)體賦予不同的權(quán)重。有學(xué)者提出“逆概率加權(quán)”的方法［4］，即以受試者進(jìn)入處理組的概率的倒數(shù)作為權(quán)重，代入估計(jì)統(tǒng)計(jì)量的公式中去計(jì)算。也就是說，構(gòu)建二值處理變量的處理水平關(guān)于協(xié)變量的多重Logistic回歸模型，求出各個(gè)體在二值處理變量上的預(yù)測(cè)值（即概率，也被稱為傾向性評(píng)分）。當(dāng)個(gè)體屬于處理組時(shí)，求出的概率P保持不變；當(dāng)個(gè)體屬于對(duì)照組時(shí)，求出的概率P需進(jìn)行簡(jiǎn)單變換，即修改為Q=1-P。與這些個(gè)體對(duì)應(yīng)的權(quán)重分別為1/P與1/Q。

2 估算方法

2.1 估算方法簡(jiǎn)介

SAS/STAT的proc causaltrt過程中介紹了6種與處理變量的因果效應(yīng)統(tǒng)計(jì)量估計(jì)有關(guān)的方法［1，4］。通過擬合處理變量T或結(jié)果變量Y或兩者的模型來調(diào)整混雜變量的影響。在model語句中指定結(jié)果變量Y，在psmodel語句中指定處理變量T，通過在proc causaltrt語句中使用method=選項(xiàng)指定估算方法（若未使用method=選項(xiàng)，系統(tǒng)將使用默認(rèn)估算方法）。

6種估算方法包括：①逆概率加權(quán)（IPW）估算法；②具有比率調(diào)整的逆概率加權(quán)（IPWR）估算法；③具有比率和尺度調(diào)整的逆概率加權(quán)（IPWS）估算法；④增強(qiáng)的逆概率加權(quán)（AIPW）估算法；⑤逆概率加權(quán)回歸調(diào)整（IPWREG）估算法；⑥回歸調(diào)整（REGADJ）估算法。在使用前5種估算法時(shí)，使用者必須在psmodel語句內(nèi)為處理變量指定模型；而在使用第四種和第五種估算法時(shí)，使用者還必須在model語句內(nèi)為結(jié)果變量指定模型；在使用第六種估算法時(shí)，使用者只需要在model語句內(nèi)為結(jié)果變量指定模型。

2.2 估算方法的要求

使用者可在proc causaltrt語句的method=選項(xiàng)中指定估算方法，并且，還可在model語句和psmodel語句中提供特定的附加建模信息，如表1所示。

表1 估算方法的要求Table 1 Requirements for estimation methods

2.3 具體的估算方法

因估算方法比較復(fù)雜，此處僅給出上述提及的6種估算方法中的前3種估算方法的計(jì)算公式，其他估算方法的計(jì)算公式見文獻(xiàn)［1，4］。

Proc causaltrt過程使用psmodel語句中指定的效應(yīng)將Logistic回歸擬合到傾向性評(píng)分模型。然后通過取估計(jì)的傾向性評(píng)分的倒數(shù)來計(jì)算逆概率權(quán)重。假設(shè)使用者有觀測(cè)值（yi，ti，xi），其中，i=1，2，…，n。傾向性評(píng)分模型的參數(shù)估計(jì)值由給出，傾向性評(píng)分的預(yù)測(cè)值計(jì)算公式見式（4）。

Proc causaltrt過程實(shí)施的三種加權(quán)方法通過求解式（5）獲得潛在結(jié)果均值的無偏估計(jì)。

對(duì) 于 μ =(μ0，μ1)，其 中 ，Sipw，i的 計(jì) 算 公 式見式（6）。

這些估算方法在（η0，η1）的值上有所不同。每種方法的（η0，η1）選擇和相應(yīng)的潛在結(jié)果平均值估計(jì)的計(jì)算公式非常復(fù)雜，因篇幅所限，此處從略。

3 實(shí)例與SAS實(shí)現(xiàn)

3.1 問題與數(shù)據(jù)結(jié)構(gòu)

【例1】假設(shè)在一項(xiàng)臨床試驗(yàn)研究中，486名受試者有患2型糖尿病的風(fēng)險(xiǎn)，他們可以選擇自己想要接受的兩種預(yù)防藥物中的某一種，研究者不對(duì)受試者的任何協(xié)變量進(jìn)行控制。在這項(xiàng)研究中，關(guān)注的結(jié)果是受試者是否在5年內(nèi)發(fā)展為2型糖尿?。?］。假設(shè)數(shù)據(jù)集Drugs包含以下變量，Age：試驗(yàn)開始時(shí)受試者的年齡（歲）；BMI：試驗(yàn)開始時(shí)受試者的體重指數(shù)；Diabetes2：受試者是否發(fā)展為2型糖尿病的指標(biāo)，其取值為是和否；Drug：處理變量，其取值為Drug_A（對(duì)照藥）和Drug_X（試驗(yàn)藥）；Gender：性別，取值為男性和女性。試分析處理變量Drug對(duì)于結(jié)果變量Diabetes2的平均處理效應(yīng)（ATE）。Drugs數(shù)據(jù)集中的前10個(gè)觀測(cè)結(jié)果見表2。

表2 某種藥物預(yù)防2型糖尿病的臨床試驗(yàn)數(shù)據(jù)結(jié)構(gòu)Table 2 Clinical trial data structure of a drug to prevent type 2 diabetes

3.2 采用傳統(tǒng)的多重Logistic回歸分析

比較兩種藥物（Drug_A與Drug_X）中的3個(gè)協(xié)變量（Age、Gender、BMI）是否具有可比性，發(fā)現(xiàn)年齡（Age）在兩組之間差異有統(tǒng)計(jì)學(xué)意義（t=5.420，P＜0.01）。

變量Diabetes2為二值結(jié)果變量，可采用二值結(jié)果變量的多重Logistic回歸分析全部協(xié)變量（包括處理變量）對(duì)二值結(jié)果變量的影響情況。設(shè)所需要的SAS程序如下：

【SAS主要輸出結(jié)果及解釋】上面這段SAS程序的主要輸出結(jié)果見表3。

表3 多重Logistic回歸模型中參數(shù)的最大似然估計(jì)結(jié)果Table 3 Maximum likelihood estimation results of parameters in the multiple Logistic regression model

以上結(jié)果表明：截距項(xiàng)和3個(gè)自變量對(duì)結(jié)果變量的影響均有統(tǒng)計(jì)學(xué)意義，可以得到多重Logistic回歸模型如下：

在式（7）顯示的多重Logistic回歸模型中，研究者最關(guān)注的是處理變量Drug對(duì)結(jié)果變量的影響，其回歸系數(shù)為0.383，它是“Drug_A”相對(duì)于“Drug_X”而計(jì)算出來的回歸系數(shù)。該回歸系數(shù)大于0，表明與使用Drug_X治療相比，接受Drug_A治療的受試者更易于患2型糖尿病。各變量?jī)?yōu)勢(shì)比的估計(jì)結(jié)果見表4。

表4 各變量?jī)?yōu)勢(shì)比的估計(jì)結(jié)果Table 4 Estimated results of odds ratio of each variable

由以上結(jié)果可知，受試者接受Drug_A治療比接受Drug_X治療更易患2型糖尿病的優(yōu)勢(shì)比為OR=2.152，95%置信區(qū)間為（1.428～3.244）。

顯然，傳統(tǒng)的Logistic回歸分析的結(jié)果中，并沒有包含變量Drug對(duì)于結(jié)果變量Diabetes2的平均處理效應(yīng)（ATE）的估計(jì)結(jié)果。

3.3 采用causaltrt過程估計(jì)平均處理效應(yīng)

3.3.1 采用逆概率加權(quán)法（IPW）進(jìn)行參數(shù)估計(jì)

設(shè)所需要的SAS程序如下：

【SAS主要輸出結(jié)果及解釋】以上SAS程序的輸出結(jié)果見表5。

表5是基于傾向性評(píng)分分析后，再基于匹配后的數(shù)據(jù)集并采用IPW法對(duì)“結(jié)果變量取值為患2型糖尿病”構(gòu)建Logistic回歸模型，得到結(jié)果：Drug_X與Drug_A的潛在結(jié)果均值（POM）分別為0.459與0.649；兩種藥物平均因果效應(yīng)ATE=-0.190，表明Drug_X在預(yù)防2型糖尿病方面比Drug_A更有效。正如95%置信區(qū)間（-0.290～-0.092）或P值（＜0.01）所示，在0.05水平上，ATE與0之間的差異有統(tǒng)計(jì)學(xué)意義。

表5 因果效應(yīng)的分析結(jié)果Table 5 Analysis results of causal effects

“ATE=-0.190”的真實(shí)含義可用通俗的語言解釋如下：因?yàn)閅=1代表受試者患2型糖尿病，Y=0代表受試者不患2型糖尿病，故在理論上可認(rèn)為，某處理組中Y的平均值是一個(gè)“隱變量”，它的取值越接近1，表明該處理組中的受試者越易患2型糖尿?。环粗嗳?。就本例而言，服用Drug_X的受試者對(duì)應(yīng)的Y的總體平均值為0.459，而服用Drug_A的受試者對(duì)應(yīng)的Y的總體平均值為0.649。說明服用Drug_X預(yù)防2型糖尿病比Drug_A更有效，使Y的總體平均值下降了0.190。

3.3.2 采用帶比率調(diào)整的逆概率加權(quán)法（IPWR）進(jìn)行參數(shù)估計(jì)

設(shè)所需要的SAS程序如下：

將“第3.3.1節(jié)”SAS程序的第1句中的方法選項(xiàng)修改為“method=IPWR”，其他內(nèi)容不變。

【SAS主要輸出結(jié)果及解釋】SAS程序的輸出結(jié)果見表6。

表6是基于傾向性評(píng)分分析后，再基于匹配后的數(shù)據(jù)集并采用IPWR法對(duì)“結(jié)果變量取值為患2型糖尿病”構(gòu)建Logistic回歸模型，得到結(jié)果：Drug_X與Drug_A的潛在結(jié)果均值（POM）分別為0.468與0.647；兩種藥物平均因果效應(yīng)ATE=-0.179，表明Drug_X在預(yù)防2型糖尿病方面比Drug_A更有效。正如95%置信區(qū)間（-0.277～-0.080）或P值（＜0.01）所示，在0.05水平上，ATE與0之間的差異有統(tǒng)計(jì)學(xué)意義。

表6 因果效應(yīng)的分析結(jié)果Table 6 Analysis results of causal effects

3.3.3 采用帶比率和尺度調(diào)整的逆概率加權(quán)法（IPWS）進(jìn)行參數(shù)估計(jì)

設(shè)所需要的SAS程序如下：

將“第3.3.1節(jié)”SAS程序的第1句中的方法選項(xiàng)修改為“method=IPWS”，其他內(nèi)容不變。

【SAS主要輸出結(jié)果及解釋】SAS程序的輸出結(jié)果見表7。

表7是基于傾向性評(píng)分分析后，再基于匹配后的數(shù)據(jù)集并采用IPWS法對(duì)“結(jié)果變量取值為患2型糖尿病”構(gòu)建Logistic回歸模型，得到結(jié)果：Drug_X與Drug_A的潛在結(jié)果均值（POM）分別為0.469與0.647；兩種藥物平均因果效應(yīng)ATE=-0.178，表明藥物X在預(yù)防2型糖尿病方面比藥物A更有效。正如95%置信區(qū)間（-0.276～-0.079）或P值（＜0.01）所示，在0.05水平上，ATE與0之間的差異有統(tǒng)計(jì)學(xué)意義。

表7 因果效應(yīng)的分析結(jié)果Table 7 Analysis results of causal effects

3.3.4 采用增廣逆概率加權(quán)法（AIPW）進(jìn)行參數(shù)估計(jì)

設(shè)所需要的SAS程序如下：

將“第3.3.1節(jié)”SAS程序的第1句中的方法選項(xiàng)修改為“method=AIPW”，其他內(nèi)容不變。

【SAS主要輸出結(jié)果及解釋】SAS程序的輸出結(jié)果見表8。

表8是基于傾向性評(píng)分分析后，再基于匹配后的數(shù)據(jù)集并采用AIPW法對(duì)“結(jié)果變量取值為患2型糖尿病”構(gòu)建Logistic回歸模型，得到結(jié)果：Drug_X與Drug_A的潛在結(jié)果均值（POM）分別為0.478與0.649；兩種藥物平均因果效應(yīng)ATE=-0.171，表明藥物X在預(yù)防2型糖尿病方面比藥物A更有效。正如95%置信區(qū)間（-0.268～-0.074）或P值（0.001）所示，在0.05水平上，ATE與0之間的差異有統(tǒng)計(jì)學(xué)意義。

表8 因果效應(yīng)的分析結(jié)果Table 8 Analysis results of causal effects

3.3.5 采用回歸調(diào)整法（REGADJ）進(jìn)行參數(shù)估計(jì)

設(shè)所需要的SAS程序如下：

將“第3.3.1節(jié)”SAS程序的第1句中的方法選項(xiàng)修改為“method=REGADJ”，其他內(nèi)容不變。

【SAS主要輸出結(jié)果及解釋】SAS程序的輸出結(jié)果見表9。

表9是基于傾向性評(píng)分分析后，再基于匹配后的數(shù)據(jù)集并采用REGADJ法對(duì)“結(jié)果變量取值為患2型糖尿病”構(gòu)建Logistic回歸模型，得到結(jié)果：Drug_X與Drug_A的潛在結(jié)果均值（POM）分別為0.473與0.647；兩種藥物平均因果效應(yīng)ATE=-0.174，表明藥物X在預(yù)防2型糖尿病方面比藥物A更有效。正如95%置信區(qū)間（-0.273～-0.075）或P值（0.001）所示，在0.05水平上，ATE與0之間的差異有統(tǒng)計(jì)學(xué)意義。

表9 因果效應(yīng)的分析結(jié)果Table 9 Analysis results of causal effects

3.3.6 采用逆概率加權(quán)回歸調(diào)整的雙重穩(wěn)健估計(jì)法（IPWREG）進(jìn)行參數(shù)估計(jì)

設(shè)所需要的SAS程序如下：

將“第3.3.1節(jié)”SAS程序的第1句中的方法選項(xiàng)修改為“method=IPWREG”，其他內(nèi)容不變。

【SAS主要輸出結(jié)果及解釋】SAS程序的輸出結(jié)果見表10。

表10是基于傾向性評(píng)分分析后，再基于匹配后的數(shù)據(jù)集并采用IPWREG法對(duì)“結(jié)果變量取值為患2型糖尿病”構(gòu)建Logistic回歸模型，得到結(jié)果：Drug_X與Drug_A的潛在結(jié)果均值（POM）分別為0.479與0.649；兩種藥物平均處理因果效應(yīng)ATE=-0.170，表明藥物X在預(yù)防2型糖尿病方面比藥物A更有效。正如 95% 置信區(qū)間（-0.269～-0.071）或 P值（0.001）所示，在0.05水平上，ATE與0之間的差異有統(tǒng)計(jì)學(xué)意義。

表10 因果效應(yīng)的分析結(jié)果Table 10 Analysis results of causal effects

【說明】因篇幅所限，以上6小節(jié)中“傾向性評(píng)分模型估計(jì)”的輸出結(jié)果均從略。

3.4 結(jié)論

由“第3.2節(jié)”的結(jié)果可知，受試者接受Drug_A治療比接受Drug_X治療，更易患2型糖尿病的優(yōu)勢(shì)比為OR=2.152，它是基于原始數(shù)據(jù)集得到的結(jié)果和結(jié)論。值得一提的是，僅采用多重Logistic回歸分析，無法獲得平均處理因果效應(yīng)的估計(jì)值。

由“第3.3節(jié)”的計(jì)算結(jié)果可知，受試者接受Drug_A治療比接受Drug_X治療，更易患2型糖尿病。基于6種估算方法，得到兩種藥物平均處理因果效應(yīng) ATE 分別為-0.190、-0.179、-0.178、-0.171、-0.174和-0.170。相對(duì)而言，本例資料采用IPW估算法可以獲得絕對(duì)值最大的平均處理因果效應(yīng)值，即ATE=-0.190。它是基于匹配后數(shù)據(jù)集得到的結(jié)果。

4 討論與小結(jié)

4.1 討論

通常在進(jìn)行多重回歸分析時(shí)，假定所有自變量是同時(shí)被觀測(cè)的、互相獨(dú)立的，完全基于各變量的觀測(cè)結(jié)果并采用統(tǒng)計(jì)學(xué)處理技術(shù)（包括異常點(diǎn)的診斷與處理、共線性診斷與處理）來篩選自變量。事實(shí)上，在實(shí)際問題中，自變量具有許多不同的特性，例如，從時(shí)間角度看，自變量的觀測(cè)有先后順序之分，而不是在同一個(gè)時(shí)間點(diǎn)或時(shí)間段上收集數(shù)據(jù)；從變量的相互關(guān)系角度看，有些自變量之間不僅不是相互獨(dú)立的，還可能存在“中介變量”［1，5］；從研究目的角度看，有些自變量是研究者特別關(guān)注的或?qū)ｉT施加的，即“處理變量”［1，6］。另外，在非隨機(jī)對(duì)照研究的資料中，處理變量?jī)伤浇M的協(xié)變量向量之間常不具有可比性，提高協(xié)變量在處理變量?jī)伤浇M之間均衡性的一個(gè)有效方法就是采用傾向性評(píng)分分析法［7］，通常直接對(duì)原始數(shù)據(jù)集進(jìn)行回歸分析是不合適的。針對(duì)具有不同特性的自變量，在進(jìn)行多重回歸分析的過程中，需要采取相應(yīng)的處置策略或添加相應(yīng)的分析內(nèi)容，以便獲得更多的、更合理的、更有價(jià)值的分析結(jié)果。

4.2 小結(jié)

本文介紹了與潛在結(jié)果均值和平均處理效應(yīng)有關(guān)的基本概念、估算方法、實(shí)例及SAS實(shí)現(xiàn)?；靖拍钌婕啊皾撛诮Y(jié)果與潛在結(jié)果均值”“平均處理效應(yīng)”“處理組中的平均處理效應(yīng)”和“逆概率權(quán)重”；估算方法有6種，它們分別是IPW估算法、IPWR估算法、IPWS估算法、AIPW估算法、IPWREG估算法和REGADJ估算法；針對(duì)一個(gè)假設(shè)的藥物臨床試驗(yàn)資料，用SAS軟件實(shí)現(xiàn)了多重Logistic回歸分析，并采用6種估算方法估計(jì)了潛在結(jié)果均值和平均處理效應(yīng)。