使用R语言做单因素方差和双因素方差分析

1.方差分析目的:

方差分析(ANOVA)主要检验不同处理间是否有显著差异，当仅有一个类别型变量(即只有一类处理)，称为单因素方差分析(one-way ANOVA)，当有两个类别型变量(即有两类处理),称为双因素方差分析(two-way ANOVA)，依次类推。单因素方差分析只能告诉处理间是否有显著差异，需要结合多重比较(一般使用TukeyHSD函数)显示哪两种处理水平有显著差异。双因素方差分析能告诉两类处理对响应变量是否有显著影响和交互影响，若无交互作用，两因素方差分析结果近似于分开做两次单因素方差分析，即只考虑两类处理的差异，可以只做单因素方差分析（可能误差很大），若有交互作用，必须做双因素方差分析显示其交互作用。

2.数据格式:

以基础安装包中ToothGrowth数据集为例,随机分配60只豚鼠，采用两种喂食方法(supp,包括橙汁: OJ和维生素C: VC), 各喂食方法有3种抗坏血酸水平（dose: 0.5; 1, 2 mg）,每种处理水平10只豚鼠，即10个重复，研究两种处理和3个水平对豚鼠牙齿长度的影响（len）。逐行运行代码查看数据及数据格式。

> attach(ToothGrowth)     ##加载数据  

> ToothGrowth               ##注意数据格式 VC和OJ下3个水平是一致的，第一列是R自动生成的编号，不是整理好的数据。

3.方差分析前提条件:

方差分析数据需符合正态分布，各组方差相等。这里使用前三十行数据进行方差齐性检验，命名为ToothGrowth1，单因素方差分析同样使用此数据集。两因素方差分析使用全部数据，即ToothGrowth。逐行运行代码。

> ToothGrowth1 <-ToothGrowth[1:30,1:3]   ##截取前30行数据

> ToothGrowth1        ##ToothGrowth1[,3] <- as.factor(ToothGrowth1[,3])###因为dose是数字需要变成因子变量，单独运行此代码，查看数据格式。

Q-Q图发现数据落在95%的置信区间内，说明满足正态性假设，方差齐性检验同样发现3组水平方差无显著不同(P=0.1048),但方差齐性检验对离群点比较敏感，可用car包outlierTest（）函数来检测离群点。这里不做举例，我目前做方差分析一般Q-Q图和bartlett.test（）检验后符合的话后面还没有出现不满足的情况。

4.代码和说明:

> library(car)     # 加载car包

> qqPlot(lm(len~dose,data=ToothGrowth1),simulate = TRUE,main = "Q-Q PLot", lables=FALSE)

> bartlett.test(len~dose,data=ToothGrowth1)   ##方差齐性检验

Bartlett test of homogeneity of variances

data:  len by dose

Bartlett's K-squared = 4.5119, df = 2, p-value = 0.1048   ##方差齐性检验发现各组方差没有显著不同，即P=0.1048，大于0.05.但是齐性检验对离群点比较敏感，可以利用car包中outlierTest（）函数来检测离群点。

      我一般是用HH包的interaction2wt（）函数，因为我交互和主效应表只放在组会的PPT里面，好看，HH包比较早，可能不好下载，需要手动安装，我也上传保存。其实也有好多其他作图方式，比如interaction.plot() 函数和gplots包中的plotmeans（）函数。我一般用R做统计分析，sigmplot作图，R 作图也在学习当中。主效应和交互效应如何理解是最开始做双因素方差分析时困扰我的问题，一般是若无交互效应，只看主效应，若有交互，不能单纯看主效应，需要结合图像来理解，或做简单效应分析。这个我单独做笔记。

HH_3.0-4.zip

5.代码汇总

attach(ToothGrowth)

ToothGrowth

ToothGrowth1 <-ToothGrowth[1:30,1:3]

fit<- aov(len~dose,data=ToothGrowth1)

ToothGrowth1[,3] <- as.factor(ToothGrowth1[,3])

ToothGrowth1

library(car)

qqPlot(lm(len~dose,data=ToothGrowth1),simulate = TRUE,main = "Q-Q PLot", lables=FALSE)

bartlett.test(len~dose,data=ToothGrowth1)

fit<- aov(len~dose,data=ToothGrowth1)

summary(fit)

TukeyHSD(fit)

fit<- aov(len~supp*dose,data=ToothGrowth)

summary(fit)

library(HH)

interaction2wt(len~supp*dose)

单因素和双因素方差分析代码，只需要把数据框的名字和列的名字改动就可以直接运行。需要注意TukeyHSD() 函数与HH包不兼容，加载HH包后，TukeyHSD()函数会自动失效。可以使用detach（）函数将HH从路径中暂时去除。我一般都是单因素方差做完之后做双因素方差分析。方差分析简单实用，原理随意参考一本生物统计书，相关操作主要参考了R语言实战第9章。统计学习，代码留存，以防下次电脑中毒再丢失。