方差分析与R实现

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

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

单因素方差分析及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,因此不能拒绝原假设,我们认为不同水平下的数据是等方差的。

单因素方差分析 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,不能拒绝原假设,我们认为不同水平下的数据是等方差的。

多重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值都很小,说明几个样本之间差异明显。

Kruskal-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,因此拒绝原假设,说明来自这三个不同的大学的雇员在管理岗位上的表现有比较显著的差异。

双因素方差分析及R实现 无交互作用的分析 例: 某商品在不同地区、不同包装的销售数据

首先为了建立数据集,引入生成因子水平的函数g1(),其调用格式为: gl(n, k, length=n*k,labels=1:n,ordered=FALSE) n是因子的水平个数;k表示每一水平上的重复次数;length=n*k表示总观测数;可通过参数labels对因子的不同水平添加标签;ordered为逻辑值,指示是否排序。

> x=c(20,12,20,10,14,22,10,20,12,6,24,14,18,18,10,16,4,8,6,18,26,22,16,20,10)
> sales=data.frame(x,A=gl(5,5),B=gl(5,1,25))
> sales$B
[1] 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 1 2 3 4 5 12 3 4 5
Levels: 1 2 3 4 5

分析前先对因素A和B作方差齐性检验,使用函数bartlett.test()

> bartlett.test(x~A,data=sales)
 
Bartlett test of homogeneity of variances
 
data: x by A
Bartlett's K-squared =0.66533, df = 4, p-value = 0.9555
 
> bartlett.test(x~B,data=sales)
 
Bartlett test of homogeneity of variances
 
data: x by B
Bartlett's K-squared =1.2046, df = 4, p-value = 0.8773

因素A和B的P值都远大于0.05的显著性水平,不能拒绝原假设,说明因素A, B的各水平是满足方差齐性的。这时再进行双因素方差分析,输入指令

> sales.aov=aov(x~A+B,data=sales)
> summary(sales.aov)
Df Sum Sq Mean Sq F valuePr(>F)
A 4 199.4 49.84 2.303 0.1032
B 4 335.4 83.84 3.874 0.0219 *
Residuals 16 346.2 21.64
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’1

检验的结论:因素B的P值=0.0219<0.05,拒绝原假设,说明销售地区对饮料的销售量有显著影响;而因素A的P值=0.1032>0.05,不能拒绝原假设,因此没有充分的理由可以说明包装方式对销售有明显影响。

有交互作用的分析 R仍然用函数aov()作双因素方差分析,只需将formula改为x~A+B+A:B或x~A*B的形式即可。 例: 不同路段和不同时段的行车时间数据

首先构造数据集,对因素A和B作方差齐性检验,利用函数bartlett.test()

> time=c(25,24,27,25,25,19,20,23,22,21,29,28,31,28,30,20,17,22,21,17,18,17,13,16,12,22,18,24,21,22)
> traffic=data.frame(time,A=gl(2,15,30),B=gl(3,5,30,labels=c("I","II","III")))
> bartlett.test(time~A,data=traffic)
 
Bartlett test of homogeneity of variances
 
data: time by A
Bartlett's K-squared =0.053302, df = 1, p-value = 0.8174
 
> bartlett.test(time~B,data=traffic)
 
Bartlett test of homogeneity of variances
 
data: time by B
Bartlett's K-squared =0.57757, df = 2, p-value = 0.7492

检验结果的P值均远大于显著性水平0.05,说明两个因素下的各水平都满足方差齐性的要求,可以进一步做方差分析。画图来观察一下数据的特点,首先是箱线图。

> op=par(mfrow=c(1,2)) #分割图形区域
> plot(time~A+B,data=traffic)
Hit <Return> tosee next plot:

从图形上单独观察时段和路段对行车时间的影响,可以发现因素的不同水平还是有明显差别的。为了考察因素间的交互作用是否存在,利用函数interaction.plot()绘制交互效应图:

interaction.plot(x.factor, trace.factor,response, fun = mean,type = c("l","p", "b", "o", "c"), legend = TRUE,trace.label =deparse(substitute(trace.factor)),fixed = FALSE,xlab =deparse(substitute(x.factor)),ylab = ylabel,ylim = range(cells, na.rm =TRUE),lty = nc:1, col = 1, pch =c(1:9, 0, letters),xpd = NULL, leg.bg =par("bg"), leg.bty = "n", xtick = FALSE, xaxt = par("xaxt"),axes = TRUE,...) x.factor表示横轴的因子 trace.factor表示分类绘图的因子 response是数值向量,要输入响应变量 fun表示汇总数据的方式,默认为计算每个因子水平下的均值 type指定图形类型 legend是逻辑值,指示是否生成图例 trace.label给出图例中的标签。

> attach(traffic)
> interaction.plot(A,B,time,legend=F)
> interaction.plot(B,A,time,legend=F)

曲线均没有相交,所以可以初步判断两个因素之间应该没有交互作用。用方差分析进行确认:

> traf.aov=aov(time~A*B,data=traffic)
> summary(traf.aov)
Df Sum Sq Mean Sq F value Pr(>F)
A 1 313.63 313.63 84.766 2.41e-09 ***
B 2 261.60 130.80 35.351 7.02e-08 ***
A:B 2 6.67 3.33 0.901 0.42
Residuals 24 88.80 3.70
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

根据检验结果的P值作判断:引素A时段和B路段对行车时间有显著影响;而交互作用A:B的P值=0.42>0.05 ,因此不能拒绝原假设H0,说明两个因素间没有明显的交互效应。

协方差分析及R实现 为了提高试验的精确性和准确性,我们对除研究因素以外的一切条件都需要采取有效措施严加控制,使它们在因素的不同水平间尽量保持一致,这叫做试验控制。但当我们进行试验设计时,即使做出很大努力控制,也经常会碰到试验个体的初始条件不同的情况,如果不考虑这些因素有可能导致结果失真。如果考虑这些不可控的因素,这种方差分析就叫做协方差分析,其是将回归分析和方差分析结合在一起的方法。它的基本原理如下:将一些对响应变量Y有影响的变量X(未知或难以控制的因素)看作协变量,建立响应变量Y随X变化的线性回归分析,从Y的总的平方和中扣除X对Y的回归平方和,对残差平方和作进一步分解后再进行方差分析。 例: 施用3种肥料的苹果产量

> Weight_Initial=c(15,13,11,12,12,16,14,17,17,16,18,18,21,22,19,18,22,24,20,23,25,27,30,32)
> Weight_Increment=c(85,83,65,76,80,91,84,90,97,90,100,95,103,106,99,94,89,91,83,95,100,102,105,110)
> feed=gl(3,8,24)
> data_feed=data.frame(Weight_Initial,Weight_Increment,feed)
> library(HH)
> m=ancova(Weight_Increment~Weight_Initial+feed,data=data_feed)
> summary(m)
Df Sum Sq Mean Sq F value Pr(>F)
Weight_Initial 1 1621.1 1621.1 142.44 1.50e-10
feed 2 707.2 353.6 31.07 7.32e-07
Residuals 20 227.6 11.4
 
Weight_Initial ***
feed ***
Residuals
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

协方差分析的P值非常小,说明结果非常显著,应该拒绝原假设,认为各因素在不同水平下的试验结果有显著差别,即三种肥料对苹果产量有很大的影响。

原文发布于微信公众号 - 大数据挖掘DT数据分析(datadw)

原文发表时间:2015-08-18

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏CDA数据分析师

测试:你是否具备企业的数据挖掘能力?

1.某超市研究销售纪录数据后发现,买啤酒的人很大概率也会购买尿布,这种属于数据挖掘的哪类问题? A. 关联规则发现 B. 聚类 C. 分类 D. 自然语言处理 ...

1836
来自专栏大数据文摘

机器学习单挑数学界:最新算法仲裁数列之美(附论文)

它揭示了表面看似无关的数学领域之间的深层联系,是数学界的伟大奇观之一。而这也指出了数学之美的另一个组成部分:数学模式必须在某种角度上是有趣的。

774
来自专栏机器之心

教程 | 如何从TensorFlow转入PyTorch

选自Medium 作者:Illarion Khlestov 机器之心编译 参与:李泽南、黄小天 当我第一次尝试学习 PyTorch 时,没几天就放弃了。和 Te...

46115
来自专栏ATYUN订阅号

如何为地图数据使用tSNE聚类

在本文中,我会展示如何在经纬度坐标对上使用tSNE来创建地图数据的一维表示。这种表示有助于开发新的地图搜索算法。这对于诸如“这个经纬度坐标是新泽西或者纽约的吗?...

773
来自专栏人工智能头条

深度学习一种变相的马尔可夫链吗?

1684
来自专栏AI研习社

博客 | Word2Vec 学习心得

好嘛博主食言了。不过本文没什么干货,主要是前后看了大概一个星期,反复去读源码和解读文章,终于感觉这东西不那么云山雾罩了。同时也发现网上很多材料有点扯淡,99% ...

612
来自专栏杨熹的专栏

机器学习&人工智能博文链接汇总

? 争取每天更新 ? 126 ? ---- 蜗牛的历程: [入门问题] [机器学习] [聊天机器人] [好玩儿的人工智能应用实例] [Tensor...

3406
来自专栏智能算法

机器人算法专题介绍

算法 算法(Algorithm)是指解题方案的准确而完整的描述,是一系列解决问题的清晰指令,算法代表着用系统的方法描述解决问题的策略机制。也就是说,能够对一定规...

3336
来自专栏机器之心

谷歌微软等科技巨头数据科学面试107道真题:你能答出多少?

选自Learndatasci 机器之心编译 参与:李泽南 来自 Glassdoor 的最新数据可以告诉我们各大科技公司最近在招聘面试时最喜欢向候选人提什么问题。...

2517
来自专栏懒人开发

(8.3)James Stewart Calculus 5th Edition:Applications to Physics and Engineering

any point in a liquid the pressure is the same in all directions 液体中任何点在任何方向收到的...

764

扫描关注云+社区