方差分析
"克服懒惰，坚持更新！"
提到方差分析（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.3412
27
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，一个因变量len
len 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
扫一扫
1232
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
