方差分析

魔法师LQ

1.问题呈现

下表给出了某种化工产品在3中浓度,4中温度水平下得率的数据,试检验各因素及交互作用对产品得率的影响是否显著。

浓度/% 10℃ 24℃ 38℃ 52℃
2 14,10 11,11 13,9 10,12
4 9,7 10,8 7,11 6,10
6 5,11 13,14 12,13 14,10

2.问题分析和实验设计

2.1 问题分析

由分析可知,此问题属于有交互作用的(可重复)双因素方差分析,由原理我们可以知道,误差的组成有以下关系:

其中$a, b$分别代表因素$\mathrm{A, B}$的水平数,$r$表示重复试验的次数,$X_{ijk}$表示实验值,$X_{\dots}$表示所有实验值的总平均值,有:

$X_{ij.}$表示组合水平$\left( \mathrm{A_i, B_j}\right)$上的$r$次实验值的算数平均值,有:

$X_{i..}$表示$\mathrm{A_i}$水平时实验值的算数平均值:

$X_{.j.}$表示$\mathrm{B_j}$水平时实验值的算数平均值:

随后根据方差分析表进行计算:

方差分析表

Source(误差源) Df(自由度) SS(平方和) MS(均方) F(F统计量) Signif.(显著性)
$\mathrm{A}$ a-1 $SS_A$ $MS_A$ $MS_A/MS_{within}$
$\mathrm{B}$ b-1 $SS_B$ $MS_B$ $MS_B/MS_{within}$
$\mathrm{A\times B}$ (a-1)(b-1) $SS_{A\times B}$ $MS_{A\times B}$ $MS_{A\times B}/MS_{within}$
$\mathrm{Within}$ ab(r-1) $SS_{within}$ $MS_{within}$
$\mathrm{Total}$ abr-1 $SS_{Total}$

计算方法见公式(1)

2.2 问题计算

由分析我们可以对问题直接进行求解:

首先进行各个均值的计算:

浓度/% 10℃ 24℃ 38℃ 52℃ $m_A$
2 14,10(12.00) 11,11(11.00) 13,9(11.00) 10,12(11.00) 11.25
4 9,7(8.00) 10,8(9.00) 7,11(9.00) 6,10(8.00) 8.50
6 5,11(8.00) 13,14(13.50) 12,13(12.50) 14,10(12.00) 11.50
$m_B$ 9.33 11.17 10.83 10.33 10.42

进而有如下计算:

  • 计算总离差平方和$SS_{Total}$
  • 计算$SS_A$, $df_A$和$MS_A$
  • 计算$SS_B$, $df_B$和$MS_B$
  • 计算误差平方和$SS_{within}$, $df_{within}$和$MS_{within}$
  • 计算$SS_{A\times B}$, $df_{AB}$和$MS_{AB}$

随后进行F检验,有:

进而计算出F统计量,以及P值,从而在一定显著性水平下,对原假设和备择假设进行是否拒绝的判断。以上就是我们对于问题的原理和计算过程,下面开始进行相关的实验设计和编程实现。

2.3 实验设计

2.3.1 实验目的

实现对于给定问题但不限于给定问题和其他同等问题的计算机自动处理和方差分析。

2.3.2 实验准备

实验平台:

  • Windows 10 Professional 64-bit (10.0, Build 16299), 8192MB RAM。
  • Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz (8 CPUs), ~2.8GHz。

实验软件:

  • RStudio Version 1.1.456
  • R3.5.1
  • Microsoft Excel

3 基于Excel的双因素有交互作用的方差分析

3.1 数据预处理

根据问题,我们将数据以Excel表格的形式进行录入,并保存为Excel工作簿(*.xlsx)格式。经过整理可以得到汇总后的工作表:

接下来使用Excel自带的数据分析工具对原始数据进行有交互的双因素方差分析。

3.2 关键步骤

  1. 点击Excel的左上角的文件,转到侧边栏后点击最下方的选项

  2. 进入选项-加载项-分析工具库, 点击确认开启

  3. 回到工作区,点击数据-数据分析

  4. 进入选择可重复双因素分析

  5. 选择输入区域,设置每一样本行数为2,设置$\alpha(\underline A)$为0.05或者0.01等

  6. 最后点击确定

3.3 分析结果

得到分析结果如下:

和我们分析和计算的结果完全一致。

由此我们可以得出:

从方差分析表中可以看出,行(样本)因素及浓度的F统计量为4.09,大于在$\alpha=0.05$显著性条件下的F临界值,说明这时应该拒绝原假设,且此时得到的P值为0.04<0.05,说明浓度对产品得率有显著影响。

由此同时,列和交互因素的F值均小于临界值,而P值较大(大于0.5),说明其对试验结果的影响不显著。

3 基于R内置函数的有交互作用的双因素方差分析

R有内置的用于方差分析的函数aov()可用于本任务。

1
> ?aov

查看使用说明:

功能上是通过调用lm,即提供了lm函数(线性模型进行拟合函数)的封装,从而对每一层进行方差分析。

3.1 数据预处理

我们进行的是有交互作用的双因素方差分析,原始的二维数据表并不适合我们使用R对其进行方差分析。这里首先需要我们手动对数据进行预处理,转换成三列,前两列代表两个因素,即浓度和温度,最后一列代表观测值,如图所示:

即将原始数据表每个点所对应的数据集作为转换后的数据帧(或叫数据框)前两列,例如第一列是第一个元素的值,第二列表示第二个因素的值,相当于是二维的坐标值,记录值则处于最后一列。

接着,将预处理后的数据保存为方便R进行读取和处理的.csv格式,并命名为product_data.csv

3.2 数据加载

使用read.csv()来加载product_data.csv:

3.3 查看数据

class(product_data)结果表明和我们预想的一样,数据的类型是data-frame类型的。使用head()查看数据的前几行,使用summary()产看数据的基本信息,包含每一列的最大值,最小值,平均值,中位数,1/4位点和3/4位点。

3.4 方差分析

使用内置aov()函数,输入如下:

1
> mX_model <- aov(data = product_data, formula = result~concentration*temperature)

其中~符号左侧表示因变量,右侧表示自变量,从而构成了一个formula。其中*标志两因素之间有交互作用,若要改成没有交互作用的双因素方差分析,可以改写成:

1
> mX_model <- aov(data = product_data, formula = result~concentration*temperature)

将结果保存到mX_model变量中。

接着输入summary(mX_model),即可查看分析结果。

但是很明显,即使是最容易验证的Df,程序也算错了。

原因在于,read.csv,将我们的观测值全部按照数值类型进行读取

而实际上我们的因素concentration(浓度)temperature(温度)均应该设置成Factor类型。操作和结果如下:

1
2
3
4
5
> product_data$concentration <- as.factor(product_data$concentration)
> product_data$temperature <- as.factor(product_data$temperature)
> sapply(product_data, class)
concentration temperature result
"factor" "factor" "integer"

接着重新执行之前的步骤:

可以看出和之前的结果有冥想的不同,显示了factor类型的水平以及记录的个数。

紧接着执行最关键的一步:

1
> mX_model2 <- aov(data = product_data, result~concentration*temperature)

即可得到有交互作用的双因素方差分析结果。

3.5 分析结果

1
2
3
4
5
6
7
8
> summary(mX_model2)
Df Sum Sq Mean Sq F value Pr(>F)
concentration 2 44.33 22.167 4.092 0.0442 *
temperature 3 11.50 3.833 0.708 0.5657
concentration:temperature 6 27.00 4.500 0.831 0.5684
Residuals 12 65.00 5.417
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

得到正确的分析结果,从结果可以看出,在显著性条件为0.05下,浓度对产品得率有显著影响,而温度,交互作用对于得率的影响并不显著。

为了使结果更加直观地呈现,下面进行数据的可视化:

1
2
> par(mar=c(3, 4, 4, 0.5))
> with(product_data, interaction.plot(concentration, temperature, result, type = "b", col = c("green", "blue", "black", "red"), pch = c(16, 9), main="Interaction Plot for Concentration and Temperature", ylab="Rate of Production"))

得到两个因素交互(分别处于不同水平)的散点图:

同时,还可以做出依照对结果有显著影响的浓度因素分组的结果分布箱线图:

1
> boxplot(result~concentration*temperature, data=product_data, col=(c("gold", "green", "blue")), ylab="Rate of Prodution", main="Different Concentration")

可以直观地看出箱线图的中位数总体随着浓度的升高而升高,照应了我们的分析结果。

4 基于R自定义函数的有交互作用的双因素方差分析

通过之前的实验我们以简便的方式均得到了实验分析的正确结果,此外还有了一些新的启发。但是我们也发现,无论还是Excel工具还是R内置函数,都有很大的局限性已经不灵活,例如:

  • Excel完全是黑盒,我们并不能看到内部运行的机;
  • R内置函数则在处理数据,特别是输入时,不够自然;
  • 且两者都有一些我们想看到的信息,却没能显示出来,比如R中没有显示总离差平方和,Excel中没有基于我们限制性水平提示等。

因此,下面我们基于R语言,实现自定义方差分析函数。

4.1 数据表读取和转换

首先我们希望原始数据的输入应该是自然的,也就是问题呈现给我们的样子:

所以我们希望R就以这种方式自动读取,过程和结果如下:

读取后,显然R不能直接进行处理,因此我们希望实现函数能够实现将原始的数据表自动转换成章节3.1当中的数据帧(框),因此编写编写下面自定义函数:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 功能:将原始的数据表转换为dataframe: 
#
# 参数:
# raw_tbl: 输入数据表
# 举例:(约定重复试验观察值使用“,”号来分隔)
# 浓度 10℃ 24℃ 38℃ 52℃
# 2 14,10 11,11 13,9 10,12
# 4 9,7 10,8 7,11 6,10
# 6 5,11 13,14 12,13 14,10
# 表示重复次数为2,因素A为浓度,水平有(2,4,6), 因素B为温度,水平有(10, 24, 38, 52)
# r: 重复次数
# (r表示重复次数,这里强制人为输入,并且程序也自动进行检测,假如有不一致则抛出异常)
#
# 输出:
# 表格转换成data-frame:
# 转换后的dataframe行数为a*b*r(其中a,b分别表示两个因素的水平,r表示重复次数),
# 第一列为因素A,第二列为因素B,第三列为观察值。

核心代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
tb2df <- function(raw_tbl, r){
# 最终要转化成的能够处理的dataframe, 初始化为空的
data_frame <- data.frame("A"= character(0), "B"=character(0), "R"=character(0), stringsAsFactors=FALSE)
#names(data_frame) <- c("A", "B", "R") # 第一列为因素A, 随后是B,最后是观察值

fA <- factor(raw_tbl[[1]]) # 第一列为因素A
level_A <- levels(fA)
fB <- factor(names(raw_tbl[-1])) # 行名称为因素B(除了第一列为空)
level_B <- levels(fB)

num_fA <- length(levels(fA))
num_fB <- length(levels(fB))

# 先列后行
for ( j in 2:(num_fB+1)){ # 加1为了跳过第一列

df_B <- level_B[j-1] # 因素B的水平,作为新dataframe的第二列
#prt("level B=", df_B)
for (i in 1:num_fA){
df_A <- level_A[i] # 因素A的水平,作为新dataframe的第一列
#prt("level A=", df_A)

#print(raw_tbl[[j]][i])
c <- unlist(strsplit(as.character(raw_tbl[[j]][i]), split = ","))
#print(c)
r_ <- length(c)
# prt("len=", r_)
if (r_!=r){return(NA)} # 处理不一致情况
for (k in c){
#print(k)
# 每次创建一个一行的dataframe
data_row <- list(df_A, df_B, k)
#print(data_row)
# 每次向总data_frame中插入一行
# 参考https://stackoverflow.com/questions/28467068/add-row-to-dataframe
data_frame[nrow(data_frame)+1, ] <- data_row
}
}
}

data_frame[, 1] <- as.factor(data_frame[, 1]) # 将第一列转换成因子类型
data_frame[, 2] <- as.factor(data_frame[, 2]) # 将第二列转换成因子类型
data_frame[, 3] <- as.numeric(data_frame[, 3]) # 将第二列转换成数值类型

data_frame
}

经过转化我们得到结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
> tb2df(table_raw, 2)
A B R
1 2 X10. 14
2 2 X10. 10
3 4 X10. 9
4 4 X10. 7
5 6 X10. 5
6 6 X10. 11
7 2 X24. 11
8 2 X24. 11
9 4 X24. 10
10 4 X24. 8
11 6 X24. 13
12 6 X24. 14
13 2 X38. 13
14 2 X38. 9
15 4 X38. 7
16 4 X38. 11
17 6 X38. 12
18 6 X38. 13
19 2 X52. 10
20 2 X52. 12
21 4 X52. 6
22 4 X52. 10
23 6 X52. 14
24 6 X52. 10

前两列是因素,最后一列是观察值。

4.2 标准函数实现

依据原理我们使用R对双因素方差分析函数进行实现。

4.2.1 函数接口

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 功能:进行双因素有交互方差分析

# 输入: 经过转化后的dataframe
# 参数:
# - data_frame:输入的datafram
# - r:表示重复的次数
# - alpha, 显著性水平,默认是0.01,可以设定成其他
# 输出:
# - tbl: 方差分析表
# - 其他:诸如各种平均值
anova_std <- function(data_frame, r, alpha=0.01){
...
return(tbl)
}

4.2.2 均值计算

1
2
3
4
mean_T <- mean(data_frame[, 3]) # 总的均值
mean_A <- tapply(data_frame[, 3], data_frame[, 1], mean)
mean_B <- tapply(data_frame[, 3], data_frame[, 2], mean)
mean_AB <- tapply(data_frame[, 3], data_frame[, 1]:data_frame[, 2], mean) # 组合水平

4.2.3 离差平方和

1
2
3
4
5
SSE <- sum(tapply(data_frame[, 3], data_frame[, 1]:data_frame[, 2], function(x) var(x)*(length(x)-1)))  # 误差平方和, within
SSA <- r*num_fB*sum(tapply(data_frame[, 3], data_frame[, 1], function(x) (mean(x)-mean_T)^2))
SSB <- r*num_fA*sum(tapply(data_frame[, 3], data_frame[, 2], function(x) (mean(x)-mean_T)^2))
SST <- sum(sapply(data_frame[, 3], function(x) (x-mean_T)^2))
SSAB <- SST - SSA - SSB - SSE # 也可以根据定义是来算SSAB

4.2.4 自由度

1
2
3
4
5
df_within <- (r-1)*num_fA*num_fB
df_A <- num_fA - 1
df_B <- num_fB - 1
df_T <- r*num_fA*num_fB-1
df_AB <- (num_fA-1)*(num_fB-1)

4.2.5 均方

1
2
3
4
MSE <- SSE/df_within
MSA <- SSA/df_A
MSB <- SSB/df_B
MSAB <- SSAB/df_AB

4.2.6 F统计量

1
2
3
4
# F
FA <- MSA/MSE
FB <- MSB/MSE
FAB <- MSAB/MSE

4.2.7 P-value

1
2
3
4
pr_A <- pf(q=FA, df1=df_A, df2=df_within, lower.tail = FALSE)
pr_B <- pf(q=FB, df1=df_B, df2=df_within, lower.tail = FALSE)
pr_AB <- pf(q=FAB, df1=df_AB, df2=df_within, lower.tail = FALSE)
pr <- c(pr_A, pr_B, pr_AB, NA, NA)

4.2.8 临界值

1
2
3
4
5
# 求临界值 q代表,quantile
qf_A <- qf(1-alpha, df1=df_A, df2=df_within)
qf_B <- qf(1-alpha, df1=df_B, df2=df_within)
qf_AB <- qf(1-alpha, df1=df_AB, df2=df_within)
qf_ <- c(qf_A, qf_B, qf_AB, NA, NA)

4.2.9 显著性水平

为了让显著性更能直观表现,定义下面函数:

1
2
3
4
5
6
7
8
9
10
11
# 辅助函数
# 依据输出显著性水平标志,"*"越多越显著
signif_code <- function(pr){
if (is.na(pr) | as.numeric(pr)<0 | as.numeric(pr) >1 ) {return("NA")}

else if (as.numeric(pr)<0.001) return("***")
else if (as.numeric(pr)<0.01) return("**")
else if (as.numeric(pr)<0.05) return("*")
else if (as.numeric(pr)<0.1) return(".")
else if (as.numeric(pr)<1) return(" ")
}

完整代码见附录A

4.3 简化公式实现

除了标准的计算外,我们还能根据推导,采用简化后的公式来进行运算。

具体做法如下表,先计算多种试验值(或称观察值)之和整理如下:

接着根据下列公式进行计算:

和标准做法相比,简化公式避免了计算平均值,理论上来讲比标准做法复杂度低一些。

下面我们进行简化公式的设计和实现。

4.3.1 函数接口

设计和标准做法统一的接口。

1
2
3
4
5
6
# 方法2. 简化公式
# 参数:见标准做法
anova_simplify <- function(data_frame, r, alpha=0.01){
...
return(tbl)
}

4.3.2 试验值之和

1
2
3
4
5
SUM_AB <- tapply(data_frame[, 3], data_frame[, 1]:data_frame[, 2], sum)
SUM_A <- tapply(data_frame[, 3], data_frame[, 1], sum)
SUM_B <- tapply(data_frame[, 3], data_frame[, 2], sum)
SUM_T <- sum(tapply(data_frame[, 3], data_frame[, 2], sum)) # 所有项算数总和
# SUM_T <- sum(tapply(data_frame[, 3], data_frame[, 1], sum)) # 计算结果一样

4.3.3 试验值平方和

1
2
3
4
SS_Ai <- sum(sapply(SUM_A, function(x) x^2))
SS_Bj <- sum(sapply(SUM_B, function(x) x^2))
SS_ABij <- sum(sapply(SUM_AB, function(x) x^2))
SS_ijk <- sum(sapply(data_frame[, 3], function(x) x^2)) # 总的所以观测值的平方和

4.3.4 离差平方和

1
2
3
4
5
SST <- SS_ijk - SUM_T^2/num_n 
SSA <- SS_Ai/(num_fB*r)-SUM_T^2/num_n
SSB <- SS_Bj/(num_fA*r)-SUM_T^2/num_n
SSE <- SS_ijk - SS_ABij/r
SSAB <- SST - SSA -SSB-SSE

后面的运算和假设检验和标准做法相同。

4.3 运行结果

运行结果如下:

两种计算结果完全一致,且和我们之前的分析以及实验结果完全一致。

对于两种算法运行效率,由于运行时间受多种因素控制,且问题的数据规模非常小,所以结果并没有显著差异。

4.4 更多例子

这里有一个参考文献[2]中的例子,有两种洗衣粉,以及三种不同的水温,每种情况下重复四次实验,试检验各因素对于洗涤效果(洗出的污渍量)的显著性影响。

读入和转换:

方差分析和结果:

1
2
3
4
5
6
7
8
9
10
> test_df <- tb2df(test_tbl, 4)
> anova_std(test_df, 4)
差异源 平方和SS 自由度df 均方MS F统计量 P值 显著性 临界值(α=0.01)
1 A 20.16667 1 20.166667 9.810811 5.758440e-03 ** 8.285420
2 B 200.33333 2 100.166667 48.729730 5.439849e-08 *** 6.012905
3 A:B 16.33333 2 8.166667 3.972973 3.722434e-02 * 6.012905
4 Residuals 37.00000 18 2.055556 NA NA NA NA
5 Total 273.83333 23 NA NA NA NA NA
[1] "---"
[1] "注意:显著性标志码(Signif. codes): [0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1]"

使用anova_simplify得到相同的分析结果,这里不再重复。

5.总结和展望

本文主要围绕双因素有交互方差分析并以某化工产品得率收浓度和温度影响的问题为例,主要围绕两方面对问题进行了解决,第一方面是依据原理对问题进行了分析和手工计算。

第二个大的方面是针对问题进行并以多种方式开展了计算机自动处理的实验:其中实验的第一和第二部分是分别基于Excel自带数据分析工具,以及R语言内置方差分析函数对问题进行了解决,并对结果进行了分析和可视化,得出了问题中浓度对于得率的影响显著,而温度以及交互对于实验结果影响不显著;实验的第三部分,实现了两种自定义的方差分析函数,一种是标准公式的计算方法,另一种是简化公式的计算方法,此外针对之前方法存在的输入数据不自然,输出结果部分关键指标缺失,显示不友好等不足进行了改进。结果表明两种方法均能够更自然地输入原始数据并进行处理,并得到了相同的且用户友好的正确结果。对比发现实际两种方法在处理小规模数据上性能并无明显差异。

对于未来的工作,虽然文中实现了双因素的有交互方差分析,但是对于更简单的单因素方差分析以及双因素无交互方差分析等没有自定义实现,此外可以就如何提高分析效率和效果,扩展性等方面对本文方法进行改进。

声明

作者保留关于此报告的一切权利。

最后更新日期:2018年10月26日

参考

[1] How do I get R to spit out the critical value for F-statistic based on ANOVA?.Stackoverflow

[2] Example of Doing Two way ANOVA.[w].Introduction to Statistics for Biology and Biostatistics.Standord.s141

[3] Swirl

[4] Add row to dataframe.Stackoverflow

[5] ANalysis+Of+VArianceV2.张凯.实验方差分析.

附录

附录A

标准版本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
# 辅助函数
# 依据输出显著性水平标志,"*"越多越显著
signif_code <- function(pr){
if (is.na(pr) | as.numeric(pr)<0 | as.numeric(pr) >1 ) {return("NA")}

else if (as.numeric(pr)<0.001) return("***")
else if (as.numeric(pr)<0.01) return("**")
else if (as.numeric(pr)<0.05) return("*")
else if (as.numeric(pr)<0.1) return(".")
else if (as.numeric(pr)<1) return(" ")
}

# 功能:进行双因素有交互方差分析
# 版本:2018-10-25 @QL @THU.sz
# 输入: 经过转化后的dataframe
# 参数:
# data_frame:输入的datafram
# r:表示重复的次数
# alpha, 显著性水平,默认是0.01,可以设定成其他
# 输出:
# tbl: 方差分析表
anova_std <- function(data_frame, r, alpha=0.01){
# 方法1.标准做法
# 1. 各种平均值
# 2. 计算各种平方和
# 3. 计算各种均方
# 4. 得到F值
# 5. 对比临界值,得到分析结果
# 6. 其他计算量,显著性水平表示
# 7. 返回计算表
num_fA <- length(levels(data_frame[[1]]))
num_fB <- length(levels(data_frame[[2]]))

mean_T <- mean(data_frame[, 3]) # 总的均值
mean_A <- tapply(data_frame[, 3], data_frame[, 1], mean)
mean_B <- tapply(data_frame[, 3], data_frame[, 2], mean)
mean_AB <- tapply(data_frame[, 3], data_frame[, 1]:data_frame[, 2], mean)

SSE <- sum(tapply(data_frame[, 3], data_frame[, 1]:data_frame[, 2], function(x) var(x)*(length(x)-1))) # 误差平方和, within
df_within <- (r-1)*num_fA*num_fB
MSE <- SSE/df_within

SSA <- r*num_fB*sum(tapply(data_frame[, 3], data_frame[, 1], function(x) (mean(x)-mean_T)^2))
df_A <- num_fA - 1
MSA <- SSA/df_A

SSB <- r*num_fA*sum(tapply(data_frame[, 3], data_frame[, 2], function(x) (mean(x)-mean_T)^2))
df_B <- num_fB - 1
MSB <- SSB/df_B

SST <- sum(sapply(data_frame[, 3], function(x) (x-mean_T)^2))
df_T <- r*num_fA*num_fB-1

SSAB <- SST - SSA - SSB - SSE # 也可以根据定义是来算SSAB
df_AB <- (num_fA-1)*(num_fB-1)
MSAB <- SSAB/df_AB

# F
FA <- MSA/MSE
FB <- MSB/MSE
FAB <- MSAB/MSE

# 返回
# 方差分析表
na <- c(c(names(data_frame)[1:2]), paste(names(data_frame)[1],names(data_frame)[2], sep = ":"), "Residuals", "Total")
sr <- c(SSA, SSB, SSAB, SSE, SST)
df <- c(df_A, df_B, df_AB, df_within, df_T)
ms <- c(MSA, MSB, MSAB, MSE, NA)
f <- c(FA, FB, FAB, NA, NA)

pr_A <- pf(q=FA, df1=df_A, df2=df_within, lower.tail = FALSE)
pr_B <- pf(q=FB, df1=df_B, df2=df_within, lower.tail = FALSE)
pr_AB <- pf(q=FAB, df1=df_AB, df2=df_within, lower.tail = FALSE)
pr <- c(pr_A, pr_B, pr_AB, NA, NA)

sig <- sapply(pr, signif_code) # 显著性标志
# 求临界值 q代表,quantile
qf_A <- qf(1-alpha, df1=df_A, df2=df_within)
qf_B <- qf(1-alpha, df1=df_B, df2=df_within)
qf_AB <- qf(1-alpha, df1=df_AB, df2=df_within)
qf_ <- c(qf_A, qf_B, qf_AB, NA, NA)


tbl <- data.frame(TYPE=na ,SS=sr, df=df, MS=ms, F=f, Pr=pr, Sig=sig, QF=qf_)

# 命名
names(tbl) <- c("差异源", "平方和SS", "自由度df", "均方MS", "F统计量", "P值","显著性", paste("临界值(α=", alpha, ")", sep=""))

# 打印到屏幕
print(tbl)
print("---")
# 显著性标志提示
print("注意:显著性标志码(Signif. codes): [0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1]")


# 返回值
return(tbl)
}

简化公式版本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

# 方法2. 简化公式
# 参数:见标准做法
anova_simplify <- function(data_frame, r, alpha=0.01){
num_fA <- length(levels(data_frame[[1]]))
num_fB <- length(levels(data_frame[[2]]))
num_n <- length(dat[, 3]) # 所有记录的个数

# 2018-10-24 原理差不多,这次就先不实现了
# 2018-10-26 决定还是实现一下
SUM_AB <- tapply(data_frame[, 3], data_frame[, 1]:data_frame[, 2], sum)
SUM_A <- tapply(data_frame[, 3], data_frame[, 1], sum)
SUM_B <- tapply(data_frame[, 3], data_frame[, 2], sum)
SUM_T <- sum(tapply(data_frame[, 3], data_frame[, 2], sum))
# SUM_T <- sum(tapply(data_frame[, 3], data_frame[, 1], sum)) # 计算结果一样

SS_Ai <- sum(sapply(SUM_A, function(x) x^2))
SS_Bj <- sum(sapply(SUM_B, function(x) x^2))
SS_ABij <- sum(sapply(SUM_AB, function(x) x^2))
SS_ijk <- sum(sapply(data_frame[, 3], function(x) x^2)) # 总的所以观测值的平方和

SST <- SS_ijk - SUM_T^2/num_n
SSA <- SS_Ai/(num_fB*r)-SUM_T^2/num_n
SSB <- SS_Bj/(num_fA*r)-SUM_T^2/num_n
SSE <- SS_ijk - SS_ABij/r
SSAB <- SST - SSA -SSB-SSE

df_within <- (r-1)*num_fA*num_fB
df_A <- num_fA - 1
df_B <- num_fB - 1
df_T <- r*num_fA*num_fB-1
df_AB <- (num_fA-1)*(num_fB-1)

MSE <- SSE/df_within
MSA <- SSA/df_A
MSB <- SSB/df_B
MSAB <- SSAB/df_AB

# F
FA <- MSA/MSE
FB <- MSB/MSE
FAB <- MSAB/MSE

# 返回
# 方差分析表
na <- c(c(names(data_frame)[1:2]), paste(names(data_frame)[1],names(data_frame)[2], sep = ":"), "Residuals", "Total")
sr <- c(SSA, SSB, SSAB, SSE, SST)
df <- c(df_A, df_B, df_AB, df_within, df_T)
ms <- c(MSA, MSB, MSAB, MSE, NA)
f <- c(FA, FB, FAB, NA, NA)

pr_A <- pf(q=FA, df1=df_A, df2=df_within, lower.tail = FALSE)
pr_B <- pf(q=FB, df1=df_B, df2=df_within, lower.tail = FALSE)
pr_AB <- pf(q=FAB, df1=df_AB, df2=df_within, lower.tail = FALSE)
pr <- c(pr_A, pr_B, pr_AB, NA, NA)

sig <- sapply(pr, signif_code) # 显著性标志
# 求临界值 q代表,quantile
qf_A <- qf(1-alpha, df1=df_A, df2=df_within)
qf_B <- qf(1-alpha, df1=df_B, df2=df_within)
qf_AB <- qf(1-alpha, df1=df_AB, df2=df_within)
qf_ <- c(qf_A, qf_B, qf_AB, NA, NA)


tbl <- data.frame(TYPE=na ,SS=sr, df=df, MS=ms, F=f, Pr=pr, Sig=sig, QF=qf_)

# 命名
names(tbl) <- c("差异源", "平方和SS", "自由度df", "均方MS", "F统计量", "P值","显著性", paste("临界值(α=", alpha, ")", sep=""))

# 打印到屏幕
print(tbl)
print("---")
# 显著性标志提示
print("注意:显著性标志码(Signif. codes): [0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1]")


# 返回值
return(tbl)
}

用于演示的函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
# 参数:
# - wh:
# 选择1:(奇数)表示标准算法,默认选项
# 选择2:(偶数)表示用简化公式计算

run <- function(wh=2){
tbl_raw <- read.csv(in_path)
summary(tbl_raw)
data_frame <- tb2df(tbl_raw, r)
summary(data_frame)

t1 <- Sys.time()
if (wh%%2==1){
anova_result <- anova_std(data_frame, r, alph)
}else{
anova_result <- anova_simplify(data_frame, r, alph)
}
t2 <- Sys.time()
print(t2-t1)

# 写到文件
# write.csv(anova_result, "../out/anova_result.csv")
# write.csv(data_frame, "../out/tbl_2_df.csv")
# write.csv(tbl_raw, "../out/tbl_raw.csv")
}
run()

附录B

Swirl教程: R Programming学习笔记

请见作者博客: R介绍