一、简介
方差分析(ANOVA)是检验多个总体均值是否相等的方法。虽然它形式上是比较总体均值,但是本质上是研究变量之间的关系。这里的变量中,自变量是分类型的,因变量是数值型的,所研究的关系是是指自变量对因变量的影响。因此,我们可以说:
方差分析是通过检验各个总体均值是否相等来判断分类型自变量对数值型自变量是否有显著影响。
在方差分析中,所要检验的对象称为因子(factor),因子的不同表现称为水平(treatment),因子的每一个水平都可以看作一个总体,每个因子水平下得到样本数据称为观测值。
比如,为了比较药物
A
,
B
,
C
A,B,C
A,B,C对治疗某疾病的疗效,我们将实验对象分成三组,分别记录服用三种药物的治疗效果,得到三组样本:
X
1
,
…
,
X
n
1
;
Y
1
,
…
,
Y
n
2
;
Z
1
,
…
,
Z
n
3
X_1,\dots,X_{n_1};\ Y_1,\dots,Y_{n_2};\ Z_1,\dots,Z_{n_3}
X1,…,Xn1; Y1,…,Yn2; Z1,…,Zn3
这个例子中,药物称为因子,
A
,
B
,
C
A,B,C
A,B,C称为该因子的水平。
由于这个实验只涉及单个因子——“药物”,我们称之为单因子实验,对应的方差分析也叫单因子方差分析。此外,如果比较不同的药物和性别对疗效的影响,这就是两因子实验,相应的方差分析叫做两因子方差分析。不难推广到多因子实验。
二、基本思想
因子的不同水平下的均值会有差异,但这种差异有可能是由抽样误差带来的,所以需要检验这种差异是否显著。虽然我们感兴趣的是均值,但我们在判断时需要借助于方差(构造检验统计量),这也是方差分析这一名称的来源。
2.1 误差分解
来自于水平内部的数据误差称为组内误差,它是由抽样的随机性造成的随机误差。显然,组内误差只含有随机误差。
来自不同水平之间的数据误差称为组间误差,它可能由抽样造成的随机误差,也可能是由因素的不同水平造成的系统误差。组间误差是随机误差和系统误差的总和。
反映全部数据误差大小的平方和称为总平方和(SST),反映组内误差大小的平方和称为组内平方和(SSE),反映组间误差大小的平方和称为组间平方和(SSA)。换句话说,SST刻画全部数据的波动程度,SSE刻画组内数据的波动程度,SSA刻画不同组均值的差异引起的波动程度。
2.2 误差分析
如果因子的不同水平对每个水平下的均值没有影响,则组间误差只有随机误差而没有系统误差。组内误差和组间误差的均方之比应该接近1;否则它们的比值就会大于1,当大到某个程度时,就认为因子的不同水平之间存在着显著差异,也即自变量(示例中的药物)对因变量(示例中的效果)有显著影响。
由以上的分析可知,在方差分析中,要研究分类型自变量对因变量的影响,在形式上就转化为了检验不同总体(因素的不同水平)的均值是否相等。
2.3 模型假定
- 每个总体都应符合正态分布
- 各个总体的方差 σ 2 \sigma^2 σ2必须相同
- 观测是独立的
三、单因子方差分析
3.1 模型
考虑一个因子 A A A,有 r r r个水平 A 1 , … , A r A_1,\dots,A_r A1,…,Ar, r ≥ 2 r\ge 2 r≥2,设在水平 A i A_i Ai下重复进行了 n i n_i ni次实验,数据是 y i 1 , y i 2 , … , y i n i y_{i1},y_{i2},\dots,y_{in_i} yi1,yi2,…,yini,根据模型的假定,这些数据之间相互独立且 y i j ∼ N ( μ i , σ 2 ) y_{ij}\sim N(\mu_i,\sigma^2) yij∼N(μi,σ2),其中 σ \sigma σ未知。我们关心的问题是 μ i \mu_i μi是否全相等,即要检验 H 0 : μ 1 = ⋯ = μ r v s . H 1 : μ i 不全相等 H_0:\mu_1=\dots=\mu_r\ vs.\ H_1:\mu_i\text{不全相等} H0:μ1=⋯=μr vs. H1:μi不全相等
令
n
=
∑
i
=
1
r
n
i
n=\sum_{i=1}^r n_i
n=∑i=1rni,即总共进行了
n
n
n次实验。令:
y
ˉ
=
1
n
∑
i
=
1
r
∑
j
=
1
n
i
y
i
j
,
y
ˉ
i
⋅
=
1
n
i
∑
j
=
1
n
i
y
i
j
,
i
=
1
,
…
,
r
\bar y = \frac{1}{n}\sum_{i=1}^r\sum_{j=1}^{n_i}y_{ij},\ \bar y_{i\cdot}=\frac 1{n_i}\sum_{j=1}^{n_i}y_{ij},i=1,\dots,r
yˉ=n1i=1∑rj=1∑niyij, yˉi⋅=ni1j=1∑niyij,i=1,…,r
这里,
y
ˉ
\bar y
yˉ表示所有观测的平均值,
y
ˉ
i
⋅
\bar y_{i\cdot}
yˉi⋅表示水平
A
i
A_i
Ai下的观测平均值。那么三个误差分别为:
S
T
2
=
∑
i
=
1
r
∑
j
=
1
n
i
(
y
i
j
−
y
ˉ
)
2
,
S
e
2
=
∑
i
=
1
r
∑
j
=
1
n
i
(
y
i
j
−
y
ˉ
i
⋅
)
2
,
S
A
2
=
∑
i
=
1
r
n
i
(
y
ˉ
i
⋅
−
y
ˉ
)
2
S_T^2=\sum_{i=1}^r\sum_{j=1}^{n_i}(y_{ij}-\bar y )^2,\ S_e^2=\sum_{i=1}^r\sum_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot} )^2,S_A^2 = \sum_{i=1}^r n_i(\bar y_{i\cdot} -\bar y)^2
ST2=i=1∑rj=1∑ni(yij−yˉ)2, Se2=i=1∑rj=1∑ni(yij−yˉi⋅)2,SA2=i=1∑rni(yˉi⋅−yˉ)2
它们满足:
S
T
2
=
S
e
2
+
S
A
2
S_T^2 = S_e^2+S_A^2
ST2=Se2+SA2
这是由于:
S
T
2
=
∑
i
=
1
r
∑
j
=
1
n
i
(
y
i
j
−
y
ˉ
i
⋅
+
y
ˉ
i
⋅
−
y
ˉ
)
2
=
∑
i
=
1
r
∑
j
=
1
n
i
(
y
i
j
−
y
ˉ
i
⋅
)
2
+
∑
i
=
1
r
∑
j
=
1
n
i
(
y
ˉ
i
⋅
−
y
ˉ
)
2
+
2
∑
i
=
1
r
(
y
ˉ
i
⋅
−
y
ˉ
)
∑
j
=
1
n
i
(
y
i
j
−
y
ˉ
i
⋅
)
=
S
e
2
+
S
A
2
S_T^2=\sum\limits_{i=1}^r\sum\limits_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot}+\bar y_{i\cdot}-\bar y )^2\\ =\sum\limits_{i=1}^r\sum\limits_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot})^2+\sum\limits_{i=1}^r\sum\limits_{j=1}^{n_i}(\bar y_{i\cdot}-\bar y )^2+2\sum\limits_{i=1}^r(\bar y_{i\cdot}-\bar y )\sum\limits_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot})\\ =S_e^2+S_A^2
ST2=i=1∑rj=1∑ni(yij−yˉi⋅+yˉi⋅−yˉ)2=i=1∑rj=1∑ni(yij−yˉi⋅)2+i=1∑rj=1∑ni(yˉi⋅−yˉ)2+2i=1∑r(yˉi⋅−yˉ)j=1∑ni(yij−yˉi⋅)=Se2+SA2
我们已经知道当组内误差和组间误差 S A 2 / S e 2 S_A^2/S_e^2 SA2/Se2大到某个程度时,就拒绝原假设,但为了确定具体的拒绝域,我们还需知道在原假设 H 0 H_0 H0成立下, ( S A 2 , S e 2 ) (S_A^2, S_e^2) (SA2,Se2)的分布。
3.2 分析
结论:考虑上述模型假设,有 S e 2 / σ 2 ∼ χ 2 ( n − r ) S_e^2/\sigma^2\sim \chi^2(n-r) Se2/σ2∼χ2(n−r)且 S e 2 S_e^2 Se2与 S A 2 S_A^2 SA2独立,如果 H 0 : μ 1 = ⋯ = μ r H_0:\mu_1=\dots=\mu_r H0:μ1=⋯=μr成立,则有 S A 2 σ 2 ∼ χ 2 ( r − 1 ) , F = S A 2 / ( r − 1 ) S e 2 / ( n − r ) ∼ F ( r − 1 , n − r ) \frac{S_A^2}{\sigma^2}\sim \chi^2(r-1),\ F=\frac{S_A^2/(r-1)}{S_e^2/(n-r)}\sim F(r-1,n-r) σ2SA2∼χ2(r−1), F=Se2/(n−r)SA2/(r−1)∼F(r−1,n−r)
证明:
由单个正态总体的抽样分布定理有
V
i
:
=
1
σ
2
∑
j
=
1
n
i
(
y
i
j
−
y
ˉ
i
⋅
)
2
∼
χ
2
(
n
i
−
1
)
V_i:=\frac{1}{\sigma^2} \sum_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot} )^2\sim \chi^2(n_i-1)
Vi:=σ21j=1∑ni(yij−yˉi⋅)2∼χ2(ni−1)
由于 y i j y_{ij} yij之间独立, 所以 V i V_i Vi相互独立。由卡方分布的可加性,我们有 S e 2 / σ 2 = ∑ i = 1 r V i ∼ χ 2 ( n − r ) S_e^2/\sigma^2=\sum_{i=1}^r V_i\sim \chi^2(n-r) Se2/σ2=i=1∑rVi∼χ2(n−r)
由于
{
V
1
,
…
,
V
r
}
\{V_1,\dots,V_r\}
{V1,…,Vr}与
{
y
ˉ
1
⋅
,
…
,
y
ˉ
r
⋅
}
\{\bar y_{1\cdot},\dots,\bar y_{r\cdot}\}
{yˉ1⋅,…,yˉr⋅}独立,而
S
A
2
S_A^2
SA2是
y
ˉ
i
⋅
\bar y_{i\cdot}
yˉi⋅的函数,所以
S
e
2
S_e^2
Se2与
S
A
2
S_A^2
SA2独立
H
0
H_0
H0成立时,
μ
1
=
⋯
=
μ
r
=
μ
\mu_1=\dots=\mu_r=\mu
μ1=⋯=μr=μ,
y
ˉ
i
⋅
∼
i
i
d
N
(
μ
,
σ
2
/
n
i
)
\bar y_{i\cdot}\stackrel{iid}{\sim} N(\mu,\sigma^2/n_i)
yˉi⋅∼iidN(μ,σ2/ni),令
x
i
=
n
i
(
y
ˉ
i
⋅
−
μ
)
/
σ
∼
i
i
d
N
(
0
,
1
)
x_i = \sqrt{n_i}(\bar y_{i\cdot}-\mu)/\sigma\stackrel{iid}\sim N(0,1)
xi=ni(yˉi⋅−μ)/σ∼iidN(0,1),此时有:
S
A
2
=
∑
i
=
1
r
n
i
(
y
ˉ
i
⋅
−
y
ˉ
)
2
σ
2
S_A^2=\frac{\sum_{i=1}^rn_i(\bar y_{i\cdot} -\bar y)^2}{\sigma^2}
SA2=σ2∑i=1rni(yˉi⋅−yˉ)2
=
∑
i
=
1
r
(
n
i
y
ˉ
i
⋅
−
n
i
∑
j
=
1
r
n
j
n
y
ˉ
j
⋅
)
2
σ
2
=\frac{\sum_{i=1}^r(\sqrt{n_i}\bar y_{i\cdot} -\sqrt{n_i}\sum_{j=1}^r \frac{n_j}{n}\bar y_{j\cdot} )^2}{\sigma^2}
=σ2∑i=1r(niyˉi⋅−ni∑j=1rnnjyˉj⋅)2
=
∑
i
=
1
r
[
n
i
(
y
ˉ
i
⋅
−
μ
)
−
n
i
∑
j
=
1
r
n
j
n
(
y
ˉ
j
⋅
−
μ
)
]
2
σ
2
=\frac{\sum_{i=1}^r[\sqrt{n_i}(\bar y_{i\cdot}-\mu) -\sqrt{n_i}\sum_{j=1}^r \frac{n_j}{n}(\bar y_{j\cdot}-\mu) ]^2}{\sigma^2}
=σ2∑i=1r[ni(yˉi⋅−μ)−ni∑j=1rnnj(yˉj⋅−μ)]2
=
∑
i
=
1
r
(
x
i
−
n
i
∑
j
=
1
r
n
j
n
x
j
)
2
=\sum\limits_{i=1}^r(x_i-\sqrt{n_i}\sum_{j=1}^r \frac{\sqrt{n_j}}{n}x_j)^2
=i=1∑r(xi−ni∑j=1rnnjxj)2
=
∑
i
=
1
r
x
i
2
−
(
∑
i
=
1
r
n
i
/
n
x
i
)
2
=\sum\limits_{i=1}^rx_i^2-(\sum\limits_{i=1}^r\sqrt{n_i/n} x_i)^2
=i=1∑rxi2−(i=1∑rni/nxi)2
=
∣
∣
x
1
:
r
∣
∣
2
−
(
α
⊤
x
1
:
r
)
2
=||x_{1{:}r}||^2-(\alpha^\top x_{1{:}r})^2
=∣∣x1:r∣∣2−(α⊤x1:r)2
其中, α = ( n 1 / n , … , n r / n ) ⊤ \alpha=(\sqrt{n_1/n},\dots,\sqrt{n_r/n})^\top α=(n1/n,…,nr/n)⊤,注意到 ∣ ∣ α ∣ ∣ = 1 ||\alpha|| = 1 ∣∣α∣∣=1,所以可构造一个正交矩阵 A A A使得 A A A的第一行为 α ⊤ \alpha^\top α⊤。
令 z 1 : r = A x 1 : r ∼ N ( 0 , I r ) z_{1{:}r}=Ax_{1{:}r}\sim N(0,I_r) z1:r=Ax1:r∼N(0,Ir),此时 ∣ ∣ z 1 : r ∣ ∣ 2 = ∣ ∣ x 1 : r ∣ ∣ 2 ||z_{1{:}r}||^2=||x_{1{:}r}||^2 ∣∣z1:r∣∣2=∣∣x1:r∣∣2, z 1 = α ⊤ x 1 : r z_1=\alpha^\top x_{1{:}r} z1=α⊤x1:r,所以 S A 2 = ∣ ∣ z 1 : r ∣ ∣ 2 − z 1 2 = ∑ i = 2 r z i 2 ∼ χ 2 ( r − 1 ) S_A^2=||z_{1{:}r}||^2-z_1^2=\sum_{i=2}^r z_i^2\sim \chi^2(r-1) SA2=∣∣z1:r∣∣2−z12=i=2∑rzi2∼χ2(r−1)
由于 S e 2 S_e^2 Se2与 S A 2 S_A^2 SA2独立,所以 F ∼ F ( r − 1 , n − r ) F\sim F(r-1,n-r) F∼F(r−1,n−r)。
2.3 方差分析表
根据2.2中的结论,我们可以利用F检验,拒绝域为 W = { F > F 1 − α ( r − 1 , n − r ) } W=\{F>F_{1-\alpha}(r-1,n-r)\} W={F>F1−α(r−1,n−r)}。可总结得到如下表格:
来源 | 自由度 | 平方和 | 均方和 | F F F值 | p p p值 |
---|---|---|---|---|---|
因子A | r − 1 r-1 r−1 | S A 2 = ∑ i = 1 r n i ( y ˉ i ⋅ − y ˉ ) 2 S_A^2=\sum_{i=1}^r n_i(\bar y_{i\cdot} -\bar y)^2 SA2=∑i=1rni(yˉi⋅−yˉ)2 | S A 2 / ( r − 1 ) S_A^2/(r-1) SA2/(r−1) | S A 2 / ( r − 1 ) S e 2 / ( n − r ) \frac{S_A^2/(r-1)}{S_e^2/(n-r)} Se2/(n−r)SA2/(r−1) | |
误差 | n − r n-r n−r | S e 2 = ∑ i = 1 r ∑ j = 1 n i ( y i j − y ˉ i ⋅ ) 2 S_e^2=\sum_{i=1}^r\sum_{j=1}^{n_i}(y_{ij}-\bar y_{i\cdot} )^2 Se2=∑i=1r∑j=1ni(yij−yˉi⋅)2 | S e 2 / ( n − r ) S_e^2/(n-r) Se2/(n−r) | ||
总和 | n − 1 n-1 n−1 | S T 2 = ∑ i = 1 I ∑ j = 1 n i ( y i j − y ˉ ) 2 S_T^2=\sum_{i=1}^I\sum_{j=1}^{n_i}(y_{ij}-\bar y )^2 ST2=∑i=1I∑j=1ni(yij−yˉ)2 | S T 2 / ( n − 1 ) S_T^2/(n-1) ST2/(n−1) |
当 p p p值小于给定的显著性水平时,拒绝原假设,当接受原假设时,后续的处理可以参考方差分析——联立区间估计
三、两因子方差分析
3.1 模型
考虑因子 A A A有 r r r个水平: A 1 , … , A r A_1,\dots,A_r A1,…,Ar,因子 B B B有 s s s个水平: B 1 , … , B s B_1,\dots,B_s B1,…,Bs。设在水平 ( A i , B j ) (A_i,B_j) (Ai,Bj)进行了 n i j n_{ij} nij次实验,实验数据是 y i j 1 , y i j 2 , … , y i j n i j y_{ij1},y_{ij2},\dots,y_{ijn_{ij}} yij1,yij2,…,yijnij,假设这些数据相互独立且 y i j k ∼ N ( μ i j , σ 2 ) y_{ijk}\sim N(\mu_{ij} ,\sigma^2) yijk∼N(μij,σ2),其中 σ \sigma σ未知。为分析各个因子对指标的影响,令 μ = 1 r s ∑ i = 1 r ∑ j = 1 s μ i j \mu = \frac{1}{rs}\sum_{i=1}^r\sum_{j=1}^s \mu_{ij} μ=rs1i=1∑rj=1∑sμij α i = 1 s ∑ j = 1 s μ i j − μ , i = 1 , … , r , \alpha_i = \frac 1 s\sum_{j=1}^s \mu_{ij}-\mu,\ i=1,\dots,r, αi=s1j=1∑sμij−μ, i=1,…,r, β j = 1 r ∑ i = 1 r μ i j − μ , j = 1 , … , s , \beta_j = \frac 1 r \sum_{i=1}^r \mu_{ij}-\mu,\ j=1,\dots,s, βj=r1i=1∑rμij−μ, j=1,…,s, λ i j = μ i j − μ − α i − β j \lambda_{ij} = \mu_{ij}-\mu-\alpha_i-\beta_j λij=μij−μ−αi−βj
其中, α i \alpha_i αi表示因子 A A A的第 i i i个水平 A i A_i Ai的主效应, β j \beta_j βj表示因子 B B B的第 j j j个水平 B j B_j Bj的主效应, λ i j \lambda_{ij} λij表示 A i A_i Ai和 B j B_j Bj下的交互作用的效应。
于是模型可以被表示为: y i j k = μ + α i + β j + λ i j + ϵ i j k , ϵ i j k ∼ i i d N ( 0 , σ 2 ) y_{ijk} = \mu+ \alpha_i+\beta_j +\lambda_{ij}+\epsilon_{ijk},\ \epsilon_{ijk}\stackrel{iid}{\sim}N(0,\sigma^2) yijk=μ+αi+βj+λij+ϵijk, ϵijk∼iidN(0,σ2)
3.2 检验
我们关心因子 A A A或者因子 B B B或者它们的交互作用 A × B A\times B A×B对指标有没有影响。对应的检验问题为: H 0 : α 1 = ⋯ = α r = 0 v s . H 1 : α i 不全为0 H_0:\alpha_1=\dots=\alpha_r=0\ vs.\ H_1:\alpha_i\text{不全为0} H0:α1=⋯=αr=0 vs. H1:αi不全为0 H 0 : β 1 = ⋯ = β s = 0 v s . H 1 : β j 不全为0 H_0:\beta_1=\dots=\beta_s=0\ vs.\ H_1:\beta_j\text{不全为0} H0:β1=⋯=βs=0 vs. H1:βj不全为0 H 0 : λ 11 = ⋯ = λ r s = 0 v s . H 1 : λ i j 不全为0 H_0:\lambda_{11}=\dots=\lambda_{rs}=0\ vs.\ H_1:\lambda_{ij}\text{不全为0} H0:λ11=⋯=λrs=0 vs. H1:λij不全为0
3.3 ANOVA
可记 n i ⋅ = ∑ j = 1 s n i j , n ⋅ j = ∑ i = 1 r n i j , n = ∑ i = 1 r ∑ j = 1 s n i j n_{i\cdot} = \sum\limits_{j=1}^s n_{ij},\ n_{\cdot j} =\sum\limits_{i=1}^r n_{ij},\ n = \sum\limits_{i=1}^r\sum\limits_{j=1}^s n_{ij} ni⋅=j=1∑snij, n⋅j=i=1∑rnij, n=i=1∑rj=1∑snij,类似的,有下列ANOVA表:
来源 | 自由度 | 平方和 | F F F值 | p p p值 |
---|---|---|---|---|
A A A | r − 1 r-1 r−1 | S A 2 = ∑ i = 1 r n i ⋅ ( y ˉ i ⋅ ⋅ − y ˉ ) 2 S_A^2=\sum\limits_{i=1}^r n_{i\cdot}(\bar y_{i\cdot\cdot} -\bar y)^2 SA2=i=1∑rni⋅(yˉi⋅⋅−yˉ)2 | S A 2 / ( r − 1 ) S e 2 ( n − r s ) \frac{S_A^2/(r-1)}{S_e^2(n-rs)} Se2(n−rs)SA2/(r−1) | |
B B B | s − 1 s-1 s−1 | S B 2 = ∑ j = 1 s n ⋅ j ( y ˉ ⋅ j ⋅ − y ˉ ) 2 S_B^2=\sum\limits_{j=1}^s n_{\cdot j}(\bar y_{\cdot j\cdot} -\bar y)^2 SB2=j=1∑sn⋅j(yˉ⋅j⋅−yˉ)2 | S B 2 / ( s − 1 ) S e 2 / ( n − r s ) \frac{S_B^2/(s-1)}{S_e^2/(n-rs)} Se2/(n−rs)SB2/(s−1) | |
A × B A\times B A×B | ( r − 1 ) ( s − 1 ) (r-1)(s-1) (r−1)(s−1) | S A B 2 = ∑ i = 1 r ∑ j = 1 s n i j ( y ˉ i j ⋅ − y ˉ ) 2 S_{AB}^2=\sum\limits_{i=1}^r\sum\limits_{j=1}^s n_{ij}(\bar y_{ij\cdot} -\bar y)^2 SAB2=i=1∑rj=1∑snij(yˉij⋅−yˉ)2 | S A B 2 / [ ( r − 1 ) ( s − 1 ) ] S e 2 / ( n − r s ) \frac{S_{AB}^2/[(r-1)(s-1)]}{S_e^2/(n-rs)} Se2/(n−rs)SAB2/[(r−1)(s−1)] | |
误 差 误差 误差 | n − r s n-rs n−rs | S e 2 = ∑ i = 1 r ∑ j = 1 s ∑ k = 1 n i j ( y i j k − y ˉ i j ⋅ ) 2 S_e^2=\sum\limits_{i=1}^r\sum_{j=1}^{s}\sum\limits_{k=1}^{n_{ij}}(y_{ijk}-\bar y_{ij\cdot})^2 Se2=i=1∑r∑j=1sk=1∑nij(yijk−yˉij⋅)2 | ||
总 和 总和 总和 | n − 1 n-1 n−1 | S T 2 = ∑ i = 1 r ∑ j = 1 s ∑ k = 1 n i j ( y i j k − y ˉ ) 2 S_T^2=\sum\limits_{i=1}^r\sum\limits_{j=1}^{s}\sum\limits_{k=1}^{n_{ij}}(y_{ijk}-\bar y )^2 ST2=i=1∑rj=1∑sk=1∑nij(yijk−yˉ)2 |