呂小康 付英濤
摘?要?隨機(jī)化檢驗(yàn)是基于實(shí)驗(yàn)中對(duì)實(shí)驗(yàn)處理的隨機(jī)化分配,通過計(jì)算所有可能分配方法的結(jié)果得出某一統(tǒng)計(jì)量的隨機(jī)化分布,并據(jù)此進(jìn)行實(shí)驗(yàn)效應(yīng)是否存在的統(tǒng)計(jì)推斷。相較基于從某一總體中進(jìn)行重復(fù)隨機(jī)抽樣而得到抽樣分布推論模式,隨機(jī)化檢驗(yàn)不需要正態(tài)總體假定,尤其適合樣本數(shù)據(jù)存在明顯離群值或小樣本情形,更適合作為隨機(jī)化實(shí)驗(yàn)的推論框架。借助免費(fèi)開源的R軟件及相關(guān)軟件包已能快速實(shí)現(xiàn)雙處理組和多處理組均值差比較及其他統(tǒng)計(jì)量比較的隨機(jī)化檢驗(yàn),但在心理統(tǒng)計(jì)教育與應(yīng)用中還需進(jìn)一步推廣。
關(guān)鍵詞?隨機(jī)化分布; 隨機(jī)化檢驗(yàn); 顯著性檢驗(yàn); 置換檢驗(yàn); 心理統(tǒng)計(jì)
分類號(hào)?B841.2
DOI: 10.16842/j.cnki.issn2095-5588.2019.05.001
1?統(tǒng)計(jì)推論中的總體模型與隨機(jī)化模型
目前心理統(tǒng)計(jì)教材中介紹的統(tǒng)計(jì)推論方法大多基于抽樣分布(sampling distribution)的框架,即將樣本數(shù)據(jù)視為從特定總體中隨機(jī)抽樣得到的一次觀測結(jié)果,并考慮對(duì)該總體進(jìn)行同樣本容量的重復(fù)抽樣,以得到所關(guān)注的樣本統(tǒng)計(jì)量的假想分布,進(jìn)而計(jì)算相關(guān)p值來判斷實(shí)驗(yàn)效應(yīng)是否存在。心理統(tǒng)計(jì)的相關(guān)教材、包括其他領(lǐng)域的諸多入門統(tǒng)計(jì)教材一般都只是將隨機(jī)化實(shí)驗(yàn)得到的數(shù)據(jù)默認(rèn)為隨機(jī)樣本處理,對(duì)數(shù)據(jù)的正態(tài)性進(jìn)行檢驗(yàn)后,直接應(yīng)用抽樣分布進(jìn)行統(tǒng)計(jì)推論。但從嚴(yán)格意義上說,這些隨機(jī)化實(shí)驗(yàn)的數(shù)據(jù)并不是在總體中抽樣得到的隨機(jī)樣本,并未直接滿足應(yīng)用抽樣分布進(jìn)行統(tǒng)計(jì)推論的條件。同時(shí),隨機(jī)化實(shí)驗(yàn)的目的通常也不是像基于隨機(jī)抽樣的研究那樣要求將樣本結(jié)論推論至總體、即追求外在效度,而是驗(yàn)證實(shí)驗(yàn)的處理效應(yīng)是否真實(shí)存在、即追求內(nèi)在效度。
事實(shí)上,統(tǒng)計(jì)學(xué)中進(jìn)行統(tǒng)計(jì)推論有兩種模型(Ludbrook, 2005):總體模型(population model)和隨機(jī)化模型(randomization model)。將隨機(jī)化實(shí)驗(yàn)的數(shù)據(jù)視為隨機(jī)樣本屬于總體模型的思想,此模型假定實(shí)驗(yàn)數(shù)據(jù)是相應(yīng)總體的隨機(jī)樣本,統(tǒng)計(jì)推論用到的抽樣分布常為正態(tài)分布或由正態(tài)分布推導(dǎo)出的其他理論分布,對(duì)正態(tài)性假定的依賴性較強(qiáng)。而隨機(jī)化模型不需要總體和樣本的假定,其思想基于實(shí)驗(yàn)中的隨機(jī)化操作,其統(tǒng)計(jì)推論無需依賴正態(tài)性前提、而是一種精確分布——隨機(jī)化分布(randomization distribution),基于此種分布進(jìn)行的顯著性檢驗(yàn)可稱為隨機(jī)化檢驗(yàn)(randomization test)。
基于隨機(jī)化分布進(jìn)行顯著性檢驗(yàn)的思想,從其起源上講并不晚于基于抽樣分布的顯著性檢驗(yàn),但囿于計(jì)算便利性上的欠缺而一直沒有得到足夠重視。隨著計(jì)算機(jī)軟件的發(fā)展,隨機(jī)化分布的計(jì)算或模擬已不是問題,因此其思想又重新得到重視與實(shí)踐。隨機(jī)化檢驗(yàn)早期常用來檢驗(yàn)t檢驗(yàn)或F檢驗(yàn)的正確性,F(xiàn)isher(1936)曾指出與隨機(jī)化方法不一致的結(jié)論(t檢驗(yàn)和F檢驗(yàn))是不合理的,因?yàn)槠渌椒ㄍǔI婕袄碚撋系慕疲鴶?shù)據(jù)的真實(shí)形態(tài)可能并不能夠滿足使用這些近似公式所需要的前提條件。而在當(dāng)代,隨機(jī)化檢驗(yàn)甚至被一些統(tǒng)計(jì)學(xué)家認(rèn)為是“隨機(jī)化實(shí)驗(yàn)情形下的金標(biāo)準(zhǔn)”(Edgington & Onghega, 2007),基于隨機(jī)化分布的檢驗(yàn)?zāi)J郊跋嚓P(guān)探討日漸增多(Basu, 2011; Dugard, 2014; Lu, Ding, & Dasgupta, 2015; Mielke, Berry, & Johnston, 2011),相關(guān)教材也不斷面世(Berry, Johnston, & Mielke, 2014; Berry, Mielke, & Johnston, 2016; Dugard, File, & Todman, 2011)。但此種推論方式目前在國內(nèi)外心理統(tǒng)計(jì)教材中仍介紹不多、在實(shí)踐研究中的應(yīng)用也較少,因此有必要加以說明與推廣。本文將介紹隨機(jī)化分布的理論假定及其實(shí)現(xiàn)方法,用免費(fèi)開源軟件R模擬隨機(jī)化實(shí)驗(yàn)數(shù)據(jù)及其隨機(jī)化分布,對(duì)比隨機(jī)化分布和抽樣分布的異同,并提出相應(yīng)的教學(xué)與應(yīng)用建議。鑒于國內(nèi)心理學(xué)界的教學(xué)與研究目前多數(shù)仍使用商業(yè)化的SPSS或SAS等軟件,對(duì)R軟件的認(rèn)識(shí)與應(yīng)用尚未普及,正文中有少數(shù)R語言命令的示范,所有圖片也使用R軟件繪制,以增進(jìn)國內(nèi)對(duì)該軟件的了解。
2?隨機(jī)化檢驗(yàn)的基本思想
隨機(jī)化檢驗(yàn)的思想其實(shí)在20世紀(jì)初期就已經(jīng)出現(xiàn)。較早的一個(gè)例子可見于統(tǒng)計(jì)學(xué)家Fisher(1935)的“女士品茶”名例。其情境大致如下。一名女士聲稱自己可以辨別奶茶里面先放的奶還是先放的茶,F(xiàn)isher想通過實(shí)驗(yàn)驗(yàn)證這一點(diǎn)。故選取8杯奶茶,隨機(jī)選取其中4杯先放奶,另外4杯先放茶,其他條件保持一致。以隨機(jī)順序讓女士品嘗后辨別出哪4杯先放了奶,哪4杯先放了茶。結(jié)果如表1。
該女士正確辨別出了先加奶的4杯中的3杯,能否據(jù)此說明該女士具有辨別能力?Fisher的推論模式如下:假設(shè)該女士沒有辨別能力,那么其判斷是完全隨機(jī)的,此時(shí)的辨別結(jié)果共有84=70種。其中(從判斷先加奶的4杯來說),0對(duì)4錯(cuò)有1種;1對(duì)3錯(cuò)有16種;2對(duì)2錯(cuò)有36種;3對(duì)1錯(cuò)有16種;4對(duì)0錯(cuò)有1種。則得出此結(jié)果的概率為p=17/70=0.24。從雙側(cè)檢驗(yàn)的角度看,此結(jié)果的p=34/70=0.486。故以0.05的顯著性水平看,不能由此認(rèn)為該女士具有辨別能力;在這個(gè)實(shí)驗(yàn)設(shè)計(jì)中,此顯著性水平下只有該女士全部辨別正確,才可以認(rèn)為她具有辨別能力。
此例中隨機(jī)化檢驗(yàn)的思想已有所體現(xiàn),即在原假設(shè)成立的條件下,根據(jù)所有可能的實(shí)驗(yàn)結(jié)果判斷出現(xiàn)實(shí)際結(jié)果的可能性。但嚴(yán)格來說,F(xiàn)isher的這種檢驗(yàn)思想應(yīng)稱為置換檢驗(yàn)(permutation test,其中permutation即排列的意思)而不是隨機(jī)化檢驗(yàn),因?yàn)檫@一實(shí)驗(yàn)中并未涉及隨機(jī)化分組。隨機(jī)化檢驗(yàn)實(shí)際上利用了置換檢驗(yàn)的計(jì)算方式進(jìn)行p值計(jì)算,但置換檢驗(yàn)本身不一定僅適用于隨機(jī)化實(shí)驗(yàn)的情形,也可適用于隨機(jī)抽樣的情形,是含義更為寬泛的一種檢驗(yàn)。換言之,隨機(jī)化檢驗(yàn)可視為置換檢驗(yàn)的一個(gè)子集(Edgington & Patrick, 2007)。不過,由于兩者的區(qū)分主要在應(yīng)用情境的區(qū)別而非計(jì)算方法的區(qū)分,在實(shí)際使用中兩種方法也常被視為是同一種方法。
此外,F(xiàn)isher(1935)在另一個(gè)配對(duì)比較問題中應(yīng)用了隨機(jī)化檢驗(yàn)的思想。問題是比較有15個(gè)觀測值的配對(duì)樣本,兩個(gè)樣本的差異值是314。在此配對(duì)設(shè)計(jì)中給被試隨機(jī)分配實(shí)驗(yàn)處理共有215種方式,在零假設(shè)條件下,所有可能的結(jié)果中有1726種情況大于等于314,即p=0.05267。而根據(jù)抽樣分布得出的t值算出的p=0.0497,兩者相差較小。
在Fisher(1925)、Geary(1927)、Eden和Yates(1933)、Pitman(1937a, 1937b, 1938)等人置換檢驗(yàn)思想的基礎(chǔ)上,Kempthorne(1955)等人開始著手發(fā)展更具一般性的隨機(jī)化分布及檢驗(yàn)的理論框架。實(shí)驗(yàn)設(shè)計(jì)中比較重要的一步是給被試分配實(shí)驗(yàn)處理,隨機(jī)化要求這一過程在實(shí)驗(yàn)的約束條件下(比如區(qū)組內(nèi))是完全隨機(jī)的,即在有限的所有可能的分配方式中隨機(jī)選取一種,每種方式被選中的概率相等。選定一個(gè)統(tǒng)計(jì)量,計(jì)算零假設(shè)條件下所有可能的分配方式下該統(tǒng)計(jì)量的值,即得到該統(tǒng)計(jì)量的隨機(jī)化分布。根據(jù)實(shí)驗(yàn)結(jié)果得出的統(tǒng)計(jì)量在隨機(jī)化分布中的相對(duì)位置和設(shè)定的顯著性水平,即可得出相應(yīng)p值并判斷是否拒絕零假設(shè)。
Kempthorne(1955)提出,隨機(jī)化分布的應(yīng)用依賴于被試—處理可加性假定(unit-treatment additivity)。其基本思想如下:每名被試的觀測值可以視為被試本身的基本量和所接受處理效應(yīng)的加和。用i表示參加實(shí)驗(yàn)的被試,i=1,2 …,N;t表示實(shí)驗(yàn)處理,t=1,2,…,T。隨機(jī)分配每名被試接受一種實(shí)驗(yàn)處理,如果被試i接受了實(shí)驗(yàn)處理t,得到的觀測值可以表示為:
在假設(shè)實(shí)驗(yàn)處理效應(yīng)相同的條件下,只需要對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行隨機(jī)化排列即可得到所有可能的結(jié)果,進(jìn)而選擇某個(gè)適合的統(tǒng)計(jì)量并計(jì)算它在每種分配方式下的取值,就可以得到該統(tǒng)計(jì)量對(duì)應(yīng)的隨機(jī)化分布。再通過計(jì)算該分布中如此次樣本觀測值這么極端、更為極端值(所謂“極端”即指偏離原假設(shè)的設(shè)定)占整個(gè)分布中所有可能取值的比例,即可得到相應(yīng)p值。這就是利用隨機(jī)化分布進(jìn)行顯著性檢驗(yàn)的基本模式。
隨機(jī)化分布的思想框架較為直觀,并不需要假想存在實(shí)驗(yàn)處理組之外的一個(gè)“(正態(tài))總體”,這對(duì)于實(shí)際研究者理解統(tǒng)計(jì)推論的思維框架是較為便利的。問題在于其計(jì)算過程比較麻煩。隨著樣本量和處理組數(shù)的增加,實(shí)驗(yàn)結(jié)果的可能性會(huì)大大增加;對(duì)于復(fù)雜問題的可能結(jié)果排列需要花費(fèi)相當(dāng)長的時(shí)間,往往超過手動(dòng)計(jì)算或公式推導(dǎo)的可能。因此隨機(jī)化檢驗(yàn)在計(jì)算機(jī)軟件興起之前并未得到太多重視。但隨著計(jì)算機(jī)的發(fā)展,隨機(jī)化結(jié)果的計(jì)算不再是問題,并且可以通過軟件對(duì)隨機(jī)化分布實(shí)現(xiàn)可視化,理論假定更少、結(jié)果更加精確的隨機(jī)化檢驗(yàn)又逐漸開始受到重視(Ludbrook & Dudley, 1998; Rubin, 1991)。
3?基于隨機(jī)化分布的統(tǒng)計(jì)推論示例:基于R的應(yīng)用
隨機(jī)化檢驗(yàn)通常用于檢驗(yàn)不同實(shí)驗(yàn)組之間的處理效應(yīng)是否真實(shí)存在,故通常不能應(yīng)用于單樣本數(shù)據(jù)的情形。這里僅就最一般意義上的雙處理組和多處理組情形做出示例。
3.1?雙處理組情形:與傳統(tǒng)t檢驗(yàn)的比較
不妨先考慮隨機(jī)分配被試到兩種實(shí)驗(yàn)處理的情況,此時(shí)的隨機(jī)化過程發(fā)生在被試之間,實(shí)驗(yàn)結(jié)果為獨(dú)立的雙樣本數(shù)據(jù)。用R模擬服從特定正態(tài)分布的隨機(jī)數(shù),樣本量為7。
此時(shí)=20.29,s2x=18.90,=15.57,s2y=21.95。能否根據(jù)上面的數(shù)據(jù)得出結(jié)論:實(shí)驗(yàn)處理X的效果大于處理Y的效果?
這顯然是一個(gè)單側(cè)檢驗(yàn)問題?;趥鹘y(tǒng)抽樣分布的統(tǒng)計(jì)推斷模式如下:將兩種實(shí)驗(yàn)處理下的結(jié)果看作是來自兩個(gè)獨(dú)立同方差的正態(tài)分布總體的簡單隨機(jī)樣本。結(jié)果可算得t=1.95,p=0.037,在0.05顯著性水平下可認(rèn)為實(shí)驗(yàn)處理X的效應(yīng)大于Y。由于這一檢驗(yàn)公式較為常見,故這里不再具體敘述相關(guān)過程。R中的命令如下:
t.test(X, Y, var.equal=TRUE, alternative="greater") # 執(zhí)行同方差前提下的雙樣本單側(cè)t檢驗(yàn)
基于隨機(jī)化分布的統(tǒng)計(jì)推斷模式如下:實(shí)驗(yàn)?zāi)康氖且獙?duì)比兩個(gè)處理組的差異,最簡單直接的統(tǒng)計(jì)量是兩個(gè)處理組的均值差X-Y。在總體方差未知的情況下,這一統(tǒng)計(jì)量的精確抽樣分布不易從公式推導(dǎo)得出,目前使用t檢驗(yàn)公式僅是一種近似公式。但是,計(jì)算3432種分配方式下的X-Y則可得到此統(tǒng)計(jì)量的精確隨機(jī)化分布(圖1(a))。此例中,分別隨機(jī)分配14名被試到兩種處理中,分配方法共有147=3432種,這一數(shù)字并不龐大,可以利用R中的combn函數(shù)窮盡所有情況,命令為:
calcu.stat<-function(x){
Sx<-sum(x)
Sy<-sum(c(X, Y))-Sx
return(mean(x)-Sy/length(Y))
} # 定義計(jì)算均值差的函數(shù)
stat<-combn(c(X, Y), length(X), calcu.stat) # 計(jì)算3432種分配方法下的均值差
通過計(jì)算這3432種分配方式后的X-Y值,再求出其值大于等于此次實(shí)驗(yàn)結(jié)果中x-y=20.29-15.57=4.72的概率,此即為隨機(jī)化分布中的p值,可求得p值約為0.042。此結(jié)果與基于抽樣分布得到的結(jié)果(0.037)相差并不大。
如要直接對(duì)上述過程進(jìn)行R語言計(jì)算,需要進(jìn)行程序代碼的編寫,可能超出R語言初學(xué)者的要求。但若能靈活使用已有的R包中的命令,則可簡化上述流程。使用perm包(初次使用時(shí)需在聯(lián)網(wǎng)狀態(tài)下用install.packages("perm")自行安裝),使用如下命令即可得出隨機(jī)化檢驗(yàn)的p值:
library(perm) # 調(diào)用perm包
permTS(X,Y,alternative="greater", method="exact.ce") # 對(duì)X、Y兩組數(shù)據(jù)進(jìn)行隨機(jī)化檢驗(yàn)
其中,alternative="greater"表示備擇假設(shè)為處理X的效應(yīng)大于處理Y的效應(yīng), method="exact.ce"表示使用精確隨機(jī)化分布進(jìn)行計(jì)算。另外,permTS命令中的TS表示two samples。結(jié)果給出的p值為0.04225。
除此之外,還可以另一種思路進(jìn)行隨機(jī)化檢驗(yàn),即計(jì)算每種排列情形下的傳統(tǒng)雙樣本t檢驗(yàn)觀測值,從而求出此t值的隨機(jī)化分布。再基于分布計(jì)算此次實(shí)驗(yàn)中的t觀測值在多大程度上偏離原假設(shè)(t=0),從而得到相應(yīng)p值。圖1(b)是覆蓋了t分布曲線的隨機(jī)化分布,p(r)表示隨機(jī)化(randomization)情形下的p值,p(t)表示傳統(tǒng)抽樣分布t檢驗(yàn)下的p值。此次實(shí)驗(yàn)得出的結(jié)果為t=1.95,在此隨機(jī)化分布中p(t≥1.95)=0.042。顯然,以X-Y作為統(tǒng)計(jì)量計(jì)算出的p值和以t值作為統(tǒng)計(jì)量得出的p值相同,這是因?yàn)楦鶕?jù)每種分配方式的結(jié)果計(jì)算出的t值和X-Y是一一對(duì)應(yīng)的。
對(duì)于上述實(shí)驗(yàn)處理數(shù)和樣本量比較少的情況,可以計(jì)算所選統(tǒng)計(jì)量的精確隨機(jī)化分布,這里所謂“精確”(exact),其實(shí)只是“完整”的意思,即窮盡了所有可能的隨機(jī)分配情況。但在實(shí)際研究中經(jīng)常會(huì)遇到實(shí)驗(yàn)處理和樣本量比較多的情況,即使有計(jì)算機(jī)的幫助,計(jì)算所有的分配方式也往往不現(xiàn)實(shí),但可以在R中運(yùn)用sample函數(shù)進(jìn)行一定次數(shù)的模擬,即在所有可能的分配方式中進(jìn)行重復(fù)抽樣,得出近似隨機(jī)化分布。用R對(duì)此例進(jìn)行10000次抽樣得出的近似隨機(jī)化分布,p值為0.039,與精確隨機(jī)化分布稍有差異,這是因?yàn)殡S機(jī)取樣所產(chǎn)生的抽樣誤差所致。
現(xiàn)考慮上文中的隨機(jī)數(shù)為配對(duì)樣本的情形。傳統(tǒng)t檢驗(yàn)?zāi)J较?,可算得p值為0.08248,在0.05的顯著性水平下沒有充分理由認(rèn)為 X的實(shí)驗(yàn)處理效應(yīng)顯著高于 Y 的實(shí)驗(yàn)處理效應(yīng)。R語言命令如下:
t.test(X, Y, paired=TRUE, alternative="greater")
如果采用隨機(jī)化檢驗(yàn),此例中由于隨機(jī)化發(fā)生在每個(gè)被試內(nèi)部,每個(gè)被試共有2種處理分配順序,故共有27=124種分配方式。通過計(jì)算每種分配方式下的配對(duì)均值差,即可得到其精確隨機(jī)化分布。再計(jì)算此分布中大于等于此次實(shí)驗(yàn)結(jié)果中的配對(duì)差值所占的比例,即可得到相應(yīng)p值。這里使用exactRankTests包中的函數(shù)來做示范。
library(exactRankTests)
perm.test(X, Y, alternative="greater", paired=TRUE)
結(jié)果給出的p值為0.09375,與配對(duì)樣本t檢驗(yàn)的結(jié)果也很接近。
3.2?多處理組情形:與傳統(tǒng)方差分析的比較
如果把上例中實(shí)驗(yàn)處理數(shù)擴(kuò)大為四組,其結(jié)果如表3。其中A、B、C和D的數(shù)據(jù)分別取自正態(tài)分布總體N1(20, 52),N2(17, 52),N3(15, 52),N4(13, 52)。其命令如下:
檢驗(yàn)多個(gè)實(shí)驗(yàn)組的均值是否具有顯著差異的常用方法是方差分析。此例中所有樣本數(shù)據(jù)來自同方差的正態(tài)總體,故適宜使用傳統(tǒng)方差分析進(jìn)行顯著性檢驗(yàn),結(jié)果為F=3.74,p=0.0246,在0.05顯著性水平下可認(rèn)為各實(shí)驗(yàn)處理的效應(yīng)并不完全相同。對(duì)應(yīng)的R語言命令如下(這里先對(duì)數(shù)據(jù)格式進(jìn)行變動(dòng),以變成軟件處理所需要的長格式數(shù)據(jù)):
dataABCD<-data.frame(A, B, C, D)
library(tidyr) # 調(diào)用tidyr包,以進(jìn)行數(shù)據(jù)操縱
ABCD<-gather(dataABCD, group, value, factor_key=TRUE) # 將數(shù)據(jù)變成長格式數(shù)據(jù),各分組信息存為一列,命令為group,各取值信息存為另一列,命令為value,factor_key=TRUE 用于確保group為因子變量、即類型變量
fit1<-aov(value~group, data=ABCD) # 進(jìn)行方差分析
summary(fit1) # 給出方差分析結(jié)果
基于隨機(jī)化分布的統(tǒng)計(jì)推斷模式如下。此例中隨機(jī)分配28名被試到4種處理中的方式共有
287×217×147×77=4.73×1014
種。此時(shí)要窮盡所有的分配方式得出精確的隨機(jī)化分布是不現(xiàn)實(shí)的,但可以通過對(duì)所有分配方式的隨機(jī)抽樣得到近似的隨機(jī)化分布。先選取F值作為統(tǒng)計(jì)量,分別計(jì)算每種分配方式下的F統(tǒng)計(jì)量即得其近似隨機(jī)化分布,再根據(jù)此分布、此次觀察結(jié)果得到的F值,即可算出相應(yīng)p值。F值得隨機(jī)化分布如圖2,由此得到的p值為0.0226,這與傳統(tǒng)F檢驗(yàn)的結(jié)果非常接近。
c=10000 # 設(shè)定模擬次數(shù)
F.stat<-numeric(c) # 構(gòu)建F統(tǒng)計(jì)量變量
F.stat[1]<-summary(aov(value~group, data=ABCD))[[1]][4][[1]][1] # 設(shè)此次觀察為第一個(gè)值
set.seed(1234) # 設(shè)定循環(huán)種子數(shù)
for(i in 2:c){
F.stat[i]<-summary(aov(sample(value)~group, data=ABCD)) [[1]] [4][[1]][1]
} # 隨機(jī)抽取其他分配方法的F值
p=mean(F.stat>=F.stat[1]) # 計(jì)算p值
對(duì)上述過程進(jìn)行編程計(jì)算稍顯麻煩,使用用coin包中的函數(shù)可簡介上述過程:
library(coin)
oneway_test(value~group, data=ABCD, distribution=approximate(B=10000))
其中B=10000即表示進(jìn)行10000次模擬。計(jì)算結(jié)果可得p值為0.02433。這與傳統(tǒng)F檢驗(yàn)的p值非常接近。實(shí)際上,當(dāng)方差分析的假定未被明顯違背時(shí),隨機(jī)化檢驗(yàn)的優(yōu)勢(shì)并不明顯。但當(dāng)數(shù)據(jù)明顯呈現(xiàn)偏態(tài)、尤其是各組方差相差較大時(shí),使用傳統(tǒng)檢驗(yàn)方法會(huì)存在統(tǒng)計(jì)上二類錯(cuò)誤的增大問題,此時(shí)使用隨機(jī)化檢驗(yàn)來做判定則可使統(tǒng)計(jì)決策更為穩(wěn)健。值得注意的是oneway_test函數(shù)使用的隨機(jī)化檢驗(yàn)方法為Pitman-Fisher法,其檢驗(yàn)統(tǒng)計(jì)量的選取與前面的模擬有所不同,具體可參見Boik(1987)。
下面考慮隨機(jī)區(qū)組設(shè)計(jì)的情形。將表2問題的表述稍作變動(dòng),假設(shè)表中的數(shù)據(jù)是7個(gè)區(qū)組接受4種實(shí)驗(yàn)處理的情況,每個(gè)區(qū)組有4名被試,每名被試隨機(jī)接受一種實(shí)驗(yàn)處理。此時(shí)就變成了一個(gè)隨機(jī)區(qū)組問題,從數(shù)據(jù)處理方法上講則是雙因素方差分析的過程。傳統(tǒng)方差分析結(jié)果為F=3.46,p=0.0382。R中的命令如下:
library(dplyr)
block<-as.factor(c(1:7)) # 生成區(qū)組標(biāo)簽
blockABCD<-data.frame(A, B, C, D, block)
newABCD<-gather(blockABCD, group, value,-block, factor_key=TRUE)
fit2<-aov(value~block+group, data=newABCD)
summary(fit2)
隨機(jī)化檢驗(yàn)的模式如下。在4種處理沒有差異的原假設(shè)下,由于對(duì)隨機(jī)化在區(qū)組內(nèi)進(jìn)行,此時(shí)隨機(jī)化分配方式共有(4×3×2×1)7=4586471424種。雖然用軟件可以一一計(jì)算各分配下的F值,但耗時(shí)較長,故仍可考慮使用隨機(jī)抽樣的方式進(jìn)行模擬。這里只用coin包中的oneway_test函數(shù)進(jìn)行示范。
library(coin)
oneway_test(value~group | block, data=newABCD, distribution=approximate(B=100000)) # | block表示指定區(qū)組名稱
結(jié)果得到的p值為0.03259,與傳統(tǒng)檢驗(yàn)的結(jié)果也很接近。
4?隨機(jī)化分布和抽樣分布的比較:樣本量與離群值的影響
對(duì)于正態(tài)性擬合良好的數(shù)據(jù),抽樣分布是隨機(jī)化分布的良好近似。在圖1(a)和圖2中,將t分布和F分布曲線覆蓋在相應(yīng)隨機(jī)化分布上,兩者均擬合良好。這是因?yàn)橛肦模擬的數(shù)據(jù)均是來自正態(tài)分布總體,數(shù)據(jù)的正態(tài)性可以得到保證。下面以兩個(gè)處理組的情況為例,考慮兩個(gè)影響數(shù)據(jù)正態(tài)性的因素——樣本量和離群值(outliers)。
此處考慮服從偏態(tài)分布和正態(tài)分布的兩種數(shù)據(jù)。對(duì)數(shù)正態(tài)分布是一種右偏態(tài)分布,對(duì)服從對(duì)數(shù)正態(tài)分布的數(shù)據(jù)取對(duì)數(shù)可以得到正態(tài)分布數(shù)據(jù)。圖3(a),(c),(e)的樣本是從對(duì)數(shù)正態(tài)分布中抽樣得到的右偏態(tài)數(shù)據(jù),樣本量分別為5,10和20。圖3(b),(d),(f)的樣本是對(duì)相應(yīng)樣本量的右偏態(tài)數(shù)據(jù)取對(duì)數(shù)得到的正態(tài)數(shù)據(jù)。
取樣本量為5的偏態(tài)樣本程序如下:
set.seed(123)
x<-rlnorm(5, 1)
set.seed(321)
y<-rlnorm(5, 0.5)
其他樣本量的情形可依此得出。圖3中可以直觀地看到不同分布形態(tài)和不同樣本量的數(shù)據(jù)對(duì)隨機(jī)化分布和t分布擬合狀況的影響??梢园l(fā)現(xiàn):(1)偏態(tài)數(shù)據(jù)和正態(tài)數(shù)據(jù)相比而言,正態(tài)數(shù)據(jù)的隨機(jī)化分布和t分布擬合得較好,這一點(diǎn)在樣本量較小的情況下體現(xiàn)得尤其明顯;(2)當(dāng)樣本量較小的時(shí)候,即便對(duì)于正態(tài)性數(shù)據(jù)而言,抽樣分布與隨機(jī)化分布的擬合狀況也是比較差的;(3)而樣本量較大時(shí),即使數(shù)據(jù)呈現(xiàn)偏態(tài),抽樣分布與隨機(jī)化分布也擬合得較好,這正是傳統(tǒng)t檢驗(yàn)在雙處理組數(shù)據(jù)中得到普遍應(yīng)用的原因之一。
另外,比較由隨機(jī)化分布和t分布得出的p值可以作為衡量擬合狀況的一個(gè)指標(biāo)。圖4是對(duì)于偏態(tài)數(shù)據(jù)而言,樣本量從5~50變化時(shí),由隨機(jī)化分布和t分布得出的p值變化。圖4(a)隨機(jī)化分布p值和t分布p值隨樣本量的變化情況,(b)是兩個(gè)p值差隨樣本量的變化情況。
從圖4可以看出,總體而言隨機(jī)化分布和t分布得出的p值是比較一致的。對(duì)于偏態(tài)數(shù)據(jù)而言,樣本量的增加可以有效地使p值減小,尤其是在樣本量增加到30之后。這一方面增大了統(tǒng)計(jì)檢驗(yàn)力,另一方面也使隨機(jī)化分布和t分布得出的p值差異逐漸減小。
下面考慮離群值的影響。以表2中的數(shù)據(jù)為例,圖5中(a),(b),(c),(d)是假設(shè)表2中X組最大值分別為28,38,48,58時(shí)的隨機(jī)化分布??梢园l(fā)現(xiàn),離群值距離均值越遠(yuǎn),隨機(jī)化分布和t分布的擬合狀況越差,根據(jù)兩者得出的p值相差越大。Ernst(2009)通過設(shè)置樣本最小值的方法探討了離群值對(duì)隨機(jī)化分布和t分布擬合狀況的影響,與此結(jié)論一致。由于t檢驗(yàn)假定其數(shù)據(jù)是取自正態(tài)分布總體的隨機(jī)樣本,因此當(dāng)不能保證樣本的隨機(jī)性和正態(tài)性時(shí),t檢驗(yàn)的結(jié)果理論上不能保證其準(zhǔn)確性。而基于隨機(jī)化分布的檢驗(yàn)不需要正態(tài)性假定,也無需隨機(jī)樣本。若其分布呈現(xiàn)偏態(tài)、樣本量較少或存在離群值時(shí),基于隨機(jī)化分布的檢驗(yàn)結(jié)果更令人信服。不過,當(dāng)樣本量較大時(shí),t檢驗(yàn)的結(jié)果對(duì)數(shù)據(jù)的正態(tài)性要求不高,與隨機(jī)化分布的結(jié)果差異不大。
5?總結(jié)與建議
隨機(jī)化檢驗(yàn)依賴于實(shí)驗(yàn)處理的隨機(jī)化分配和處理效應(yīng)的可加性假定,不適合觀測數(shù)據(jù),但較適合具有兩個(gè)或多個(gè)處理組的隨機(jī)化實(shí)驗(yàn)結(jié)果的顯著性檢驗(yàn)。在推理框架上,隨機(jī)化檢驗(yàn)并不需要假想總體的存在,這更有利于初學(xué)者理解隨機(jī)化實(shí)驗(yàn)的實(shí)際情境,更適宜作為隨機(jī)化實(shí)驗(yàn)的統(tǒng)計(jì)推論框架。同時(shí),隨機(jī)化檢驗(yàn)并不需要假定總體的分布形態(tài),因此比基于抽樣分布的假設(shè)更具靈活性,可以適用于傳統(tǒng)顯著性檢驗(yàn)不能勝任的情形。雖然目前隨機(jī)化實(shí)驗(yàn)的顯著性檢驗(yàn)常常通過抽樣分布理論給出,這是因?yàn)楹芏嗲闆r下抽樣分布理論對(duì)數(shù)據(jù)正態(tài)性的要求并不十分嚴(yán)格。隨著實(shí)驗(yàn)處理數(shù)和被試數(shù)的增多,基于抽樣分布的顯著性檢驗(yàn)與基于隨機(jī)化分布的顯著性檢驗(yàn)越來越接近。在許多情況下,兩者在實(shí)驗(yàn)處理數(shù)和被試數(shù)較少時(shí)仍有良好的近似性。因此,運(yùn)用抽樣分布對(duì)隨機(jī)化實(shí)驗(yàn)的結(jié)果進(jìn)行推論多數(shù)情況下是比較可靠的。
但在以下幾種情況中,建議基于隨機(jī)化分布進(jìn)行推論而非抽樣分布,或者至少將隨機(jī)化分布的結(jié)果作為必要參考。(1)樣本數(shù)據(jù)量過少。樣本數(shù)據(jù)量過少帶來的問題是無法保證數(shù)據(jù)的正態(tài)性,即使是從正態(tài)總體中抽取的小樣本,抽樣分布和隨機(jī)化分布的擬合狀況也很差。(2)樣本中存在較嚴(yán)重的離群值。此時(shí)或者找出產(chǎn)生離群值的原因,將離群值進(jìn)行剔除使數(shù)據(jù)符合抽樣分布條件,或者沒有合理的理由剔除離群值,此時(shí)需考慮隨機(jī)化分布方法。(3)需選取的統(tǒng)計(jì)量難以從公式推導(dǎo)得出。隨機(jī)化分布的統(tǒng)計(jì)量選取更加自由,許多統(tǒng)計(jì)量的抽樣分布難以計(jì)算,但是研究者可以根據(jù)需要選擇合適的統(tǒng)計(jì)量計(jì)算其隨機(jī)化分布。例如離群值對(duì)于抽樣分布的統(tǒng)計(jì)量往往影響較大,但對(duì)中位數(shù)、四分位數(shù)之類的統(tǒng)計(jì)量影響就較小,它們的理論分布通常難以求得。除此之外,對(duì)于計(jì)數(shù)數(shù)據(jù)和等級(jí)數(shù)據(jù)而言,常常利用一些近似估計(jì)和檢驗(yàn),可能會(huì)造成近似誤差,此時(shí)隨機(jī)化分布的應(yīng)用更應(yīng)當(dāng)?shù)玫街匾暎‥udey, Kerr, & Trumbo, 2010)。
最后需要說明的是,對(duì)于心理統(tǒng)計(jì)的教育與應(yīng)用而言,隨機(jī)化檢驗(yàn)的思想并不復(fù)雜,但計(jì)算比較繁瑣,教學(xué)過程中可以借助R之類的軟件進(jìn)行輔助教學(xué)或計(jì)算。這需要對(duì)軟件的應(yīng)用有一定了解。鑒于目前R、Python等開源統(tǒng)計(jì)軟件尚未成為國內(nèi)心理統(tǒng)計(jì)主流軟件,加強(qiáng)對(duì)此類軟件的宣傳和應(yīng)用仍是必要的。實(shí)際上,這些軟件已經(jīng)具備不亞于SPSS或SAS的統(tǒng)計(jì)分析功能,在某些方面甚至更勝一籌。促進(jìn)軟件應(yīng)用的開源化、免費(fèi)化,對(duì)于國內(nèi)的心理學(xué)教學(xué)和研究單位節(jié)省統(tǒng)計(jì)軟件方面的支出或避免軟件使用上的版權(quán)問題,也是具有積極意義的。其中,除了文中提到的perm、exactRankTests、coin等R包外,還有ez、AUtests、flip、jmuOutlier、jmuOutlier、treeperm等軟件包提供了豐富的隨機(jī)化檢驗(yàn)函數(shù)與算法,可供研究者進(jìn)一步探索利用。當(dāng)然,相較于發(fā)展更為成熟的總體模型推論,一些復(fù)雜實(shí)驗(yàn)設(shè)計(jì)的隨機(jī)化檢驗(yàn)?zāi)J饺赃€有待開發(fā),基于這一模式的統(tǒng)計(jì)功效與效應(yīng)值探討也還有待深入。
參考文獻(xiàn)
Basu, D.(2011). Randomization analysis of experimental data: The fisher randomization test. Journal of the American Statistical Association, 75(371), 305-325.
Berry, K. J., Johnston, J. E., & Mielke, Jr, P. W.(2014). A chronicle of permutation statistical methods: 1920-2000 and beyond. Cham, Switzerland: Springer.
Berry, K. J., Mielke, Jr, P. W. & Johnston, J. E.(2016). Permutation statistical methods: an integrated approach. Cham, Switzerland: Springer.
Boik, R. J.(1987). The fisher-pitman permutation test: A non-robust alternative to the normal theoryFtest when variances are heterogeneous. British Journal of Mathematical & Statistical Psychology, 40(1), 26-42.
Box, G. E. P., & Anderson, S. L.(1955). Permutation theory in the derivation of robust criteria and the study of departures from assumption. Journal of the Royal Statistical Society, 17(1), 1-34.
Dugard, P.(2014). Randomization tests: A new gold standard?Journal of Contextual Behavioral Science, 3(1), 65-68.
Dugard, P., File, P., & Todman, J.(2011). Single-case and small-n experimental designs: A practical guide to randomization tests. London: Routledge.
Eden, T., & Yates, F.(1933). On the validity of Fishersztest when applied to an actual example of non-normal data. Journal of Agricultural Science, 23(1), 6-17.
Edgington, E. S., & Onghena, P.(2007). Randomization tests. Boca Raton, FL: CRC Press.
Ernst, M. D.(2009). Teaching inference for randomized experiments. Journal of Statistics Education, 7(1).
Eudey, T. L., Kerr, J. D., & Trumbo, B. E.(2010). Using R to simulate permutation distributions for some elementary experimental designs. Journal of Statistics Education, 18(1).
Fisher, R. A.(1925). Statistical Methods for Research Workers. Edinburgh: Oliver and Boyd.
Fisher, R. A.(1935). The Design of Experiments. Edinburgh: Oliver and Boyd.
Fisher, R. A.(1936). The coefficient of racial likeness and the future of craniometry. The Journal of the Royal Anthropological Institute of Great Britain and Ireland, 66, 57-63.
Geary, R. C.(1927). Some properties of correlation and regression in a limited universe. Metron Rivista Internazionale de Statistica, 7, 83-119.
Kempthorne, O.(1955). The Randomization Theory of Experimental Inference. Journal of the American Statistical Association. 50(271), 946-967.
Lu, J., Ding, P., & Dasgupta, T.(2015). Construction of alternative hypotheses for evaluation of randomization tests with ordinal outcomes. Statistics & Probability Letters, 107(12), 348-355.
Ludbrook, J. & Dudley, H. A. F.(1998). Why permutation tests are superior tot- andF-tests in biomedical research. The American Statistician, 52(2), 127-132.
Ludbrook, J.(2005). Randomization based tests. Encyclopedia of statistics in behavioral science. Hoboken, NJ: John Wiley & Sons, Inc.
Mielke, P. W., Berry, K. J., & Johnston, J. E.(2011). Robustness without rank order statistics. Journal of Applied Statistics, 38(1), 207-214.
Pitman, E. J. G.(1937a). Significance tests which may be applied to samples from any populations. Supplement to the Journal of the Royal Statistical Society, 4(1), 119-130.
Pitman, E. J. G.(1937b). Significance tests which may be applied to samples from any populations II: The correlation coefficient test. Supplement to the Journal of the Royal Statistical Society, 4(2), 225-232.
Pitman, E. J. G.(1938). Significance tests which may be applied to samples from any populations III: The analysis of variance test. Biometrika, 29(201), 322-335.
Rubin, D. B.(1991). Practical applications of modes of statistical inference for causal effects and the critical role of the assignment mechanism. Biometrics, 47(4), 1213-1234.