軍事醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)統(tǒng)計(jì)學(xué)咨詢中心(100850) 周詩國 胡良平
若要采用Newman-Keuls檢驗(yàn)對單因素多水平設(shè)計(jì)一元定量資料進(jìn)行兩兩比較,則需要用到q臨界值;若要進(jìn)行多個總體均數(shù)比較時樣本含量或檢驗(yàn)效能的估計(jì),則通常需要獲得相應(yīng)條件下的ψ值;若要進(jìn)行多個總體率的比較時樣本含量或檢驗(yàn)效能的估計(jì),則通常需要獲得相應(yīng)條件下的λ值〔1-3〕。絕大多數(shù)統(tǒng)計(jì)學(xué)教科書或?qū)V际珍浟擞嘘P(guān)這些參數(shù)在特定條件下的數(shù)值表供讀者查閱〔4-11〕,但明確闡述它們的含義及在特定條件下如何計(jì)算出它們的數(shù)值,卻很難尋覓相應(yīng)的文獻(xiàn)資料。本文的目的就在于介紹q臨界值、ψ值和λ值的含義及其計(jì)算方法。
1.q臨界值的含義
當(dāng)單因素多水平設(shè)計(jì)一元定量資料的方差分析結(jié)果為拒絕H0,接受H1,得出“各總體均數(shù)不同或不全相同”的結(jié)論時,并不能說明各總體均數(shù)兩兩之間是否不同,為此,可在前述方差分析的基礎(chǔ)上,對均數(shù)進(jìn)一步作兩兩比較,也稱多重比較(multiple comparison)。均數(shù)間兩兩比較的方法有多種,Newman-Keuls檢驗(yàn)只是其中之一。Newman-Keuls檢驗(yàn)亦稱Student-Newman-Keuls(SNK)檢驗(yàn),簡稱q檢驗(yàn)。其檢驗(yàn)統(tǒng)計(jì)量q的計(jì)算公式為
式中,ˉXi、ˉXj分別為兩對比組的樣本均數(shù);SˉXi-ˉXj為兩對比組樣本均數(shù)差值的標(biāo)準(zhǔn)誤;MSe為方差分析的誤差均方;ni、nj分別為兩對比組的樣本含量。統(tǒng)計(jì)量q又被稱為學(xué)生化極差(studentized range)統(tǒng)計(jì)量,檢驗(yàn)拒絕域?yàn)?/p>
q≥q(a,v,α) (2)式中,q(a,v,α)為統(tǒng)計(jì)量q的臨界值,即統(tǒng)計(jì)量q所服從的分布(不妨記為q分布)的第100(1-α)百分位數(shù),α 的取值為小數(shù),且有 q≥q(a,v,α)的概率為 α,即P[q≥q(a,v,α)]=α。一般取 α =0.05 或 α =0.01。q(a,v,α)是所比較的兩個經(jīng)排序后的均數(shù)之間的跨度a(即每次比較所涉及到的處理組的組數(shù),其最小值為2,最大值為試驗(yàn)因素的水平數(shù)k)、方差分析中誤差項(xiàng)的自由度v及完全無效假設(shè)下的試驗(yàn)誤差率(EERC)α的函數(shù)。
2.q臨界值的計(jì)算方法
將正態(tài)密度函數(shù)看作一個簡單函數(shù),則q分布的概率密度函數(shù)是一個二重積分的形式,其形狀(或圖形)隨著k、v兩個參數(shù)的變化而劇烈變化。
要求臨界值q(a,v,α),則需求解概率方程P[q≥值就是直接通過數(shù)值積分法計(jì)算出來的。
目前的SAS軟件已經(jīng)提供了一個可以用來計(jì)算q臨界值的 SAS 函數(shù),即 PROBMC 函數(shù)〔12,13〕。其語法規(guī)則如下:
PROBMC(distribution,q,prob,υ,nparms,<parameters>)
式中各個參數(shù)的含義如下:
distribution:一個標(biāo)識分布類型的字符串。對PROBMC函數(shù)有效的分布共有5種,見表1。
表1 對PROBMC函數(shù)有效的5種分布及其對應(yīng)的參數(shù)值
q:當(dāng)q等于某個具體的數(shù)值時(此時,參數(shù)prob的取值必須用·代替),其表示服從參數(shù)distribution所代表的某種分布的統(tǒng)計(jì)量的值;當(dāng)q的取值用·代替時(此時,參數(shù)prob必須等于一個位于(0,1)之間的具體的數(shù)值),表示將通過PROBMC函數(shù)求取參數(shù)distribution所代表的某種分布曲線下左側(cè)概率為prob的分位數(shù)。
prob:相應(yīng)分布的概率密度曲線下隨機(jī)變量的取值位于特定分位數(shù)左側(cè)的概率。
q與prob必須有且只有一個的取值用·來代替,而PROBMC函數(shù)所返回的值正是q和prob兩個參數(shù)中取值用·來代替的那個參數(shù)的值。
υ:方差分析時誤差項(xiàng)的自由度,其計(jì)算公式為υ=(N-1)-(k-1)=N-k。式中的k表示試驗(yàn)因素的水平數(shù),N表示總樣本含量,即k個水平組的樣本含量之和。若υ的取值用·來代替,則表示自由度為無窮大(∞)。
nparms:處理組的組數(shù)。對DUNNETT1和DUNNETT2所對應(yīng)的分布,組數(shù)不包括對照組。比如有一項(xiàng)單因素5水平設(shè)計(jì),若以樣本均值最小的那個組作為對照組,采用DUNNETT法比較其他組與對照組的總體均值差異有無統(tǒng)計(jì)學(xué)意義(此時,參數(shù)distribution的取值為‘DUNNETT1'或‘DUNNETT2'),在比較樣本均值最大的那個組與對照組之間的總體均值的差異有無統(tǒng)計(jì)學(xué)意義時,nparms的取值為4;若采用SNK法對5個總體均值進(jìn)行兩兩比較,在比較樣本均值最大的那個組與樣本均值最小的那個組之間的總體均值的差異有無統(tǒng)計(jì)學(xué)意義時,nparms的取值為5。
parameters:是一個可選項(xiàng)。該選項(xiàng)為一個序列,組成了“nparms”這個參數(shù)的具體內(nèi)容。當(dāng)試驗(yàn)因素各水平組的樣本含量不等時,必須指定該選項(xiàng)。nparms的含義取決于分布的類型。若不指定這些參數(shù),則意味著假定各組的樣本含量相等。
對學(xué)生化極差,當(dāng)自由度υ≠∞且各組的樣本含量不等時,函數(shù)PROBMC無效。
對Williams檢驗(yàn),當(dāng)各組的樣本含量不等時,函數(shù)PROBMC也無效。
【例1】 設(shè)有某項(xiàng)單因素5水平設(shè)計(jì),每個水平組做4次獨(dú)立重復(fù)試驗(yàn),測量某項(xiàng)定量指標(biāo)的取值,并假定資料滿足方差分析的前提條件,取檢驗(yàn)水準(zhǔn)α=0.05,經(jīng)單因素5水平設(shè)計(jì)一元定量資料的方差分析處理,發(fā)現(xiàn)試驗(yàn)因素對試驗(yàn)結(jié)果的影響有統(tǒng)計(jì)學(xué)意義,現(xiàn)欲采用SNK檢驗(yàn)進(jìn)行5個總體均值之間的多重比較,請用SAS計(jì)算所需的q臨界值。
【分析與解答】 由題意可知,檢驗(yàn)水準(zhǔn)α=0.05,試驗(yàn)因素的水平數(shù)k=5,各水平組的樣本含量n=4,故多重比較時的自由度υ=N-k=kn-k=k(n-1)=15。所需求取的 q 臨界值為 q(5,15,0.05)、q(4,15,0.05)、q(3,15,0.05)和 q(2,15,0.05),共4 個??捎孟旅娴腟AS程序計(jì)算出上述q臨界值。
data c_q_v;
alpha=0.05;
df=15;/*用df表示自由度υ*/
do a=5 to 2 by-1;
q=probmc("range",.,1-alpha,df,a);
output;
end;
run;
ods html;
proc print data=c_q_v noobs;run;
ods html close;
quit;
計(jì)算結(jié)果如表2所示。
表2 例1所需的q臨界值
alpha df a q 0.05 15 5 4.366 99 0.05 15 4 4.075 97 0.05 15 3 3.673 38 0.05 15 2 3.014 32
1.ψ值的含義
從ψ值表中可以看出,ψ值是4個參數(shù)的函數(shù),這4個參數(shù)分別為進(jìn)行單因素多水平設(shè)計(jì)一元定量資料方差分析時的Ⅰ型錯誤概率α、Ⅱ型錯誤概率β以及檢驗(yàn)統(tǒng)計(jì)量F的分子和分母的自由度v1和v2。
此外,方差分析是以F分布為理論依據(jù)的。當(dāng)H0成立時,方差分析的檢驗(yàn)統(tǒng)計(jì)量F服從分子和分母的自由度分別為v1和v2的中心F分布。當(dāng)H0不成立時,方差分析的檢驗(yàn)統(tǒng)計(jì)量F服從分子和分母的自由度分別為v1和v2、非中心參數(shù)δ≠0的非中心F分布。中心F分布只是非中心F分布的一個特例(非中心參數(shù)δ=0)。根據(jù)非中心F分布的非中心參數(shù)值的計(jì)算方法及相應(yīng)的SAS函數(shù),筆者對非中心F分布的非中心參數(shù)值進(jìn)行了探索性的計(jì)算,并將計(jì)算結(jié)果與統(tǒng)計(jì)學(xué)教材或?qū)V薪o出的相應(yīng)條件下的ψ值進(jìn)行比較,結(jié)果表明:ψ值并不是非中心F分布的非中心參數(shù)值,而是使方差分析同時滿足下面兩個條件時的非中心F分布的非中心參數(shù)值δ除以試驗(yàn)因素的自由度v1(即檢驗(yàn)統(tǒng)計(jì)量F的分子的自由度)后開算術(shù)平方根的結(jié)果,即ψ=。兩個條件為:
(1)拒絕H0時可能犯錯誤的最大概率為α;
(2)不能拒絕H0時,若接受H0,可能犯錯誤的最大概率為β。
2.ψ值的計(jì)算方法
(1)確定單因素多水平設(shè)計(jì)一元定量資料方差分析的Ⅰ型錯誤概率α、Ⅱ型錯誤概率β以及檢驗(yàn)統(tǒng)計(jì)量F的分子和分母的自由度v1和v2;
(2)用F分布的分位數(shù)函數(shù)FINV求出與分子和分母的自由度分別為ndf和ddf、分布曲線下左側(cè)概率為1-α的中心F分布對應(yīng)的分位數(shù)FINV(1-α,ndf,ddf);
(3)用非中心F分布的非中心參數(shù)函數(shù)FNONCT求出與分子和分母的自由度分別為ndf和ddf、分位數(shù)為FINV(1-α,ndf,ddf)、分布曲線下左側(cè)概率為β的非中心F分布的非中心參數(shù)值FNONCT(FINV(1-α,ndf,ddf),ndf,ddf,β);
此結(jié)果與從統(tǒng)計(jì)學(xué)教材或?qū)V薪o出的ψ值表中查到的結(jié)果相吻合。
1.λ值的含義
從λ值表中可以看出,λ值是3個參數(shù)的函數(shù),這3個參數(shù)分別為進(jìn)行單因素多水平設(shè)計(jì)定性資料(結(jié)果為二值變量)χ2檢驗(yàn)時的Ⅰ型錯誤概率α、Ⅱ型錯誤概率β以及自由度v=k-1(k為試驗(yàn)因素的水平數(shù))。
此外,χ2檢驗(yàn)是以χ2分布為理論依據(jù)的。當(dāng)H0成立時,χ2檢驗(yàn)的檢驗(yàn)統(tǒng)計(jì)量χ2服從自由度為v=(k-1)(2-1)=k-1的中心χ2分布(結(jié)果變量為二值變量)。當(dāng)H0不成立時,χ2檢驗(yàn)的檢驗(yàn)統(tǒng)計(jì)量χ2服從自由度為v=(k-1)(2-1)=k-1、非中心參數(shù)δ≠0的非中心χ2分布(結(jié)果變量為二值變量)。中心χ2分布只是非中心χ2分布的一個特例(非中心參數(shù)δ=0)。經(jīng)過對非中心χ2分布的研究及相應(yīng)SAS函數(shù)的運(yùn)用,筆者發(fā)現(xiàn):λ值就是使χ2檢驗(yàn)同時滿足本文“ψ值的含義”一節(jié)所列出的兩個條件時的非中心χ2分布的非中心參數(shù)值δ。
2.λ值的計(jì)算方法
(1)確定結(jié)果變量為二值變量的單因素多水平設(shè)計(jì)一元定性資料χ2檢驗(yàn)的Ⅰ型錯誤概率α、Ⅱ型錯誤概率β及自由度v;
(2)根據(jù)α和v的值,用χ2分布的分位數(shù)函數(shù)CINV,計(jì)算出自由度為v、分布曲線下左側(cè)概率為1-α的中心χ2分布的分位數(shù)CINV(1-α,v);
(3)根據(jù)剛剛計(jì)算出來的中心χ2分布的分位數(shù)CINV(1-α,v)、v和β的值,用計(jì)算χ2分布的非中心參數(shù)值的函數(shù)CNONCT,計(jì)算出自由度為v、分位數(shù)為CINV(1-α,v)、分布曲線下左側(cè)概率為β的非中心χ2分布的非中心參數(shù)值δ=CNONCT(CINV(1-α,v),v,β)。該非中心參數(shù)值就是所要計(jì)算的λ值,即λ=δ。
比如:當(dāng) α =0.05,β =0.10,v=4 時,
λ = δ=CNONCT(CINV(1 - α,v),v,β)
=CNONCT(CINV(1 - 0.05,4),4,0.10)
=15.405 051 859
此結(jié)果與λ值表中查到的結(jié)果相吻合。
采用SAS函數(shù) PROBMC來計(jì)算 Newman-Keuls檢驗(yàn)用的q臨界值,既克服了查q臨界值表的不便,又提高了涉及q臨界值的SAS程序的自動化水平,還降低了相關(guān)問題的求解難度。借助SAS函數(shù)PROBMC,不但可以計(jì)算相應(yīng)條件下Newman-Keuls檢驗(yàn)的臨界值或統(tǒng)計(jì)量對應(yīng)的概率值,而且可以計(jì)算相應(yīng)條件下 其 他 四 種 檢 驗(yàn) (One-sided Dunnett、Two-sided Dunnett、Maximum Modulus、Williams)的臨界值或統(tǒng)計(jì)量所對應(yīng)的概率值。
此外,TUKEY法、DUNCAN法和 REGWQ法三種兩兩比較法用的q臨界值或尾端概率值的計(jì)算問題也一并得到了解決,因?yàn)檫@三種方法跟Newman-Keuls檢驗(yàn)法一樣,都要用到q臨界值,而且都可以通過SAS函數(shù)PROBMC來進(jìn)行計(jì)算,唯一需要注意的是,它們各自所對應(yīng)的臨界值q(a,v,α)中的參數(shù)a及α的含義有所不同。
在借助SAS函數(shù)FNONCT及FINV計(jì)算ψ值時,我們發(fā)現(xiàn):當(dāng)與檢驗(yàn)統(tǒng)計(jì)量F的分子和分母對應(yīng)的兩個自由度中至少有一個過大時,SAS程序的返回值為缺失值。這是由于前述條件下SAS函數(shù)FNONCT因計(jì)算過程不收斂所導(dǎo)致的結(jié)果。可見,用SAS函數(shù)FNONCT及FINV來計(jì)算ψ值的方法并不完美,尚有一些缺陷。不過,實(shí)際工作中幾乎不會有檢驗(yàn)統(tǒng)計(jì)量F的分子對應(yīng)的自由度超過1萬或分母對應(yīng)的自由度超過10億的情況。因此,用SAS函數(shù)FNONCT及FINV來計(jì)算ψ值,完全能夠滿足實(shí)際需要。
1.錢俊,陳平雁.假設(shè)檢驗(yàn)中計(jì)算觀察檢驗(yàn)效能的意義的探討.中國衛(wèi)生統(tǒng)計(jì),2005,22(3):133-137.
2.郭靜,徐勇勇,何大衛(wèi).多組比較樣本含量及檢驗(yàn)效能的線性算圖估計(jì).中國衛(wèi)生統(tǒng)計(jì),2002,19(2):94-95.
3.錢俊,陳平雁.樣本率多重比較方法的模擬研究.中國衛(wèi)生統(tǒng)計(jì),2009,26(2):131-134.
4.劉桂芬主編.醫(yī)學(xué)統(tǒng)計(jì)學(xué).第2版.北京:中國協(xié)和醫(yī)科大學(xué)出版社,2007:66-68,177-178,408,422.
5.中國科學(xué)院數(shù)學(xué)研究所概率統(tǒng)計(jì)室編.常用數(shù)理統(tǒng)計(jì)表.北京:科學(xué)出版社,1974:14-16,108-112.
6.郭祖超主編.醫(yī)用數(shù)理統(tǒng)計(jì)方法.第3版.北京:人民衛(wèi)生出版社,1988:899.
7.倪宗瓚主編.衛(wèi)生統(tǒng)計(jì)學(xué).第4版.北京:人民衛(wèi)生出版社,2002:251.
8.馬斌榮主編.醫(yī)學(xué)科研中的統(tǒng)計(jì)方法.第3版.北京:科學(xué)出版社,2005:192.
9.孫振球主編.醫(yī)學(xué)統(tǒng)計(jì)學(xué).北京:人民衛(wèi)生出版社,2002:526.
10.楊樹勤主編.衛(wèi)生統(tǒng)計(jì)學(xué).第3版.北京:人民衛(wèi)生出版社,1996:148-150,202,217,220.
11.Louis L,F(xiàn)ran?ois-A D.Statistical Tables,Explained and Applied.World Scientific Publishing Company,New Jersey.London.Singapore.HongKong,2002:65-71.
12.朱世武編著.SAS編程技術(shù)教程.北京:清華大學(xué)出版社,2007:519-522.
13.http://www.sfu.ca/sasdoc/sashtml/lgref/z1016947.htm.