【R<-练习】双因子方差分析
对于两因素的方差分析, 基本思想和方法与单因素的方差分析相似, 前提 条件仍然是要满足独立、正态、方差齐性. 所不同的是在双因素方差分析中, 有时会出现交互作用, 即二因素的不同水平交叉搭配对指标产生影响. 我们先 讨论无交互作用的双因素方差分析.
1 无交互作用的方差分析
在R软件中, 方差分析函数aov( )既适合于单因素方差分析, 也同样适用 于双因素方差分析, 其中方差模型公式为x~A+B, 加号表示两个因素具有 可加的. 下面用一个例子来说明.
例 题1
原来检验果汁中含铅量有三种方法A1、A2、A3, 现研究出另 一种快速检验法A4, 能否用A4代替前三种方法, 需要通过实验考察. 观察的对 象是果汁, 不同的果汁当做不同的水平: B1为苹果, B2为葡萄汁, B3为西红柿 汁, B4为苹果饮料汁, B5桔子汁, B6菠萝柠檬汁. 现进行双因素交错搭配试验, 即用四种方法同时检验每一种果汁, 其检验结果如表8.4所示. 问因素A(检验方 法)和B(果汁品种) 对果汁的含铅量是否有显著影响?
统计表
解答：
# 将数据导入到命名为juice的数据框中
juice<-data.frame( X = c(0.05, 0.46, 0.12, 0.16, 0.84, 1.30, 0.08, 0.38, 0.4, 0.10, 0.92, 1.57, 0.11, 0.43, 0.05, 0.10, 0.94, 1.10, 0.11, 0.44, 0.08, 0.03, 0.93, 1.15), A = gl(4, 6), B = gl(6, 1, 24) )
> juice
X A B
1 0.05 1 1
2 0.46 1 2
3 0.12 1 3
4 0.16 1 4
5 0.84 1 5
6 1.30 1 6
7 0.08 2 1
8 0.38 2 2
9 0.40 2 3
10 0.10 2 4
11 0.92 2 5
12 1.57 2 6
13 0.11 3 1
14 0.43 3 2
15 0.05 3 3
16 0.10 3 4
17 0.94 3 5
18 1.10 3 6
19 0.11 4 1
20 0.44 4 2
21 0.08 4 3
22 0.03 4 4
23 0.93 4 5
24 1.15 4 6
注: 这里函数gl( )用来给出因子水平
其调用格式为 gl(n, k, length=n*k, labels=1:n, ordered=FALSA)
说明: n是水平数, k是每一水平上的重复次数, length是总观测值数, ordered指明各水平是否先排序.
# 接下来作双因子方差分析
> juice.aov <- aov(X~A+B, data = juice)
> summary(juice.aov)
Df Sum Sq Mean Sq F value Pr(>F)
A 3 0.057 0.0190 1.629 0.225
B 5 4.902 0.9804 83.976 2e-10 ***
Residuals 15 0.175 0.0117
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
结论: p值说明果汁品种(因素B)对含铅量有显著影响, 而没有充分理由说明检 验方法(因素A)对含铅量有显著影响.
# 最后用函数bartlett.test( )分别对因素A和因素B作方差的齐性检验:
> bartlett.test(X~A,data = juice)
Bartlett test of homogeneity of variances
data: X by A
Bartlett's K-squared = 0.26802, df = 3, p-value = 0.9659
> bartlett.test(X~B,data = juice)
Bartlett test of homogeneity of variances
data: X by B
Bartlett's K-squared = 17.422, df = 5, p-value = 0.003766
结论: 对因素A, p值(0.966)远大于0.05, 接受原假设, 认为因素A的各水平下 的数据是等方差的; 对因素B, p值(0.003766)小于0.05, 拒绝原假设, 即认为因 素B不满足方差齐性要求
2 有交互作用的方差分析
R软件中仍用函数aov( )进行有交互作用的方差分析, 但其中的方差模型格式 为x A+B+A : B. 下面用一个例子来全面展示有交互作用方差分析过程.
例 2
有一个关于检验毒品强弱的试验, 给48只老鼠注射I、II、III三 种毒药(因素A), 同时有A、B、C、D 4种治疗方案(因素B), 这样的试验在每一 种因素组合下都重复四次测试老鼠的存活时间, 数据如表8.5所示. 试分析毒药 和治疗方案以及它们的交互作用对老鼠存活时间有无显著影响.
统计表2
首先以数据框形式输入数据, 并用函数plot( )作图.
# 将统计表的数据输入到命名为rats的数据框中
> rats<-data.frame( Time=c(0.31, 0.45, 0.46, 0.43, 0.82, 1.10, 0.88, 0.72, 0.43, 0.45, 0.63, 0.76, 0.45, 0.71, 0.66, 0.62, 0.38, 0.29, 0.40, 0.23, 0.92, 0.61, 0.49, 1.24, 0.44, 0.35, 0.31, 0.40, 0.56, 1.02, 0.71, 0.38, 0.22, 0.21, 0.18, 0.23, 0.30, 0.37, 0.38, 0.29, 0.23, 0.25, 0.24, 0.22, 0.30, 0.36, 0.31, 0.33),
+ toxicant=gl(3,16,48,labels = c("I","II","III")),
+ cure=gl(4,4,48,labels = c("A","B","C","D")))
> op <- par(mfrow=c(1,2)) #更改图形显示的界面排版
> plot(Time~toxicant + cure, data = rats)
Hit <Return> to see next plot:
显示出的统计图：
图2
图2显示两因素的各水平均有较大差异存在.
下面再用函数interaction.plot( )作出交互效应图, 以考查因素之间交互作用是否存在, R程序为
> with(rats,interaction.plot(toxicant, cure, Time, trace.label = "cure"))
> with(rats,interaction.plot(cure, toxicant, Time, trace.label = "toxicant"))
交互效应图如下所示：
图3
两图中的曲线并没有明显的相交情况出现, 因此我们初步认为两个因素没有交互作用
尽管如此, 由于实验误差的存在, 我们用方差分析函数aov( )对此进行确 认, 其中方差模型格式为x~A*�B, 或A+B+A : B, 表示不仅考虑因素A、B各 自的效应, 还考虑两者的交互效应. 若仅考虑A与B的交互效应则方差模型格式为A : B.
> rats.aov <- aov(Time~toxicant+cure+toxicant:cure,data = rats)
> summary(rats.aov)
Df Sum Sq Mean Sq F value Pr(>F)
toxicant 2 1.0356 0.5178 23.225 3.33e-07 ***
cure 3 0.9146 0.3049 13.674 4.13e-06 ***
toxicant:cure 6 0.2478 0.0413 1.853 0.116
Residuals 36 0.8026 0.0223
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
根据p值知, 因素Toxicant和Cure对Time的影响是高度显著的, 而交互作用对Time的影响却是不显著的.
再 进 一 步 使 用 前 面 的Bartlett检 验 因 素Toxicant和Cure下的数据是否满足方差齐性的要求, R程序如下.
> bartlett.test(Time~toxicant,data=rats)
Bartlett test of homogeneity of variances
data: Time by toxicant
Bartlett's K-squared = 25.806, df = 2, p-value = 2.49e-06
> bartlett.test(Time~cure,data = rats)
Bartlett test of homogeneity of variances
data: Time by cure
Bartlett's K-squared = 13.055, df = 3, p-value = 0.004519
其中各p值均小于0.05表明在0.05显著性水平下两因素下的方差 不满足齐性的要求, 这与图2是一致的.
