统计建模与R软件第七章习题答案(方差分析)
Ex7.1
ddi(1)
>lamp<-data.frame(X=c(115,116,98,83,103,107,118,116,73,89,85,97),A=factor(rep(1:3,c(4,4,4))))
> lamp.aov<-aov(X~A,data=lamp);summary(lamp.aov)
Df Sum Sq Mean Sq F value Pr(>F)
A 2 1304 652.0 4.923 0.0359 *
Residuals 9 1192 132.4
P值小于0.05，有显著差异。
(2)
对甲的区间估计：
> a<-c(115,116,98,83)
> t.test(a)
One Sample t-test
data: a
t = 13.1341, df = 3, p-value = 0.0009534
生活大爆炸第3季alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
78.04264 127.95736
sample estimates:
小角色mean of x
103
或者用这个命令更简单：
>attach(lamp)
> t.test(X[A==1])
乙的均值估计为111，95%置信区间为99.59932, 122.40068。
丙的均值估计为86，95%置信区间为70.08777, 101.91223。
(3)多重检验：
> attach(lamp)
P值不做调整：
> st(X,A,hod = "none")
Pairwi comparisons using t tests with pooled SD
data: X and A
1 2
2 0.351 -
3 0.066 0.013
P值进行Holm调整：
P value adjustment method: nonemicroblog
> st(X,A,hod = "holm",data)
Pairwi comparisons using t tests with pooled SD
data: X and A
1 2
2 0.35 -
3 0.13 0.04
P value adjustment method: holm
不论采取哪种方法，都可看出乙和丙有显著差异。
Ex7.2
(1)
>lamp<-data.frame(X=c(20,18,18,17,15,16,13,18,22,17,26,19,26,28,23,25,24,25,18,22,27,24,12,14),A=factor(rep(1:4,c(10,6,6,2))))
> lamp.aov<-aov(X ~ A, data=lamp);summary(lamp.aov)
Df Sum Sq Mean Sq F value Pr(>F)
A 3 351.7 117.24 15.11 2.28e-05 ***magnifier
Residuals 20 155.2 7.76
P值小于0.05，可认为四个厂生产的产品的变化率有显著差异。
(2)
> attach(lamp)
P值不做调整：
> st(X,A,hod = "none")
Pairwi comparisons using t tests with pooled SD
data: X and A
1 2 3
2 8.0e-05 - -
党在我心中的演讲稿3 0.00053 0.47666 -
4 0.05490 6.1e-05 0.00020
P value adjustment method: none
P值进行Holm调整：
> st(X,A,hod = "holm")
Pairwi comparisons using t tests with pooled SD
data: X and A
1 2 3
2 0.00040 - -
3 0.00158 0.47666 -
4 0.10979 0.00036 0.00079
P value adjustment method: holm
由此可得，除了A1和A4，A2和A3这两组的差异不显著外，其他组合的差异都很显著。
Ex7.3
>lamp1<-data.frame(X=c(30,27,35,35,29,33,32,36,26,41,33,31,43,45,53,44,51,53,54,37,47,57,48,42,82,66,66,86,56,52 ,76,83,72,73,59,53),A=factor(rep(1:3,c(12,12,12))))
>attach(lamp1)
正态性检验：
> st(X[A==1])
Shapiro-Wilk normality test
data: X[A == 1]
W = 0.9731, p-value = 0.9407
> st(X[A==2])
Shapiro-Wilk normality test
data: X[A == 2]
W = 0.9708, p-value = 0.9193
> st(X[A==3])
Shapiro-Wilk normality test
data: X[A == 3]
W = 0.9371, p-value = 0.4613
数据在三种水平下均是正态的。
方差齐性检验：
> st(X~A,data=lamp1)
Bartlett test of homogeneity of variances
jino
data: X by A
Bartlett's K-squared = 12.139, df = 2, p-value = 0.002312
P值小于0.05，认为各组方差不等。
Ex7.4
>lamp<-data.frame(X=c(2.79,2.69,3.11,3.47,1.77,2.44,2.83,2.52,3.83,3.15,4.70,3.97,2.03,2.87,3.65,5.09,5.41,3.47,4.92,4.07,2.18,3.13,3.77,4.26), g=factor(rep(1:3,c(8,8,8))))
先进行正态性和方差齐性检验以选择使用方差分析aov()还是KW检验st()。
正态性检验：
> attach(lamp)
> st(X[g==1])
Shapiro-Wilk normality test
data: X[g == 1]
W = 0.9659, p-value = 0.8638
> st(X[g==2])
Shapiro-Wilk normality test
data: X[g == 2]
W = 0.983, p-value = 0.9763
> st(X[g==3])
梯子英文 Shapiro-Wilk normality test
data: X[g == 3]
W = 0.99, p-value = 0.9951
三组数据都服从正态分布。
方差齐性检验：
> st(X~g,data=lamp)
宣布的意思 Bartlett test of homogeneity of variances
data: X by g
Bartlett's K-squared = 3.4559, df = 2, p-value = 0.1776
p值大于0.05，可认为三组方差齐。因此选用方差分析aov()或者KW检验st()均可。
> st(X~g,data=lamp)
Kruskal-Wallis rank sum test
data: X by g
Kruskal-Wallis chi-squared = 7.9322, df = 2, p-value = 0.01895
> lamp.aov<-aov(X~g,data=lamp);summary(lamp.aov)
Df Sum Sq Mean Sq F value Pr(>F)
g 2 6.437 3.218 4.284 0.0275 *
Residuals 21 15.776 0.751
两种方法得到的p值均小于0.05，可认为3种不同处理的诱导作用不同。
Ex7.5
根据题意，适用配伍组设计的Friedman秩和检验。
难受的英文
>lamp<-data.frame(X=c(23.1,57.6,10.5,23.6,11.9,54.6,21.0,20.3,22.7,53.2,9.7,19.6,13.8,47.1,13.6,23.6,22.5,53.7,10.8,21.1,13.7,39.2,13.7,16.3,22.6,53.1,8.3,21.6,13.3,37.0,14.8,14.8), g=gl(4,8),b=gl(8,1,32)) #其中g表示group，b表示block。
> st(X~g|b,data=lamp)
Friedman rank sum test
data: X and g and b
Friedman chi-squared = 6.45, df = 3, p-value = 0.09166
P值大于0.05，尚不能拒绝原假设。
Ex7.6
(1)>qua<-data.frame(x=c(4.6,4.3,6.1,6.5,6.8,6.4,6.3,6.7,3.4,3.8,4.0,3.8,4.7,4.3,3.9,3.5,6.5,7.0),a=gl(3,6,18), b=gl(3,2,18))
> qua.aov<-aov(x~a+b+a:b,data=qua);summary(qua.aov)
Df Sum Sq Mean Sq F value Pr(>F)
a 2 3.974 1.987 26.69 0.000164 ***
b 2 4.441 2.221 29.83 0.000107 ***
a:b 4 21.159 5.290 71.06 8.34e-07 ***
Residuals 9 0.670 0.074
两种因素以及其交互作用对产品质量的影响都很显著。
(2)最优条件为A3和B3组合。
> t.test(c(6.5,7.0))
One Sample t-test
data: c(6.5, 7)
t = 27, df = 1, p-value = 0.02357
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
3.573449 9.926551
sample estimates:
mean of x
6.75
点估计为6.75，区间估计为 3.573449，9.926551
(3)双因素方差分析的多重比较？不会...
Ex7.7
正交试验的方差分析。
L9(3^4)正交表。
>pro<-data.frame(Y=c(62.925,57.075,51.6,55.05,58.05,56.55,63.225,50.7,54.45),A=gl(3,3),B=gl(3,1,9),C=factor(c(1,2,3,2,3,1,3,1,2)))
> pro.aov<-aov(Y~A+B+C,data=pro);summary(pro.aov)
Df Sum Sq Mean Sq F value Pr(>F)
A 2 1.76 0.88 0.022 0.978
B 2 65.86 32.93 0.836 0.545
C 2 6.66 3.33 0.085 0.922
Residuals 2 78.78 39.39
结果显示这三个因素对水稻产量的影响均不明显。
K[i,j]<-mean(rate$Y[rate[j]="
> for (j in 1:3)
for (i in 1:3)
K[i,j]<-mean(rate$Y[rate[j]==i])
> K
A B C
1 41 47 45
2 48 55 57
3 61 48 48
B取水平2，A取水平3，C取水平2.
Ex7.8
表示不懂交互作用表。who knows?
Ex7.9
有重复试验的方差分析。
>out<-data.frame(Y=c(1.5,1.7,1.3,1.5,1.0,1.2,1.0,1.0,2.5,2.2,3.2,2.0,2.5,2.5,1.5,2.8,1.5,1.8,1.7,1.5,1.0,2.5,1.3,1.5,1.8,1.5,1.8,2.2,1.9,2.6,2.3,2.0),A=gl(2,16),B=gl(2,8,32),C=gl(2,4,32))
> out.aov<-aov(Y~A+B+C+A:B+A:C+B:C,data=out);summary(out.aov)
Df Sum Sq Mean Sq F value Pr(>F)
A 1 0.008 0.008 0.051 0.8224
B 1 4.728 4.728 31.142 8.37e-06 ***
C 1 0.038 0.038 0.249 0.6221
A:B 1 1.015 1.015 6.688 0.0159 *
A:C 1 0.428 0.428 2.818 0.1057
B:C 1 0.263 0.263 1.731 0.2002
Residuals 25 3.795 0.152
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
B的影响最显著，A和B的交互作用影响也很显著。
计算A和B的交互作用：
> ab<-function(x,y){
+ n<-length(x);z<-rep(0,n)
+ for (i in 1:n)
+ if (x[i]==y[i]){z[i]<-1} el{z[i]<-2}
+ factor(z)
+ }
> out$AB<-ab(out$A,out$B)
计算各水平的均值：
> K<-matrix(0,nrow=2,ncol=4,dimnames=list(1:2,c("A","B","C","AB")))
> for (j in 2:5)
+ for (i in 1:2)
+ K[i,j-1]<-mean(out$Y[out[j]==i])
> K
A B C AB
1 1.83750 1.43750 1.85625 1.64375
2 1.80625 2.20625 1.78750 2.00000
由于值越小越好，因此B应用水平1，AB也用水平1，因此A也得用水平1，C用水平2.
