解法并不单一，下列方法带有璇子个人的偏好，因此仅供参考。如有错误，欢迎在评论区斧正！
8.1 三个工厂生产同一种零件，现从各厂产品中分别抽取4 件产品作检测，其检测强度如表8-7 所示。
（1）对数据作方差分析，判断三个厂的产品的零件强度是否有显著差异；
#导入数据
> factory <- as.factor(c(rep("甲",4),rep("乙",4),rep("丙",4)))
> strength <- c(115,116,98,83,103,107,118,116,73,89,85,97)
> data8.1 <- data.frame(factory,strength)
#各组情况
aggregate(strength, by=list(factory), FUN=mean)
Group.1 x
1 丙 86
2 甲 103
3 乙 111
> aggregate(strength, by=list(factory), FUN=sd)
Group.1 x
1 丙 10.000000
2 甲 15.684387
3 乙 7.164728 #初步可见乙工厂的零件强度最好、质量最稳定
#组间差异
> fit8.1 <- aov(strength ~ factory)
> summary(fit8.1)
Df Sum Sq Mean Sq F value Pr(>F)
factory 2 1304 652.0 4.923 0.0359 *
Residuals 9 1192 132.4
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
由单因素方差分析结果对应P值<0.05，拒绝原假设，即认为三个厂的产品的零件强度的确存在显著差异。
（2）求每个工厂生产产品零件强度的均值，作出相应的区间估计(α= 0.05）；
#各组均值
aggregate(strength, by=list(factory), FUN=mean)
Group.1 x
1 丙 86
2 甲 103
3 乙 111
#置信区间
confint <- function(x,sigma=-1,alpha=0.5){
+ n <- length(x)
+ m <- mean(x)
+ if(sigma >= 0){
+ deta <- sigma/sqrt(n)*qnorm(1-alpha/2);df <- n
+ }
+ else{
+ deta <- sd(x)/sqrt(n)*qt(1-alpha/2,n-1);df <- n-1
+ }
+ data.frame(mean=m, df=df, a=m-deta, b=m+deta)
+ }
>
> confint1 <- confint(data8.1[c(1:4),]$strength)
> confint2 <- confint(data8.1[c(5:8),]$strength)
> confint3 <- confint(data8.1[c(9:12),]$strength)
confint1
mean df a b
1 103 3 97.00157 108.9984
> confint2
mean df a b
1 111 3 108.2599 113.7401
> confint3
mean df a b
1 86 3 82.17554 89.82446
为便于观察，下面通过gplots包可视化各组对应均值和置信区间。
library(gplots)
> plotmeans(strength ~ factory, xlab = "factory", ylab = "strength", main="Mean plot \n with 95% CI")
（3）对数据作多重检验。
TukeyHSD(fit8.1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = strength ~ factory)
$factory
diff lwr upr p adj
甲-丙 17 -5.720515 39.72051 0.1471244
乙-丙 25 2.279485 47.72051 0.0323077
乙-甲 8 -14.720515 30.72051 0.6050019
> par(las=2)
> plot(TukeyHSD(fit8.1))
从结果看，乙工厂生产的零件强度比甲、丙都强，甲工厂的零件强度次之。同时，结合前面的分析，乙工厂生产的零件质量也最稳定。不难得出乙工厂在生产该零件方面最优的结论。
8.2 有四种产品，Ai,i=1,2,3 分别为国内甲、乙、丙三个工厂生产的产品，A4 为国外同类产品，现从各厂分别取10,6,6 和2 个产品做300 小时连续磨损老化试验，得变化率如表8-8 所示，假定各厂产品试验变化率服从等方差的正态分布。
（1）试问四个厂生产的产品的变化率是否有显著差异？
#数据录入
> product <- as.factor(c(rep("A1",10),rep("A2",6),rep("A3",6),rep("A4",2)))
> ratio <- 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)
> data8.2 <- data.frame(product,ratio)
> #各组情况
> aggregate(ratio, by=list(product), FUN=mean)
Group.1 x
1 A1 17.50000
2 A2 24.50000
3 A3 23.33333
4 A4 13.00000
> aggregate(ratio, by=list(product), FUN=sd)
Group.1 x
1 A1 2.549510
2 A2 3.146427
3 A3 3.076795
4 A4 1.414214
> #检验组间差异
> fit8.2 <- aov(ratio ~ product)
> summary(fit8.2)
Df Sum Sq Mean Sq F value Pr(>F)
product 3 346.0 115.33 14.66 2.79e-05 ***
Residuals 20 157.3 7.87
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> #可视化
> library(gplots)
> plotmeans(ratio ~ product, xlab = "product", ylab = "ratio", main="Mean plot \n with 95% CI")
与第一问类似，由上述结果P值远小于0.05得：四个厂生产的产品的变化率的确有显著差异（这里一个工厂对应一个产品类型）。可视化结果如下：
（2）若有差异，请做进一步的检验。i）国内产品与国外产品有无显著差异？ii）国内各厂家的产品有无显著差异？
> TukeyHSD(fit8.2) #多重比较
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = ratio ~ product)
$product
diff lwr upr p adj
A2-A1 7.000000 2.946103 11.053897 0.0005414
A3-A1 5.833333 1.779436 9.887231 0.0033963
A4-A1 -4.500000 -10.580846 1.580846 0.1964846
A3-A2 -1.166667 -5.699062 3.365728 0.8877803
A4-A2 -11.500000 -17.909774 -5.090226 0.0003529
A4-A3 -10.333333 -16.743108 -3.923559 0.0011245
#可视化
> par(las=2)
> plot(TukeyHSD(fit8.2))
由图可见：
（1）除A1，A4（国外产品）与A2、A3产品的变化率有显著差异，并且国外产品（A4）的变化率比国内产品的变化率都低。
（2）同时，就国内厂家的产品来说，甲与乙、丙的差异显著，而乙丙两家的产品无明显差异。
8.3 某单位在大白鼠营养试验中，随机将大白鼠分为三组，测得每组12 只大白鼠尿中氨氮的排出量X(mg/6d)，数据由表8-9 所示，试对该资料做正态性检验和方差齐性检验。
#数据录入
> group <- as.factor(c(rep("g1",12),rep("g2",12),rep("g3",12)))
> output <- 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)
> data8.3 <- data.frame(group,output)
>
> #正态性检验
> library(car)
> qqPlot(lm(output ~ group), simulate=T, main = "Q-Q plot")
[1] 28 30
由图，数据落在95%的置信区间范围内，说明满足正态性假设。
#方差齐性检验
> bartlett.test(output ~ group)
Bartlett test of homogeneity of variances
data: output by group
Bartlett's K-squared = 12.139, df = 2, p-value = 0.002312
Bartlett检验表明三组的方差有显著不同（P<0.05）
> outlierTest(aov(output ~ group))
No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferroni p
28 2.410369 0.021857 0.78686
同时，从输出结果看，该数据中没有离群点（P<1）
总结：因数据仅满足了正态性检验，不满足方差齐性检验，故并不能用ANOVA模型拟合出较好的结果
8.4 以小白鼠为对象研究正常肝核糖核酸（RNA）对癌细胞的生物作用，试验分别为对照组（生理盐水），水层RNA 组和酚层RNA 组，分别用此三种不同处理诱导肝癌细胞的果糖二磷酸酯酶（FDP 酶）活力，数据如表7.28 所示，问三种不同处理的诱导作用是否相同？
#数据录入
> method <- as.factor(c(rep("对照组",8),rep("水层RNA",8),rep("酚层RNA",8)))
> result <- 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)
> data8.4 <- data.frame(method,result)
#各组情况
> aggregate(result, by=list(method), FUN=mean)
Group.1 x
1 对照组 2.70250
2 酚层RNA 3.90125
3 水层RNA 3.66125
> aggregate(result, by=list(method), FUN=sd)
Group.1 x
1 对照组 0.5001357
2 酚层RNA 1.0164425
3 水层RNA 0.9850807
> #检验组间差异
> fit8.4 <- aov(result ~ method)
> summary(fit8.4)
Df Sum Sq Mean Sq F value Pr(>F)
method 2 6.437 3.218 4.284 0.0275 *
Residuals 21 15.776 0.751
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
P值<0.05，故拒绝原假设，认为三种不同处理的诱导作用不同
TukeyHSD(fit8.4)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = result ~ method)
$method
diff lwr upr p adj
酚层RNA-对照组 1.19875 0.1064201 2.2910799 0.0298641
水层RNA-对照组 0.95875 -0.1335799 2.0510799 0.0922432
水层RNA-酚层RNA -0.24000 -1.3323299 0.8523299 0.8456673
具体来看，水层和酚层两者间、水层与对照组的诱导作用无显著差异，而酚组与对照组的差异明显
8.5 为研究人们在催眠状态下对各种情绪的反应十分有差异，选取了8 个受试者，在催眠状态下，要求每人按任意次序作出恐惧、愉快、忧虑和平静4 中反映。表7.29 给出了各受试者在处于这4 中清醒状态下皮肤的电位变化值，试在(\alpha=0.05)下，检验受试者在催眠状态下对这4 种情绪的反应力是否有显著差异。
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
Friedman检验对应P值大于0.05，故不拒绝原假设H0，认为受试者在催眠状态下对这4 种情绪的反应力没有显著差异。
8.6 为了提高化工厂的产品质量，需要寻求最优反应温度与反应压力的配合，为此选择如下水平：
A：反应温度（0C） 60 70 80
B：反应压力（kg） 2 2.5 3
在每个AiBj 条件下做两次试验，其产量如表8-12 所示.
（1）对数据作方差分析（应考虑交互作用）；
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.441 2.221 29.83 0.000107 ***
B 2 3.974 1.987 26.69 0.000164 ***
A:B 4 21.159 5.290 71.06 8.34e-07 ***
Residuals 9 0.670 0.074
---
Signif. codes:
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
可以看出，因素A和B，及其交互因素对应P值<<0.05，即对产品质量有显著影响
（2）求最优条件下平均产量的点估计和区间估计；
> tapply(product$X, product$A, mean)
1 2 3
5.150000 4.533333 5.750000
> tapply(product$X, product$B, mean)
1 2 3
5.783333 4.666667 4.983333
> 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.058759 9.141241
sample estimates:
mean of x
6.6
即最优条件下平均产量的点估计为6.6，置信度为95%的区间估计为[4.058759 9.141241]
（3）对AiBj 条件下平均产量作多重比较。
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.6166667 -1.0564835 -0.1768499 0.0089220
3-1 0.6000000 0.1601832 1.0398168 0.0104492
3-2 1.2166667 0.7768499 1.6564835 0.0000777
$B
diff lwr upr p adj
2-1 -1.1166667 -1.5564835 -0.6768499 0.0001517
3-1 -0.8000000 -1.2398168 -0.3601832 0.0017204
3-2 0.3166667 -0.1231501 0.7564835 0.1653096
$`A:B`
diff lwr upr p adj
2:1-1:1 1.85 0.7706087 2.9293913 0.0015185
3:1-1:1 2.15 1.0706087 3.2293913 0.0004844
1:2-1:1 2.05 0.9706087 3.1293913 0.0006999
2:2-1:1 -0.85 -1.9293913 0.2293913 0.1557229
3:2-1:1 -0.55 -1.6293913 0.5293913 0.5677071
1:3-1:1 0.05 -1.0293913 1.1293913 0.9999999
2:3-1:1 -0.75 -1.8293913 0.3293913 0.2502167
3:3-1:1 2.30 1.2206087 3.3793913 0.0002854
3:1-2:1 0.30 -0.7793913 1.3793913 0.9600405
因篇幅原因，交互影响后面的输出结果略
8.7 某良种繁殖场为了提高水稻产量，指定试验的因素如表7.31 所示。试选择L9(34)正交表安排试验，假定相应的产量为（单位：kg/100m2）：
62.925, 57.075, 51.6, 55.05, 58.05, 56.55, 63.225, 50.7, 54.45。
试对试验结果进行方差分析，并给出一组较好的种植条件
设三个因素为A(品种),B(密度),C(施肥量)，采用L9(34)L9(34)正交表安排试验，将ABC放到前三列，得到试验结果表
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.200 56.550 56.125
> tapply(rice$Y, rice$B, mean)
1 2 3
60.400 55.275 54.200
> tapply(rice$Y, rice$C, mean)
1 2 3
56.725 55.525 57.625
> # 可见，A取1，B取1，C取3时，达到最优条件。即最有种植条件为，品种：窄叶青8号；密度：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.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
发现ABC对结果的影响均不显著，说明后续还需要对各因素的选择进行调整或增加实验数据
扫一扫
1241
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
