R学习-方差分析

单因素方差分析
题解 
单因素四水平试验
四种不同配方材料A1,A2,A3,A4生产出的原件，寿命如下表，问：四种不同配方下原件的使用寿命有无显著差异
lamp<-data.frame(
X=c(1600, 1610, 1650, 1680, 1700, 1700, 1780, 1500, 1640,
1400, 1700, 1750, 1640, 1550, 1600, 1620, 1640, 1600,
1740, 1800, 1510, 1520, 1530, 1570, 1640, 1600),
A=factor(c(rep(1,7),rep(2,5), rep(3,8), rep(4,6)))
)
lamp.aov<-aov(X ~ A, data=lamp)
summary(lamp.aov)
例 小白鼠在接种了3种不同菌型的伤寒杆菌后的存活天数如表所示，判断小白鼠被注射三种菌型后的平均存活天数有无显著差异
mouse<-data.frame(
X=c( 2, 4, 3, 2, 4, 7, 7, 2, 2, 5, 4, 5, 6, 8, 5, 10, 7,
12, 12, 6, 6, 7, 11, 6, 6, 7, 9, 5, 5, 10, 6, 3, 10),
A=factor(c(rep(1,11),rep(2,10), rep(3,12)))
)
mouse.aov<-aov(X ~ A, data=mouse)
7.1.4 均值的多重比较
多重t检验：缺点，如果因素的水平较多，多次重复使用t检验会增加犯第一类错误的概率
P值的修正 
p.adjust()
调整方法有
Bonferroni "bonferroni"
Holm (1979) "holm"
Hochberg (1988) "hochberg"
Hommel (1988) "hommel"
Benjamini & Hochberg (1995) "BH"
Benjamini & Yekutieli (2001) "BY"
均值的多重比较
pairwise.t.test()
attach(mouse)
mu<-c(mean(X[A==1]), mean(X[A==2]), mean(X[A==3])); mu
pairwise.t.test(X, A, p.adjust.method = "none")
调整后的p值
pairwise.t.test(X, A, p.adjust.method = "holm")
pairwise.t.test(X, A, p.adjust.method = "bonferroni")
plot(mouse$X~mouse$A)
方差的齐次检验
方差分析的三个条件
1 可加性 假设模型是线性可加模型，每个处理效应与随机误差是可以叠加的
2 独立正态性 试验误差应当服从正态分布，而其相互独立
3 方差齐性，不同处理间的方差是一致的。
误差的正态性检验 shapiro.test()
attach(lamp)
shapiro.test(X[A==1])
shapiro.test(X[A==2])
shapiro.test(X[A==3])
shapiro.test(X[A==4])
方差齐性检验 Baerlett检验
bartlett.test()
bartlett.test(X~A, data=lamp)
bartlett.test(lamp$X, lamp$A)
对两个以上样本进行比较的Kruskal-Wallis 秩和检验
kruskal.test()
kruskal.test(formula, data, subset, na.action, ...)
例子 为了比较属于同一类的四种不同食谱的营养效果，将25只老鼠随机的氛围4组，每组分别为8只，4只，7只和6只，采用食谱甲乙丙丁喂养，假设其他条件均值相同，12周后测得体重增加量如表所示，在alpha=0.05水平上，检验各食谱的营养效果是否有显著差异。
food<-data.frame(
x=c(164, 190, 203, 205, 206, 214, 228, 257,
185, 197, 201, 231,
187, 212, 215, 220, 248, 265, 281,
202, 204, 207, 227, 230, 276),
g=factor(rep(1:4, c(8,4,7,6)))
)
kruskal.test(x~g, data=food)
另一种写法
kruskal.test(food$x, food$g)
或
A<-c(164, 190, 203, 205, 206, 214, 228, 257)
B<-c(185, 197, 201, 231)
C<-c(187, 212, 215, 220, 248, 265, 281)
D<-c(202, 204, 207, 227, 230, 276)
kruskal.test(list(A,B,C,D))
方差分析的过程
1 正态性检验
attach(food)
shapiro.test(x[g==1])
shapiro.test(x[g==2])
shapiro.test(x[g==3])
shapiro.test(x[g==4])
2 方差齐性检验
bartlett.test(x~g, data=food)
aovtable<-aov(x~g, data=food)
summary(aovtable)
Friedman秩和检验
在配伍组设计中，如果总体不能满足正态性和方差齐性的要求，可以采用Friedman秩和检验
friedman.test()
friedman.test(y, groups, blocks, ...)
friedman.test(formula, data, subset, na.action, ...)
24只小鼠按不同窝氛围八个区组，再把每个区组中的观察单位随机分配到3种不同的饲料组，喂养一定时间之后，测得小鼠肝中的铁含量，结果如表，所示，试分析不同饲料的小鼠肝中铁含量是否不同。
X<-matrix(
c(1.00, 1.01, 1.13, 1.14, 1.70, 2.01, 2.23, 2.63,
0.96, 1.23, 1.54, 1.96, 2.94, 3.68, 5.59, 6.96,
2.07, 3.72, 4.50, 4.90, 6.00, 6.84, 8.23, 10.33),
ncol=3, dimnames=list(1:8, c("A", "B", "C"))
)
friedman.test(X)
双因素方差分析
在一个农业试验中，考虑四种不同的品种A1,A2,A3,A4和三种不同的施肥方法B1,B2,B3得到的产量数据表。分析种子与施肥对产量有无显著影响。
1 不考虑交互作用时
agriculture<-data.frame(
Y=c(325, 292, 316, 317, 310, 318,
310, 320, 318, 330, 370, 365),
A=gl(4,3),
B=gl(3,1,12)
)
agriculture.aov <- aov(Y ~ A+B, data=agriculture)
summary(agriculture.aov)
考虑交互作用
分析树种与地理位置对松树生长的影响。A表示树种，B表示地区，对同一水平组合进行了五次测量，对结果进行方差分析。
tree<-data.frame(
Y=c(23, 25, 21, 14, 15, 20, 17, 11, 26, 21,
16, 19, 13, 16, 24, 20, 21, 18, 27, 24,
28, 30, 19, 17, 22, 26, 24, 21, 25, 26,
19, 18, 19, 20, 25, 26, 26, 28, 29, 23,
18, 15, 23, 18, 10, 21, 25, 12, 12, 22,
19, 23, 22, 14, 13, 22, 13, 12, 22, 19),
A=gl(3,20,60),
B=gl(4,5,60)
)
tree.aov <- aov(Y ~ A+B+A:B, data=tree)
summary（tree.aov）
双因素方差分析的方差齐性检验
attach(tree)
正态性检验
shapiro.test(Y[A==1])
shapiro.test(Y[A==2])
shapiro.test(Y[A==3])
shapiro.test(Y)
shapiro.test(Y)
shapiro.test(Y)
shapiro.test(Y)
方差齐性检验
bartlett.test(Y~A, data=tree)
正交试验设计与方差分析
正交试验设计分析
为提高某种化学产品的转化率，考虑三个有关因素，反应温度A，反应时间B，和用碱量C，各因素选三个水平，
rate<-data.frame(
A=gl(3,3),
B=gl(3,1,9),
C=factor(c(1,2,3,2,3,1,3,1,2)),
Y=c(31, 54, 38, 53, 49, 42, 57, 62, 64)
)
K<-matrix(0, nrow=3, ncol=3, dimnames=list(1:3, c("A","B","C")))
> for (j in 1:3)
for (i in 1:3)
K<-mean(rate$Y[rate[j]==i])
K
plot(as.vector(K), axes=F, xlab="Level", ylab="Rate")
xmark<-c(NA,"A1","A2","A3","B1","B2","B3","C1","C2","C3",NA)
axis(1,0:10,labels=xmark)
axis(2,4*10:16)
axis(3,0:10,labels=xmark)
axis(4,4*10:16)
lines(K[,"A"]); lines(4:6, K[,"B"]); lines(7:9,K[,"C"])
正交试验的方差分析
假设没有交互作用
rate.aov<-aov(Y~A+B+C, data=rate)
summary(rate.aov)
有交互作用的实验
在梳棉机上纺纱，为了提高质量，选了三个因素，每个因素两个水平，三个因素之间可能有交互作用，设计实验。并做方差分析
cotton<-data.frame(
Y=c(0.30, 0.35, 0.20, 0.30, 0.15, 0.50, 0.15, 0.40),
A=gl(2,4), B=gl(2,2,8), C=gl(2,1,8)
)
cotton.aov<-aov(Y~A+B+C+A:B+A:C+B:C, data=cotton)
summary(cotton.aov)
结果表明，A，A:B, B:C的F值都很小，P值很大，可将这三个因素去掉。
cotton.new<-aov(Y~B+C+A:C, data=cotton)
研究四种药物对淡色库蚊的杀灭作用，每种药物取三水平，采用L9(34)正交表，不考虑交互作用，相同实验条件下均做四次实验，检验四种药物对淡色库蚊杀灭作用有无差别。
mosquito<-data.frame(
A=gl(3, 12), B=gl(3,4,36),
C=factor(rep(c(1,2,3,2,3,1,3,1,2),rep(4,9))),
D=factor(rep(c(1,2,3,3,1,2,2,3,1),rep(4,9))),
Y=c( 9.41, 7.19, 10.73, 3.73, 11.91, 11.85, 11.00, 11.72,
10.67, 10.70, 10.91, 10.18, 3.87, 3.18, 3.80, 4.85,
4.20, 5.72, 4.58, 3.71, 4.29, 3.89, 3.88, 4.71,
7.62, 7.01, 6.83, 7.41, 7.79, 7.38, 7.56, 6.28,
8.09, 8.17, 8.14, 7.49)
)
mosquito.aov<-aov(Y~A+B+C+D, data=mosquito)
source("anova.tab.R"); anova.tab(mosquito.aov)