一、从哑变量到方差分析

1.1 哑变量的引入

在前面的学习中，我们已经掌握了线性回归模型的基本形式，即通过线性组合的方式，用一个或多个自变量去解释响应变量的变化。当自变量是连续变量时，线性回归能够直接刻画变量之间的线性关系。然而，在许多实际研究中，自变量可以不是连续变量，它可以是离散变量/属性变量，比如病情的描述可以是“良好”，“中等”，“较差”，“恶化”。如果遇到这样的属性变量，该如何运用到线性回归中呢？

此时，我们仍然可以利用线性回归的框架，只需要引入哑变量(dummy variable)来将其编码为可计算的数值即可。哑变量是通过0-1编码来表征变量所处的状态或属性，从而能够在拟合时使用某一个状态来参与建模。

具体来说，如果一个属性变量有kk个状态，那么我们向模型中添加k−1k-1个哑变量x(1),⋯ ,x(k−1)x_{(1)},\cdots,x_{(k-1)}，这些哑变量的取值是：

x(i)={1,处于状态i0,其他状态x_{(i)}=\begin{cases} 1,&处于状态i\\ 0,&其他状态 \end{cases}

对于第kk个状态，我们只需要将所有哑变量设置为0，即x(1)=⋯=x(k−1)=0x_{(1)}=\cdots=x_{(k-1)}=0。

此时，我们的模型可以写作：

Y=β0+β1x(1)+⋯+βqx(q)+eY=\beta_0+\beta_1x_{(1)}+\cdots+\beta_qx_{(q)}+e

例如，在状态jj下，E(Y∣状态j)=β0+βjE(Y|状态j)=\beta_0+\beta_j，在状态kk下，E(Y∣状态k)=β0E(Y|状态k)=\beta_0。这样一来就分别表示了属性变量不同状态下的观测。

除了上面这种取法，还有一种取法，就是在上述取法的基础上，对状态kk取-1，即：

x(i)={1,处于状态i−1,处于状态k0,其他状态x_{(i)}=\begin{cases} 1,&处于状态i\\ -1,&处于状态k\\ 0,&其他状态 \end{cases}

此时，在状态jj下，E(Y∣状态j)=β0+βjE(Y|状态j)=\beta_0+\beta_j不变，但是在状态kk下：

E(Y∣状态k)=β0−(β1+⋯+βk−1)E(Y|状态k)=\beta_0-(\beta_1+\cdots+\beta_{k-1})

此时：

β0=∑j=1kE(Y∣状态j)/k\beta_0=\sum_{j=1}^kE(Y|状态j)/k

β0\beta_0就代表了这个属性变量的平均效应。

1.2 方差分析的概念

我们再回到1.1小节所提到的模型：

Y=β0+β1x(1)+⋯+βqx(q)+eY=\beta_0+\beta_1x_{(1)}+\cdots+\beta_qx_{(q)}+e

这个模型虽然是线性回归，但是这里的回归系数还是原来的含义吗？

在一般的线性回归模型当中，回归系数刻画了当自变量变化一个单位时响应变量的变化量。而在这个模型中你会发现，每一个回归系数只刻画了某个状态的贡献。也就是说，因为每一个样本只可能处于某一个状态，所以当样本值输入到模型当中时，除了常数项β0\beta_0之外，只有一个回归系数能够参与运算。再进一步说，我们是以某个状态的均值为基准，考察其他状态均值相对于基准状态均值的差异，而对回归系数的假设检验则是考察状态之间均值的差异是否显著。

既然是考虑各个状态的均值差异是否显著，那么除了通过回归模型检验回归系数，还可以对响应变量的方差（波动）进行分解。也就是说，响应变量的方差可以被分解为：

总方差=状态/组别之间方差+状态/组别内部方差总方差=状态/组别之间方差+状态/组别内部方差

对于同一个状态/分组，其内部存在波动；不同状态/组别的均值之间也存在波动；如果组别之间的波动远远大于组别内部的波动，就能说明每个组别是显著不同的。这就是方差分析的概念。

方差分析(Analysis of Variance, ANOVA)模型是一类特殊的线性回归模型，它的响应变量往往是连续的，而自变量则是属性变量的状态。方差分析并不是与回归分析平行的另一种统计方法，而是线性回归模型的一个特例，它关注的核心问题是：不同组别的均值差异是否显著。

方差分析在许多研究中都有应用，例如不同肥料是否导致作物产量的差异、不同教学方法是否影响学生成绩等，我们关心的问题在于不同处理组或条件下的均值是否不同。

你或许听说过协方差分析(Analysis of Covariance, ANCOVA)。协方差分析就是在属性变量的基础上再添加一些连续变量进行混合分析，即：

Y=β0+β1x(1)+⋯+βqx(q)+Zγ+eY=\beta_0+\beta_1x_{(1)}+\cdots+\beta_qx_{(q)}+Z\gamma+e

其中γ=(γ1,⋯ ,γp)T\gamma=(\gamma_1,\cdots,\gamma_p)^T代表连续变量的回归系数，ZZ为设计矩阵。协方差分析本质上就是在考虑连续变量效应（这些连续变量称为协变量）的基础上进行方差分解，也是非常常见的模型。但是我们并不进行展开，这是因为协方差分析理论较为落后，其主要原因是受到算力的限制，被迫发展出便于手工计算的方法。在目前的高算力时代，直接引入哑变量进行拟合已不再是难题，所以我们不再学习传统的协方差分析。

当然，后面你会看到，方差分析也是可以手工计算而不需要拟合模型的，但是它的思想非常重要，应用也不局限于方差分析本身。

二、单因素方差分析

我们先从最简单的情况开始，也就是只有一个属性变量的情况。

2.1 模型概述与拟合

假设我们仅仅研究因素A对响应值的影响，因素A有a个水平，并且在第Ai(i=1,⋯ ,a)A_i(i=1,\cdots,a)个水平下有nin_i个响应值yi1,⋯ ,yiniy_{i1},\cdots,y_{in_i}。如果使用方差分析来建立模型，那么观测值要服从正态分布，即：

yij∼N(μi,σ2),i=1,⋯ ,a,j=1,⋯ ,niy_{ij}\sim N(\mu_i,\sigma^2),i=1,\cdots,a,j=1,\cdots,n_i

其中μi\mu_i表示第ii个水平的均值。因此，我们可以给出一种可能的模型结构：

yij=μi+eij,eij∼N(0,σ2)y_{ij}=\mu_i+e_{ij},e_{ij}\sim N(0,\sigma^2)

这种结构我们称之为均值模型。在均值模型下，我们很容易写出估计值：

μ^i=∑j=1niyij/ni=y‾i⋅\hat\mu_i=\sum_{j=1}^{n_i}y_{ij}/n_i=\overline y_{i\cdot}\\

也就是说，每个水平均值的估计值就是这个水平下所有观测值的平均值。

均值模型对于单因素方差分析来说已经足够，但不便于把所有方差分析模型都统一起来。于是，我们使用一个更普遍地写法，称之为主效应模型，即：

yij=μi+eij=μ+αi+eij,μ=∑i=1aniμi/Ny_{ij}=\mu_i+e_{ij}=\mu+\alpha_i+e_{ij},\mu=\sum_{i=1}^an_i\mu_i/N

其中μ\mu称为主效应，是因素A各个水平带来的平均效应；αi\alpha_i称为额外效应，是由因素A的不同水平带来的效应。也就是说，我们认为只要存在因素A，就应该有一个平均效应μ\mu，然后根据样本所处不同组别进行额外效应的调整。

如果我们把上述分量形式写成矩阵形式（其中EnE_n代表(1,⋯ ,1)1×nT(1,\cdots,1)^T_{1\times n}）：

(y11⋮y1n1y21⋮ya1⋮yana)=(En1En10⋯0En20En2⋯0⋮⋮⋮⋱⋮Ena00⋯Ena)(μα1α2⋮αa)+(e11⋮e1n1e21⋮ea1⋮eana)\begin{pmatrix}y_{11}\\ \vdots\\ y_{1n_1}\\ y_{21}\\ \vdots\\ y_{a1}\\ \vdots\\ y_{an_a}\end{pmatrix}= \begin{pmatrix}E_{n_1}&E_{n_1}&0&\cdots&0\\ E_{n_2}&0&E_{n_2}&\cdots&0\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ E_{n_a}&0&0&\cdots&E_{n_a}\end{pmatrix} \begin{pmatrix}\mu\\ \alpha_1\\ \alpha_2\\ \vdots\\ \alpha_a\end{pmatrix}+ \begin{pmatrix}e_{11}\\ \vdots\\ e_{1n_1}\\ e_{21}\\ \vdots\\ e_{a1}\\ \vdots\\ e_{an_a}\end{pmatrix}

你就会发现这就是一个线性回归模型，只不过没使用哑变量。另外，由于每一次观测都只属于某一个因素水平，所以这个模型的设计矩阵XX是精确设计的，这与一般的线性回归模型有很大不同。

有了这个表达式并不够，你会发现XN×(a+1)X_{N\times(a+1)}的秩为a，即设计矩阵不是列满秩的，所以根据之前所学，这个模型存在多个解。那么为了使得模型有解，我们需要给额外效应加上约束条件：

∑i=1aniαi=0\sum_{i=1}^an_i\alpha_i=0

这个约束条件其实非常符合逻辑。在统计学中我们知道，一组样本x1,⋯ ,xnx_1,\cdots,x_n的离均差之和为0，即：

∑i=1n(xi−x‾)=0\sum_{i=1}^n(x_i-\overline x)=0

这里也一样，因为我们定义了因素A的平均效应/主效应，那么因素A的各个水平的额外效应（相当于离均差）的总和应当为0，所以这个条件非常自然，同时也使得模型有解。

因此，我们可以写出最终的主效应模型：

{Y=Xβ+e,β=(μ,α1,⋯ ,αa)TLβ=0,L=(0,n1,⋯ ,nα)\begin{cases} Y=X\beta+e,&\beta=(\mu,\alpha_1,\cdots,\alpha_a)^T\\ L\beta=0,&L=(0,n_1,\cdots,n_\alpha) \end{cases}

如果你学习了之前的内容，你会感到非常熟悉，这就是有约束的线性回归模型，并且设计矩阵被精确设计，不存在共线性等一系列问题。根据正规方程不难解得：

μ^=y‾⋅⋅,α^i=y‾i⋅−y‾⋅⋅\hat\mu=\overline y_{\cdot\cdot},\hat\alpha_{i}=\overline y_{i\cdot}-\overline y_{\cdot\cdot}

同时，由于Lβ=0L\beta=0是边界条件，所以上述估计是无偏估计。

这里的y‾⋅⋅\overline y_{\cdot\cdot}代表总体均值（也就是对行列所有数据求平均），y‾i⋅\overline y_{i\cdot}代表对第ii行/水平的所有列求平均。

2.2 模型假设检验

我们已经知道，方差分析模型要检验的就是不同水平之间是否存在显著差异。那么对于单因素方差分析来说，其检验的目标是：

H0:α1=⋯=αaH_0:\alpha_1=\cdots=\alpha_a

仿照线性回归模型的检验方法，我们构造组内平方和，即每个观测值减去其所属水平的均值的平方：

ESS=∑i=1a∑j=1ni(yij−μ^−α^i)2=∑i=1a∑j=1ni(yij−y‾i⋅)2ESS=\sum_{i=1}^a\sum_{j=1}^{n_i}(y_{ij}-\hat\mu-\hat\alpha_i)^2=\sum_{i=1}^a\sum_{j=1}^{n_i}(y_{ij}-\overline y_{i\cdot})^2

同理，我们构造总平方和，即每个观测值减去总体均值的平方：

TSS=∑i=1a∑j=1ni(yij−y‾⋅⋅)2TSS=\sum_{i=1}^a\sum_{j=1}^{n_i}(y_{ij}-\overline y_{\cdot\cdot})^2

用总平方和减去组内平方和，就可以构造组间平方和，或者等价于用每个水平的均值减去总体均值的平方，再赋予每组的样本数：

TSS−ESS=∑i=1ani(y‾i⋅−y‾⋅⋅)2≡SSATSS-ESS=\sum_{i=1}^an_i(\overline y_{i\cdot}-\overline y_{\cdot\cdot})^2\equiv SSA

根据我们之前的讨论，如果因素A的组间平方和SSA显著有别于组内平方和ESS，那么说明因素A的不同水平之间有显著差异。类似于线性回归模型，我们仍然可以构造FF分布：

F=SSA/(a−1)ESS/(N−a)∼Fa−1,N−aF=\frac{SSA/(a-1)}{ESS/(N-a)}\sim F_{a-1,N-a}

如果通过检验，说明因素A的不同水平之间是有差异的。

为了记录计算过程，便于手工计算，Fisher提出了方差分析表。我们重新定义分子分母为两个均方和

MSA=SSA/(a−1)MSE=ESS/(N−a)\begin{aligned} &MSA=SSA/(a-1)\\ &MSE=ESS/(N-a) \end{aligned}

那么可以构造如下方差分析表：

来源 自由度 平方和 均方和 F值 显著性标识
因子A a−1a-1 SSA MSA F=MSA/MSE ***
误差E N−aN-a ESS MSE
总和T N−1N-1 TSS

基于这个表就可以很清晰的看到方差分析的过程，这也是很多方差分析软件返回的结果形式。值得一提的是，MSEMSE就是线性回归模型中的σ^2\hat\sigma^2，是对σ2\sigma^2的无偏估计。

如果你学习过之前的内容，那么你或许对“方差分析表”一词有点熟悉。没错，在第4讲的内容当中，我们对线性回归模型进行假设检验时用到了方差分析表，当时的输出是：

> anova(nullmod, lmod)
Analysis of Variance Table

Model 1: Species ~ 1
Model 2: Species ~ Area + Elevation + Nearest + Scruz + Adjacent
  Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
1     29 381081                                  
2     24  89231  5    291850 15.699 6.838e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

你可以看到这个命令返回了一个方差分析表，这个表和上面的表大同小异。那么为什么线性回归模型检验本质是方差分析呢？

1.2小节中我们提到，方差分析本质也就是检验每个水平前面的回归系数是否显著不为0，那么在线性回归中原理也是一样的：回归平方和相当于组间平方和，残差平方和相当于组内平方和。这也是为什么我们可以把anova()函数用于线性回归模型。

2.3 多重比较

如果原假设被拒绝，那么说明总体上来说因素A的不同水平之间存在差异，但是我们并不清楚每两个水平之间是否真的存在差异。就像线性回归模型，我们只知道一定有回归系数显著不为0，但是不能保证每一个回归系数都显著不为0，因此要对每一个回归系数进行检验。

在方差分析中，这种检验被称为多重比较(Multiple Comparison)，即我们要检验的目标是αi=αj(i≠j)\alpha_i=\alpha_j(i\ne j)，或者等价于要检验αi−αj(i≠j)\alpha_i-\alpha_j(i\ne j)的置信区间是否包含0点。因此，我们要设法求出置信区间。下面给出三种可行的方法：

  1. Bonferroni同时置信区间
  2. Scheffe同时置信区间
  3. TukeyHSD

前面两种方法我们在之前的讲解中已经遇到，不过我们需要在方差分析的版块中进行一定的变化；第三种是方差分析特有的方法，当然也是更流行、更推荐的方法。

为了模型的通用性，这里我们引入对照的概念，即∑i=1aciαi\sum_{i=1}^ac_i\alpha_i，其中∑i=1aci=0\sum_{i=1}^ac_i=0。对于方差分析来说，多重比较的对照是遍历αi−αj\alpha_i-\alpha_j，对照总个数是m=a(a−1)/2m=a(a-1)/2。例如，如果要检验α1−α2\alpha_1-\alpha_2，那么其对应的对照是c1=(1−10⋯0)c_1=\begin{pmatrix}1&-1&0&\cdots&0\end{pmatrix}。

Bonferroni同时置信区间。对于mm个对照∑i=1aci(k)αi(k=1,⋯ ,m)\sum_{i=1}^ac_i^{(k)}\alpha_i(k=1,\cdots,m)来说，其1−α1-\alpha的Bonferroni同时置信区间为：

∑i=1aci(k)y‾i⋅±tN−a(α2m)σ^2∑i=1a[ci(k)]2n,k=1,⋯ ,m\sum_{i=1}^ac_i^{(k)}\overline y_{i\cdot}\pm t_{N-a}(\frac{\alpha}{2m})\hat\sigma^2\sqrt{\sum_{i=1}^a\frac{[c_i^{(k)}]^2}{n}},k=1,\cdots,m

对于方差分析来说，Bonferroni同时置信区间为：

(y‾i⋅−y‾j⋅)±tN−a(α2m)σ^21ni+1nj,1≤i≠j≤a(\overline y_{i\cdot}-\overline y_{j\cdot})\pm t_{N-a}(\frac{\alpha}{2m})\hat\sigma^2\sqrt{\frac{1}{n_i}+\frac{1}{n_j}},1\le i\ne j\le a

Scheffe同时置信区间。无论对照个数，所有对照的1−α1-\alpha的Scheffe同时置信区间为：

∑i=1aciy‾i⋅±σ^(a−1)Fa−1,N−a(α)∑i=1aci2ni,1≤i≠j≤a\sum_{i=1}^ac_i\overline y_{i\cdot}\pm\hat\sigma\sqrt{(a-1)F_{a-1,N-a}(\alpha)\sum_{i=1}^a\frac{c_i^2}{n_i}},1\le i\ne j\le a

对于方差分析来说，Scheffe同时置信区间为：

(y‾i⋅−y‾j⋅)±σ^2(a−1)Fa−1,N−a(α)(1ni+1nj),1≤i≠j≤a(\overline y_{i\cdot}-\overline y_{j\cdot})\pm \hat\sigma^2\sqrt{(a-1)F_{a-1,N-a}(\alpha)(\frac{1}{n_i}+\frac{1}{n_j})},1\le i\ne j\le a

TukeyHSD(Tukey Honest Significant Difference)。TukeyHSD是由Tukey提出来的检验方法，它没有采用熟知的tt分布或FF分布，而是提出了一个新的分布，称为学生化极差分布。设一组独立样本样本Z1,⋯ ,Zn∼N(0,1)Z_1,\cdots,Z_n\sim N(0,1)，且mW2∼χm2mW^2\sim\chi_m^2，那么定义：

Q={max⁡1≤i≤nZi−min⁡1≤i≤nZi}W∼Qn,mQ=\frac{\{\max_{1\le i\le n}Z_i-\min_{1\le i\le n}Z_i\}}{W}\sim Q_{n,m}

因此，如果设yi∼N(μi,σ2),i=1,⋯ ,ny_i\sim N(\mu_i,\sigma^2),i=1,\cdots,n，且mσ^2/σ2∼χm2m\hat\sigma^2/\sigma^2\sim\chi_m^2，那么Tukey的1−α1-\alpha同时置信区间应当为：

(yi−yj)±σ^Qn,m(α)(y_i-y_j)\pm \hat\sigma Q_{n,m}(\alpha)

在单因素方差分析中，TukeyHSD只能用于平衡数据，即n1=⋯=na=nn_1=\cdots=n_a=n，此时y‾i⋅∼N(μ+αi,σ2/n),i=1,⋯ ,a\overline y_{i\cdot}\sim N(\mu+\alpha_i,\sigma^2/n),i=1,\cdots,a，(N−a)σ^2/σ2∼χN−a2(N-a)\hat\sigma^2/\sigma^2\sim\chi_{N-a}^2，那么所有αi−αj\alpha_i-\alpha_j的1−α1-\alpha同时置信区间为：

(y‾i⋅−y‾j⋅)±σ^nQa,N−a(α)(\overline y_{i\cdot}-\overline y_{j\cdot})\pm \frac{\hat\sigma}{\sqrt n}Q_{a,N-a}(\alpha)

当然，我们可以推广到所有对照：

∑i=1aciy‾i⋅±Qa,N−a(α)σ^2n∑i=1a∣ci∣\sum_{i=1}^ac_i\overline y_{i\cdot}\pm Q_{a,N-a}(\alpha)\frac{\hat\sigma}{2\sqrt n}\sum_{i=1}^a|c_i|

三、双因素方差分析

上一节我们处理了单因素的方差分析，现在我们将其推广到两个因素。不过其中的原理和方法都是类似的，我们可以仿造上面的内容衍生对应的结论。

3.1 无交互双因素方差分析

无交互，指的是两个因素之间互不干扰，对观测值的影响相互独立。在无交互模型中，我们考虑因素A有a个水平，因素B有b个水平。根据主效应模型，我们定义总体的主效应μ\mu，因素A的额外效应αi\alpha_i，因素B的额外效应βj\beta_j。我们认为效应之间是可加的，且由于不存在交互，每一组效应(αi,βj)(\alpha_i,\beta_j)之下只进行一次测定。

模型的基本假定为：yij∼N(μij,σ2)y_{ij}\sim N(\mu_{ij},\sigma^2)，eij∼N(0,σ2)e_{ij}\sim N(0,\sigma^2)，μ=∑i=1a∑j=1bμij/(ab)\mu=\sum_{i=1}^a\sum_{j=1}^b\mu_{ij}/(ab)。

于是模型可以写作：

yij=μ+αi+βj+eij,i=1,⋯ ,a,j=1,⋯ ,by_{ij}=\mu+\alpha_i+\beta_j+e_{ij},i=1,\cdots,a,j=1,\cdots,b

容易想到设计矩阵XN×(a+b+1)X_{N\times (a+b+1)}的秩是a+b−1a+b-1，解不唯一，于是添加约束：

∑i=1aαi=0,∑j=1bβj=0\sum_{i=1}^a\alpha_i=0,\sum_{j=1}^b\beta_j=0

于是模型转化为矩阵形式：

{Y=Xβ+e,β=(μ,α1⋯ ,αa,β1,⋯ ,βb)T(L1L2)β=0,L1=(0EaT),L2=(0EbT)\begin{cases} Y=X\beta+e,&\beta=(\mu,\alpha_1\cdots,\alpha_a,\beta_1,\cdots,\beta_b)^T\\ \begin{pmatrix}L_1&L_2\end{pmatrix}\beta=0,&L_1=\begin{pmatrix}0&E^T_a\end{pmatrix},L_2=\begin{pmatrix}0&E^T_b\end{pmatrix} \end{cases}

求解可得：

μ^=y‾⋅⋅,α^i=y‾i⋅−y‾⋅⋅,β^j=y‾⋅j−y‾⋅⋅\hat\mu=\overline y_{\cdot\cdot},\hat\alpha_i=\overline y_{i\cdot}-\overline y_{\cdot\cdot},\hat\beta_j=\overline y_{\cdot j}-\overline y_{\cdot\cdot}

得到估计值后就要进行假设检验。无交互双因素方差分析的原假设有两个：α1=⋯=αa\alpha_1=\cdots=\alpha_a，β1=⋯=βj\beta_1=\cdots=\beta_j，因此总平方和的来源是因素A、因素B、残差，于是我们可以定义方差分析表的几个指标：

TSS=SSA+SSB+ESSTSS=∑i=1a∑j=1b(yij−y‾⋅⋅)2ESS=∑i=1a∑j=1b(yij−y‾i⋅−y‾⋅j+y‾⋅⋅)2SSA=∑i=1a∑j=1b(y‾i⋅−y‾⋅⋅)2SSB=∑i=1a∑j=1b(y‾⋅j−y‾⋅⋅)2\begin{aligned} TSS&=SSA+SSB+ESS\\ TSS&=\sum_{i=1}^a\sum_{j=1}^b(y_{ij}-\overline y_{\cdot\cdot})^2\\ ESS&=\sum_{i=1}^a\sum_{j=1}^b(y_{ij}-\overline y_{i\cdot}-\overline y_{\cdot j}+\overline y_{\cdot\cdot})^2\\ SSA&=\sum_{i=1}^a\sum_{j=1}^b(\overline y_{i\cdot}-\overline y_{\cdot\cdot})^2\\ SSB&=\sum_{i=1}^a\sum_{j=1}^b(\overline y_{\cdot j}-\overline y_{\cdot\cdot})^2 \end{aligned}

无交互双因素方差分析的方差分析表为：

来源 自由度 平方和 均方和 F值 显著性
因子A a−1a-1 SSA MSA MSA/MSE ***
因子B b−1b-1 SSB MSB MSB/MSE ***
误差E (a−1)(b−1)(a-1)(b-1) ESS MSE
总和T ab−1ab-1 TSS

类似地，我们仿造单因素方差分析的结论，给出同时置信区间。

任意mm个αi−αi′\alpha_i-\alpha_{i^{'}}和βj−βj′\beta_j-\beta_{j^{'}}的Bonferroni1−α1-\alpha同时置信区间分别为：

(y‾i⋅−y‾i′⋅)±t(a−1)(b−1)(α2m)σ^2b(y‾⋅j−y‾⋅j′)±t(a−1)(b−1)(α2m)σ^2a\begin{aligned} (\overline y_{i\cdot}-\overline y_{i^{'}\cdot})&\pm t_{(a-1)(b-1)}(\frac{\alpha}{2m})\hat\sigma\sqrt{\frac{2}{b}}\\ (\overline y_{\cdot j}-\overline y_{\cdot j^{'}})&\pm t_{(a-1)(b-1)}(\frac{\alpha}{2m})\hat\sigma\sqrt{\frac{2}{a}} \end{aligned}

所有αi−αi′\alpha_i-\alpha_{i^{'}}和βj−βj′\beta_j-\beta_{j^{'}}的Scheffe1−α1-\alpha同时置信区间分别为：

(y‾i⋅−y‾i′⋅)±σ^(a−1)2bF(a−1)(b−1)(α)(y‾⋅j−y‾⋅j′)±σ^(b−1)2aF(a−1)(b−1)(α)\begin{aligned} (\overline y_{i\cdot}-\overline y_{i^{'}\cdot})&\pm \hat\sigma\sqrt{(a-1)\frac{2}{b}F_{(a-1)(b-1)}(\alpha)}\\ (\overline y_{\cdot j}-\overline y_{\cdot j^{'}})&\pm \hat\sigma\sqrt{(b-1)\frac{2}{a}F_{(a-1)(b-1)}(\alpha)} \end{aligned}

所有αi−αi′\alpha_i-\alpha_{i^{'}}和βj−βj′\beta_j-\beta_{j^{'}}的Tukey1−α1-\alpha同时置信区间分别为：

(y‾i⋅−y‾i′⋅)±Qa,(a−1)(b−1)(α)σ^b(y‾⋅j−y‾⋅j′)±Qb,(a−1)(b−1)(α)σ^a\begin{aligned} (\overline y_{i\cdot}-\overline y_{i^{'}\cdot})&\pm Q_{a,(a-1)(b-1)}(\alpha)\frac{\hat\sigma}{\sqrt b}\\ (\overline y_{\cdot j}-\overline y_{\cdot j^{'}})&\pm Q_{b,(a-1)(b-1)}(\alpha)\frac{\hat\sigma}{\sqrt a} \end{aligned}

3.2 有交互双因素方差分析

3.1小节介绍了无交互下的双因素方差分析，这一小节我们来看存在交互项的情况。

存在交互项，意味着在(αi,βj)(\alpha_i,\beta_j)水平下还有二者交互产生的作用rijr_{ij}。因此，如果要拟合交互项，那么每一个水平组合下就应当不止测量1次。于是模型可以写作：

yijk=μ+αi+βj+rij+eijk,k=1,⋯ ,nijy_{ijk}=\mu+\alpha_i+\beta_j+r_{ij}+e_{ijk},k=1,\cdots,n_{ij}

模型中的kk表示了在(αi,βj)(\alpha_i,\beta_j)水平下测量的次数nijn_{ij}。

为什么存在交互项就必须要重复测量？

这是因为，如果每一个水平组合下只有一个观测，那么所谓的“组内方差”就为0，这样就无法分离误差项，也就无法解释主效应和额外效应了。就好像只用一个点做回归分析，毫无意义。

为了简化，我们假定是平衡数据，即所有nijn_{ij}均为常数cc，那么设计矩阵Xabc×(ab+a+b+1)X_{abc\times (ab+a+b+1)}的秩应当为abab（根据正规方程，最后只有abab个独立的方程），因此我们需要a+b+1a+b+1个约束条件，即：

∑i=1aαi=0∑j=1bβj=0∑j=1brij=0,1≤i≤a∑i=1arij=0,1≤j≤b\begin{aligned} &\sum_{i=1}^a\alpha_i=0\\ &\sum_{j=1}^b\beta_j=0\\ &\sum_{j=1}^br_{ij}=0,1\le i\le a\\ &\sum_{i=1}^ar_{ij}=0,1\le j\le b\\ \end{aligned}

根据正规方程可以解得：

μ^=y‾⋅⋅⋅α^i=y‾i⋅⋅−y‾⋅⋅⋅β^j=y‾⋅j⋅−y‾⋅⋅⋅r^ij=y‾ij⋅−y‾i⋅⋅−y‾⋅j⋅+y‾⋅⋅⋅\begin{aligned} \hat\mu&=\overline y_{\cdot\cdot\cdot}\\ \hat\alpha_i&=\overline y_{i\cdot\cdot}-\overline y_{\cdot\cdot\cdot}\\ \hat\beta_j&=\overline y_{\cdot j\cdot}-\overline y_{\cdot\cdot\cdot}\\ \hat r_{ij}&=\overline y_{ij\cdot}-\overline y_{i\cdot\cdot}-\overline y_{\cdot j\cdot}+\overline y_{\cdot\cdot\cdot} \end{aligned}

同理，假设检验的统计量为：

TSS=SSA+SSB+SSAB+SSESSA=∑i=1a∑j=1b∑k=1c(y‾i⋅⋅−y‾⋅⋅⋅)2SSB=∑i=1a∑j=1b∑k=1c(y‾⋅j⋅−y‾⋅⋅⋅)2SSAB=∑i=1a∑j=1b∑k=1c(y‾ij⋅−y‾i⋅⋅−y‾⋅j⋅+y‾⋅⋅⋅)2ESS=∑i=1a∑j=1b∑k=1c(yijk−y‾ij⋅)2\begin{aligned} TSS&=SSA+SSB+SSAB+SSE\\ SSA&=\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c(\overline y_{i\cdot\cdot}-\overline y_{\cdot\cdot\cdot})^2\\ SSB&=\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c(\overline y_{\cdot j\cdot}-\overline y_{\cdot\cdot\cdot})^2\\ SSAB&=\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c(\overline y_{ij\cdot}-\overline y_{i\cdot\cdot}-\overline y_{\cdot j\cdot}+\overline y_{\cdot\cdot\cdot})^2\\ ESS&=\sum_{i=1}^a\sum_{j=1}^b\sum_{k=1}^c(y_{ijk}-\overline y_{ij\cdot})^2 \end{aligned}

方差分析表为：

来源 自由度 平方和 均方和 F值 显著性
因子A a−1a-1 SSA MSA MSA/MSE ***
因子B b−1b-1 SSB MSB MSB/MSE ***
交互项 (a−1)(b−1)(a-1)(b-1) SSAB MSAB MSAB/MSE ***
误差E ab(c−1)ab(c-1) ESS MSE
总和T abc−1abc-1 TSS

四、方差分析的诊断

前面我们介绍了方差分析的模型构建、假设检验、多重比较等，但是一直没有诊断模型的假设。如果不满足假设，那么上述所有的检验统计量就都不服从FF分布，那么方差分析的结果就不可靠。所以，有必要对方差分析的前提假设进行诊断。对于方差分析来说，我们主要检验残差的正态性和方差齐性。这里我们以单因素方差分析为例进行说明。

4.1 正态性诊断

回顾单因素方差分析的模型表达式为yij=μ+αi+eijy_{ij}=\mu+\alpha_i+e_{ij}，其中i=1,⋯ ,ai=1,\cdots,a，j=1,⋯ ,nij=1,\cdots,n_i，令残差的估计值为e^ij=yij−y‾i⋅\hat e_{ij}=y_{ij}-\overline y_{i\cdot}，那么有：

E(e^ij)=0Var(e^ij)=ni−1niσ2Ee^ije^i′j′={0,i≠i′−σ2ni,i=i′,j≠j′\begin{aligned} &E(\hat e_{ij})=0\\ &Var(\hat e_{ij})=\frac{n_i-1}{n_i}\sigma^2\\ &E\hat e_{ij}\hat e_{i^{'}j^{'}}= \begin{cases} 0,&i\ne i^{'}\\ -\frac{\sigma^2}{n_i},&i=i^{'},j\ne j^{'} \end{cases} \end{aligned}

从上面结论中可以看到，同一水平下残差方差相等但不独立，不同水平下残差方差不等但相互独立。于是做如下线性变换：

Zil=ll+1(1l∑j=1le^ij−e^i,j+1)=ll+1(1l∑j=1lyij−yi,j+1)l=1,⋯ ,ni−1;i=1,⋯ ,a\begin{aligned} Z_{il}&=\sqrt{\frac{l}{l+1}}(\frac{1}{l}\sum_{j=1}^l\hat e_{ij}-\hat e_{i,j+1})\\ &=\sqrt{\frac{l}{l+1}}(\frac{1}{l}\sum_{j=1}^l y_{ij}-y_{i,j+1})\\ &l=1,\cdots,n_i-1;i=1,\cdots,a \end{aligned}

此时我们把N=∑i=1aniN=\sum_{i=1}^an_i个残差变为了N−aN-a个ZilZ_{il}，且满足EZil=0EZ_{il}=0，Var(Zil)=σ2Var(Z_{il})=\sigma^2，Cov(Zil,Zi′l′)=0Cov(Z_{il},Z_{i^{'}l^{'}})=0，其中i≠i′i\ne i^{'}，l≠l′l\ne l^{'}。只需要把ZilZ_{il}看作从N(0,σ2)N(0,\sigma^2)总体中抽出的一组独立样本，使用通常检验残差正态分布的方法做检验即可。

4.2 方差齐性诊断

FF检验对方差齐性很敏感，因此必须要做方差齐性检验。对于单因素方差分析，yij=μ+αi+eijy_{ij}=\mu+\alpha_i+e_{ij}，其中i=1,⋯ ,ai=1,\cdots,a，j=1,⋯ ,nij=1,\cdots,n_i，eij∼N(0,σi2)e_{ij}\sim N(0,\sigma^2_i)且相互独立，那么方差齐性检验的原假设是σ12=⋯=σa2\sigma_1^2=\cdots=\sigma^2_a。下面介绍四种方法：

  • Levene检验
  • Hartley检验
  • Cochran C检验
  • Bartlett检验

Levene检验。该检验只能用于平衡数据，即n1=⋯=na=nn_1=\cdots=n_a=n。令lij=e^ij2l_{ij}=\hat e_{ij}^2，e^ij=yij−y‾i⋅\hat e_{ij}=y_{ij}-\overline y_{i\cdot}，那么在原假设成立的情况下有：

L=a(n−1)a−1∑i=1a∑j=1n(l‾i⋅−l‾⋅⋅)2∑i=1a∑j=1n(lij−l‾i⋅)2∼Fa−1,a(n−1)L=\frac{a(n-1)}{a-1}\frac{\sum_{i=1}^a\sum_{j=1}^n(\overline l_{i\cdot}-\overline l_{\cdot\cdot})^2}{\sum_{i=1}^a\sum_{j=1}^n(l_{ij}-\overline l_{i\cdot})^2}\sim F_{a-1,a(n-1)}

Hartley检验(最大F比检验)。令SSEi=∑j=1ni(yij−y‾i⋅)2SSE_{i}=\sum_{j=1}^{n_i}(y_{ij}-\overline y_{i\cdot})^2，由于在正态分布下有SSEi∼σ2χni−12SSE_{i}\sim\sigma^2\chi^2_{n_i-1}，因此令MSEi=SSEi/(ni−1)MSE_i=SSE_i/(n_i-1)，那么在原假设成立的情况下有：

Fmax=max⁡1≤i≤aMSEimin⁡1≤i≤aMSEi∼Fa,n−1F_{max}=\frac{\max_{1\le i\le a MSE_i}}{\min_{1\le i\le a MSE_i}}\sim F_{a,n-1}

Cochran C检验。该检验只能用于平衡数据，即n1=⋯=na=nn_1=\cdots=n_a=n。与Hartley检验类似，在原假设成立下：

G=max⁡1≤i≤aMSEi∑i=1aMSEi∼Ga,n−1G=\frac{\max_{1\le i\le a}MSE_i}{\sum_{i=1}^aMSE_i}\sim G_{a,n-1}

查询临界值表即可。

Bartlett检验。需要进行如下的定义：

MSE=1N−a∑i=1a∑j=1ni(yij−y‾i⋅)2MSEi=∑j=1ni(yij−y‾i⋅)2/(ni−1)C=1+13(a−1)(∑i=1a1ni−1−1N−a)B=1C[(N−a)ln⁡MSE−∑i=1a(ni−1)ln⁡MSEi]\begin{aligned} MSE&=\frac{1}{N-a}\sum_{i=1}^a\sum_{j=1}^{n_i}(y_{ij}-\overline y_{i\cdot})^2\\ MSE_i&=\sum_{j=1}^{n_i}(y_{ij}-\overline y_{i\cdot})^2/(n_i-1)\\ C&=1+\frac{1}{3(a-1)}(\sum_{i=1}^a\frac{1}{n_i-1}-\frac{1}{N-a})\\ B&=\frac{1}{C}\left[(N-a)\ln MSE-\sum_{i=1}^a(n_i-1)\ln MSE_i\right] \end{aligned}

Bartlett证明，在大样本下（一般ni≥5n_i\ge5），B∼χa−12B\sim\chi_{a-1}^2。

当然，针对样本小于5的情况，Box提出了修正的Bartlett检验。此时定义：

B~=f2BCf1(A−BC)f1=a−1f2=a+1(C−1)2A=f22−C+2/f2BC同上\begin{aligned} \widetilde B&=\frac{f_2BC}{f_1(A-BC)}\\ f_1&=a-1\\ f_2&=\frac{a+1}{(C-1)^2}\\ A&=\frac{f_2}{2-C+2/f_2}\\ &BC同上 \end{aligned}

Box证明了B~∼Ff1,f2\widetilde B\sim F_{f_1,f_2}。如果f2f_2不是整数，可以通过FF分布的分位数表进行内插法得到结果。

本讲的内容到这里就结束了，但是要注意，方差分析远远不止这些，本讲只能算是引入。如果要使用更高级的方差分析，还需要后续知识的补充。