一團糟的數據做方差分析,幸好還有Wilcox的robust ANOVA(R統計專用)

伍教練在倫敦大學學院完成的碩士課題是Bilingual Advantage In Visual Selective Attention Task Under Various Perceptual Loads with European and Chinese samples，主要用到混合設計方差分析(mixed ANOVA)，對比雙語人士和單語人士在4種不同知覺負荷下的答題表現，探測雙語認知優勢。

在北大的時候我學過至少4次統計學，每次都是似懂非懂，糊裡糊塗通過考試之後就忘光了。在UCL攻讀語言發育方向，統計學的要求比較高，足足上了兩個學期，要上機考試還要寫2800字的案例分析論文。最近跟語言學另一個方向的同學溝通，才知道他們壓根沒這門頭痛的課！不過這次（希望是最後一次）學統計的經歷是最深刻也最特別的，不光是用英語學習，而且還覺得很有用。寫統計課論文是對已在學術期刊上公開發表的論文進行挑刺，為此我翻遍了英國統計大神Andy Field的Discovering Statistics Using IBM SPSS，還找了很多文獻研讀，鑽牛角尖的過程居然第二次發現統計學很好玩（第一次是在高中的時候用Pascal語言編了個簡陋的統計程序）。

在分析自己課題的數據時，更是要吃透統計學了。即使是已發表的研究，有相當比例的統計分析是有漏洞的，這也是統計老師讓我們挑刺的原因。但是，剛學完一些統計基礎的研究生，水平不會高過那些發論文的博士、教授，難道他們不知道自己的統計方法有問題嗎？現在我明白了，他們是「臣妾做不到」，因為實際收集到的數據往往是一團糟的，做同樣的測試有人認真做有人瞎做，有天才也有低能兒，最後你發現數據有缺失的，有高或低得離譜的，更麻煩的是在不大的樣本裡，有些數據就是不符合正態分布，或者缺乏方差齊性、圓性等，導致經典的統計方法如ANOVA無法使用。其中最麻煩的是正態分布，數據不符合就是不能做常規的ANOVA，連做Greenhouse-Geisser修正的機會都沒有，而混合設計方差分析你還無法轉做非參數檢驗。

我在之前研究相同問題的研究生論文(Lukic, 2017) 裡見到作者「霸王硬上弓」，不知道從哪兒引用了一篇論文說「個別次數的數據不符合正態分布也可以做ANOVA，因為ANOVA本身很強大」 (Schmider, Ziegler, Danay, Beyer, & Bühner, 2010)。Lukic引用的這個論文估計救了很多人，Google Scholar顯示被引用了541次。Lukic也盡力了，後面核心數據實在太亂，不敢用ANOVA了，就把數據縱向做Friedman檢驗（一維ANOVA的非參數檢驗替代方法），橫向逐層做Mann-Whitney檢驗（非參數檢驗），我都看傻眼了。還有一位參與科研的本科生直接選擇無視，用ANOVA分析出全套結果……研究生導師建議嘗試對數轉換，但Andy Field說這不一定能解決問題，還可能產生新的問題。

其實這個問題很經典，Andy Field在他的教材裡頭就多次寫道，如果每次他回答非正態分布數據做混合設計方差分析的問題假如收取1英鎊，那麼他就有錢買一套新的架子鼓了(Field, 2013)。他的回答是：

There are robust methods that can be used based on bootstrapping (Wilcox, 2012). They can’t be done directly in SPSS but they can be implemented in R, and are explained in the sister textbook for that package.

美國統計大神Rand Wilcox有一套基於自展法的強力的辦法做ANOVA，但是要用到R語言，一種學習曲線很陡峭的方法，而我學了這麼多年統計學都是用SPSS。沒辦法，只好根據Andy Field在Discovering statistics Using R這本書中的建議下載安裝R，後來自學了YouTube上一套教程又改用更方便的R Studio，都是開源免費的，後者在用之前還得先安裝好R。

R是執行命令行的，沒有SPSS那麼傻瓜，開始用起來需要適應。而且要用R做某方面的統計，還得額外下載專門的擴展功能包，還得每次都重新調用。按照書上一步步操作就好，但跟無數電腦軟體問題一樣，到了實際操作的時候就會出現各種意外，出錯導致執行不下去。Wilcox的統計大法要安裝WRS這個包，但老是安裝失敗，原來2018年6月改包已經升級到WRS2，連統計函數的自變量格式也發生了重大變化。最關鍵的是WRS要求分析的數據不能是long結構，而要改成wide結構，而WRS2又要求數據回到long結構。Andy Field的Discovering statistics Using R（2012版）和Rand Wilcox的Introduction to robust estimation and hypothesis testing（2012版）兩本能下載到的書都是介紹WRS一代，能說清楚WRS2的Introduction to robust estimation and hypothesis testing（2017版）要買或者在UCL的教育學院圖書館借。

我本來打算去那個世界第一教育學院借書的，最後還是宅著解決了問題。原來，Wilcox一直在自己的主頁更新他的大法：

https://dornsife.usc.edu/labs/rwilcox/software/

用R的常規package安裝方法不成功（聯想到WRS2的升級，我懷疑大神跟R社區的人鬧翻了），我只能下載他的txt函數集，用source(file.choose())直接導入——這些括號誤導我把文件名和地址填進去，折騰半天才知道不用填，直接ctrl+回車執行然後彈窗選文件。終於可以統計了，我用到的函數如下：

setwd("C:/R")  # 設置默認工作檯文件夾。R裡頭只認/，正好跟Windows默認的文件地址\相反，這也是坑人。

hw<-read.csv("hw.csv", header = TRUE)  # 讀入數據文件。建議把Excel和SPSS另存為csv或dat。

english=bw2list(hw,4,c(14,17,20,23))  # 用bw2list定義一個數據矩陣english，其中4是我的數據文件的第4列，組間因素（兩個語種），c()裡頭是組內因素的4個層次（知覺載荷的4個水平）。

bwtrim(2,4,english)  #  2*4混合設計的強力ANOVA。結果如下：

$`Qa`

[1] 0.01028409

 

$Qa.p.value

          [,1]

[1,] 0.9202307

 

$Qb

[1] 32.37664

 

$Qb.p.value

             [,1]

[1,] 1.578525e-07

 

$Qab

[1] 0.6168458

 

$Qab.p.value

          [,1]

[1,] 0.6128597

WRS一代自由度欠奉。$`Qa`是組間因素主效應，$Qb是組內因素主效應，$Qab是互動效應。Andy Field書中的tsplit函數和bwtrim完全等價。如果是用WRS2，那麼要寫成bwtrim(error ~ speakerstatus*level, id = participantID, data = hw)，不用定義數據矩陣了，直接用數據列計算，但數據要搞成long形式：error是要算的數據，speakerstatus、level是指定數據屬於哪一組哪一層，對應唯一id，在Excel上乾坤大挪移一下就好。

                      value df1     df2 p.value

speakerstatus        0.0103   1 20.0732  0.9202

level               32.3766   3 18.2990  0.0000

speakerstatus:level  0.6168   3 18.2990  0.6129

這樣的顯示結果就比較清晰，比SPSS還清楚，還提供自由度。其他WRS2的函數用法就沒有研究了，我懶得去借那本Introduction to robust estimation and hypothesis testing（2017版），你要搞的話可以等Library Genesis更新。WRS的一代二代計算結果應該是一樣的，就是精確位數有點出入，我用兩代方法計算過所有二維ANOVA，確保無誤。以下還是介紹WRS一代。

ESmainMCP(2,4,english,tr=0.2,nboot=100,SEED=T)  # 計算主效應的effect size，多於2組的因素會兩兩給出一個effect size，不像SPSS總是報1個，請自己取捨。我請教了北大的統計高手仍搞不明白，就報一個範圍拉倒。

$`Factor.A`

     Level Level  Effect.Size

[1,]     1     2 9.517426e-05

 

$Factor.B

     Level Level Effect.Size

[1,]     1     2   0.1151448

[2,]     1     3   0.8777893

[3,]     1     4   0.9323554

[4,]     2     3   0.7370218

[5,]     2     4   0.7749252

[6,]     3     4   0.5531625

esImcp(2,4,english,tr=0.2,nboot=100,SEED=T)  # 互動效應的effect size，同上。

mcp2atm(2,4,english, tr=0.2,alpha=0.05,grp=NA,op=F)  # 做組內因素的post hoc兩兩比較，結果看上去很恐怖，完全不是SPSS那種一目了然的表格。Andy Field在他的R書539頁有詳細解讀，簡而言之就是兩種看法，1是看test值是否大過crit值，是就significant；2是看psihat的上下限是否跨越0。1和2有時結果不一致，一個說顯著一個說不顯著，我一般就看2（比較符合普通ANOVA的結果）。$Factor.B就是兩兩比較，一行對應下面相應矩陣的一列，1和-1比較，誰是誰請自己琢磨，就是那幾個組合。

$`Factor.A`$test

     con.num      test     crit       se       df

[1,]       1 0.1365107 2.013117 4.654467 45.81346

 

$`Factor.A`$psihat

     con.num    psihat  ci.lower ci.upper   p.value

[1,]       1 0.6353846 -8.734601 10.00537 0.8920157

 

 

$Factor.B

$Factor.B$`n`

[1] 21 21 21 21 21 21 21 21

 

$Factor.B$test

     con.num       test     crit       se       df

[1,]       1 -0.8093546 2.753126 1.715517 42.43072

[2,]       2 -6.8477007 2.813179 2.604917 28.14531

[3,]       3 -8.8196922 2.819383 3.924609 27.34946

[4,]       4 -6.5736440 2.840275 2.502300 24.97138

[5,]       5 -8.6137220 2.833324 3.857262 25.71528

[6,]       6 -3.8772795 2.777507 4.326785 34.90722

 

$Factor.B$psihat

     con.num     psihat   ci.lower   ci.upper      p.value

[1,]       1  -1.388462  -6.111495   3.334572 4.228268e-01

[2,]       2 -17.837692 -25.165791 -10.509593 1.879895e-07

[3,]       3 -34.613846 -45.678822 -23.548870 1.749087e-09

[4,]       4 -16.449231 -23.556452  -9.342010 6.954278e-07

[5,]       5 -33.225385 -44.154260 -22.296510 4.713615e-09

[6,]       6 -16.776154 -28.793830  -4.758477 4.459276e-04

 

 

$Factor.AB

$Factor.AB$`n`

[1] 21 21 21 21 21 21 21 21

 

$Factor.AB$test

     con.num       test     crit       se       df

[1,]       1 -0.5582529 2.753126 1.715517 42.43072

[2,]       2  0.4521036 2.813179 2.604917 28.14531

[3,]       3  0.7628393 2.819383 3.924609 27.34946

[4,]       4  0.8533687 2.840275 2.502300 24.97138

[5,]       5  1.0244412 2.833324 3.857262 25.71528

[6,]       6  0.4197468 2.777507 4.326785 34.90722

 

$Factor.AB$psihat

     con.num     psihat   ci.lower  ci.upper   p.value

[1,]       1 -0.9576923  -5.680726  3.765341 0.5796058

[2,]       2  1.1776923  -6.150407  8.505791 0.6546596

[3,]       3  2.9938462  -8.071130 14.058822 0.4520892

[4,]       4  2.1353846  -4.971836  9.242605 0.4015668

[5,]       5  3.9515385  -6.977336 14.880413 0.3151679

[6,]       6  1.8161538 -10.201523 13.833830 0.6772428

 

 

$All.Tests

[1] NA

 

$conA

     [,1]

[1,]    1

[2,]    1

[3,]    1

[4,]    1

[5,]   -1

[6,]   -1

[7,]   -1

[8,]   -1

 

$conB

     [,1] [,2] [,3] [,4] [,5] [,6]

[1,]    1    1    1    0    0    0

[2,]   -1    0    0    1    1    0

[3,]    0   -1    0   -1    0    1

[4,]    0    0   -1    0   -1   -1

[5,]    1    1    1    0    0    0

[6,]   -1    0    0    1    1    0

[7,]    0   -1    0   -1    0    1

[8,]    0    0   -1    0   -1   -1

 

$conAB

     [,1] [,2] [,3] [,4] [,5] [,6]

[1,]    1    1    1    0    0    0

[2,]   -1    0    0    1    1    0

[3,]    0   -1    0   -1    0    1

[4,]    0    0   -1    0   -1   -1

[5,]   -1   -1   -1    0    0    0

[6,]    1    0    0   -1   -1    0

[7,]    0    1    0    1    0   -1

[8,]    0    0    1    0    1    1

 

z=bw2list(hw, 4, c(15,18,21,24,25,26,27,28))  # 為三維強力ANOVA定義一個2*2*4的數據矩陣z，結構為between-by-within-by-within。Andy Field說自己也搞不太清楚3-way ANOVA，我就更不懂了，但看到Rand Wilcox提供了函數就照做，而麻煩的是在定義數據矩陣的時候居然用了和二維一樣的定義函數bw2list，這就讓人困惑了，因為另一種三維結構between-by-between-by-within設計是用bbw2list的，自變量的規定是很清楚的。Rand Wilcox的書這個部分我沒看太懂，幸好Google到網上一份義大利的暑期培訓班講義，裡頭講到bw2list的三維用法：c(15,18,21,24,25,26,27,28)之中，前四列是第一層，後四列是第二層。如果是2*3*4結構，我就只能投降了。

bwwtrim(2,2,4,z)  # 三維強力ANOVA。

$`Qa`

         [,1]

[1,] 5.011362

 

$Qa.p.value

           [,1]

[1,] 0.02540106

 

$Qb

         [,1]

[1,] 48.11958

 

$Qb.p.value

             [,1]

[1,] 7.196466e-12

 

$Qc

          [,1]

[1,] 0.8477969

 

$Qc.p.value

          [,1]

[1,] 0.4678504

 

$Qab

           [,1]

[1,] 0.09024845

 

$Qab.p.value

         [,1]

[1,] 0.763924

 

$Qac

          [,1]

[1,] 0.3643552

 

$Qac.p.value

          [,1]

[1,] 0.7787607

 

$Qbc

         [,1]

[1,] 2.869379

 

$Qbc.p.value

           [,1]

[1,] 0.03550038

 

$Qabc

          [,1]

[1,] 0.5556118

 

$Qabc.p.value

          [,1]

[1,] 0.6444565

看到了，這次abc三個因素兩兩互動、三個一起互動都算出來了。Post hoc和effect size就不知道怎麼搞了，就是follow up其中的顯著差別便可。

yuend(hw$PL1dD, hw$PL1dS, tr = .2, alpha =0.05)  # 獨立樣本t檢驗，提供自由度。

yuenv2(hw$PL1dD, hw$PL1dS,tr=0.2,alpha=0.05)  # t檢驗的effect size

不知道robust t-test和非參數檢驗Mann-Whitney誰更有power，反正我就全套用Wilcox的方法了。

估計不少讀研的同學可能也會遇到類似的問題，我推薦英美兩位統計大神近幾年的新成果。可能統計高手在Matlab、SAS上面有更多的手段。伍教練並非統計高手，只是出於實際研究的需要自學了如何在R上實現Rand Wilcox的robust ANOVA，一個接一個解決技術問題，比較有快感（實際上走過的彎路比本文所述多了好幾倍，就不贅述了）。其實這套方法的統計結果跟普通的ANOVA是差不多的，但勝在沒有違規操作，寫到論文中更理直氣壯。

 

 

References

Field, A. (2013). Discovering statistics using IBM SPSS statistics. Sage.

Field, A., Miles, J., Field, Z. (2012). Discovering statistics Using R. Sage.

Lukic, S. (2017, unpublished MSc manuscript). Eyes full of language: examining perceptual capacity in bilingual speakers.

Schmider, E., Ziegler, M., Danay, E., Beyer, L., & Bühner, M. (2010). Is it really robust?. Methodology.

Wilcox, R. R. (2012). Introduction to robust estimation and hypothesis testing. Academic press.

文/伍君儀，透析英語使用法創始人、暢銷書系列《把你的英語用起來！》《把你的詞彙用起來！》作者、北京大學心理學碩士、倫敦大學學院語言科學（語言發育方向）研究生

透析法微信群：請加微信號DIALYSIS55，添加時請標明「兒童」、「打卡」、「報名」，然後拉你入相應的群。

透析法QQ群：59634714或71929542

透析英語新浪微博：@伍君儀透析英語

讀懂原著，聽懂新聞，詞彙量直奔10,000！「把你的英語用起來」VIP訓練營課程價值只需1個月，您就能從「學英語」飛躍到「用英語」，通過手機平板實現全英化的生活方式，讓您讀懂英文原著，聽懂英語新聞，詞彙量穩步擴充到10,000以上，英語考試一路通關。課程內容課程有效期為1年，每周開課，每次課2小時，總共100小時以上，「4+1」結構滾動開課。

必修課

聽力·短讀

口語·寫作

詞彙·長讀

特色課

打造孩子的英語母語

選修課

主題包括全英化教育、雙語智力開發、英語健康養生等。

每次課程保留精華錄音，伍教練每日布置作業並一對一指導，24/7教學答疑及操作示範，口語圈寫作圈實戰。

報名方法點擊「閱讀原文」，或者向本公眾號發送報名，又或者長按二維碼訪問透析法官方微店：