统计建模与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.