15.8 灯丝的配料方案选优

   某灯泡厂用四种不同配料方案制成的灯丝生产了四批灯泡，在每批灯泡中随机地抽取若干个，测得其使用寿命（单位：小时），所得数据如表 $15.5$。

试问：这四种灯丝生产的灯泡其使用寿命有无显著差异？倘若使用寿命无显著差异，我们就可以从中选一种既经济又方便的配料方案；如果有显著差异，则希望选一种较优的方案，以便提高灯泡的使用寿命。

MATLAB程序

x = [1600 1610  1650  1680  1700  1720  1800,...

   1580 1640  1640  1700  1750,...

   1460 1550  1600  1620  1640  1740  1660  1820,...

    1510  1520  1530  1570  1680  1600];

g = [ones(1,7),2*ones(1,5),3*ones(1,8),4*ones(1,6)];

[P,ANOVATAB,STATS] = anova1(x,g)

得到

P =

   0.1229

ANOVATAB =

   'Source'    'SS'             'df'    'MS'             'F'         'Prob>F'

   'Groups'    [4.4361e+004]    [ 3]    [1.4787e+004]    [2.1494]    [0.1229]

   'Error'     [1.5135e+005]    [22]    [6.8796e+003]          []          []

   'Total'     [1.9571e+005]    [25]               []          []          []

STATS =

   gnames: {4x1 cell}

        n: [7 5 8 6]

   source: 'anova1'

    means: [1680 1662 1.6363e+003 1.5683e+003]

       df: 22

        s: 82.9433

R语言程序

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)

得到

           Df Sum Sq Mean Sq F value Pr(>F)

A            3  49212   16404   2.166  0.121

Residuals   22 166622    7574

> boxplot(X~A,data=lamp,col="red")

得到

以下转自网络，以供参考

例 小白鼠在接种了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)