【数据分析 R语言实战】学习笔记 第八章 单因素方差分析与R实现

方差分析泛应用于商业、经济、医学、农业等诸多领域的数量分析研究中。例如商业广告宣传方面，广告效果可能会受广告式、地区规模、播放时段、播放频率等多个因素的影响，通过方差分析研究众多因素中，哪些是主要的以及如何产生影响等。而在经济管理中，方差分析常用于分析变量之间的关系，如人民币汇率对股票收益率的影响、存贷款利率对债券市场的影响，等等。

协方差是在方差分析的基础上，综合回归分析的方法，研究如何调节协变量对因变量的影响效应，从而更加有效地分析实验处理效应的一种统计技术。

8.1单因素方差分析及R实现

(1)正态性检验

对数据的正态性，利用Shapiro-Wilk正态检验方法(W检验)，它通常用于样本容量n≤50时，检验样本是否符合正态分布。

R中，函数shapiro.test()提供了W统计量和相应P值，所以可以直接使用P值作为判断标准，其调用格式为shapiro.test(x)，参数x即所要检验的数据集，它是长度在35000之间的向量。

例:

某银行规定VIP客户的月均账户余额要达到100万元，并以此作为比较各分行业绩的一项指标。这里分行即因子，账户余额是所要检验的指标，先从三个分行中，分别随机抽取7个VIP客户的账户。为了用单因素方差分析判断三个分行此项业绩指标是否相同，首先对二个分行的账户余额分别进行正态检验。

> x1=c(103,101,98,110,105,100,106) > x2=c(113,107,108,116,114,110,115) > x3=c(82,92,84,86,84,90,88) > shapiro.test(x1) Shapiro-Wilk normality test data: x1 W = 0.97777, p-value =0.948 > shapiro.test(x2) Shapiro-Wilk normality test data: x2 W = 0.91887, p-value =0.4607 > shapiro.test(x3) Shapiro-Wilk normality test data: x3 W = 0.95473, p-value =0.7724

P值均大于显著性水平a=0.05,因此不能拒绝原假设，说明数据在因子A的三个水平下都

是来自正态分布的。

(2)方差齐性检验

方差分析的另一个假设：方差齐性，需要检验不同水平卜的数据方差是否相等。R中最常用的Bartlett检验,bartlett.test()调用格式为

bartlett.test(x，g…)

其中，参数X是数据向量或列表(list) ; g是因子向量，如果X是列表则忽略g.当使用数据集时，也通过formula调用函数:

bartlett.test(formala, data, subset，na.action…)

formula是形如lhs一rhs的方差分析公式;data指明数据集:subset是可选项，可以用来指定观测值的一个子集用于分析:na.action表示遇到缺失值时应当采取的行为。

续上例:

> x=c(x1,x2,x3)

> account=data.frame(x,A=factor(rep(1:3,each=7)))

> bartlett.test(x~A,data=account)

        Bartlett test of homogeneity of variances

data:  x by A

Bartlett's K-squared = 0.13625, df = 2, p-value = 0.9341

由于P值远远大于显著性水平a=0.05，因此不能拒绝原假设，我们认为不同水平下的数据是等方差的。

8.1.2单因素方差分析

R中的函数aov()用于方差分析的计算，其调用格式为:

aov(formula, data = NULL, projections =FALSE, qr = TRUE,contrasts = NULL, ...)

其中的参数formula表示方差分析的公式，在单因素方差分析中即为x~A ; data表示做方差分析的数据框:projections为逻辑值，表示是否返回预测结果:qr同样是逻辑值，表示是否返回QR分解结果，默认为TRUE; contrasts是公式中的一些因子的对比列表。通过函数summary()可列出方差分析表的详细结果。

上面的例子已经对数据的正态性和方差齐性做了检验，接F来就可以进行方差分析:

> a.aov=aov(x~A,data=account)

> summary(a.aov)

            Df Sum Sq Mean Sq F value   Pr(>F)    

A            2   2315    1158   82.68 8.46e-10 ***

Residuals   18    252      14                     

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

> plot(account$x~account$A)

Levene检验

Levene检验，它既可以用于正态分布的数据，也可用于非正态分布的数据或分布不明的数据，具有比较稳健的特点，检验效果也比较理想。

R的程序包car中提供了Levene检验的函数levene.test()

> library(car)

> levene.test(account$x,account$A)

Levene's Test for Homogeneity of Variance (center = median)

      Df F value Pr(>F)

group  2  0.0426 0.9584

      18 

由于p值大于a=0.05，不能拒绝原假设，我们认为不同水平下的数据是等方差的。

8.1.3多重t检验

单因素方差分析是从总体的角度上说明各效应的均值之间存在显著差异，但具体哪些水平下的均值存在较人差异无从得知，所以我们要对每一对样本均值进行一一比较，即要进行均值的多重比较。

> p.adjust.methods

[1] "holm"       "hochberg"   "hommel"     "bonferroni" "BH"        

[6] "BY"         "fdr"        "none"      

> attach(account)

> pairwise.t.test(x,A,p.adjust.method="bonferroni")

        Pairwise comparisons using t tests with pooled SD 

data:  x and A 

  1       2      

2 0.0013  -      

3 3.9e-07 6.5e-10

P value adjustment method: bonferroni 

经过修正后的p值比原来会增大很多，这在一定程度上克服了多重t检验增加犯第一类错误的

概率的缺点。从检验结果来看，样本两两之问t检验的p值都很小，说明几个样本之间差异明显。

8.1.4Kruskal-Wallis秩和检验

R内置函数kruskal.test()可以完成Kruskal-Wallis秩和检验，使用如下:

kruskal.test(x, ...) kruskal.test(x, g, ...) kruskal.test(formula, data, subset,na.action, ...)

例:

某制造商雇用了来自三所本地大学的雇员作为管理人员。最近，公司的人事部门已经收集信息并考核了年度工作成绩。从三所大学来的雇员中随机地抽取了三个独立样本，样本量分别为7、6, 7，数据如表所示。制造商想知道来自这三所不同的大学的雇员在管理岗位上的表现是否有所不同，我们通过Kruskal-Wallis秩和检验来得到结论。

> data=data.frame(x=c(25,70,60,85,95,90,80,60,20,30,15,40,35,50,70,60,80,90,70,75),g=factor(rep(1:3,c(7,6,7))))

> kruskal.test(x~g, data=data)

        Kruskal-Wallis rank sum test

data:  x by g

Kruskal-Wallis chi-squared = 8.9839, df = 2, p-value = 0.0112

检验的结果为P=0.0112<0.05，因此拒绝原假设，说明来自这三个不同的大学的雇员在管理岗位上的表现有比较显著的差异。