胡純嚴(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)行列聯(lián)表資料的獨(dú)立性分析和/或廣義線性回歸模型[1]分析中,常會(huì)涉及“擬合優(yōu)度檢驗(yàn)”和“失擬檢驗(yàn)”等內(nèi)容[1-4]。本文將介紹有關(guān)概念、各種檢驗(yàn)方法的計(jì)算公式以及基于SAS軟件[5]實(shí)現(xiàn)計(jì)算的方法。
所謂擬合優(yōu)度檢驗(yàn),就是評(píng)價(jià)基于某種假設(shè)或回歸模型推算出按順序排列的一組期望頻數(shù)與相應(yīng)觀察頻數(shù)之間一致性的一種檢驗(yàn)方法[1,5]。擬合優(yōu)度檢驗(yàn)常涉及以下兩個(gè)統(tǒng)計(jì)學(xué)問題:其一,檢驗(yàn)?zāi)撤N假設(shè)是否成立。例如:在列聯(lián)表資料分析中可能會(huì)提出“某些屬性變量之間存在獨(dú)立性”的假設(shè),此時(shí)的檢驗(yàn)被稱為“獨(dú)立性擬合優(yōu)度檢驗(yàn)”;又例如:在展示一組單變量資料的頻數(shù)分布時(shí),可能會(huì)提出“該組資料服從正態(tài)分布、Poisson分布或二項(xiàng)分布”的假設(shè),此時(shí)的檢驗(yàn)被稱為“分布擬合優(yōu)度檢驗(yàn)”。其二,檢驗(yàn)所擬合的某個(gè)回歸模型對(duì)所描述資料的變化趨勢(shì)是否吻合,此時(shí)的檢驗(yàn)被稱為“回歸模型擬合優(yōu)度檢驗(yàn)”。
從功能上來劃分,擬合優(yōu)度檢驗(yàn)可分為四種:①獨(dú)立性擬合優(yōu)度檢驗(yàn)(常用于列聯(lián)表資料的獨(dú)立性檢驗(yàn));②線性擬合優(yōu)度檢驗(yàn)(常用于定量因變量的線性回歸分析中,以判斷是否需要采取更復(fù)雜的方法或非線性回歸模型);③分布擬合優(yōu)度檢驗(yàn)(以判斷給定的資料是否服從某種特定的概率分布,如正態(tài)分布、Poisson分布或二項(xiàng)分布);④回歸模型擬合優(yōu)度檢驗(yàn)(以判斷當(dāng)前所構(gòu)建的回歸模型與給定資料的吻合程度)”。從檢驗(yàn)統(tǒng)計(jì)量所服從的分布來劃分,擬合優(yōu)度檢驗(yàn)可分為三種:①F分布(常用于線性回歸模型的分析之中);②經(jīng)驗(yàn)分布(即柯爾莫哥洛夫檢驗(yàn));③χ2分布(常用于除線性回歸模型分析之外的其他場(chǎng)合)。從構(gòu)造檢驗(yàn)統(tǒng)計(jì)量的方法來劃分,擬合優(yōu)度檢驗(yàn)可分為三種:①方差分析法(常用于線性回歸模型的分析);②似然比法(常用于除線性回歸模型分析之外的其他場(chǎng)合);③觀測(cè)頻數(shù)與理論頻數(shù)之偏差的平方和法(常用于除線性回歸模型分析之外的其他場(chǎng)合)。這三種劃分方法在內(nèi)容上存在交叉和重疊。因篇幅所限,本文僅介紹基于“χ2檢驗(yàn)統(tǒng)計(jì)量”的擬合優(yōu)度檢驗(yàn)方法[5-9]。
所謂失擬檢驗(yàn),就是對(duì)回歸函數(shù)是否為自變量的線性函數(shù)進(jìn)行的檢驗(yàn)[2]。失擬檢驗(yàn)最早出現(xiàn)在簡(jiǎn)單線性回歸模型的分析過程中,用以檢驗(yàn)采用簡(jiǎn)單線性回歸模型是否可以取代曲線回歸模型。若失擬檢驗(yàn)的結(jié)果為P>0.05,表明用簡(jiǎn)單線性回歸模型可以取代曲線回歸模型;反之亦然。在SAS/STAT的REG過程中,在進(jìn)行多重線性回歸分析時(shí),也可以進(jìn)行失擬檢驗(yàn)[5]。值得注意的是:在線性回歸分析中所做的失擬檢驗(yàn)是F檢驗(yàn)(也叫做方差分析),而不是χ2檢驗(yàn)。在SAS/STAT的LOGISTIC過程中,將“Hosmer-Lemeshow檢驗(yàn)”列入“擬合優(yōu)度檢驗(yàn)(Goodness of fit tests)”標(biāo)題之下,但在指定“選項(xiàng)(lackfit=)”時(shí),卻又將其歸類于“失擬檢驗(yàn)”之中;然而,在SAS/STAT的PROBIT過程中,將“失擬檢驗(yàn)”與“擬合優(yōu)度檢驗(yàn)”視為同一種檢驗(yàn)[5]。正因如此,在下文中,統(tǒng)一稱為“擬合優(yōu)度檢驗(yàn)”,而不特別區(qū)分何時(shí)為“擬合優(yōu)度檢驗(yàn)”,何時(shí)為“失擬檢驗(yàn)”。
Pearson’s擬合優(yōu)度檢驗(yàn)可以應(yīng)用于下列三種場(chǎng)合[1-5]:①列聯(lián)表資料獨(dú)立性假定的擬合優(yōu)度檢驗(yàn);②單變量頻數(shù)分布的分布擬合優(yōu)度檢驗(yàn);③除線性回歸模型分析之外的其他回歸模型的擬合優(yōu)度檢驗(yàn)。
2.2.1 獨(dú)立性擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量
在進(jìn)行R×C列聯(lián)表(包括2×2列聯(lián)表)資料的獨(dú)立性檢驗(yàn)時(shí),Pearson’s擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量的形式見式(1):
在式(1)中,χ2服從自由度為df=(R-1)(C-1)的χ2分布;R和C分別代表二維列聯(lián)表的“行數(shù)”與“列數(shù)”。
2.2.2 分布擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量
在推斷一組單變量定量資料服從何種頻數(shù)分布時(shí),Pearson’s擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量的形式見式(2):
在式(2)中,χ2服從自由度為df=k-r-1的χ2分布;k為資料中不同數(shù)據(jù)個(gè)數(shù)或分組數(shù);r為假定概率分布函數(shù)中待估參數(shù)的個(gè)數(shù);fi為定量變量的第i觀測(cè)值(需要由小到大排序,在實(shí)際使用中,應(yīng)是第i個(gè)組的組中值)所對(duì)應(yīng)的觀測(cè)頻數(shù);Pi為與fi對(duì)應(yīng)的頻率;n為總樣本含量;nPi為基于設(shè)定概率分布推算出來的理論頻數(shù)。
值得注意的是:各nPi應(yīng)大于等于5,否則,需要將相鄰的數(shù)據(jù)或組進(jìn)行合并,以滿足前述的前提條件。
2.2.3 回歸模型擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量
2.2.3.1 在SAS/STAT的GENMOD過程中的定義
在SAS/STAT的GENMOD過程中,Pearson’s擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量的定義見式(3):
在式(3)中,χ2服從自由度為df=k-r的χ2分布;k為資料中不同數(shù)據(jù)個(gè)數(shù)或分組數(shù);fi、wi、yi、μi和V(μi)分別代表第i個(gè)組的“頻數(shù)”“權(quán)重”“因變量的取值”“yi的期望值”和“yi的方差”;r為所構(gòu)建的回歸模型中待估參數(shù)的個(gè)數(shù)(包括截距項(xiàng))。值得注意的是:廣義線性模型涵蓋面很寬泛,其因變量通常涉及以下常見概率分布,例如:正態(tài)分布、Poisson分布、伽瑪分布、逆高斯分布、二項(xiàng)分布、多項(xiàng)分布、負(fù)二項(xiàng)分布、零膨脹Poisson分布和零膨脹負(fù)二項(xiàng)分布等;依據(jù)不同的概率分布,yi、μi和V(μi)將取不同的值或計(jì)算公式,限于篇幅,不再贅述。
2.2.3.2 在SAS/STAT的LOGISTIC過程中的定義
在SAS/STAT的LOGISTIC過程中,Pearson’s擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量的定義見式(4):
2.2.3.3 在SAS/STAT的PROBIT過程中的定義
在SAS/STAT的PROBIT過程中,Pearson’s擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量的定義見式(5):
偏差或稱似然比擬合優(yōu)度檢驗(yàn)的應(yīng)用場(chǎng)合較少,目前僅應(yīng)用于回歸模型擬合優(yōu)度檢驗(yàn)和檢驗(yàn)整個(gè)回歸模型中全部參數(shù)是否都為0。
3.2.1 在SAS/STAT的GENMOD過程中的定義
在SAS/STAT的GENMOD過程中[5],偏差擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量的定義見式(6):
在式(6)中,Y與μ分別代表因變量的向量與因變量預(yù)測(cè)均值向量;l(Y,Y)和l(Y,μ)分別代表含有全部自變量的回歸模型的似然函數(shù)的自然對(duì)數(shù)與含有部分自變量的回歸模型的似然函數(shù)的自然對(duì)數(shù);D(Y,μ)服從自由度為df=n-q的χ2分布,這里,n和q分別代表總樣本含量與所構(gòu)建的回歸模型中待估參數(shù)的個(gè)數(shù)(包括截距項(xiàng))。
3.2.2 在SAS/STAT的LOGISTIC過程中的定義
在SAS/STAT的LOGISTIC 過程中[5],偏差擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量的定義見式(7):
3.2.3 在SAS/STAT的PROBIT過程中的定義
在SAS/STAT的PROBIT過程中[5],偏差擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量的定義見式(8):
前面介紹的Pearson’s擬合優(yōu)度檢驗(yàn)與偏差擬合優(yōu)度檢驗(yàn)都是基于數(shù)據(jù)分組之上的,當(dāng)資料中存在連續(xù)性變量時(shí),就不可避免地會(huì)出現(xiàn)某些組中的樣本含量很少,甚至只有1~2個(gè)樣本,這種現(xiàn)象被稱為“稀疏資料”。事實(shí)上,前面提及的兩種擬合優(yōu)度檢驗(yàn)要求各組樣本含量應(yīng)大于等于5,否則,其檢驗(yàn)結(jié)果不可靠。針對(duì)“稀疏資料(某些組中沒有重復(fù)數(shù)據(jù))”,Hosmer和Lemeshow于2000年提出了一種新的擬合優(yōu)度檢驗(yàn)方法,該擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量近似服從χ2分布[5]。
二值因變量條件下Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量見式(9):
多值因變量條件下Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)統(tǒng)計(jì)量[5]見式(10):
盡管Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)對(duì)稀疏資料具有一定的校正能力,但它仍然需要將資料分組后才能進(jìn)行計(jì)算。當(dāng)資料中“稀疏情形”非常嚴(yán)重時(shí),其檢驗(yàn)效能將會(huì)明顯下降。為此,統(tǒng)計(jì)學(xué)家提出了另外6種稀疏資料擬合優(yōu)度檢驗(yàn)方法。迄今為止,這些方法僅適用于因變量為二值變量的logistic回歸模型之中[5]。設(shè)j代表事件出現(xiàn)的預(yù)測(cè)概率,代表擬合模型的協(xié)方差矩陣。
5.2.1 信息矩陣檢驗(yàn)
Orme于1988年將White于1982年提出的“一般錯(cuò)誤指定檢驗(yàn)”應(yīng)用于二項(xiàng)反應(yīng)資料。設(shè)Xj為第j個(gè)觀測(cè)的設(shè)計(jì)向量(即由全部協(xié)變量形成的向量),于是,XjX'j就是信息矩陣。通過信息矩陣XjX'j的右上角矩陣,將設(shè)計(jì)向量Xj展開成為一個(gè)新向量,記為Vech(Xj)。通過將Xj尺度化,通過將Vech(Xj)尺度化。使用殘差ej=Yj-j的二項(xiàng)形式創(chuàng)建新反應(yīng)變量的值Y*,其表達(dá)式見式(11):
構(gòu)建Y*關(guān)于展開后的協(xié)變量的多重線性回歸模型,求出該模型的回歸平方和,這個(gè)平方和服從自由度為df的χ2分布。這里,df等于尺度化后的Vech(Xj)中未退化的協(xié)變量的數(shù)目?;诖朔椒óa(chǎn)生的檢驗(yàn)結(jié)果在擬合優(yōu)度檢驗(yàn)表中的標(biāo)簽為“信息矩陣”。
5.2.2 信息矩陣對(duì)角線檢驗(yàn)
Kuss于2002年修改了上述“信息矩陣檢驗(yàn)”方法,即僅通過信息矩陣XjX'j對(duì)角線上的元素,將設(shè)計(jì)向量Xj展開成為一個(gè)新向量,記為Vech(Xj)?;诖朔椒óa(chǎn)生的檢驗(yàn)結(jié)果在擬合優(yōu)度檢驗(yàn)表中的標(biāo)簽為“信息矩陣對(duì)角”。
5.2.3 Osius-Rojek檢驗(yàn)
Osius和Rojek于1992年用固定格漸近方法推導(dǎo)出 Pearson’sχ2統(tǒng)計(jì)量的均值和方差,均值為m,即子總體的個(gè)數(shù)(即“組數(shù)”);方差的計(jì)算方法見式(12):
基于式(13)可得到Osius-Rojek檢驗(yàn)統(tǒng)計(jì)量見式(14):
在式(14)中,B的定義見式(13);為服從自由度為1的χ2分布的檢驗(yàn)統(tǒng)計(jì)量。
5.2.4 未加權(quán)的殘差平方和檢驗(yàn)
Hosmer等于1997年在Copas于1989年提出的基本公式的基礎(chǔ)上,構(gòu)造出未加權(quán)的殘差平方和檢驗(yàn)統(tǒng)計(jì)量,見式(15):
在式(15)中,σ為回歸的殘差標(biāo)準(zhǔn)差,其他符號(hào)的含義分別如下:
5.2.5 Spiegelhalter檢驗(yàn)
Spiegelhalter于1986年基于Brier評(píng)分推導(dǎo)出一個(gè)檢驗(yàn)統(tǒng)計(jì)量,見式(18):
在式(18)中,各符號(hào)的含義如下:
5.2.6 Stukel檢驗(yàn)
Stukel于1988年提出了如下方法,在現(xiàn)有模型中增加兩個(gè)新變量,然后檢驗(yàn)這兩個(gè)新變量是否為無統(tǒng)計(jì)學(xué)意義的變量。對(duì)于二值反應(yīng)變量而言,設(shè)ηj=X'jβ,I(·)代表指示變量,增加的兩個(gè)新變量如下:
然后,應(yīng)用“評(píng)分檢驗(yàn)(說明:此檢驗(yàn)常被用于檢驗(yàn)?zāi)P椭胁糠只貧w系數(shù)是否為0)”檢驗(yàn)這個(gè)大模型(指增加兩個(gè)新變量后的模型)與已擬合的模型(即未增加新變量的模型)之間差異是否有統(tǒng)計(jì)學(xué)意義。若檢驗(yàn)的結(jié)果為“P>0.05”,說明原先已擬合的模型對(duì)資料的擬合效果滿意。
【例1】為探討新生兒乙肝病毒宮內(nèi)感染的危險(xiǎn)因素,某醫(yī)科大學(xué)流行病學(xué)教研室于2004年-2005年收集了145例HBsAg陽性孕婦及其145例新生兒作為研究對(duì)象,檢測(cè)新生兒出生24小時(shí)內(nèi)血清中HBV DNA,以判斷是否發(fā)生了乙肝病毒宮內(nèi)感染,新生兒血清HBV DNA陽性者為發(fā)生宮內(nèi)感染,陰性者為未發(fā)生宮內(nèi)感染。同時(shí)檢測(cè)孕婦血清中HBVDNA(陰性=0,陽性=1)、HbeAg(陰性=0,陽性=1),所得資料見表1[10],試分析新生兒乙肝宮內(nèi)感染的影響因素。
表1 145例孕婦及其新生兒血清HBV相關(guān)指標(biāo)檢測(cè)結(jié)果
【分析與解答】設(shè)所需要的SAS程序如下:
以上輸出的是“擬合優(yōu)度檢驗(yàn)”結(jié)果,前兩行是最常用的方法,后六行是稀疏資料擬合優(yōu)度檢驗(yàn)結(jié)果,P均>0.05,說明所構(gòu)建的回歸模型對(duì)資料的擬合效果好。
以上輸出的是“模型擬合統(tǒng)計(jì)量”,它們只能相對(duì)反映模型對(duì)資料的擬合效果,因?yàn)樗鼈儾皇羌僭O(shè)檢驗(yàn)的結(jié)果,只有在有兩個(gè)或兩個(gè)以上模型時(shí),才可根據(jù)模型擬合統(tǒng)計(jì)量的數(shù)值,判斷哪個(gè)模型對(duì)資料的擬合效果相對(duì)更好(例如,各種“R2”值越大越好,AIC、AICC和SC的值越小越好)。
以上輸出的是“Hosmer和Lemeshow擬合優(yōu)度檢驗(yàn)”結(jié)果,因P=0.9085>0.05,說明所構(gòu)建的回歸模型對(duì)資料的擬合效果較好,不存在“失擬”問題。
以上輸出的是“檢驗(yàn)全局原假設(shè):BETA=0”的結(jié)果,三種檢驗(yàn)方法所得結(jié)果和結(jié)論基本一致,即認(rèn)為回歸模型中的回歸系數(shù)不全為0,所構(gòu)建的回歸模型具有一定的實(shí)用價(jià)值。
以上輸出的是“優(yōu)比估計(jì)”結(jié)果,X2的95%Wald置信區(qū)間包括“1”,表明X2的兩個(gè)水平對(duì)結(jié)果的影響接近相等;X1的,其95%Wald置信區(qū)間為[1.694,31.393],未包括“1”,表明X1的兩個(gè)水平之間對(duì)結(jié)果的影響是不同的,即“孕婦血中HBV DNA(X1)”為陽性條件下“新生兒出現(xiàn)乙肝病毒宮內(nèi)感染”的概率是“孕婦血中HBV DNA(X1)”為陰性條件下“新生兒出現(xiàn)乙肝病毒宮內(nèi)感染”概率的7.293倍。
以上輸出的是“用于Hosmer和Lemeshow檢驗(yàn)的分區(qū)”結(jié)果,實(shí)際上就是將全部資料劃分成4個(gè)組。
以上輸出的是“最大似然估計(jì)分析”結(jié)果,此結(jié)果表明:X2的回歸系數(shù)與0之間的差別無統(tǒng)計(jì)學(xué)意義,表明“孕婦血中HbeAg(X2)”是否為“陽性”對(duì)“新生兒是否出現(xiàn)乙肝病毒宮內(nèi)感染”的影響可以忽略不計(jì)。要想獲得更簡(jiǎn)潔的回歸模型,可以淘汰掉自變量X2后重新擬合logistic回歸模型。
【統(tǒng)計(jì)結(jié)論與專業(yè)結(jié)論】由多種擬合優(yōu)度檢驗(yàn)結(jié)果可知,所擬合的回歸模型對(duì)資料的擬合效果滿意;由回歸系數(shù)的檢驗(yàn)結(jié)果可知,孕婦血中HBV DNA(X1)是新生兒乙肝宮內(nèi)感染的主要影響因素。
擬合優(yōu)度檢驗(yàn)與失擬檢驗(yàn)在本質(zhì)上是相同的,都是用來評(píng)價(jià)某種假設(shè)(例如獨(dú)立性假設(shè)或關(guān)于資料來自某種分布的假設(shè))是否成立或所擬合的回歸模型與資料的吻合程度。其微小差別在于:擬合優(yōu)度檢驗(yàn)一般要求將資料整理成一維頻數(shù)表或高維頻數(shù)表資料形式,當(dāng)檢驗(yàn)結(jié)果為P>0.05時(shí),表明所做的假設(shè)成立或所擬合的回歸模型與資料的變化趨勢(shì)是吻合的;當(dāng)檢驗(yàn)結(jié)果為P≤0.05時(shí),表明所做的假設(shè)可能不成立或所擬合的回歸模型與資料的變化趨勢(shì)不夠吻合;然而,失擬檢驗(yàn)一般要求在自變量或自變量向量取特定值的條件下應(yīng)做多次重復(fù)試驗(yàn),采用失擬檢驗(yàn)的目的是揭示資料中所包含的“曲線變化的成分”是否足夠多。當(dāng)檢驗(yàn)結(jié)果為P>0.05時(shí),表明所擬合的回歸模型與資料的變化趨勢(shì)是吻合的;當(dāng)檢驗(yàn)結(jié)果為P≤0.05時(shí),表明所擬合的回歸模型與資料的變化趨勢(shì)不夠吻合,分析者需要構(gòu)建更復(fù)雜的回歸模型或非線性回歸模型。
本文介紹了9種擬合優(yōu)度檢驗(yàn)方法,它們的主要應(yīng)用場(chǎng)合是評(píng)價(jià)廣義線性回歸模型(包括logistic回歸模型)對(duì)資料的擬合效果。其中,偏差擬合優(yōu)度檢驗(yàn)和Pearson’s擬合優(yōu)度檢驗(yàn)最常用;Hosmer-Lemeshow擬合優(yōu)度檢驗(yàn)對(duì)稀疏資料具有一定的校正作用;當(dāng)資料的稀疏程度十分嚴(yán)重時(shí),需要使用專門用于稀疏資料的6種擬合優(yōu)度檢驗(yàn)。本文結(jié)合一個(gè)實(shí)例并借助SAS軟件實(shí)現(xiàn)了全部9種擬合優(yōu)度檢驗(yàn),對(duì)SAS輸出結(jié)果做了解釋,給出了統(tǒng)計(jì)結(jié)論和專業(yè)結(jié)論。