R语言对正交实验结果（含交互作用）进行极差分析与方差分析实例

题目

某工厂为了提高某产品的收率，根据经验和分析，认为反应温度A、反应时间B、碱用量C和催化剂种类D可能对产品的收率造成较大的影响,并考虑交互作用AB，AC。用正交表L8(27)安排试验，试验方案及结果如下表所示，试用直观分析和方差分析来分析结果。

解答

1 数据预处理

读取数据，对各因子列数据进行因子化处理，并将最后一列y对应值赋值给response

mydata = read.csv("data2.csv",fileEncoding = "UTF-8-BOM")
head(mydata)

attach(mydata)
for (i in 3:ncol(mydata)-1) {
  mydata[,i] <- as.factor(mydata[,i])
}

response = mydata[,ncol(mydata)]

> mydata = read.csv("data2.csv",fileEncoding = "UTF-8-BOM")
> head(mydata)
  试验号 A B A.B C A.C 空列 D 收率
1      1 1 1   1 1   1    1 1   65
2      2 1 1   1 2   2    2 2   74
3      3 1 2   2 1   1    2 2   71
4      4 1 2   2 2   2    1 1   73
5      5 2 1   2 1   2    1 2   70
6      6 2 1   2 2   1    2 1   73
> 
> #attach(mydata)
> for (i in 3:ncol(mydata)-1) {
+   mydata[,i] <- as.factor(mydata[,i])
+ }
> 
> response = mydata[,ncol(mydata)]
> 

2 直观分析法

分别计算每一个因子对应的极差（各水平下最大均值-最小均值），将计算结果存于数据框K中

k.max = c(NULL)#存最大效应值
k.min = c(NULL)#存最小效应值
k.maxlevel = c(NULL)#最大效应值对应水平
J = 1
for (i in 3:ncol(mydata)-1) {
  KK = tapply(response, mydata[,i], mean)
  k.max[J] = each(max)(KK)
  k.min[J] = each(min)(KK)
  k.maxlevel[J] = which(KK==k.max[J])
  J=J+1
}
k.mean = k.max - k.min#计算极差
factors = colnames(mydata)[3:ncol(mydata)-1]
K <- data.frame(factors,k.mean, k.maxlevel, k.max, k.min) ##生成数据框
K

> k.max = c(NULL)#存最大效应值
> k.min = c(NULL)#存最小效应值
> k.maxlevel = c(NULL)#最大效应值对应水平
> J = 1
> for (i in 3:ncol(mydata)-1) {
+   KK = tapply(response, mydata[,i], mean)
+   k.max[J] = each(max)(KK)
+   k.min[J] = each(min)(KK)
+   k.maxlevel[J] = which(KK==k.max[J])
+   J=J+1
+ }
> k.mean = k.max - k.min#计算极差
> factors = colnames(mydata)[3:ncol(mydata)-1]
> K <- data.frame(factors,k.mean, k.maxlevel, k.max, k.min) ##生成数据框
> K
  factors k.mean k.maxlevel k.max k.min
1       A   2.75          1 70.75 68.00
2       B   2.25          1 70.50 68.25
3     A.B   4.75          2 71.75 67.00
4       C   4.75          2 71.75 67.00
5     A.C   0.75          2 69.75 69.00
6    空列   1.25          2 70.00 68.75
7       D   2.25          2 70.50 68.25
>

由计算结果可知，若只考虑极差的效果，则各因素对响应值的效应主次排序为：
$$C=A\times B>A>B=D>空列>A\times C$$
绘制因子各水平作用效应均值图：

#只有当正交实验中各因子为两水平时才可以调用以下代码绘制图形
library(DoE.base)
library(FrF2)
lm.model = lm(response~A+B+A.B+C+A.C+空列+D,data=mydata)
MEPlot(lm.model, abbrev = 2, response