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