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

王琪，胡良平

針對高維列聯(lián)表（表中涉及到的定性變量的個數(shù) k ≥3）資料的分析，前四期的講座分別介紹了結果變量為二值變量的高維列聯(lián)表和結果變量為多值有序變量的高維列聯(lián)表資料的統(tǒng)計分析與SAS 軟件實現(xiàn)，本文將詳細介紹結果變量為多值名義變量的高維列聯(lián)表資料及其用 SAS 軟件實現(xiàn)統(tǒng)計分析的內(nèi)容。

分析結果變量為多值名義變量的高維列聯(lián)表資料時，可以選用的統(tǒng)計分析方法有CMH χ2檢驗、對數(shù)線性模型和擴展的多重 logistic 回歸分析，本文重點介紹前兩種方法。

1 CMH χ2 檢驗[1]

在高維列聯(lián)表中，當結果變量為多值名義變量時，選用一般關聯(lián)統(tǒng)計量，也就是 FREQ過程 CMH 檢驗輸出結果中的第三項，其自由度 ν=(R-1) (C-1)。它的具體計算公式可以參見本刊 2013年第 3 期，在計算該統(tǒng)計量時，行的評分陣 Rh與計算行平均得分統(tǒng)計量的行評分矩陣相同，即Rh=[IR－1,–JR－1]；列的評分陣 Ch可以被類似地定義為：

其中 IC－1是秩為C－1的單位陣，JC－1是元素均為1的(C－1)×1的列向量。在一般關聯(lián)統(tǒng)計量中，行的評分陣 Rh與列的評分陣 Ch都是由 FREQ 過程內(nèi)部產(chǎn)生的。

一般關聯(lián)統(tǒng)計量不要求原因變量或結果變量是有序的。無論原因變量是多值有序變量還是名義變量，只要結果變量是多值名義的，都可以采用該統(tǒng)計量。與之相對應的原假設和備擇假設分別為：

H0：每層中原因變量和結果變量之間不存在關聯(lián)；

H1：至少有一層，原因變量和結果變量之間存在某種關聯(lián)。

當僅有一層時，該 CMH 統(tǒng)計量與Pearson χ2統(tǒng)計量的關系為：

其中 n 為總例數(shù)；當有多層時，該統(tǒng)計量為層修正的Pearson χ2統(tǒng)計量。當然，相似的校正也能夠通過對各層Pearson χ2統(tǒng)計量求和而得到，但是這種校正方法需要每層的樣本含量都要足夠大，而 CMH 統(tǒng)計量僅僅需要總的樣本含量比較大。

需要說明的是，當各層的各比較組間的趨勢方向一致時，CMH 方法比較有效；當各層的各比較組間的趨勢方向不一致時，CMH 方法則不容易檢出差別，此時應單獨考察各層或采用其他方法。

2 對數(shù)線性模型

結果變量為多值名義變量的高維列聯(lián)表資料同樣也可以采用對數(shù)線性模型進行分析，其基本的原理與分析結果變量為二值變量的高維列聯(lián)表資料基本相同[2]，詳見本刊 2013年第 2 期，本文僅結合實例向讀者介紹如何通過SAS 軟件使用 CMH χ2檢驗和對數(shù)線性模型處理結果變量為多值名義變量的高維列聯(lián)表資料。對軟件實現(xiàn)和結果解釋作進一步說明。

【例1】調(diào)查某中醫(yī)院一日內(nèi)醫(yī)生開出的針對甲、乙兩種疾病的處方情況，結果見表1，試對數(shù)據(jù)進行分析。

表1 不同疾病、不同性別患者的藥物使用情況

分析與解答：此表中含有三個定性變量，分別為疾病類型、患者性別、藥物種類，結果變量為藥物種類（多值名義變量），數(shù)據(jù)以列聯(lián)表的形式呈現(xiàn)，因此該表被稱作“結果變量為多值名義變量的三維列聯(lián)表”。對于該高維列聯(lián)表資料，分析目的是考察不同疾病、不同性別的患者所用藥物種類頻數(shù)構成有無差別，分析時通?？蛇x用 CMH χ2檢驗、對數(shù)線性模型或擴展的多重 logistic 回歸分析。本期重點介紹前兩種。

⑴ CMH χ2檢驗：CMH χ2檢驗可以控制外側影響相對次要的變量而得到最內(nèi)側兩個變量之間的關系。如果我們想控制住性別的影響，分析疾病類型與藥物種類之間的關系，此時可以使用 CMH χ2檢驗。SAS 程序如下，設程序名為li1.sas。

data a1; /*1*/do a=1 to 2;do b=1 to 2;do c=1 to 3;input f @@;output;end; end; end;cards;52 7 234 8 123 19 418 11 3;run;ods html;proc freq data=a1; /*2*/tables b*a*c/cmh;weight f;run;ods html close;

程序說明：本例第一步建立數(shù)據(jù)集 a1，a 表示疾病類型，a=1 表示甲，a=2 表示乙；b 表示患者性別，b=1 表示男性，b=2 表示女性；c 表示藥物種類，c=1 表示中藥，c=2 表示西藥，c=3 表示復方藥；變量 f 表示頻數(shù)。多值名義資料的CMH 檢驗仍然采用 FREQ 過程，在tables語句中依次列出性別、疾病類型和藥物種類，列在第一位的變量是需要控制的原因變量，列在第二位的變量是想要考察的原因變量，列在第三位的變量是結果變量。Tables 語句中的選項 cmh 指定輸出CMH 統(tǒng)計量。ods html 語句則要求以 HTML（網(wǎng)頁）格式輸出結果。

SAS 程序運行結果：首先輸出的是兩個二維列聯(lián)表（略），這是按性別分層以后，疾病類型和藥物種類所形成的2×3 列聯(lián)表，其中包括頻數(shù)、百分比、行百分比和列百分比。

FREQ 過程

輸出結果的第二部分是三個 CMH 統(tǒng)計量，這些統(tǒng)計量都是在使用 table 方法打分的基礎上計算的。本例中結果變量是多值名義的，因此無論原因變量是二值的、多值有序的或多值名義的，都應該使用一般關聯(lián)統(tǒng)計量。此處自由度ν=2=19.0157，P＜0.0001?？偟臉颖竞繛?82 例。

統(tǒng)計學與專業(yè)結論：本例CMH 檢驗中的P＜0.05，拒絕原假設 H0，接受 H1，說明在控制了性別因素的基礎上，疾病類型與藥物種類之間存在一定的關聯(lián)。結合列聯(lián)表中的行百分數(shù)可以看到，相對來說，甲病患者中使用中藥的人數(shù)比例要大于乙病患者，也就是說甲病患者的中藥使用率比乙病患者高。

⑵對數(shù)線性模型：若想要考察疾病類型、性別和藥物種類三者之間的關系，可以考慮采用對數(shù)線性模型。SAS 程序如下，設程序名為li2.sas。

data a2; /*1*/do a=1 to 2;do b=1 to 2;do c=1 to 3;input f @@;output;end; end; end;cards;52 7 234 8 123 19 418 11 3;run;ods html;proc catmod data=a2; /*2*/weight f;model a*b*c=_response_;loglin a|b|c;run;proc catmod data=a2; /*3*/weight f;model a*b*c=_response_;loglin a|b|c @2;run;proc catmod data=a2; /*4*/weight f;model a*b*c=_response_;loglin a b c a*c;run;ods html close;

程序說明：本例第一步建立數(shù)據(jù)集 a2，同 li1.sas。然后第二步調(diào)用 catmod 過程，采用飽和模型進行對數(shù)線性模型分析，“model a*b*c=_response_;”語句中等號左端的a*b*c 指明要分析的變量，等號右端的_response_ 表示擬合對數(shù)線性模型；“l(fā)oglin a|b|c;”語句中的a|b|c 表示擬合所有的主效應與交互效應，即飽和模型，該語句等同于“l(fā)oglin a b c a*b a*c b*c a*b*c”。第三步調(diào)用 catmod 過程，采用一階交互效應模型進行對數(shù)線性模型分析，該模型包括主效應和所有的一階交互效應（即在飽和模型中去掉最后一項，它是二階交互效應項），“l(fā)oglin a|b|c @2;”語句表示擬合一階交互效應模型。第四步調(diào)用 catmod 過程，擬合最優(yōu)模型，語句“l(fā)oglin a b c a*c;”中包含了主效應和a*c 一階交互效應。

SAS 程序運行結果：

⑴飽和模型：

Maximum likelihood analysis of variance

以上是方差分析表，使用最大似然法進行分析。這里三個因素之間的二階交互作用 a*b*c 經(jīng)檢驗沒有統(tǒng)計學意義（P ＞0.05），一階交互作用中僅 a*c 有統(tǒng)計學意義（P＜0.05）。由于這里建立的是飽和模型，已經(jīng)沒有剩余的自由度分配給似然比檢驗，所以并沒有關于模型擬合情況的似然比檢驗的結果。

Analysis of maximum likelihood estimates

上表給出了模型中參數(shù)的估計值以及對其進行假設檢驗的結果。這部分結果與方差分析表的結果一致。參數(shù)估計值默認將每個效應的最后一類作為參照類，其他各個水平通過與參照類相比來分析效應大小。如疾病類型與藥物種類的交互效應（a*c）在兩個變量取值都為1 時的參數(shù)估計值為0.5010（P=0.0006），與0的差別有統(tǒng)計學意義，表示甲病組（a=1）與乙病組（a=2）相比，使用中藥的傾向性存在差別，甲病患者比乙病患者更傾向于使用中藥。由于飽和模型無法進行假設檢驗，且二階交互效應項也無統(tǒng)計學意義，可以考慮首先將此項從模型中刪除，以簡化模型。

⑵一階交互效應模型：

Maximum likelihood analysis of variance

一階交互效應模型僅列出方差分析表，結果提示，疾病類型是患者所用藥物種類的影響因素。似然比檢驗的結果為χ2=1.42，P=0.4905（此結果表明該模型對資料的擬合效果較好）?？紤]到模型中還包含有無統(tǒng)計學意義的項，需要進一步優(yōu)化。

⑶最優(yōu)模型的篩選：

Maximum likelihood analysis of variance

由于在一階交互效應模型擬合過程中，a*b 項和b*c項無統(tǒng)計學意義，我們將其剔除出模型。最終的模型包括的一階交互作用僅有 a*c，方差分析的結果表明在新模型中a*c 項的作用仍然是有統(tǒng)計學意義的。似然比檢驗的結果χ2=1.46，P=0.9176 ＞0.05，說明模型對資料的擬合效果非常好，可以認為該模型不僅成立而且已達到最優(yōu)程度。盡管因素 a的假設檢驗結果無統(tǒng)計學意義，但由于對數(shù)線性模型是嵌套式結構，只要模型中有高階項，所有比它低階的項都必須保留在模型中。

Analysis of maximum likelihood estimates

在參數(shù)估計的結果中，疾病類型和患者性別有兩個水平，自由度均為1，所以參數(shù)估計部分對應一行結果；而藥物種類有三個水平，自由度為2，因此參數(shù)估計的結果有兩行，分別代表中藥相對于復方藥、西藥相對于復方藥對相應網(wǎng)格中理論頻數(shù)的對數(shù)值的影響情況[前者差別有統(tǒng)計學意義（P＜0.0001）、后者差別無統(tǒng)計學意義（P=0.2922）]；疾病類型和藥物種類兩個因素交互作用項的參數(shù)估計結果也有兩行，其結果解釋參見下面的專業(yè)結論部分。

專業(yè)結論：由于疾病類型和藥物種類的交互作用項有統(tǒng)計學意義，同時結合實際資料可知，在甲病和乙病患者中，三種藥物使用的構成比不同，在甲病患者中，使用中藥所占的比例要明顯高于乙病患者；而在乙病患者中，使用中藥和使用西藥的比例沒有明顯的區(qū)別。

[1]Hu LP.Medical statistics-analysis of quantitative and qualitative data applying the triple-type theory.Beijing: People’s Military Medical Press, 2009:376-393.(in Chinese)胡良平.醫(yī)學統(tǒng)計學-運用三型理論分析定量與定性資料.北京:人民軍醫(yī)出版社, 2009:376-393.

[2]Hu LP.Statistics facing practical scientific issues -- (2) multi-factor designs and linear model analysis.Beijing: People’s Medical Publishing House, 2012:547-566.(in Chinese)胡良平.面向問題的統(tǒng)計學——(2)多因素設計與線性模型分析.北京:人民衛(wèi)生出版社, 2012:547-566.