方差分析

"克服懒惰，坚持更新！"

提到方差分析（Analysis of Variance），简写为ANOVA，相信只要接触过统计学或者有过科研经历的小伙伴们对此不会陌生。之前我更多的是使用SPSS来操作，那么怎么用R语言来实现呢？
首先，我们先来看一下方差分析的前提假设：

样本数据独立
每组数据的总体服从正态分布
每组数据方差齐性

我的第一篇博客介绍了T检验，其前提假设也是以上三条，事实上，二者在某些情况下是等价的，比如独立样本T检验（两组数据）和方差分析。下面重点介绍单因素方差分析和多因素方差分析。

单因素方差分析

先来看一下数据，这里选择R自带的数据集PlantGrowth，weight是植物的重量，group是不同处理，包括(trt1，trt2)和空白对照组(ctrl)。

> My_data <- PlantGrowth
> head(My_data)weight group
1   4.17  ctrl
2   5.58  ctrl
3   5.18  ctrl
4   6.11  ctrl
5   4.50  ctrl
6   4.61  ctrl

老规矩，在进行方差分析之前先要进行正态分布检验和方差齐性检验。首先是正态分布检验，把数据分为两列，构建函数，对3个水平下的每组数据都做一次正态分布检验。


> #正态分布检验
> My_data1 <- split(My_data[,1], My_data[,2])
> unlist(lapply(My_data1, function(x){shapiro.test(x)$p.value}))ctrl      trt1      trt2
0.7474734 0.4519440 0.5642519

结果很明显，P值全部大于0.05，满足正态分布，可以进行下一步了！但是代码中的lapply，function又是什么东西？？？

除了lapply，还有apply、tapply、sapply、mapply，这些都是一个家族的，功能也都相似，他们都是向量化函数。至于function，就是构建或使用的函数。

简单介绍一下吧。比如apply()，举个栗子：

> M <- matrix(c(1:3, 4:6), nrow = 2)#设置一个矩阵M
> M[,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> apply(M,1,sum)
#对一个数组按行或者按列进行加和，1代表行，2代表列
[1]  9 12#得到一个矩阵或向量

如果你想加上function，其实多此一举而已：


> apply(M,1,function(x){sum(x)})
[1]  9 12

lapply(x，function …)
lapply的返回值是和一个和x有相同的长度的list对象，这个list对象中的每个元素是将函数function应用到X的每一个元素。
在这个”普拉爱”家族里，sapply()可能是使用最为频繁的了，如果把lapply换成sapply，就不需要unlist()了，其输出格式是较为友好的矩阵格式。

My_data1 <- split(My_data[,1], My_data[,2])
sapply(My_data1, function(x){shapiro.test(x)$p.value
})ctrl      trt1      trt2
0.7474734 0.4519440 0.5642519

啰嗦了这么多，下面开始方差齐性检验！

> library(car)
> My_data <- PlantGrowth
> leveneTest(weight~group, data = My_data)
Levene's Test for Homogeneity of Variance (center = median)Df F value Pr(>F)
group  2  1.1192 0.341227

P值为0.3412>0.05，顺利通关！

接下来就是方差分析了，稍等，如果你还是不放心的话，可以再检查一下是否有离群点，只需要一行代码：

> outlierTest(lm(weight~group, data=My_data))
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:rstudent unadjusted p-value Bonferroni p
17 2.537341           0.017509      0.52526

恩，没有离群点。

终于可以开始方差分析了！一行代码而已~~~

> summary(aov(weight~group, data = My_data))Df Sum Sq Mean Sq F value Pr(>F)
group        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

答案是一个小星星，说明p<0.05，意味着不同处理之间植物的重量有显著差异。一个不错的结果~

双因素方差分析

双因素方差分析分为无交互作用和有交互作用两种。
无交互作用的双因素与单因素的区别就是多了一个因素。举栗子，先看数据，这里选用的是R自带的数据集ToothGrowth。下面的分析就不再做正态分布检验和方差齐性检验了，假设都通过了~

> My_data <- ToothGrowth
> head(My_data)#两个自变量supp和dose，一个因变量lenlen supp dose
1  4.2   VC  0.5
2 11.5   VC  0.5
3  7.3   VC  0.5
4  5.8   VC  0.5
5  6.4   VC  0.5
6 10.0   VC  0.5

相信聪明的你已经猜到代码是什么了：

> summary(aov(len~supp+dose, data = My_data))Df Sum Sq Mean Sq F value   Pr(>F)
supp         1  205.4   205.4   11.45   0.0013 **
dose         1 2224.3  2224.3  123.99 6.31e-16 ***
Residuals   57 1022.6    17.9
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

结果一目了然，就不做解释了。

再来看考虑交互作用的方差分析怎么做：

> summary(aov(len ~ supp * dose, data = My_data))Df Sum Sq Mean Sq F value   Pr(>F)
supp         1  205.4   205.4  12.317 0.000894 ***
dose         1 2224.3  2224.3 133.415  < 2e-16 ***
supp:dose    1   88.9    88.9   5.333 0.024631 *
Residuals   56  933.6    16.7
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

没想到竟然这么简单，就是把+换成了*。如果你以为大功告成了，那么恭喜你做错了。为什么？？？
来看看数据结构：

> str(My_data)
'data.frame': 60 obs. of  3 variables:$ len : num  4.2 11.5 7.3 5.8 6.4 10 11.2 11.2 5.2 7 ...$ supp: Factor w/ 2 levels "OJ","VC": 2 2 2 2 2 2 2 2 2 2 ...$ dose: num  0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 ...

dose不是Factor，竟然是num…，而我们分析的是什么，是双因素！
所以重新开始吧，首先要把dose转换成factor：

> My_data$dose <- factor(My_data$dose)
> summary(aov(len~supp+dose, data = My_data))Df Sum Sq Mean Sq F value   Pr(>F)
supp         1  205.4   205.4   14.02 0.000429 ***
dose         2 2426.4  1213.2   82.81  < 2e-16 ***
Residuals   56  820.4    14.7
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> summary(aov(len~supp*dose, data = My_data))Df Sum Sq Mean Sq F value   Pr(>F)
supp         1  205.4   205.4  15.572 0.000231 ***
dose         2 2426.4  1213.2  92.000  < 2e-16 ***
supp:dose    2  108.3    54.2   4.107 0.021860 *
Residuals   54  712.1    13.2
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

可以与上边的两个结果对比一下，不一样吧。结果不解释了。一定不要忘记加载数据之后要先看一下数据的结构str()。

学会了双因素，多因素也就可以照葫芦画瓢了。除了这两种之外，还有方差不齐性的数据（并且每组样本数也不相同），那么可以使用oneway.test()。

> oneway.test(weight~group, data = My_data, var.equal = F)One-way analysis of means (not assuming equal variances)data:  weight and group
F = 5.181, num df = 2.000, denom df = 17.128, p-value = 0.01739

可以和上面的对比一下，P值一个是0.0159，一个是0.01739。这说明在方差齐性和样本数相同的情况下结果差别不大，否则可能就很大了，在此就不验证了。
另外，还有一种是有协方差的数据，代码结构如下

summary(aov(因变量~协方差+自变量,data = 数据集名称))
#协方差在前，自变量在后

—— END

坚持就是胜利~~

菜鸟学R语言（方差分析）相关推荐

菜鸟学R语言（PCA）
主成分分析(PCA) 理论介绍 PCA的全拼是Principal Component Analysis,翻译成中文就是主成分分析.在探索性数据分析中应用广泛,能更好地可视化在具有许多变量的数据集中出现 ...
菜鸟学R语言（回归分析1）
学会用R做回归分析在统计学中,回归分析(regression analysis)指的是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法.回归分析按照涉及的变量的多少,分为一元回归和多元回 ...
“跟着菜鸟一起学R语言” 现已更名为“数据志”
大家好,我的公众号"跟着菜鸟一起学R语言" 现已更名为"数据志",欢迎大家关注,谢谢.
r语言中矩阵QR分解_从零开始学R语言Day4|向量、矩阵和数组
从零开始学R语言Day4|向量.矩阵和数组 1.1向量 1.1.1向量在Day2中我们提及过用和c()函数来构建向量,具体实例如下. 我们还可以采用vector("类型",长度) ...
菜鸟学C语言-环境搭建
菜鸟学C语言-环境搭建本人由于工作需要,需要用C写一套程序,虽然以前看过一两本C的书,但是时隔几年早已忘得一干二净,所以相当于是需要从头开始学习,在此记录一下我得学习成果环境所需工具 window ...
R语言方差分析的注意事项
本文首发于公众号:医学和生信笔记,完美观看体验请至公众号查看本文. 文章目录均衡设计和非均衡设计方差分析的3种类型示例 one-way anova two-way anova 协方差分析 R语言 ...
入门必学 | R语言数值型、字符型及因子型数据之间的差异与转换
字符型.数值型及因子型数据之间的转换数据类型的基本知识不同数据类型之间的差异数值型与字符型或因子型绘图时的差异数值型与因子型和字符型的模型构建时的差异三种数据类型之间进行转换常用的三 ...
一个菜鸟学习R语言的历程（一）
菜鸟也有梦想! 写在前面:本人的确是一只菜鸟,大学没学过编程,对于C语言和java也只是耳闻,R语言和Python是大学毕业后才知道(卑微).本科只用过SPSS和EXCEL分析数据,目前马上研二,才认 ...
R语言方差分析ANOVA
自己整理编写的R语言常用数据分析模型的模板,原文件为Rmd格式,直接复制粘贴过来,作为个人学习笔记保存和分享.部分参考薛毅的<统计建模与R软件>和<R语言实战> I. 单因素方 ...
入门必学 | R语言数据的独立性，正态性及方差齐性检验
参数分析的三大前提检验检验数据独立性的方法 Chisq检验 Fisher检验 Cochran-Mantel-Haenszel检验检验数据正态性的方法 shapiro.test函数 qqnorm函数 ...

菜鸟学R语言（方差分析）

方差分析

单因素方差分析

双因素方差分析

坚持就是胜利~~

菜鸟学R语言（方差分析）相关推荐

最新文章

热门文章