基于R实现统计中的检验方法---方差分析

栏目: R语言 · 发布时间: 5年前

作者：徐涛，19年应届毕业生，专注于珊瑚礁研究，喜欢用R各种清洗数据。

知乎：

https://www.zhihu.com/people/parkson-19/posts

前言

方差分析（均数的显著性检验）

方差分析(Analysis of Variance，简称ANOVA)，又称“变异系数分析”，是R.A.Fisher发明的，用于两个及两个以上样本均数差别的显著性检验。

1.适用条件

方差分析数据必须满足3个条件即独立、正态，方差齐次性。

---------------------#正态性检验#---------------------
#Shapiro-Wilk检验用来检验数据是否符合正态分布，类似于线性回归的方法一样，是检验其于回归曲线的残差。
#该方法作者推荐在样本量很小的时候使用，比如N<20。但是也有作者推荐在大数据集上使用。该作者将这种修改
#后的方法运用在R语言的stats包中的shapiro.test函数中。该检验原假设为H0：数据集符合正态分布。

#R语言中的Shapiro-Wilk检验（W检验），shapiro.test(x)该函数只有一个参数即数据集x。x可以是数值型向量，允许存在NA
#但是非丢失数据需要在3-5000内。
##example##
#11个随机抽取的样本的体重(斤)数据为：148,154,158,160,161,162,166,170,182,195,236
k<-c(148,154,158,160,161,162,166,170,182,195,236)
shapiro.test(k)
        Shapiro-Wilk normality test

data:  k
W = 0.78881, p-value = 0.006704
##p=0.0067<-0.05，拒绝原假设，该数据不符合正态分布。


---------------------#方差齐性检验#---------------------
#R中最常用的方差齐性检验（bartlett检验），bartlett.test()调用格式为bartlett.test()调用格式为
#bartlett.test(x,g...)其中，参数x是数据向量或列表(list);g是因子向量，如果x是列表则忽略g，当使用数据集时，也通过formula调用函数
##example##
#以multcomp包中的cholesterol数据集为例
library(multcomp)
head(cholesterol)
   trt response
1 1time   3.8612
2 1time  10.3868
3 1time   5.9059
4 1time   3.0609
5 1time   7.7204
6 1time   2.7139
bartlett.test(response~trt,data=cholesterol)
 Bartlett test of homogeneity of variances

data:  response by trt
Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653
#Bartlett检验表明五组的方差并没有显著不同（p=0.97）

2.分类（按因素分）

单因素方差分析、二因素方差分析和多因素方差分析。（单因素方差分析主要对因素的主效应进行分析，而二因素及多因素方差分析除了多组效应分析外，还可对因素间的交互作用进行分析。分析交互作用时，观测值必须有重复）

3.R实例

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

aov(formula,data=NULL,projections=FALSE,qr=TRUE,contrasts=NULL,...)
#其中的参数formula表示方差分析的公式，在单因素方差分析中即为X~A；data表示做方差分析的数据框；
projections为逻辑值，表示是否返回预测结果；qr同样是逻辑值，表示是否返回QR分解结果，默认为TRUE；
contrasts是公式中的一些因子的对比列表。通过summary()可列出分析表的详细结果。
#详情参见R语言实战（第2版）page:201.

单因素方差分析

---------------------#单因素方差分析#---------------------
#单因素方差分析比较的是分类因子定义的两个或多个组别中的因变量均值。以multcomp包中cholesterol数据集为例。

##正态性检验
library(multcomp)
head(cholesterol)
str(cholesterol)
shapiro.test(cholesterol$response)
      Shapiro-Wilk normality test

data:  cholesterol$response
W = 0.97722, p-value = 0.4417
#结果表明，p=0.4417>0.05，接受原假设，说明数据在五个水平下都是正态分布的。

##方差齐性检验
bartlett.test(response~trt,data=cholesterol)
 Bartlett test of homogeneity of variances

data:  response by trt
Bartlett's K-squared = 0.57975, df = 4, p-value = 0.9653
#结果表明，p=0.9653，接受原假设，说明数据在不同水平下是等方差的。

##单因素方差分析
fit<-aov(response~trt,data=cholesterol)
summary(fit)
            Df Sum Sq Mean Sq F value   Pr(>F)    
trt          4 1351.4   337.8   32.43 9.82e-13 ***
Residuals   45  468.8    10.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#结果表明，5种治疗方式间的效果差异非常显著。

##多重比较（多重T检验）
#以上单因素方差分析可以得出，五组之间治疗效果有显著差异，但是我们不知道哪种疗法与其他疗法的不同，多重比较可以解决这个问题。
#R语言中TukeyHSD()函数提供了对各组均值差异的成对检验。
TukeyHSD(fit)
  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = response ~ trt, data = cholesterol)

$trt
                  diff        lwr       upr     p adj
2times-1time   3.44300 -0.6582817  7.544282 0.1380949
4times-1time   6.59281  2.4915283 10.694092 0.0003542
drugD-1time    9.57920  5.4779183 13.680482 0.0000003
drugE-1time   15.16555 11.0642683 19.266832 0.0000000
4times-2times  3.14981 -0.9514717  7.251092 0.2050382
drugD-2times   6.13620  2.0349183 10.237482 0.0009611
drugE-2times  11.72255  7.6212683 15.823832 0.0000000
drugD-4times   2.98639 -1.1148917  7.087672 0.2512446
drugE-4times   8.57274  4.4714583 12.674022 0.0000037
drugE-drugD    5.58635  1.4850683  9.687632 0.0030633
#结果表明，1time和2times的均值差异不显著(p=0.138)，而1time和4times间的差异非常显著(p<0.01)

一般图像给的效果会更加直观，R中multcomp包中的glht()函数提供了多重均值比较更为全面的方法。

par(mar=c(5,4,6,2))
tuk<-glht(fit,linfct=mcp(trt="Tukey"))
plot(cld(tuk,level=0.05),col="lightgreen")##cld函数中的level选项设置了使用的显著水平0.05.
#图像中，有相同字母的组说明均值差异不显著。

基于R实现统计中的检验方法---方差分析

备注：

1）在单因素方差分析前需要对数据进行正态性和方差齐性检验。除了以上提及的检验方法，R中还可以使用Q-Q图来检验正态性假设。#详情参见R语言实战（第2版）page:206-207.

2）单因素方差分析的扩展为单因素协方差分析，包含一个或多个定量的协变量。（详情参见R语言实战（第2版）page:208.）

二因素方差分析（双因素方差分析）

---------------------#二因素方差分析#---------------------
#以基础安装包中的ToothGrowth数据集为例，60只豚鼠，自变量包括两种投喂方法（橙汁或VC）,不同喂食方法中抗坏血酸含量有
#三种水平，牙齿长度为因变量，每种处理都被分配10只豚鼠。
head(ToothGrowth)
   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
ToothGrowth$dose<-as.factor(ToothGrowth$dose)#将剂量改为factor，而不是数值型。

#正态性检验#
shapiro.test(ToothGrowth$len)
       Shapiro-Wilk normality test

data:  ToothGrowth$len
W = 0.96743, p-value = 0.1091
#结果显示，p>0.05，接受原假设，数据正态分布。

#方差齐性检验#
#需要检验方差在两个不同分类标准下方差是否等值。
#检验牙齿长度在不同喂食方法水平下方差是否等值
bartlett.test(len~supp,data=ToothGrowth)
        Bartlett test of homogeneity of variances

data:  len by supp
Bartlett's K-squared = 1.4217, df = 1, p-value = 0.2331
#结果显示，p>0.05，接受原假设，在不同喂食方法水平上，牙齿长度方差等值。

#检验牙齿长度在不同剂量水平下方差是否等值
bartlett.test(len~dose,data=ToothGrowth)
        Bartlett test of homogeneity of variances

data:  len by dose
Bartlett's K-squared = 0.66547, df = 2, p-value = 0.717
#结果显示，p>0.717，接受原假设，在不同剂量水平上，牙齿长度方差等值。

#双因素方差分析#
fitt<-aov(len~supp*dose,data=ToothGrowth)
summary(fitt)
            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
#结果显示，不同喂食方法以及不同剂量对豚鼠牙齿生长影响显著，且喂食方法和不同剂量之间的交互作用非常显著。