胡純嚴 ,胡良平 ,2*
(1.軍事科學院研究生院,北京 100850;2.世界中醫藥學會聯合會臨床科研統計學專業委員會,北京 100029*通信作者:胡良平,E-mail:lphu927@163.com)
在橫斷面設計的二維列聯表資料或稱為R×C列聯表資料中,根據兩變量是否為有序變量,可以歸納出三種不同的二維列聯表資料:①雙向無序R×C列聯表資料;②結果變量為有序變量R×C列聯表資料;③雙向有序且屬性不同R×C列聯表資料。對其進行統計分析的方法分別為“獨立性假設檢驗”“秩和檢驗”和“相關分析”,它們可以被概括成一種統計分析方法,即CMH χ2檢驗。本文介紹廣義CMH χ2檢驗統計量及其三種變形,并基于SAS軟件實現統計計算。
【例1】文獻[1]中有一個雙向無序R×C表資料(說明:不考慮“時間”的有序性),見表1,試分析“條目”與“時間”之間是否存在關聯性。

表1 基層精防醫護人員K6評定結果
【例2】文獻[2]給出了如下臨床資料:比較三種對肥厚性鼻炎(HR)治療方法的效果。將162例HR患者分為三組,A組(n=56)行下鼻甲骨黏骨膜下切除術,B組(n=43)行下鼻甲部分切除術,C組(n=63)行下鼻甲黏膜下微波熱凝術。治療效果見表2。試對三組治療效果進行比較。

表2 三組下鼻甲手術治療的臨床效果
【例3】在文獻[3]中,為探討新冠肺炎疫情期間居民急性應激障礙(ASD)癥狀的檢出情況及影響因素,得到有效問卷16 048份,有效問卷回收率為72.83%。在該研究中,因素“媒體暴露情況(時間)”與結果變量(出現ASD癥狀的程度)之間的數量關系如表3所示。試分析“媒體暴露情況(時間)”與“ASD癥狀的程度”之間是否存在線性相關關系。

表3 媒體暴露情況(時間)與ASD癥狀檢出情況的調查結果
“CMH”是“Cochran-Mantel-Haenszel”的縮寫,即采用三位作者姓名的第一個字母來命名該分析方法 ,實際上,Cochran、Mantel、Haenszel、Mantel、Birch、Landis、Heyman和Koch都對該分析方法做出過貢獻[4]。
一般來說,CMHχ2檢驗統計量是為分析高維列聯表資料而構造出來的。將高維列聯表資料按某一個或多個因素進行分層,每層都應該是規模相同的R×C列聯表資料。下面以表4形式表示第h層的R×C表,h=1、2、…、q。q為層數,R為行數,C為列數。

表4 第h層R×C列聯表的列表格式
廣義CMHχ2檢驗統計量[4]定義如下:

需注意的是,當各層間效應方向不一致時,CMHχ2檢驗統計量的檢驗功效很低。
2.3.1 概述
式(1)是以矩陣形式呈現的計算公式,很不直觀。由于它試圖達到高維列聯表資料的多種不同分析目的(例如獨立性分析、相關性分析、差異性分析),故其隱含諸多前提條件,包括“是否存在有序變量”“如何給有序變量賦值”。基于不同的前提條件,就會產生出不同的統計計算公式。為了直觀起見,現假定只有一層,即一個二維R×C列聯表資料,基于不同的前提條件,將式(1)轉變成三種具體的、直觀的檢驗統計量。
2.3.2 與備擇假設為“非零相關”對應的檢驗統計量
若擬分析的資料如本文表3(即雙向有序R×C列聯表資料)所示,此時式(1)就簡化成下式:

2.3.3 與備擇假設為“行平均評分不同”對應的檢驗統計量
若擬分析的資料如本文表2(即結果變量為有序變量R×C列聯表資料)所示,此時式(1)就簡化成下面兩個式子之一:

在式(3)中,“∝”代表“呈正比”;當給R×C列聯表的列變量(即結果變量)賦值或評分時,若按“表評分”方式評分,就采用式(3)(即單因素R水平設計一元定量資料方差分析)計算[4,7];否則[包括“rank(秩)評分”“Ridit評分”和“修正的Ridit評分”],就采用式(4)(即單因素R水平設計一元定量資料Kruskal-Wallis’s H秩和檢驗)計算[4,7]。在式(4)中,ni與Mi分別代表第i組中的“樣本含量”與“秩和”。
在式(3)中,可對F統計量進行變換,見下式:

通過對式(5)變形,可以得到下式[8-9]:

式(6)的含義是:采用單因素R水平設計一元定量資料方差分析得到的檢驗統計量F值后,乘以組間自由度df組間所得到的結果,近似服從自由度為df組間的χ2分布。
2.3.4 與備擇假設為“一般關聯性”對應的檢驗統計量
若擬分析的資料如本文表1(即雙向無序R×C列聯表資料)所示,此時式(1)就簡化成下式:


【例4】沿用例3的資料,試將其分別視為“雙向有序R×C列聯表資料”“結果變量為有序變量R×C列聯表資料”和“雙向無序R×C列聯表資料”,同時對其進行CMHχ2檢驗。基于“表評分”算法所需要的SAS程序如下:


以上是基于“表評分”算法所得CMH χ2檢驗的三種計算結果,分別對應三種不同的備擇假設:第一種,“非零相關(即行、列兩有序變量之間具有非零的相關關系)”,=6.635,df=1,P=0.0100;第二種,“行評分均值不同(即各行上有序結果變量的評分值的平均數不等)”,=33.7083,d f=4,P<0.0001;第三種,“常規關聯(即行、列兩無序變量之間存在關聯性)”,=36.6435,df=8,P<0.0001。
以上三種檢驗都得出P<0.01,故接受備擇假設。對本例而言,可以認為:①“媒體暴露情況(時間)”與“ASD程度”之間存在“非零相關”,具體地說,隨著媒體暴露時間延長,出現重度ASD癥狀的比例呈非線性上升趨勢;②五種不同媒體暴露時間所對應的平均ASD值是不相等的,除第2個時間段之外,其他4個時間段都隨著“暴露時間延長”,平均ASD值變大(即ASD癥狀程度逐漸加重);③媒體暴露情況(時間)與ASD程度之間存在關聯性(前面的“①”和“②”可以清楚體現出兩變量之間是如何關聯的)。
在前面的SAS程序中,在“tables語句”中出現了選項“SCORES=table”,其含義是用表格里規定的方式給有序變量評分。這里的“表格”,其實是指SAS程序中如何給“行變量”或“列變量”賦值。“DO A=1 TO 5;”就是讓“行變量 A”依次取 1、2、3、4、5分;同理,“DO B=1 TO 3;”就是讓“列變量B”依次取1、2、3分。當然,也可以給變量A或B賦任意的幾個由小到大的數值,例如“DO A=0.1,1.2,13.5,21.8,34.6;”。
在SAS/STAT的FREQ過程中,還有另外三種評分方法,分別為“rank評分法”“ridit評分法”和“修正的ridit評分法”,對應的選項依次為“SCORES=rank”“SCORES=ridit”和“SCORES=modridit”。
對“常規關聯”的計算結果而言,四種評分方法所得計算結果是完全相同的,因為計算公式中不涉及“評分”;對“非零相關”的計算結果而言,“表評分法”(采用Pearson’s相關分析)與另外三種評分法(其結果是一樣的,采用Spearman’s秩相關分析)計算結果可能相差很大、甚至結論相反;對“行評分均值不同”的計算結果而言,“表評分法”(采用單因素R水平設計一元定量資料方差分析)與另外三種評分法(其結果是一樣的,采用單因素R水平設計一元定量資料Kruskal-Wallis’s H秩和檢驗)計算結果稍有差別。例如,例4“基于秩評分”的結果如下:

【說明】“ridit評分法”和“修正的ridit評分法”輸出結果與“基于秩評分”方法輸出結果相同,此處從略。
注意:“非零相關”的結果與前面“基于表評分”的結果相反;四種評分法對應的“常規關聯”結果都是χ2=36.6435、P<0.0001;四種評分法對應的“行評分均值不同”有兩個不同的計算結果,分別為“χ2=33.7083[采用方差分析計算再基于式(6)轉換的結果]”與“χ2=31.2716(采用 Kruskal-Wallis’s H秩和檢驗計算的結果)”,但P值都小于0.0001。
在給“scores=”選項指定“table”或其他三種非參數評分法(即rank、ridit和modridit)時,為什么“非零相關”有兩個完全相反的結果?原因在于:前者采用的是“Pearson’s相關分析”,而后者采用的是“Spearman’s秩相關分析”。就例3資料而言,分別基于“Pearson’s相關分析”與“Spearman’s秩相關分析”得到“r=0.02033”與“rs=0.00568”。依據式(2)可算出對應的CMH χ2值分別為6.632368(與前面相應的計算結果6.6350很接近)和0.517715(與前面相應的計算結果0.5174很接近)。
在“常規關聯”的分析中,將有序變量視為無序變量或名義變量,故給“scores=”選項指定四種內容(即table、rank、ridit和modridit)時,所得結果完全相同。就例3資料而言,。它在本質上就是=36.6458(即Pearson’s擬合優度χ2值),由式(7)可知:

在備擇假設為“行評分均值不同”時,SAS/STAT中FREQ過程將此時的資料視為“單因素R水平設計一元定量資料”。若采用“表評分(即scores=table)”,則SAS/STAT中FREQ過程采用“單因素方差分析”;若采用其他三種評分,SAS/STAT中FREQ過程采用Kruskal-Wallis’s秩和檢驗。事實上,就是把結果變量B視為“定量變量”,因此,對結果變量B的賦值方法不同,將直接影響R×C列聯表資料各行上結果變量平均值的計算結果。四種評分方法對應結果變量B的3個評分值如下:

若直接采用SAS/STAT中ANOVA過程計算,得到F=8.44、df=4,基于式(6)可算得≈ 8.44 × 4=33.76,與“scores=table”對應的“行評分均值不同”時的“=33.7083”比較接近;若直接采用SAS/STAT中NPAR1WAY過程計算,得到H=31.2716,與“scores=rank”對應的“行評分均值不同”時的“=31.2716”完全相同。
值得一提的是,在分析二維列聯表資料時,三種非參數評分法輸出的三種檢驗結果相同;但若是分析高維列聯表資料,情況就不盡相同了。因篇幅所限,此處從略。
本文呈現了三種R×C列聯表資料的實例及其三種相應的統計分析方法,這些方法可概括成廣義CMH χ2檢驗;基于SAS/STAT中FREQ過程對R×C列聯表資料進行分析時,可同時輸出三種統計分析結果;本文還詳細揭示了選取不同的“評分”方法,將會影響“非零相關”與“行評分均值不同”計算結果的真實原因。