R语言对完全随机分组实验、拉丁方实验及正交实验进行方差分析（例题，过程+代码）

第一题

题目

研究5种不同的配料（A、B、C、D、E）对某一化学过程反应时间的效应。每批新材料仅够进行5次试验。每次试验大约需要1.5小时，所以一天只能做5次试验。实验者决定用拉丁方来进行实验，使用日期和批次效应可以系统地控制。分析不同配料对反应时间的效应有无显著差异。

解答

读取数据并进行预处理

mydata = read.csv("data.csv")
head(mydata)time <- factor(time)
batch <- factor(batch)
ingredients <- factor(ingredients)

> mydata = read.csv("data.csv")
> head(mydata)time batch ingredients response
1    1     1           A        8
2    1     2           C       11
3    1     3           B        4
4    1     4           D        6
5    1     5           E        4
6    2     1           B        7
>
>
> time <- factor(time)
> batch <- factor(batch)
> ingredients <- factor(ingredients)

读取数据，将时间(time)、批次(batch)和配料(ingredients)进行因子化处理

对数据进行正态性检验

对全体响应(response)进行正态性分布，采用w检验

shapiro.test(response)#检验全体response

> shapiro.test(response)#检验全体responseShapiro-Wilk normality testdata:  response
W = 0.95056, p-value = 0.2581

检验原假设是全体响应满足正态性分布，由w检验可知p>0.05p>0.05p>0.05，接受原假设，response服从正态分布

对数据进行方差齐性检验

采用bartlett方法对不同配方下响应进行方差齐性检验

bartlett.test(response ~ ingredients, data=mydata)

> bartlett.test(response ~ ingredients, data=mydata)Bartlett test of homogeneity of variancesdata:  response by ingredients
Bartlett's K-squared = 1.5544, df = 4, p-value = 0.817

检验原假设是不同配料下的各组响应满足方差齐性，由检验结果可知p>0.05p>0.05p>0.05，接受原假设，不同配料下的各组响应满足方差齐性

对数据进行方差分析

使用aov()函数对数据进行方差分析;检验对象为响应(response),自变量为时间(time)、批次(batch)和配方(ingredients)

mydata.aov.blk = aov(response ~ ingredients + time + batch)
summary(mydata.aov.blk)

> mydata.aov.blk = aov(response ~ ingredients + time + batch)
> summary(mydata.aov.blk)Df Sum Sq Mean Sq F value   Pr(>F)
ingredients  4 141.44   35.36  11.309 0.000488 ***
time         4  12.24    3.06   0.979 0.455014
batch        4  15.44    3.86   1.235 0.347618
Residuals   12  37.52    3.13
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>

结论

由方差分析结果可知，对于原假设H0:不同配料对反应时间的效应无显著差异，而方差分析结果为原假设出现的概率p=0.000488，p<0.05p=0.000488，p<0.05p=0.000488，p<0.05，故不同配料对反应时间的效应有显著差异。同理可知，使用日期和批次效应对反应时间的效应影响不大。

第二题

题目

测量一个化学过程的产量，其中用5批原材料、5种酸性浓度、5种持续时间，以及5种催化剂浓度。分析该实验数据（α=0.05），并写出结论。

解答

读取数据并进行预处理

mydata = read.csv("data.csv")
head(mydata)#attach(mydata)
acid_concentration <- factor(acid_concentration)
batch <- factor(batch)
time <- factor(time)
catalyst_concentration <- factor(catalyst_concentration)

> mydata = read.csv("data.csv")
> head(mydata)acid_concentration batch time catalyst_concentration response
1                  1     1    A                      α       26
2                  1     2    B                      γ       18
3                  1     3    C                      ε       20
4                  1     4    D                      β       15
5                  1     5    E                      δ       10
6                  2     1    B                      β       16
>
> #attach(mydata)
> acid_concentration <- factor(acid_concentration)
> batch <- factor(batch)
> time <- factor(time)
> catalyst_concentration <- factor(catalyst_concentration)

读取数据获得表头，将酸性浓度(acid_concentration)、批次(batch)、时间(time)和催化剂浓度(ingredients)进行因子化处理

对数据进行正态性检验

对全体响应(response)进行正态性分布，做Q-Q图进行检验

qqnorm(response, main = "Normal Q-Q Plot",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles",plot.it = TRUE, datax = FALSE)qqline(response, datax = FALSE, distribution = qnorm,probs = c(0.25, 0.75), qtype = 7,col="red", lwd=2)#lwd=2设置线的宽度

由Q-Q图可知，各数据点散落分布在直线附近，故认为response服从正态分布

对数据进行方差齐性检验

采用bartlett方法对各因子不同水平下的响应进行方差齐性检验

bartlett.test(response ~ acid_concentration, data=mydata)
bartlett.test(response ~ batch, data=mydata)
bartlett.test(response ~ time, data=mydata)
bartlett.test(response ~ catalyst_concentration, data=mydata)

> bartlett.test(response ~ acid_concentration, data=mydata)Bartlett test of homogeneity of variancesdata:  response by acid_concentration
Bartlett's K-squared = 3.5251, df = 4, p-value = 0.4741> bartlett.test(response ~ batch, data=mydata)Bartlett test of homogeneity of variancesdata:  response by batch
Bartlett's K-squared = 1.1667, df = 4, p-value = 0.8835> bartlett.test(response ~ time, data=mydata)Bartlett test of homogeneity of variancesdata:  response by time
Bartlett's K-squared = 0.4853, df = 4, p-value = 0.9749> bartlett.test(response ~ catalyst_concentration, data=mydata)Bartlett test of homogeneity of variancesdata:  response by catalyst_concentration
Bartlett's K-squared = 2.9004, df = 4, p-value = 0.5746

由检验结果可知：各因子不同水平下，p>0.05p>0.05p>0.05，接受原假设，即各因子不同水平下的响应均满足方差齐性

对数据进行方差分析

使用aov()函数对数据进行方差分析;检验对象为响应(response),自变量为时间(time)、批次(batch)和配方(ingredients)

mydata.aov = aov(response ~ acid_concentration + batch + time + catalyst_concentration)
#mydata.aov = summary(mydata.aov)anova.table<-function(fm){tab<-summary(fm)k<-length(tab[[1]])-2temp<-c(sum(tab[[1]][,1]),sum(tab[[1]][,2]),rep(NA,k))tab[[1]]["total",]<-temptab
}
anova.table(mydata.aov)

> mydata.aov = aov(response ~ acid_concentration + batch + time + catalyst_concentration)
> #mydata.aov = summary(mydata.aov)
>
>
> anova.table<-function(fm){
+   tab<-summary(fm)
+   k<-length(tab[[1]])-2
+   temp<-c(sum(tab[[1]][,1]),sum(tab[[1]][,2]),rep(NA,k))
+   tab[[1]]["total",]<-temp
+   tab
+ }
> anova.table(mydata.aov)Df Sum Sq Mean Sq F value   Pr(>F)
acid_concentration      4   24.4    6.10   1.043 0.442543
batch                   4   10.0    2.50   0.427 0.785447
time                    4  342.8   85.70  14.650 0.000941 ***
catalyst_concentration  4   12.0    3.00   0.513 0.728900
Residuals               8   46.8    5.85
total                  24  436.0
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>

结论

由方差分析结果可知，反应持续时间对产量无显著效应的概率p=0.000941，p<0.05p=0.000941，p<0.05p=0.000941，p<0.05，故不同反应持续时间对产量的效应有显著差异。同理可知，酸性浓度、原料批次和催化剂浓度对产量的效应影响不大。

第三题

题目

在对工艺的研究过程中，发现有3个因子很重要，它们是因子A（反应温度），因子B（浓度），因子C（压力）。对每个因子都设定了高低两个水平。我们希望考察这3个因子中，哪些因子的主效应及交互效应是显著的，其具体取值及实验安排如下：

解答

读取数据并进行预处理

mydata = read.csv("data2.csv")
head(mydata)#attach(mydata)
temperature <- as.factor(temperature)
concentration <- as.factor(concentration)
pressure <- as.factor(pressure)

> mydata = read.csv("data2.csv")
> head(mydata)temperature concentration pressure yield
1         160            20        5    60
2         160            20        5    58
3         180            20        5    72
4         180            20        5    73
5         160            40        5    54
6         160            40        5    51
>
> #attach(mydata)
> temperature <- as.factor(temperature)
> concentration <- as.factor(concentration)
> pressure <- as.factor(pressure)

读取数据获得表头，将反应温度(temperature)、浓度(concentration)和压力(pressure)进行因子化处理

对数据进行正态性检验

对全体收率（yield）进行正态性分布，做Q-Q图进行检验

qqnorm(yield, main = "Normal Q-Q Plot",xlab = "Theoretical Quantiles", ylab = "Sample Quantiles",plot.it = TRUE, datax = FALSE)qqline(yield, datax = FALSE, distribution = qnorm,probs = c(0.25, 0.75), qtype = 7,col="red", lwd=2)#lwd=2设置线的宽度

由Q-Q图可知，各数据点散落分布在直线附近，故认为yield服从正态分布

对数据进行方差分析

使用aov()函数对数据进行方差分析;检验对象为收率(yield),考虑因子间的相互作用并制成方差分析表

mydata.aov = aov(yield ~ temperature*concentration*pressure)anova.table<-function(fm){tab<-summary(fm)k<-length(tab[[1]])-2temp<-c(sum(tab[[1]][,1]),sum(tab[[1]][,2]),rep(NA,k))tab[[1]]["total",]<-temptab
}
anova.table(mydata.aov)

> mydata.aov = aov(yield ~ temperature*concentration*pressure)
>
> anova.table<-function(fm){
+   tab<-summary(fm)
+   k<-length(tab[[1]])-2
+   temp<-c(sum(tab[[1]][,1]),sum(tab[[1]][,2]),rep(NA,k))
+   tab[[1]]["total",]<-temp
+   tab
+ }
> anova.table(mydata.aov)Df Sum Sq Mean Sq F value   Pr(>F)
temperature                         1 2328.1  2328.1 760.184 3.23e-09 ***
concentration                       1   68.1    68.1  22.224  0.00151 **
pressure                            1   10.6    10.6   3.449  0.10037
temperature:concentration           1   27.6    27.6   9.000  0.01707 *
temperature:pressure                1  280.6   280.6  91.612 1.18e-05 ***
concentration:pressure              1    0.1     0.1   0.020  0.88994
temperature:concentration:pressure  1    0.6     0.6   0.184  0.67954
Residuals                           8   24.5     3.1
total                              15 2739.9
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
>

结论

由方差分析结果中的p值可知，反应温度、浓度、反应温度和浓度的交互作用与反应温度和压力的交互作用对收率的效应有显著影响，且效应由强到弱分别为：
反应温度>反应温度和压力的交互作用>浓度>反应温度和浓度的交互作用反应温度>反应温度和压力的交互作用>浓度>反应温度和浓度的交互作用反应温度>反应温度和压力的交互作用>浓度>反应温度和浓度的交互作用
而压力、浓度与压力之间的交互作用以及反应温度、浓度和压力三者之间的交互作用均对产率无显著影响,在方差分析过程中可视为误差列。