【R-数据科学】方差分析（ANOVA）在R语言中如何实现 ...

引言方差分析，ANOVA，Analysis of Variances，又叫变异数分析，主要功能是“检验两个或多个样本的平均数的差异是否显著”。

其本质是： 分析样本的组间变异程度（方差） ... 首发于数据分析：Python+R+Stata无障碍写文章登录/注册引言方差分析，ANOVA，AnalysisofVariances，又叫变异数分析，主要功能是“检验两个或多个样本的平均数的差异是否显著”。

其本质是：分析样本的组间变异程度（方差）和组内变异程度相比，是否足够大。

如果够大，证明这个分组因素，对于样本均值的影响是显著的。

F检验。

ANOVA由F检验实现，因此某些场合也被称为F检验。

1方差分析那个故事一提到方差分析，我就想起自己一个不堪回首的经历。

那是我刚到北大读硕士的第一还是第二年，一个学医的中学同学在301医院进修，问我可以用Stata做方差分析不？（他研究的是很典型的临床医学问题，关注的是某些手术方式是否效果更好之类）我当时一愣，方差分析什么鬼，和t检验听上去有点关系？感觉好low啊。

他们医学不用跑回归的嘛？我们用的统计方法怎么差那么远？不管了，反正我们计量经济学肯定牛——我们控制了很多变量，我们用多期面板数据，我们搞因果推断，我们更靠谱。

然后我自己就百度了下，匆匆忙忙看了下关于方差分析的介绍，看了Stata的官方视频，了解个大概。

坦白说，当时我没弄透彻什么是方差分析，Stata里面跑出来的结果也说得不是很清楚。

总之，当时的我没能帮上那个忙。

之后很长一段时间，我都在纠结这个事情。

我在想，方差分析属于统计学的内容。

统计学我是学过的，概率论丢硬币嘛，统计学就讲假设检验和参数估计，最后会提到OLS。

至于方差分析，它肯定跑不出概率论和数理统计的范畴！但是。

。

但问题是，我就不会做方差分析？直到有一天，我翻到了一本忘了哪个国内教授写的《统计学》一书。

里面提到了卡方检验，方差分析，于是就认真看了一遍。

发现这本书说的很详细，也提供具体的案例分析。

于是我才想明白，这个方差分析简直就是我自己的知识盲区——之前上课用的统计学用的教材，只是简单略过方差分析；计量经济学里面会用到t检验，但是几乎也不会提到方差分析。

好吧，既然这个东西那么基础，想必也会很重要。

反正难的我不懂多，不如先搞熟基础。

于是就认真把这个方差分析，照着教材学了两三遍，用Excel、Stata、R和Python，都试着实现了一遍。

2ANOVA在四大统计软件中的实现Stata还是最简单高效。

就两个命令，一个叫oneway，一个叫anova。

oneway用于做单因素方差分析，anova则用于做多因素方差分析。

Excel也很方便（以2019为例）：【数据】-【数据分析】-【单因素方差分析/多因素方差分析】（界面如下截图）R：aov()函数。

相对Stata啰嗦很多，但过程也细致严谨很多。

Python：statsmodels包中的anova_lm()函数。

最啰嗦，当然也就最严谨了。

Excel2019的ANOVA界面2单因素方差分析这里，我们用统计软件之王R（我封给它的）来实现吧。

这里我们使用的是R自带的Iris数据集，分析不同的花的亚种，其花瓣的特征如长度，是否存在差异。

其中Species将鸢（yuan1）尾花分为三种，其它4个指标分别表示鸢尾花的花的特征。

比如Sepal.Length表示花瓣的长度。

这里Species成为分组变量，也就是所谓的“因素”（相当于回归分析中的x）；花瓣的长度为结果变量（相当于回归分析中的y），是比较组间均值差异的变量指标。

>data(iris)#导入R自带的Iris数据 >aov1aov1 Call: aov(formula=Sepal.Length~Species,data=iris) Terms: SpeciesResiduals SumofSquares63.2121338.95620 Deg.ofFreedom2147 Residualstandarderror:0.5147894 Estimatedeffectsmaybeunbalanced >summary(aov1) DfSumSqMeanSqFvaluePr(>F) Species263.2131.606119.3<2e-16*** Residuals14738.960.265 --- Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1结果解释：p值远小于0.001（经济学最小一般是0.01），说明鸢尾花的品种这个因素，对鸢尾花花瓣的长度，有显著的影响；或通俗点说，不同种类的鸢尾花的花瓣长度（均值）显著不同（至少有两种显著不同）。

为了进一步得到更详细的两两比较，我们可以使用R中的TukeyHSD()函数实现。

#使用TukeyHSD()函数进行组间花瓣长度均值的，两两差异研究。

>tukeytukey=as.data.frame(tukey$Species) >tukey$pair=rownames(tukey) #PlotpairwiseTukeyHSDcomparisonsandcolorbysignificancelevel >ggplot(tukey,aes(colour=cut(`padj`,c(0,0.01,0.05,1), label=c("p<0.01","p<0.05","Non-Sig"))))+ theme_bw(base_size=16)+ geom_hline(yintercept=0,lty="11",colour="grey30",size=1)+ geom_errorbar(aes(pair,ymin=lwr,ymax=upr),width=0.2,size=1)+ geom_point(aes(pair,diff),size=2)+ labs(colour="")+ theme(axis.text.x=element_text(size=14))结果显示，三种差异均显著不为0。

3多因素方差分析这里我们使用R自带的ToothGrowth数据集。

该数据关于在两因子交叉的分类组别中，豚鼠牙齿生长的情况。

具体来说，实验中通过将60只豚鼠随机分为6组；分组方法基于喂食方法的不同。

实验中有两种基础的喂食方法，橙汁或维生素C。

两种基础喂食方法中，加入三种不同水平（0.5mg、1mg或2mg）的抗坏血酸。

最终共2大组x3种抗坏血酸水平=6小组，每小组10只豚鼠。

最后，牙齿长度为因变量。

>data(ToothGrowth)#导入数据 #将抗坏血酸浓度数据转化为因子变量 >ToothGrowth$dosesummary(ToothGrowth) lensuppdose Min.:4.20OJ:30d0.5:20 1stQu.:13.07VC:30d1:20 Median:19.25d2:20 Mean:18.81 3rdQu.:25.27 Max.:33.90 >str(ToothGrowth) 'data.frame': 60obs.of3variables: $len:num4.211.57.35.86.41011.211.25.27... $supp:Factorw/2levels"OJ","VC":2222222222... $dose:Factorw/3levels"d0.5","d1","d2":1111111111... >aov2aov2 Call: aov(formula=len~dose+supp,data=ToothGrowth) Terms: dosesuppResiduals SumofSquares2426.434205.350820.425 Deg.ofFreedom2156 Residualstandarderror:3.82759 Estimatedeffectsmaybeunbalanced >summary(aov2) DfSumSqMeanSqFvaluePr(>F) dose22426.41213.282.81<2e-16*** supp1205.4205.414.020.000429*** Residuals56820.414.7 --- Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1结果说明：dose和supp两个因素系数的p值都小于0.001，表示两种因素分组的实验下，豚鼠牙齿生长速度显著不同。

换句话说，两种因素，都对结果产生影响。

4有交互作用的多因素方差分析以上展示的是，无交互作用的多因素方差分析。

下面尝试分析两种因素交互作用（有点像OLS中的交互项）下，结果是否显著不同。

>aov3aov3 Call: aov(formula=len~dose*supp,data=ToothGrowth) Terms: dosesuppdose:suppResiduals SumofSquares2426.434205.350108.319712.106 Deg.ofFreedom21254 Residualstandarderror:3.631411 Estimatedeffectsmaybeunbalanced >summary(aov3) DfSumSqMeanSqFvaluePr(>F) dose22426.41213.292.000<2e-16*** supp1205.4205.415.5720.000231*** dose:supp2108.354.24.1070.021860* Residuals54712.113.2 --- Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1结果显示，两种因素的交互作用dose:suup在0.05的水平上，对结果有显著影响。

另外，dose和supp两个因素对结果也分别在0.01的水平上对结果产生影响。

最后还有个问题。

t检验也是做分组均值差异检验，那它和方差分析有什么区别？事实上，t检验是两组间均值差异比较，ANOVA是多组比较。

所以说，t检验属于方差分析，是方差分析最基础的应用。

好了，经济学的实证论文中也经常用的t检验，相当于也经常用方差分析了。

顺便推荐一波本人长期自学R语言的书籍（有涉及编程入门，有偏数学统计、有偏计量）：-----全文结束-----编辑于2022-03-2114:15方差分析R（编程语言）Stata​赞同76​​3条评论​分享​喜欢​收藏​申请转载​文章被以下专栏收录数据分析：Python+R+StataPython+R+Stata+Excel