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

本文分享自微信公众号 - 机器学习与统计学(tjxj666)

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2015-05-27

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏生信小驿站

R语言之可视化(21)令人眼前一亮的颜色包

26920
来自专栏生信小驿站

R语言之可视化(22)绘制堆积条形图

经过这张图,我们可以初步得到的信息是:(1)T1到T4各个分期的患者总数(2)T1期男性患者的数目,T1女性患者的数目(3)其他分期男性或者女性的患者数目。

26920
来自专栏nummy

奇怪的编码问题

今天使用R爬取数据的时候发现一个奇怪的问题,我将每个属性的数据先保存在vector中,然后再合并到data.frame中时,发现打印names时数据正常显示中文...

8930
来自专栏生信小驿站

R语言第二章数据处理(9)数据合并

=========================================

16420
来自专栏机器之心

Python与PHP的对决:谁是工程师最喜欢和最讨厌的语言?

为了弄清楚雇主对哪些编程技能最感兴趣,Hired 研究了求职者在到六周内收到的面试邀请数量。如下图显示,谷歌的 Go 语言是雇主最需要的编程语言技能,可能因为这...

9830
来自专栏Python数据科学

经验分享 | 如何写好数据分析师简历?

我们要确定怎么样简历是一份好数据分析师简历呢?那我们就要涉及到如何评价一个好数据分析师?一般来说,优秀的数据分析师有着很好的表达能力,能通过在二分钟对自己工作内...

70630
来自专栏yw的数据分析

R语言读入数据库的中英名词互译测试并计分脚本(考试用)

    1. 分子生物学中英文.csv,输入文件,两列,以tab键分隔的txt文本,没有列名

6710
来自专栏零维领域

线性代数--MIT18.06(三)

,我们依然可以使用矩阵消元的形式来求解,只不过要比我们之前提到的矩阵消元多做一些消元而已,这就是Gauss-Jordan法。

11020
来自专栏Rude3Knife的后端开发专栏

【C#】SM2C多云平台安全数据库应用

下文为论文中的程序实现进行了英文描述,在此做备份,该应用上传在Github上https://github.com/qqxx6661/SMC_Yang,欢迎for...

9730
来自专栏生信小驿站

R语言之可视化(23)高亮某一元素

总结:假如需要高亮ggplot2中的某一元素时,首先需要新建一列,然后修改新建列中需要高亮的部分即可

11230

扫码关注云+社区

领取腾讯云代金券

年度创作总结 领取年终奖励