统计建模与R软件第七章习题…

原文地址：统计建模与R软件第七章习题答案(方差分析) 作者：

Ex7.1
(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
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值不做调整：
> pairwise.t.test(X,A,p.adjust.method = "none")

        Pairwise 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: none
> pairwise.t.test(X,A,p.adjust.method = "holm",data)

        Pairwise 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 ***
Residuals   20  155.2    7.76
P值小于0.05，可认为四个厂生产的产品的变化率有显著差异。
(2)
> attach(lamp)
P值不做调整：
> pairwise.t.test(X,A,p.adjust.method = "none")

        Pairwise 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调整：
> pairwise.t.test(X,A,p.adjust.method = "holm")

        Pairwise 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)
正态性检验：
> shapiro.test(X[A==1])

        Shapiro-Wilk normality test

data:  X[A == 1]
W = 0.9731, p-value = 0.9407

> shapiro.test(X[A==2])

        Shapiro-Wilk normality test

data:  X[A == 2]
W = 0.9708, p-value = 0.9193

> shapiro.test(X[A==3])

        Shapiro-Wilk normality test

data:  X[A == 3]
W = 0.9371, p-value = 0.4613
数据在三种水平下均是正态的。
方差齐性检验：
> bartlett.test(X~A,data=lamp1)

        Bartlett test of homogeneity of variances

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检验kruskal.test()。
正态性检验：
> attach(lamp)
> shapiro.test(X[g==1])

        Shapiro-Wilk normality test

data:  X[g == 1]
W = 0.9659, p-value = 0.8638

> shapiro.test(X[g==2])

        Shapiro-Wilk normality test

data:  X[g == 2]
W = 0.983, p-value = 0.9763

> shapiro.test(X[g==3])

        Shapiro-Wilk normality test

data:  X[g == 3]
W = 0.99, p-value = 0.9951
三组数据都服从正态分布。
方差齐性检验：
> bartlett.test(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检验kruskal.test()均可。
> kruskal.test(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。
> friedman.test(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} else{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.