代码改变世界

『原创』统计建模与R软件-第七章 方差分析

2016-01-20 11:16  Digging4  阅读(8967)  评论(0编辑  收藏  举报

摘要: 本文由digging4发表于:http://www.cnblogs.com/digging4/p/5091536.html

统计建模与R软件-第七章 方差分析


7.1 三个工厂生产同一种零件,现从各厂产品中分别抽取4件产品作检测,其检测强度如表7.25所示。

表7.25 产品检测数据

工厂 零件强度
115 116 98 83
103 107 118 116
73 89 85 97

(1)对数据作方差分析,判断三个厂的产品的零件强度是否有显著差异

(2)求每个工厂生产产品零件强度的均值,作出相应的区间估计(\(\alpha=0.05\)

(3)对数据作多重检验。

df <- data.frame(X = c(115, 116, 98, 83, 103, 107, 118, 116, 73, 89, 85, 97), 
    A = gl(3, 4))

fit.aov <- aov(X ~ A, data = df)
summary(fit.aov)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## A            2   1304     652    4.92  0.036 *
## Residuals    9   1192     132                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 由于P值为0.036,小于0.05,因此认为因素A(三个厂)对零件强度有显著差异。
plot(df$X ~ df$A)

# 从箱线图也可以看出三个工厂的零件强度有显著差异。

# 对三个厂的零件强度求均值及其区间估计
mean.1 <- mean(df$X[df$A == 1])
mean.1.t <- t.test(df$X[df$A == 1], conf.level = 0.95)
mean.2 <- mean(df$X[df$A == 2])
mean.2.t <- t.test(df$X[df$A == 2], conf.level = 0.95)
mean.3 <- mean(df$X[df$A == 3])
mean.3.t <- t.test(df$X[df$A == 3], conf.level = 0.95)

reulst <- data.frame(A = c(1, 2, 3), mean = c(mean.1, mean.2, mean.3), t.test.down = c(mean.1.t[4]$conf.int[1], 
    mean.2.t[4]$conf.int[1], mean.3.t[4]$conf.int[1]), t.test.up = c(mean.1.t[4]$conf.int[2], 
    mean.2.t[4]$conf.int[2], mean.3.t[4]$conf.int[2]))
reulst
##   A mean t.test.down t.test.up
## 1 1  103       78.04     128.0
## 2 2  111       99.60     122.4
## 3 3   86       70.09     101.9

# 对数据做多重检验
pairwise.t.test(df$X, df$A, p.adjust.method = "none")
## 
## 	Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$X and df$A 
## 
##   1     2    
## 2 0.351 -    
## 3 0.066 0.013
## 
## P value adjustment method: none
# 可以看到1和2的p值为0.351,可认为两水平无差异。
# 1和3的P值0.066有差异,不显著;2和3的P值0.013,差异显著。
# 结论为工厂丙的零件强度与其他两个厂有显著差异,甲乙两厂无差异。

7.2有四种产品,\(A_i,i=1,2,3\)分别为国内甲、乙、丙三个工厂生产的产品,\(A_4\)国外同类产品,现从各厂分别取10,6,6和2个产品做300小时连续磨损老化试验,得变化率如表7.26所示,假定各厂产品试验变化率服从等方差的正态分布。

表7.26 磨损老化试验数据

产品 变化率
\(A_1\) 20 18 19 17 15 16 13 18 22 17
\(A_2\) 26 19 26 28 23 25
\(A_3\) 24 25 18 22 27 24
\(A_4\) 12 14

(1)试问四个厂生产的产品的变化率是否有显著差异?

(2)若有差异,请做进一步的检验。i)国内产品与国外产品有无显著差异?ii)国内各厂家的产品有无显著差异?

df <- data.frame(X = c(20, 18, 19, 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))))
fit.aov <- aov(X ~ A, data = df)
summary(fit.aov)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## A            3    346   115.3    14.7 2.8e-05 ***
## Residuals   20    157     7.9                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# P值为2.8e-05,可认为四个厂生产的产品的变化率有显著差异
# 国内产品与国外产品有无显著差异? 将国内产品看做一类,国外产品看做另一类
df2 <- data.frame(X = c(20, 18, 19, 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:2, c(22, 2))))
fit.aov2 <- aov(X ~ A, data = df2)
summary(fit.aov2)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## A            1    117   117.3    6.69  0.017 *
## Residuals   22    386    17.5                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# p值为0.017,可认为国内产品与国外产品有显著差异

# 国内各厂家的产品有无显著差异?
fit.aov3 <- aov(X ~ A, data = df[df$A != 4, ])
summary(fit.aov3)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## A            2    229   114.3      14 0.00018 ***
## Residuals   19    155     8.2                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# p值为0.00018,可认为国内各厂家的产品有显著差异

7.3某单位在大白鼠营养试验中,随机将大白鼠分为三组,测得每组12只大白鼠尿中氨氮的排出量\(X(mg/6d)\),数据由表7.27所示,试对该资料做正态性检验和方差齐性检验。

表7.27 白鼠尿中氨氮检测数据

白鼠 大白鼠营养试验中各组大鼠尿氨氮排出量\(X(mg/6d)\)
第一组  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
rat <- 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))))
# 做正态性检验
shapiro.test(rat$X[rat$A == 1])
## 
## 	Shapiro-Wilk normality test
## 
## data:  rat$X[rat$A == 1] 
## W = 0.9731, p-value = 0.9407
shapiro.test(rat$X[rat$A == 2])
## 
## 	Shapiro-Wilk normality test
## 
## data:  rat$X[rat$A == 2] 
## W = 0.9708, p-value = 0.9193
shapiro.test(rat$X[rat$A == 3])
## 
## 	Shapiro-Wilk normality test
## 
## data:  rat$X[rat$A == 3] 
## W = 0.9371, p-value = 0.4613
# 得到三组的p值均大于0.05,可认为数据在三个水平下均是正态分布

# 方差齐性检验,是检验数据在不同水平下方差是否相同
bartlett.test(X ~ A, data = rat)  #bartlett.test(rat$x,rat$A) 效果相同
## 
## 	Bartlett test of homogeneity of variances
## 
## data:  X by A 
## Bartlett's K-squared = 12.14, df = 2, p-value = 0.002312
# p-value = 0.002312<0.05,因此可认为各组数据方差是不等的。

7.4 以小白鼠为对象研究正常肝核糖核酸(RNA)对癌细胞的生物作用,试验分别为对照组(生理盐水),水层RNA组和酚层RNA组,分别用此三种不同处理诱导肝癌细胞的果糖二磷酸酯酶(FDP酶)活力,数据如表7.28所示,问三种不同处理的诱导作用是否相同?

表7.28 三种不同处理的诱导结果

处理方法 诱导结果
对照组 2.79 2.69 3.11 3.47 1.77 2.44 2.83 2.52
水层RNA  3.83 3.15 4.70 3.97 2.03 2.87 3.65 5.09
酚层RNA  5.41 3.47 4.92 4.07 2.18 3.13 3.77 4.26
rat <- 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.7, 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), A = factor(rep(1:3, c(8, 8, 8))))
fit.aov <- aov(X ~ A, data = rat)
summary(fit.aov)
##             Df Sum Sq Mean Sq F value Pr(>F)  
## A            2   6.44    3.22    4.28  0.028 *
## Residuals   21  15.78    0.75                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 可认为三组的作用不相同
plot(rat$X ~ rat$A)

# 从图中可以可以看出,第一组显著低于第二和第三组

7.5 为研究人们在催眠状态下对各种情绪的反应十分有差异,选取了8个受试者,在催眠状态下,要求每人按任意次序作出恐惧、愉快、忧虑和平静4中反映。表7.29给出了各受试者在处于这4中清醒状态下皮肤的电位变化值,试在\(\alpha=0.05\)下,检验受试者在催眠状态下对这4种情绪的反应力是否有显著差异。

表7.29:4种情绪状态下皮肤的电位变化值(单位: \(mV\)

情绪状态 受试者
1 2 3 4 5 6 7 8
恐惧 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
X = c(23.1, 57.6, 10.5, 23.6, 11.9, 54.6, 21, 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, 14.8, 14.8)
g <- gl(4, 8)  #分为三组,每组8个
b <- gl(8, 1, 32)  #每组内又有8个人,有顺序
friedman.test(X, g, b)
## 
## 	Friedman rank sum test
## 
## data:  X, g and b 
## Friedman chi-squared = 6.45, df = 3, p-value = 0.09166
# p-value = 0.09166>0.05,接受原假设,可认为四种情绪的反应力无显著差异

7.6 为了提高化工厂的产品质量,需要寻求最优反应温度与反应压力的配合,为此选择如下水平:

A:反应温度(\(^0C\)) 60 70 80

B:反应压力(公斤) 2 2.5 3

在每个\(A_iB_j\)条件下做两次试验,其产量如表7.30所示

(1)对数据作方差分析(应考虑交互作用)

(2)求最优条件下平均产量的点估计和区间估计

(3)对\(A_iB_j\)条件下平均产量作多重比较。

表7.30 试验数据

\(A_1\) \(A_2\) \(A_3\)
\(B_1\) 4.6 4.3 6.1 6.5 6.8 6.4
\(B_2\) 6.3 6.7 3.4 3.8 4.0 3.8
\(B_3\) 4.7 4.3 3.9 3.5 6.5 7.0
product <- 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, 3.8, 4.7, 4.3, 3.9, 3.5, 6.5, 7), A = gl(3, 2, 18), B = gl(3, 6, 18))
# 首先应该对数据做不同水平下的正态性检验,和各水平的方差齐性检验
# 对数据作方差分析(应考虑交互作用)
product.aov <- aov(X ~ A + B + A:B, data = product)
summary(product.aov)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## A            2   4.44    2.22    29.8 0.00011 ***
## B            2   3.97    1.99    26.7 0.00016 ***
## A:B          4  21.16    5.29    71.1 8.3e-07 ***
## Residuals    9   0.67    0.07                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 可以看出,因素A和B,及其交互因素对产品质量有显著影响

# 求最优条件下平均产量的点估计和区间估计 计算两个因素在不同条件下的均值
tapply(product$X, product$A, mean)
##     1     2     3 
## 5.150 4.533 5.750
tapply(product$X, product$B, mean)
##     1     2     3 
## 5.783 4.667 4.983
# 可以看到在A3和B1条件,为最优条件

product.best <- c(6.8, 6.4)
# 均值的点估计为
mean(product.best)
## [1] 6.6
# 均值的区间估计为
t.test(product.best, alternative = c("two.sided"))
## 
## 	One Sample t-test
## 
## data:  product.best 
## t = 33, df = 1, p-value = 0.01929
## alternative hypothesis: true mean is not equal to 0 
## 95 percent confidence interval:
##  4.059 9.141 
## sample estimates:
## mean of x 
##       6.6

# 对$A_iB_j$条件下平均产量作多重比较
TukeyHSD(product.aov)
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = X ~ A + B + A:B, data = product)
## 
## $A
##        diff     lwr     upr  p adj
## 2-1 -0.6167 -1.0565 -0.1768 0.0089
## 3-1  0.6000  0.1602  1.0398 0.0104
## 3-2  1.2167  0.7768  1.6565 0.0001
## 
## $B
##        diff     lwr     upr  p adj
## 2-1 -1.1167 -1.5565 -0.6768 0.0002
## 3-1 -0.8000 -1.2398 -0.3602 0.0017
## 3-2  0.3167 -0.1232  0.7565 0.1653
## 
## $`A:B`
##          diff     lwr     upr  p adj
## 2:1-1:1  1.85  0.7706  2.9294 0.0015
## 3:1-1:1  2.15  1.0706  3.2294 0.0005
## 1:2-1:1  2.05  0.9706  3.1294 0.0007
## 2:2-1:1 -0.85 -1.9294  0.2294 0.1557
## 3:2-1:1 -0.55 -1.6294  0.5294 0.5677
## 1:3-1:1  0.05 -1.0294  1.1294 1.0000
## 2:3-1:1 -0.75 -1.8294  0.3294 0.2502
## 3:3-1:1  2.30  1.2206  3.3794 0.0003
## 3:1-2:1  0.30 -0.7794  1.3794 0.9600
## 1:2-2:1  0.20 -0.8794  1.2794 0.9965
## 2:2-2:1 -2.70 -3.7794 -1.6206 0.0001
## 3:2-2:1 -2.40 -3.4794 -1.3206 0.0002
## 1:3-2:1 -1.80 -2.8794 -0.7206 0.0019
## 2:3-2:1 -2.60 -3.6794 -1.5206 0.0001
## 3:3-2:1  0.45 -0.6294  1.5294 0.7612
## 1:2-3:1 -0.10 -1.1794  0.9794 1.0000
## 2:2-3:1 -3.00 -4.0794 -1.9206 0.0000
## 3:2-3:1 -2.70 -3.7794 -1.6206 0.0001
## 1:3-3:1 -2.10 -3.1794 -1.0206 0.0006
## 2:3-3:1 -2.90 -3.9794 -1.8206 0.0000
## 3:3-3:1  0.15 -0.9294  1.2294 0.9995
## 2:2-1:2 -2.90 -3.9794 -1.8206 0.0000
## 3:2-1:2 -2.60 -3.6794 -1.5206 0.0001
## 1:3-1:2 -2.00 -3.0794 -0.9206 0.0008
## 2:3-1:2 -2.80 -3.8794 -1.7206 0.0001
## 3:3-1:2  0.25 -0.8294  1.3294 0.9857
## 3:2-2:2  0.30 -0.7794  1.3794 0.9600
## 1:3-2:2  0.90 -0.1794  1.9794 0.1219
## 2:3-2:2  0.10 -0.9794  1.1794 1.0000
## 3:3-2:2  3.15  2.0706  4.2294 0.0000
## 1:3-3:2  0.60 -0.4794  1.6794 0.4737
## 2:3-3:2 -0.20 -1.2794  0.8794 0.9965
## 3:3-3:2  2.85  1.7706  3.9294 0.0001
## 2:3-1:3 -0.80 -1.8794  0.2794 0.1980
## 3:3-1:3  2.25  1.1706  3.3294 0.0003
## 3:3-2:3  3.05  1.9706  4.1294 0.0000
# 3-2 0.3167 -0.1232 0.7565 0.1653 表示因素B的3和2差异不显著 A:B
# 的比较类似,p值大于0.05的均可以认为扎伊不显著

7.7 某良种繁殖场为了提高水稻产量,指定试验的因素如表7.31所示。试选择\(L_9(3^4)\)正交表安排试验,假定相应的产量为(单位:\(kg/100m^2\))62.925,57.075,51.6,55.05,58.05,56.55,63.225,50.7, 54.45。试对试验结果进行方差分析,并给出一组较好的种植条件

表7.31 水稻的试验因素水平表

因素 水平
1 2 3
品种 窄叶青8号 南二矮5号 珍珠矮11号
密度 4.50棵/100\(m^2\) 3.75棵/100\(m^2\) 3.00棵/100\(m^2\)
施肥量 0.75\(kg/100m^2\) 0.375\(kg/100m^2\) 1.125\(kg/100m^2\)

设三个因素为A(品种),B(密度),C(施肥量),采用\(L_9(3^4)\)正交表安排试验,将ABC放到前三列,得到试验结果表

试验号 A(品种) B(密度) C(施肥量) 产量
1 1 1 1 62.925
2 1 2 2 57.075
3 1 3 3 51.6
4 2 1 2 55.05
5 2 2 3 58.05
6 2 3 1 56.55
7 3 1 3 63.225
8 3 2 1 50.7
9 3 3 2 54.45
rice <- 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(62.925, 57.075, 51.6, 55.05, 58.05, 56.55, 63.225, 50.7, 
    54.45))
# 计算三个因素下产量的均值
tapply(rice$Y, rice$A, mean)
##     1     2     3 
## 57.20 56.55 56.12
tapply(rice$Y, rice$B, mean)
##     1     2     3 
## 60.40 55.27 54.20
tapply(rice$Y, rice$C, mean)
##     1     2     3 
## 56.73 55.52 57.62
#
# 可见,A取2,B取1,C取3时,达到最优条件。即最有种植条件为,品种:南二矮5号;密度:4.50棵/100$m^2$;施肥量:1.125$kg/100m^2$
# 做无交互作用方差分析
rice.aov <- aov(Y ~ A + B + C, data = rice)
summary(rice.aov)
##             Df Sum Sq Mean Sq F value Pr(>F)
## A            2    1.8     0.9    0.02   0.98
## B            2   65.9    32.9    0.84   0.54
## C            2    6.7     3.3    0.08   0.92
## Residuals    2   78.8    39.4
# 发现ABC对结果的影响均不显著

7.8 某单位研究四种因素对钉螺产卵数(Y)的影响,制定试验的因素如表7.32所示,试选择\(L_8(2^7)\)正交表安排试验,假定相应的钉螺产卵数为(单位:个),86,95,91,94,91,96,83,88,试对试验结果进行方差分析,并给出一组较好的灭螺方案(考虑有交互作用)。

表7.32 钉螺产卵影响试验因素的水平表

因素 水平
1 2
温度(A) 5\(^0C\) 10\(^0C\)
含氧量(B) 0.5 5.0
含水量(C) 10% 30%
pH值(D) 6.0 8.0

设四个因素为温度(A),含氧量(B),含水量(C),pH值(D),采用\(L_8(2^7)\)两列间交互作用的正交表安排试验,第1列放A,第2列放B,第3列放A×B,第4列放C,第5列A×C,第6列B×C,第7列D,得到试验结果表

试验号 温度(A) 含氧量(B) A×B 含水量(C) A×C B×C pH值(D) 钉螺产卵数
1 1 1 1 1 1 1 1 86
2 1 1 1 2 2 2 2 95
3 1 2 2 1 1 2 2 91
4 1 2 2 2 2 1 1 94
5 2 1 2 1 2 1 2 91
6 2 1 2 2 1 2 1 96
7 2 2 1 1 2 2 1 83
8 2 2 1 2 1 1 2 88
insect <- data.frame(A = gl(2, 4), B = gl(2, 2, 8), C = gl(2, 1, 8), D = factor(c(1, 
    2, 2, 1, 2, 1, 1, 2)), Y = c(86, 95, 91, 94, 91, 96, 83, 88))

# 计算四个个因素下钉螺产卵数的均值
tapply(insect$Y, insect$A, mean)
##    1    2 
## 91.5 89.5
tapply(insect$Y, insect$B, mean)
##  1  2 
## 92 89
tapply(insect$Y, insect$C, mean)
##     1     2 
## 87.75 93.25
tapply(insect$Y, insect$D, mean)
##     1     2 
## 89.75 91.25

#
# 温度(A)取水平2,含氧量(B)取水平2,含水量(C)取水平1,pH值(D)取水平1,此时钉螺产卵数最少

# 做有交互作用方差分析
insect.aov <- aov(Y ~ A + B + A:B + C + A:C + B:C + D, data = insect)
summary(insect.aov)
##             Df Sum Sq Mean Sq
## A            1    8.0     8.0
## B            1   18.0    18.0
## C            1   60.5    60.5
## D            1    4.5     4.5
## A:B          1   50.0    50.0
## A:C          1    0.5     0.5
## B:C          1    4.5     4.5
# 发现ABCD及其交互作用,对结果的影响均不显著,去掉Sum
# Sq最大的三个因素B,C,A:B

insect.aov <- aov(Y ~ A + A:C + B:C + D, data = insect)
summary(insect.aov)
##             Df Sum Sq Mean Sq F value Pr(>F)
## A            1    8.0     8.0    0.16   0.76
## D            1    4.5     4.5    0.09   0.81
## A:C          2   61.0    30.5    0.61   0.67
## C:B          2   22.5    11.3    0.23   0.83
## Residuals    1   50.0    50.0
# 可看出各因素及其交互因素对钉螺产卵数无显著影响

7.9 某工厂为了提高零件内孔研磨工序质量进行工艺的参数优选试验,考察孔的锥度值,希望其越小越好,在试验中考察因子的水平表7.33,试选择\(L_8(2^7)\)正交表安排试验,其表头设计如表7.34所示,在每一条件下加工了四个零件,测量其锥度,试验结果如表7.35所示,试对试验结果进行方差分析,并给出一组较好的工艺参数指标。

表7.33 因子水平表

因素 水平
1 2
研孔工艺设备(A) 通用夹具 专用夹具
生铁研圈材质(B) 特殊铸铁 一般灰铸铁
留研量(mm)(C) 0.01 0.015

表7.34 试验结果

表头设计 A B --- C --- --- ---
列号 1 2 3 4 5 6 7

表7.35 试验结果

试验号 试验值
1 1.5 1.7 1.3 1.5
2 1.0 1.2 1.0 1.0
3 2.5 2.2 3.2 2.0
4 2.5 2.5 1.5 2.8
5 1.5 1.8 1.7 1.5
6 1.0 2.5 1.3 1.5
7 1.8 1.5 1.8 2.2
8 1.9 2.6 2.3 2.0

选择\(L_8(2^7)\)正交表安排试验,ABC分别放到124列

试验号 A B C 试验值
1 1 1 1 1.5 1.7 1.3 1.5
2 1 1 2 1.0 1.2 1.0 1.0
3 1 2 1 2.5 2.2 3.2 2.0
4 1 2 2 2.5 2.5 1.5 2.8
5 2 1 1 1.5 1.8 1.7 1.5
6 2 1 2 1.0 2.5 1.3 1.5
7 2 2 1 1.8 1.5 1.8 2.2
8 2 2 2 1.9 2.6 2.3 2.0
# 有重复试验的方差分析
hole <- data.frame(A = gl(2, 4, 32), B = gl(2, 2, 32), C = gl(2, 1, 32), Y = c(1.5, 
    1, 2.5, 2.5, 1.5, 1, 1.8, 1.9, 1.7, 1.2, 2.2, 2.5, 1.8, 2.5, 1.5, 2.6, 1.3, 
    1, 3.2, 1.5, 1.7, 1.3, 1.8, 2.3, 1.5, 1, 2, 2.8, 1.5, 1.5, 2.2, 2))
hole.aov <- aov(Y ~ A + B + C, data = hole)
summary(hole.aov)
##             Df Sum Sq Mean Sq F value  Pr(>F)    
## A            1   0.01    0.01    0.04    0.84    
## B            1   4.73    4.73   24.06 3.6e-05 ***
## C            1   0.04    0.04    0.19    0.66    
## Residuals   28   5.50    0.20                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# 可以看出因素B对实验结果影响显著 计算三个因素下钉螺产卵数的均值
tapply(hole$Y, hole$A, mean)
##     1     2 
## 1.837 1.806
tapply(hole$Y, hole$B, mean)
##     1     2 
## 1.438 2.206
tapply(hole$Y, hole$C, mean)
##     1     2 
## 1.856 1.788

# 因此最优的工艺参数指标为A1,B2,C1