實(shí)驗(yàn)數(shù)據(jù)的隨機(jī)化檢驗(yàn)及R語言實(shí)現(xiàn)

呂小康 付英濤

摘?要?隨機(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)槠渌椒ㄍǔＩ婕袄碚撋系慕疲鴶?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)Ｊ郊跋嚓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)Ｊ较?，可算得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)Ｊ饺赃€有待開發(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.