伍教练在伦敦大学学院完成的硕士课题是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.

文/伍君仪，透析英语使用法创始人、畅销书系列《把你的英语用起来！》《把你的词汇用起来！》作者、北京大学心理学硕士、伦敦大学学院语言科学（语言发育方向）研究生

读懂原著，听懂新闻，词汇量直奔10,000！“把你的英语用起来”VIP视频直播训练营每周开课，透析法创始人伍教练真人主讲

【课程价值】

只需1个月，就从“学英语”质变到“用英语”！通过手机平板实现全英化生活方式，让您读懂英文原著，词汇量稳步扩充到10000以上，英语考试一路通关。

【课程内容】课程有效期为1年，每周讲课2小时，“4+1”结构滚动开课，8小时听说读写必修课+每月更新的丰富选修课。QQ群内24/7教学答疑及操作示范，每日口语写作打卡实战，伍教练指导完成作业。

【报名方法】

登录http://dialysis.taobao.com进入“伍君仪透析英语”淘宝店拍下课程，或者添加伍教练微信号auguwu报名