R语言数据分析与挖掘(第五章):方差分析(2)——多因素方差分析
R语言数据分析与挖掘(第五章):方差分析(2)——多因素方差分析
在实际应用中，更多出现的是包含多因素的试验和处理。多因素试验与双因素试验背后的基本思想是一致的。与单因素方差分析不同，在双因素方差分析中因素间可能会有交互作用。假设有两个因素A和B，因素A和B没有交互作用指的是A的水平值不取决于B的水平值，反之亦然。对于有交互作用的因素，我们不可孤立地看待这些因素。对于双因素的情形，一般从图像上看，没有交互作用的因素水平图表现为两条不相交的线段，而有交互作用的因素水平图为两相交的线段。例如，下图显示的是在研究年龄和性别对身高是否有显著作用过程中，因素年龄与性别之间的交互作用。从图像上看，两曲线没有明显相交，据此可以推测二者间不存在相互作用。当然，要判定是否存在或者不存在交互作用，还需要根据相应的统计量来分析。
如果因素A和B间不存在交互作用，那么可以对因素A和B各自进行独立分析，在后续的分析中去除不显著的因素。如果方差分析的结果显示因素A和B间存在交互作用，这时候要对数据进行进一步的分析，具体包括：
在因素A的某个水平下，因素B对响应变量的作用。
在因素B的某个水平下，因素A对响应变量的作用。
在所有因素（A，B）的组合中，哪两组的差异最大。
需要指出的是，多因素的情况与双因素处理方法类似，只不过分析的因素会增多。
下面以不同剂量下老鼠妊娠重量的差异性分析为例进行介绍。
在R语言中，实现双因素方差分析的函数与实现单因素方差分析的函数一致，可以实现aov()和anova()函数，不同之处在于模型公式的设定，双因素方差分析的模型公式应设定为"X~A+B"或"X~A*B"的形式，后者表示考虑因素A与B的交互作用。
下面利用包multcomp中的litter数据集进行操作演练，该数据集是关于老鼠妊娠重量与剂量反应的研究数据，共74个行观测值和4个变量，列变量包括，剂量等级，共有4个等级，0,5,50和500，妊娠时间，作为协变量共有4个水平，分别是21.5、22、22.5和23；妊娠动物的个数作为协变量；整个妊娠过程中出生的小鼠平均重量，数据的基本情况如下：
> library(multcomp)
> data(litter)
> dim(litter)
[1] 74 4
> head(litter)
dose weight gesttime number
1 0 28.05 22.5 15
2 0 33.33 22.5 14
3 0 36.37 22.0 14
4 0 35.52 22.0 13
5 0 36.77 21.5 15
6 0 29.60 23.0 5
> summary(litter)
dose weight gesttime number
0 :20 Min. :19.22 Min. :21.50 Min. : 5.00
5 :19 1st Qu.:27.77 1st Qu.:21.50 1st Qu.:12.00
50 :18 Median :30.76 Median :22.00 Median :14.00
500:17 Mean :30.33 Mean :22.09 Mean :13.43
3rd Qu.:33.30 3rd Qu.:22.50 3rd Qu.:15.00
Max. :38.75 Max. :23.00 Max. :17.00
> table(litter$gesttime)
21.5 22 22.5 23
20 24 27 3
> aggregate(litter$weight,by=list(litter$dose),FUN=mean)
Group.1 x
1 0 32.30850
2 5 29.30842
3 50 29.86611
4 500 29.64647
> aggregate(litter$weight,by=list(litter$dose),FUN=sd)
Group.1 x
1 0 2.695119
2 5 5.092352
3 50 3.762529
4 500 5.404372
> par(mfrow=c(1,2))
> boxplot(weight~dose,data=litter,col=2:5,main="按变量dose分组")
> boxplot(weight~gesttime,data=litter,col=2:5,main="按变量gesttime分组")上述代码表示:利用函数summary()对数据集进行描述性统计，发现变量dose的4个不同水平的数量不相等，分别为20、19、18和17;需要注意的是，原始数据中变量gesttime中存储的数据类型为数值型，故利用函数table()统计该变量4个水平的数据量，发现妊娠时间为21.5、22、22.5和23的老鼠数量分别为20、24、27和30;利用函数aggregate()对变量weight进行分组统计，并计算每-组的均值和方差， 分组依据为变量dose的4个水平，并根据两个协变量的不同分组绘制变量weight的箱线图，结果下图所示。
下面对数据进行方差齐性检验，由于原始变量gesttime的数据类型为数值型，在进行方差齐性检验的时候需要将其转化成因子，具体操作如下：
> bartlett.test(weight~as.factor(gesttime),data=litter)
Bartlett test of homogeneity of variances
data: weight by as.factor(gesttime)
Bartlett's K-squared = 0.26778, df = 3, p-value = 0.966
> bartlett.test(weight~dose,data=litter)
Bartlett test of homogeneity of variances
data: weight by dose
Bartlett's K-squared = 9.6497, df = 3, p-value = 0.02179输出结果显示：变量gesttime的值为0.966，大于给定的显著性水平0.05，故不能拒绝原假设，即认为不同水平下是等方差的，然而dose的p值为0.02179，没有通过检验。接下来对数据进行方差分析，分别按照不考虑gesttime和dose的交互和考虑其交互进行分析。
> fit1<-aov(weight~gesttime+dose,data=litter)
> summary(fit1)
Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.049 0.00597 **
dose 3 137.1 45.71 2.739 0.04988 *
Residuals 69 1151.3 16.69
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1在不考虑交互作用的情况下，协变量gesttime和dose的检验p值均小于给定的显著水平，故拒绝原假设，即认为变量gesttime和dose对变量weight的影响均显著的。
> fit2<-aov(weight~gesttime*dose,data=litter)
> summary(fit2)
Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.289 0.00537 **
dose 3 137.1 45.71 2.821 0.04556 *
gesttime:dose 3 81.9 27.29 1.684 0.17889
Residuals 66 1069.4 16.20
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1上述代码表示:在考虑交互作用的情况下进行方差分析，结果显示，交互作用项的p值为0.17889，大于给定的显著性水平0.05,故不认为交互项对交量weigh的影响是显著的，该交互项可以删除。
交互作用的效果还可以进行可视化展示图，利用程辑包HH中的函数intrecion2wt()即可，其函数的用法与函数aov()类似，具体操作代码如下:
install.packages("HH")
library(HH)
interaction2wt(weight~gesttime*dose,data=litter)结果如下图所示。主要看图形的第一、四象限的曲线是否存在明显相交的情况，若存在，则说明两因素间的交互作用显著，否则认为不显著，本图中第一、四象限的曲线有一定的相交，说明在后续的方差分析中需要添加交互项，但是根据前面的分析结果，交互项并没有通过检验，即交互项的作用并不明显。
本文分享自 MedBioInfoCloud 微信公众号，前往查看
如有侵权，请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体分享计划 ，欢迎热爱写作的你一起参与！