如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料XXII. 結(jié)果變量為二值變量的高維列聯(lián)表資料的統(tǒng)計分析與SAS軟件實現(xiàn)（二）

鮑曉蕾，胡良平

?

如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料XXII. 結(jié)果變量為二值變量的高維列聯(lián)表資料的統(tǒng)計分析與SAS軟件實現(xiàn)（二）

鮑曉蕾，胡良平

100850 北京，軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心

編者按

生物統(tǒng)計學(xué)是生物學(xué)領(lǐng)域科學(xué)研究和實際工作中必不可少的工具，在分子生物學(xué)迅速發(fā)展的今天，生物統(tǒng)計學(xué)更顯示出了它的重要性。實驗設(shè)計與數(shù)據(jù)統(tǒng)計分析是現(xiàn)代生物學(xué)的基石，是生物學(xué)研究者檢驗假說、尋找模式、建立生物學(xué)理論的有利工具，也是生物學(xué)研究者探索微觀和宏觀生物世界的必備基礎(chǔ)知識。對于每天甚至是每時每刻涌現(xiàn)的大量的、以天文數(shù)字計量的分子遺傳數(shù)據(jù)，必須借助統(tǒng)計學(xué)知識加以分析處理，才能從中獲得有意義的信息?！吧锒鄻有詳?shù)據(jù)分析”是開展生物多樣性研究的一個重要方面，數(shù)據(jù)分析能力的高低極大地影響著我們對各種生態(tài)學(xué)現(xiàn)象認(rèn)識的深度和廣度。現(xiàn)在，電子計算機(jī)的普及使得生物統(tǒng)計分析過程大大簡化，生物統(tǒng)計分析軟件包的普及將生物統(tǒng)計學(xué)從統(tǒng)計學(xué)家的書本里解放了出來，簡化了生物統(tǒng)計分析過程，使之成為生物學(xué)研究者的常用工具。本刊特邀軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計學(xué)咨詢中心主任胡良平教授，以“如何用SAS軟件正確分析生物醫(yī)學(xué)科研資料”為題，撰寫系列統(tǒng)計學(xué)講座，希望該系列講座能對生物醫(yī)學(xué)科研工作者有所幫助。

當(dāng)觀測結(jié)果是定性資料時，人們習(xí)慣將資料整理成列聯(lián)表形式。比如“2′2 列聯(lián)表資料”、“R′C 列聯(lián)表資料”和“高維列聯(lián)表資料”等。所謂高維列聯(lián)表，也就是表中涉及到的定性變量的個數(shù)≥ 3。對于高維列聯(lián)表資料，根據(jù)結(jié)果變量的性質(zhì)可將高維列聯(lián)表分為以下三類：一是結(jié)果變量為二值變量的高維列聯(lián)表；二是結(jié)果變量為多值有序變量的高維列聯(lián)表；三是結(jié)果變量為多值名義變量的高維列聯(lián)表。在上一期中，我們已經(jīng)介紹了用CMH2檢驗和加權(quán)2檢驗處理結(jié)果變量為二值變量的高維列聯(lián)表資料，本文將繼續(xù)介紹結(jié)果變量為二值變量的高維列聯(lián)表資料的其他分析方法及 SAS 實現(xiàn)。

處理不同的高維列聯(lián)表資料所用的統(tǒng)計分析方法不盡相同，對于結(jié)果變量為二值變量的高維列聯(lián)表資料可使用 CMH2檢驗、多重 logistic 回歸分析、對數(shù)線性模型等統(tǒng)計分析方法。本文重點介紹對數(shù)線性模型和多重 logistic 回歸分析如何處理結(jié)果變量為二值變量的高維列聯(lián)表資料。

1 對數(shù)線性模型

在上一期中，我們介紹了用 CMH2檢驗和加權(quán)2檢驗來處理結(jié)果變量為二值變量的高維列聯(lián)表資料，這兩種方法的本質(zhì)都是控制其中一個原因變量對結(jié)果的影響，從而分析另一個原因變量的作用。這樣的做法往往不夠準(zhǔn)確，特別是當(dāng)被合并的變量是一個混雜因素，即與留下的變量之間不獨立時，極易得出錯誤的結(jié)論。那么，當(dāng)同時分析各變量的主效應(yīng)及變量間的交互效應(yīng)時，對數(shù)線性模型便顯示出其巨大的潛力。

對數(shù)線性模型是分析高維列聯(lián)表資料行之有效的方法，最先由 Yule、Bartlett 利用 Yule（1900 年）定義的交叉乘積比分析三維交互作用，然后在 Kullback（1968 年）引入方差分析的思想發(fā)展而來[1]。

對數(shù)線性模型的基本思想：把各分組變量（包括自變量和因變量）水平組合下期望（理論頻數(shù)）的自然對數(shù)表示為各分組變量及其交互作用項和定性結(jié)果變量的線性函數(shù)，通過迭代計算求得模型中參數(shù)的估計值，進(jìn)而運用方差分析的思想檢驗各主效應(yīng)和交互作用的效應(yīng)大小[1]。

主效應(yīng)與交互效應(yīng)：對數(shù)線性模型主要考察各分類變量間的交互作用（關(guān)聯(lián)性），對主效應(yīng)的分析相當(dāng)于單變量分析。交互效應(yīng)按其因素多少可分為兩變量間交互效應(yīng)和多變量間交互效應(yīng)，它們依次又被稱為一階交互效應(yīng)、二階交互效應(yīng)……，以此類推。

對數(shù)線性模型參照線性模型和方差分析的思想，將列聯(lián)表的總效應(yīng)（期望頻數(shù)）取自然對數(shù)后分解成主效應(yīng)和各因素間的交互效應(yīng)。對數(shù)線性模型雖與線性模型和方差分析類似，但又有其自身特點。簡單起見，以三維列聯(lián)表為例介紹對數(shù)線性模型的構(gòu)建原理，對于四維乃至更高維度的列聯(lián)表，不失一般性。

設(shè) 3 個分類變量分別為 A、B、C，則包含所有交互效應(yīng)的飽和模型為：

只包含 A 與 B、A 與 C 交互效應(yīng)的譜系模型為：

表 1 對數(shù)線性模型參數(shù)意義與個數(shù)

2 多重 logistic 回歸分析

Logistic 回歸分析是 1970 年 Cox 提出的[1]。設(shè)（y = 1│X）（簡記為）表示暴露因素為 X 時個體發(fā)病的概率。稱發(fā)病的概率與未發(fā)病的概率 1-之比為“優(yōu)勢（odds）”。對作 logit 變換，logit() 定義為優(yōu)勢之對數(shù)（log odds），即：

因此，多重 logistic 回歸方程的線性表達(dá)式定義為：

logit() =+βx+ βx+ … +βx

logit()與各因素間呈線性關(guān)系，β為 logistic 回歸的回歸系數(shù)，表示在其他變量都固定的情況下，對即 logit() 影響的大小。取β的反對數(shù)可得優(yōu)勢比= exp(β)，表示其他變量不變的情況下，每變化一個單位時所引起的優(yōu)勢比的自然對數(shù)改變量。當(dāng)某種疾病的發(fā)病率或死亡率很低時，≈（即相對危險度）[1]。

同多重線性回歸分析一樣，當(dāng)比較暴露因素對因變量相對貢獻(xiàn)大小時，由于各自變量取值單位不同，也不能用回歸系數(shù)的大小直接進(jìn)行比較，而需用標(biāo)準(zhǔn)化回歸系數(shù)來比較。建立 logistic 回歸方程的過程也就是求常數(shù)項及各回歸系數(shù)β的過程。

Logistic 回歸分析按照因變量（即反應(yīng)變量）的類型可分為：因變量為二值變量的 logistic 回歸分析、因變量為多值有序變量的 logistic 回歸分析、因變量為多值無序（或名義）變量的 logistic 回歸分析；logistic 回歸按研究設(shè)計類型可分為非條件（成組設(shè)計）logistic 回歸分析和條件（匹配設(shè)計）logistic 回歸分析。

對于 logistic 回歸方程的擬合一般利用 SAS 中l(wèi)ogistic 過程完成，一些較復(fù)雜的logistic 回歸分析則需要通過調(diào)用 SAS 中的 catmod 過程或 phreg 過程完成。本文僅介紹調(diào)用 SAS 中的 logistic 過程完成一般的 logistic 回歸分析。

下面仍沿用上一期的實例，向讀者介紹如何通過 SAS 軟件使用對數(shù)線性模型和多重logistic 回歸模型處理結(jié)果變量為二值變量的高維列聯(lián)表資料。

【例 1】 為調(diào)查眼科門診患者的用藥依從性及其影響因素，并提出相應(yīng)對策，提高眼科門診藥房的藥學(xué)服務(wù)水平，在 5 家醫(yī)院的眼科門診患者中，隨機(jī)發(fā)放 610 份調(diào)查問卷，進(jìn)行匯總，數(shù)據(jù)見表 2。

表 2 不同告知者、不同告知程度間患者的依從情況

⑴對數(shù)線性模型

SAS 程序如下，程序名為li3.sas。

data a1;do a=1 to 2;do b=1 to 3;do c=1 to 2;input f @@;output;end; end; end; cards;240 24 202 60 31 16 188 19 238 60 47 21;run;ods html;proc catmod;weight f;model a*b*c=_response_/noglsnoparm noresponse ml;loglin a|b|c;run;ods html close;

程序說明：a 表示告知者，a = 1 表示醫(yī)生，a = 2 表示藥師；b 表示告知程度，b = 1 表示詳細(xì)告知，b = 2 表示告知，b = 3 表示沒告知；c 表示是否依從，c = 1 表示依從，c = 2 表示不依從。過程步調(diào)用 catmod 過程，若需輸出對數(shù)線性模型中的各參數(shù)的估計值和誤差，只需去掉 model 語句中的 noparm 選項。ml 選項表示采用最大似然法進(jìn)行參數(shù)估計，nogls 表示不采用加權(quán)最小二乘法計算。Noresponse 表示不輸出反應(yīng)矩陣。Loglin 定義線性模型的效應(yīng)項。

SAS 程序運行結(jié)果：

Maximum likelihood analysis of variance

SourceDFChi-SquarePr > ChiSq a1 0.450.5014 b2173.20< .0001 a*b2 5.420.0665 c1241.62< .0001 a*c1 0.290.5929 b*c2 42.45< .0001 a*b*c2 0.210.8988 Likelihood Ratio0..

結(jié)果顯示，b 與 c 的交互作用存在統(tǒng)計學(xué)意義（< 0.0001），說明患者被告知的程度不同，其依從性不同。而因素 a（告知者）對結(jié)果的影響沒有統(tǒng)計學(xué)意義，此時把不同告知者的數(shù)據(jù)進(jìn)行簡單合并后進(jìn)行分析是可以接受的。若 a*b*c 這個高階交互作用有統(tǒng)計學(xué)意義，則說明不同告知者、不同告知程度下患者的依從性有差別，此時不適合用壓縮法消除“告知者”的影響而研究“告知程度”對患者依從性的作用。

專業(yè)結(jié)論：患者被告知的程度不同，依從性不同（實際資料顯示：被詳細(xì)告知者依從性高）；而不同的告知者對患者依從性的影響沒有統(tǒng)計學(xué)意義。

⑵多重 logistic 回歸分析

程序如下，設(shè)程序名為 li4.sas。

data a2;do a=1 to 2;do b=1 to 3;do c=1 to 2;input f @@;output;end; end; end; cards;240 24 202 60 31 16 188 19 238 60 47 21;run;ods html;proc logistic;class b/param=ref;freq f; model c=a b a*b/selection=forward;oddsratio b;run;proc logistic;class b/param=ref;freq f;model c=a b a*b/selection=backward;oddsratio b;run;proc logistic;class b/param=ref;freq f;model c=a b a*b/selection=stepwise;oddsratio b;run;ods html close;

程序說明：數(shù)據(jù)步與程序 li3.sas 一樣，不再重復(fù)。調(diào)用 logistic 過程進(jìn)行結(jié)果變量為二值變量的多重 logistic 回歸分析。Class 語句指定對變量 b 賦啞變量，param = ref表示以其中的一個水平為基準(zhǔn)建立啞變量，默認(rèn)以輸入的最后一個水平為基準(zhǔn)，本例則以變量b 的最后一個水平（即“未告知”）為基準(zhǔn)建立啞變量。語句freq f 指定f 變量為頻數(shù)變量。通過model 語句進(jìn)行建模，等號前面表示因變量，等號后面表示自變量，其中a*b 表示變量a 和b 的交互作用，斜杠后面的selection 選項表示進(jìn)行變量篩選。程序中的三個logistic 過程分別表示用forward（前進(jìn)法）、backward（后退法）和stepwise（逐步篩選法）進(jìn)行變量篩選，運用多種變量篩選方法的目的是為了使篩選結(jié)果更可信。第一個logistic 過程中selection = forward 表示用前進(jìn)法進(jìn)行變量篩選，其中forward 法默認(rèn)的入選變量進(jìn)入方程的檢驗水準(zhǔn)為0.05，若想自定義該水準(zhǔn)，可在該語句后加入選項sle = xx 來指定；第二個logistic 過程中selection = backward 表示用后退法進(jìn)行變量篩選，其默認(rèn)的從方程中剔除變量的檢驗水準(zhǔn)為0.05，若想自定義該水準(zhǔn)，可在該語句后加入選項sls = xx 來指定；第三個logistic 過程中selection = stepwise 表示用逐步篩選法進(jìn)行變量篩選，其默認(rèn)的入選水準(zhǔn)和排除水準(zhǔn)分別為0.05，若想自定義該值，可在該語句后加入sle = xx sls = xx 來指定，前者表示入選水準(zhǔn)，后者表示剔除水準(zhǔn)，里面的“xx”實際使用時都應(yīng)給予具體的數(shù)值，如0.3、0.1、0.03、0.01 等。oddsratio b 表示輸出變量b 任意兩個水平間的優(yōu)勢比，若不寫該語句，程序?qū)⒆詣虞敵鲎兞縝 各水平相對于基準(zhǔn)水平的優(yōu)勢比。

SAS 程序運行結(jié)果：

Class level information

ClassValueDesign variables b110 201 300

以上結(jié)果為通過class 語句對變量b 賦啞變量的結(jié)果，以第三水平（即“未告知”）為基準(zhǔn)，其賦值方式與寫法b1 = 0; b2 = 0; if b = 1 then b1 = 1; else if b = 2 then b2 = 1;等價。

Summary of forward selection

StepEffect enteredDFNumber inScore Chi-SquarePr > ChiSq 1b2146.0932< .0001

以上為前進(jìn)法（forward）篩選變量的結(jié)果，只有變量b入選方程（< 0.0001）。

Summary of backward elimination

StepEffect removedDFNumber inWald Chi-SquarePr > ChiSq 1a*b220.21330.8988 2a110.54780.4592

以上為后退法（backward）篩選變量的結(jié)果，因素a 以及a 和b 的交互作用，a*b 都被排除在方程外（值分別為0.4592 和0.8988），說明只有變量b 保留在方程中，該結(jié)果與前進(jìn)法的篩選結(jié)果一致。

Summary of stepwise selection

StepEffectDFNumber inScore Chi-SquareWald Chi-SquarePr > ChiSq EnteredRemoved 1b 2146.0932 < .0001

以上為逐步篩選法（stepwise）篩選變量的結(jié)果，只有因素b 最終保留在方程中（< 0.0001），該結(jié)果與前進(jìn)法、后退法的篩選結(jié)果一致。綜合以上三種變量篩選法的結(jié)果，可以認(rèn)為變量b 對結(jié)果的影響有統(tǒng)計學(xué)意義。

Type 3 analysis of effects

EffectDFWald Chi-SquarePr > ChiSq b242.6689< .0001

以上為對變量b 進(jìn)行整體假設(shè)檢驗的結(jié)果，wald2= 42.6689，< 0.0001，說明變量b 對結(jié)果的影響有統(tǒng)計學(xué)意義。

Analysis of maximum likelihood estimates

Parameter DFEstimateStandarderrorWaldChi-SquarePr > ChiSq Intercept 10.74580.199613.95810.0002 b111.55210.255836.8130< .0001 b210.55350.22466.07220.0137

以上是關(guān)于截距項以及變量b 產(chǎn)生的兩個啞變量的參數(shù)估計及假設(shè)檢驗的結(jié)果。截距項的參數(shù)估計值為0.7458，= 0.0002 < 0.05；啞變量b1的參數(shù)估計值為1.5521，< 0.0001；啞變量b2的參數(shù)估計值為0.5535，= 0.0137 < 0.05。據(jù)此可以寫出該模型的表達(dá)式為：

Odds Ratio Estimates

EffectPoint estimate95% Wald confidence limits b 1 vs 34.7212.8607.795 b 2 vs 31.7391.1202.701

以上為程序默認(rèn)的優(yōu)勢比的輸出結(jié)果，給出了告知程度的1 水平（詳細(xì)告知）和2 水平（告知）相對于3 水平（未告知）的優(yōu)勢比。1 水平（詳細(xì)告知）相對于3 水平（未告知）的優(yōu)勢比OR= 4.721，其95% 的置信區(qū)間為（2.860，7.795），說明被詳細(xì)告知患者的依從性是未被告知患者依從性的4.721 倍。2 水平（告知）相對于3 水平（未告知）的優(yōu)勢比OR= 1.739，其95% 置信區(qū)間為（1.120，2.701），說明被告知患者的依從性是未被告知患者依從性的1.739 倍。

Wald confidence interval for odds ratios

LabelEstimate95% Confidence limits b 1 vs 22.7141.8703.941 b 1 vs 34.7212.8607.795 b 2 vs 31.7391.1202.701

以上為語句oddsratio b 對應(yīng)的輸出結(jié)果。可以看到，除了給出變量b（告知程度）1 水平（詳細(xì)告知）和2 水平（告知）相對于3 水平（未告知）的優(yōu)勢比外，還給出了1 水平（詳細(xì)告知）相對于2 水平（告知）的優(yōu)勢比OR，比值為2.714，其95% 置信區(qū)間為（1.870，3.941），說明被詳細(xì)告知患者的依從性是被告知患者依從性的2.714 倍。

專業(yè)結(jié)論：根據(jù)logistic 回歸分析的結(jié)果，告知者對患者依從性沒有影響，告知程度對患者的依從性有影響（< 0.0001）。OR= 2.714，其95% 置信區(qū)間為（0.352，0.567），說明被詳細(xì)告知患者的依從性是被告知患者的2.714 倍；OR= 4.721，其95% 置信區(qū)間為（2.860，7.795），說明被詳細(xì)告知患者的依從性是未被告知患者依從性的4.721 倍；OR= 1.739，其95% 置信區(qū)間為（1.120，2.701），說明被告知患者的依從性是未被告知患者依從性的1.739 倍。

[1] Hu LP. Medical statistics-analysis of quantitative and qualitative data applying the triple-type theory. Beijing: People’s Military Medical Press, 2009:325-351. (in Chinese)

胡良平. 醫(yī)學(xué)統(tǒng)計學(xué)-運用三型理論分析定量與定性資料. 北京: 人民軍醫(yī)出版社, 2009:325-351.

[2] Hu LP. Medical statistics-analysis of regression models using the triple-type theory. Beijing: People’s Military Medical Press, 2010:203- 220. (in Chinese)

胡良平. 醫(yī)學(xué)統(tǒng)計學(xué)-運用三型理論進(jìn)行現(xiàn)代回歸分析. 北京: 人民軍醫(yī)出版社, 2010:203-220.

10.3969/cmba.j.issn.1673-713X.2013.02.016

胡良平，Email：lphu812@sina.com