关于用R语言进行方差分析的乐子二则

我们都知道R语言是非常强的统计工具，今天想着用R语言来做数理统计作业，却被摆了两道，哈哈。

方差分析

要讨论这两则乐子，首先要知道什么是方差分析。

单变量方差分析

单因素方差分析可以看成基础统计中两样本\(t\)检验的一个推广， 要比较试验观测值的某个因变量（称为“指标”）按照一个分组变量（称为“因素”）分组后， 各组的因变量均值有无显著差异。

设因素A有p个不同的水平\(A_1,A_2,\cdots,A_p\),在每个水平下，总体\(X_i\sim N(\mu_i,\sigma^2)\),需要检验的是这\(p\)个样本的平均值有没有显著差异，即： \[ H_0:\mu_1=\mu_2=\cdots=\mu_p \] 一般地，我们在每个水平中抽取\(r\)个样本，形成如下的样本表：

\(A_1\)	\(x_{11}\)	\(x_{12}\)	\(\cdots\)	\(x_{1r}\)
\(A_2\)	\(x_{21}\)	\(x_{22}\)	\(\cdots\)	\(x_{2r}\)
\(\cdots\)			\(\cdots\)
\(A_P\)	\(x_{p1}\)	\(x_{p2}\)	\(\cdots\)	\(x_{pr}\)

【例】：用R语言对下列数据进行方差分析：

A	B	C
58	58	48
64	69	57
55	71	59
66	64	47
67	68	49

【解】先把表格转换成纵列形式

grp	y
A	58
A	64
A	55
A	66
A	67
B	58
B	69
B	71
B	64
B	68
C	48
C	57
C	59
C	47
C	49

将上表保存为R1.csv，在”我的文档“中。在R软件中输入：

data <- read.csv("R1.csv")
data

可以读入数据，并展示数据。

然后建立方差分析模型，在R软件中输入：

mod <- aov(y~grp,data)
summary(mod)

即可输出分析结果，如下：

            Df Sum Sq Mean Sq F value  Pr(>F)   
grp          2    520  260.00   9.176 0.00382 **
Residuals   12    340   28.33                   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我们看到grp的\(P_r\)值有\(0.00382<0.05\)，所以可以认为各组之间有显著差异。

双变量方差分析

在实际任务中，影响实验结果的因素可能不止一个，这时就要用到双变量方差分析（或者正交实验）。在两个因素的实验中，不但每一个因素会对结果起作用，而且两个因素联合起来往往也会对结果起作用。如果你不能理解这个结论，不妨设想这样一个例子：工人在工厂工作时，对工人效率的影响因素有工人和机器。如果工人专业技能高、机器运行状况好，那么效率很自然地就会比较高。但是一个资质平平的工人如果对某个状况一般的机器非常熟悉，以至于把这个机器什么时候会出现什么误差都记得滚瓜烂熟，那么把他俩结合起来，也有可能获得很高的效率。

为了进行双变量方差分析，进行双因素重复性实验，得到的数据表往往是这样给出的：

这时，我就美美地开始做题了！

第一道题

我一看这不简单？于是美美列了个表：

   A B   Y
1  1 1  95
2  2 1  95
3  3 1 106
4  4 1  98
5  5 1 102
6  6 1 112
7  7 1 105
8  8 1  95
9  1 2  95
10 2 2  94
11 3 2 105
12 4 2  97
13 5 2  98
14 6 2 112
15 7 2 103
16 8 2  92
17 1 3  89
18 2 3  88
19 3 3  87
20 4 3  95
21 5 3  97
22 6 3 101
23 7 3  97
24 8 3  90
25 1 4  83
26 2 4  84
27 3 4  90
28 4 4  90
29 5 4  88
30 6 4  94
31 7 4  88
32 8 4  80

读取数据：

X=read.csv("R.csv")

建立模型并分析：

mod<-aov(Y~A+B,X)
summary(mod)

结果：

            Df Sum Sq Mean Sq F value   Pr(>F)    
A            1   74.0    74.0   2.385    0.133    
B            1  926.4   926.4  29.857 6.99e-06 ***
Residuals   29  899.8    31.0                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

其中：

项目	含义
Df	自由度
Sum Sq	平方和
Mean Sq	均方
F value	F值
Pr(>F)	p值

哈哈，聪明的你一定看出问题了吧，那就是这俩的自由度都是\(1\)，这是怎么会事呢？

原来是因为我的两个标签\(A\)和\(B\)都用了数字，而R语言里的数据类型有数字、字符等。它把我们的1，2，3，...当成了一个连续变量，于是就出错了。当然这个问题可以用ac.factor()函数解决，但是我直接把标签换成字母也能解决。

修改了以后的输出为：

            Df Sum Sq Mean Sq F value   Pr(>F)    
A            7  747.5   106.8   12.78 2.62e-06 ***
B            3  977.3   325.8   39.00 9.12e-09 ***
Residuals   21  175.4     8.4                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这个结果是正确的。我们可以发现\(A\)和\(B\)都对结果有不可忽视的影响。

第二道题

三位操作工分别在四台不同的机器上操作了3天，其产量如下表所示：

于是我想当然哈，以为三次实验，加一下就好了，于是列了这么个表：

   A B  Y
1  a x 47
2  b x 51
3  c x 48
4  d x 60
5  a y 54
6  b y 45
7  c y 51
8  d y 48
9  a z 55
10 b z 63
11 c z 54
12 d z 51

建立模型：

mod <- aov(Y~A+B+A:B,X)
summary(mod)

输出结果：

            Df Sum Sq Mean Sq
A            3   8.25    2.75
B            2  81.50   40.75
A:B          6 220.50   36.75

根据\(F\)值，我们可以看出....欸我\(F\)值呢？咋没了？

这是因为，在”加和“（或者求平均值）的时候，就顺带消除了\(A\times B\)的影响。这是因为\(A\times B\)的离差平方和为： \[ S_{A\times B}=\sum_{i=1}^p\sum_{j=1}^q\sum_{k=1}^r (\bar{x}_{ij\cdot} -\bar{x}_{i\cdot \cdot } - \bar{x}_{\cdot j \cdot } +\bar{x})^2 \] 对于多次重复实验，应该把它们分拆开，每次都单独算。列表如下：

   A B  Y
1  a x 15
2  b x 17
3  c x 15
4  d x 18
5  a y 19
6  b y 15
7  c y 18
8  d y 15
9  a z 16
10 b z 19
11 c z 18
12 d z 17
13 a x 15
14 b x 17
15 c x 17
16 d x 20
17 a y 19
18 b y 15
19 c y 17
20 d y 16
21 a z 18
22 b z 22
23 c z 18
24 d z 17
25 a x 17
26 b x 17
27 c x 16
28 d x 22
29 a y 16
30 b y 15
31 c y 16
32 d y 17
33 a z 21
34 b z 22
35 c z 18
36 d z 17

建立模型：

mod <- aov(Y~A*B,X)
summary(mod)

输出结果：

            Df Sum Sq Mean Sq F value   Pr(>F)    
A            3   2.75   0.917   0.532 0.664528    
B            2  27.17  13.583   7.887 0.002330 ** 
A:B          6  73.50  12.250   7.113 0.000192 ***
Residuals   24  41.33   1.722                     
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

我们可以看到，效率主要和工人的水平以及工人和机器的交互有关。

哈哈，我是铸币。