LRA系列的第三篇,介绍方差分析。协方差分析我作为单独一篇抽出来,因为我觉得其十分有用而且深刻。

手敲\(\LaTeX\)难免出现纰漏,有任何疑似错误或者不清楚的地方请直接在下方评论区留言,谢谢各位读者。

方差分析一个简单的理解,就是当自变量的取值离散化后的线性回归模型。我们将这时的自变量称为因子(Factor),而将因子的取值成为水平(level)。单因子模型中的因子的水平,或者多因子模型中因子水平的组合(笛卡儿积)成为一个处理(treatment).

单因子模型

模型假设

\[Y_{ij}=\mu_i+\varepsilon_{ij}\]

其中\(i=1,\ldots,r\)为levels, \(j=1,\ldots,n_i\)\(i\)组内的索引,我们通常假设\(\varepsilon_{ij}\sim \mathcal{N}(0,\sigma^2)\)。也就是说,各组内分别服从\(\mathcal{N}(\mu_i,\sigma^2)\)的正态分布。这个方差是组间共享的。

值得一提的是,所有多因子模型都可以转化为单因子模型,只需将所有的因子水平的组合考虑进来。但这样的模型会引入太多的参数,我们将在下一节讨论这一问题。

模型求解

同样的,我们的最小二乘解优化的目标为\(Q=\sum_{i=1}^{r}\sum_{j=1}^{n_i}(Y_{ij}-\mu_i)^2\),求导立得其最小二乘解

\[\widehat{\mu_i}=\bar{Y}_{i\cdot}=\sum_{j=1}^{n_i}Y_{ij} / n_i\]

此时\(\mathrm{SSE}=\sum_{i=1}^{n}\sum_{j=1}^{n_i}(Y_{ij}-\bar{Y}_{i\cdot})^2\)\(\sigma^2\)的无偏估计\(\mathrm{MSE}=\mathrm{SSE}/(n_T-r)\),其中\(n_T=\sum_{i=1}^{r}n_i\)为所有样本的数量。\(\mathrm{MSE}\)的另一种写法为

\[\mathrm{MSE}=\frac{1}{\sum_{i=1}^{r}n_i-1}\sum_{i=1}^{r}(n_i-1)s_i^2,\quad s_i^2\overset{\Delta}{=}\frac{1}{n_i-1}\sum_{j=1}^{n_i}(Y_{ij}-\bar{Y}_{i\cdot})^2\]

如果读者之前接触过双样本的检验,应当对这一形式并不陌生。

同时,我们知道有关系\(\bar{Y}_{i\cdot}\sim\mathcal{N}(\mu_i,\sigma^2/n_i)\),由于\(\sigma^2\)未知,推断时我们使用下面的关系

\[\frac{\bar{Y}_{i\cdot}-\mu_i}{\sqrt{\mathrm{SSE}/n_i}}\sim t_{n_T-r}\]

方差分析

定义\(\widehat{Y}_{\cdot\cdot}=\frac{1}{n_T}\sum_{i=1}^{r}\sum_{j=1}^{n_i}Y_{ij}\),我们有如下分解

\[Y_{ij}-\widehat{Y}_{\cdot\cdot}=(Y_{ij}-\widehat{Y}_{i\cdot})+(\widehat{Y}_{i\cdot}-\widehat{Y}_{\cdot\cdot})\] \[\sum_{i=1}^{r}\sum_{j=1}^{n_i}(Y_{ij}-\widehat{Y}_{\cdot\cdot})^2=\sum_{i=1}^{r}\sum_{j=1}^{n_i}(Y_{ij}-\widehat{Y}_{i\cdot})^2+\sum_{i=1}^{r}n_j(\widehat{Y}_{i\cdot}-\widehat{Y}_{\cdot\cdot})\] \[\mathrm{SSTO}:(n_T-1)=\mathrm{SSE}:(n_T-r)+\mathrm{SSTR}:(r-1)\]

Source SS df MS F
Model SSTR \(r-1\) SSTR/\((r-1)\) MSTR/MSE
Error SSE \(n_T-r\) SSE/\((n_T-r)\)
Total SSTO \(n_T-1\)

对于\(H_0:\mu_1=\cdots=\mu_r\),使用\(F\)检验

\[F^\ast=\frac{\mathrm{MSTR}}{\mathrm{MSE}}\sim F_{r-1,n_T-r}\]

和多元正态的联系

对于单因子模型\(Y_{ij}=\mu_i+\varepsilon_{ij}\),很容易将其建模为多元线性回归模型的形式(无截距项)

\[Y_i=\mu_1X_1+\mu_2X_2+\cdots+\mu_rX_r+\varepsilon_{ij}\],其中当数据点为水平\(k\)\(X_k=1, X_j=0(j\neq k)\)

若使用有截距项的形式,则相当于现在的线性模型有\(r+1\)个参数,但该模型只有\(r\)个levels即\(r\)个自由度,因此这\(r+1\)个参数必然有内部的约束。上面的无截距项给出的是\(\mu=0\)的约束。一般的,还有两种约束是我们经常采用的。

  • \(\sum_{i=1}^{r} \mu_i=0\),也就是将组与组之间的平均水平提到\(\mu\)中,此时\(\mu_i\)表示的是组水平相对于平均水平的差值。
  • \(\mu_r=0\),相当于将第\(r\)个水平取为基准,此时的\(\mu\)为第\(r\)组的绝对水平,而\(\mu_i\)则是第\(i\)组对第\(r\)组的相对水平。

作为练习,读者可以考虑下在这样的模型里\(X_i\)应当如何定义。更详细的代码和讨论在 协方差分析可以见到。

组间效应的分析

\(\mu_i\)

\(1-\alpha\) 置信区间为

\[\widehat{Y}_{i\cdot}\pm t_{1-\alpha/2,n_T-r}s\{\widehat{Y}_{i\cdot}\}\]

其中 \(s^2\{\widehat{Y}_{i\cdot}\}=\mathrm{MSE}/n_i\),这是很好理解的,因为所有的数据点的方差都是\(\sigma^2\),因此均值的方差就是\(\sigma^2/n_i\),使用\(\mathrm{MSE}\)代替就可以构造\(t\)置信区间。

\(\mu_i-\mu_j\)

类似于双样本的检验,其\(1-\alpha\) 置信区间为

\[(\widehat{Y}_{i\cdot}-\widehat{Y}_{j\cdot})\pm t_{1-\alpha/2,n_T-r}\sqrt{\mathrm{MSE}(\frac{1}{n_i}+\frac{1}{n_j})}\]

\(\theta\overset{\Delta}{=}\sum_{i=1}^{r}c_i\mu_i\)

事实上,使用\(\widehat{\theta}=\sum_{i=1}^{r}c_i\widehat{Y}_{i\cdot}\),有

\[\mathbb{E}(\widehat{\theta})=\sum_{i=1}^{r}c_i\mu_i=\theta\] \[s^2\{\widehat{\theta}\}=\sum_{i=1}^{r}c_i^2\frac{\sigma^2}{n_i}=\sigma^2\sum_{i=1}^{r}\frac{c_i^2}{n_i}\]

contrast

上一节中,当\(\sum_{i=1}^{n}c_i=0\)时,我们称这样的\(\theta\)为一个contrast。讨论这个是因为,当一个contrast的置信区间包括\(0\)时,意味着有一个水平是多余的,即不需要那么多参数就可以表现这个模型。典型的contrast\(\mu_1-\mu_2\), \(\mu_1-(\mu_2+\mu_3)/2\)等。其检验和上一节完全一样。

Tukey Method

Tukey Method 针对的是 所有的 \(\mu_i-\mu_j\)给出一个最窄的联合置信区间,而且,当\(n_i\equiv n\)时,能给出精确\(1-\alpha\) 置信区间。

对于\(Z_1,\ldots,Z_r\sim\mathcal{N}(0,1), \nu s^2\sim\chi_\nu^2\)且彼此独立,此时我们可以导出studentized range这个分布,并且有 \[\mathbb{P}(\frac{\max_i Z_i-\min_i Z_i}{s}\leq q_{1-\alpha,r,\nu})=1-\alpha\] 利用这一背景可以帮助我们构造一个更加精确的置信区间

考察\(n_i\equiv n\)的单因子模型,有

\[\sqrt{n}(\widehat{Y}_{i\cdot}-\mu_i)/\sigma\sim\mathcal{N}(0,1)\] \[(n_T-r)\mathrm{MSE}/\sigma^2\sim\chi_{n_T-r}^2\]

使用上面的结论,有

\[\mathbb{P}(\frac{\max_{i,j}\lvert(\widehat{Y}_{i\cdot}-\widehat{Y}_{j\cdot})-(\mu_i-\mu_j)\rvert }{\sqrt{\mathrm{MSE}/n}}\leq q_{1-\alpha,r,\nu})=1-\alpha\]

因此,Tukey 给出的\(\mu_i-\mu_j\)联合置信区间为\((\widehat{Y}_{i\cdot}-\widehat{Y}_{j\cdot})+Bs\{(\widehat{Y}_{i\cdot}-\widehat{Y}_{j\cdot})\}\),其中

\[B=q_{1-\alpha,r,n_T-r}/\sqrt{2}\] \[s\{(\widehat{Y}_{i\cdot}-\widehat{Y}_{j\cdot})\}=\sqrt{\mathrm{MSE}(\frac{1}{n_i}+\frac{1}{n_j})}\]

Scheffe Method

这是一个更加万能的构造contrast的置信区间的方法(由于其万能性也意味着通常区间会更长),如果读者了解过多元线性回归分析,应当对这一形式并不陌生。事实上,它给出了所有contrast精确的(\(1-\alpha\))联合置信区间。

\[\mathbb{P}\left(\frac{(\sum_{i=1}^{r}c_i(\widehat{Y}_{i\cdot}-\mu_i))^2}{(r-1)s^2\{\sum_{i=1}^{r}c_i\widehat{Y}_{i\cdot}\}}\leq F_{1-\alpha, r-1, n_T-r},\,\forall \sum_{i=1}^{r}c_i=0\right)=1-\alpha\]

因此,对于contrast\(\sum_{i=1}^{r}c_i\mu_i\)的置信区间为\((\sum_{i=1}^{r}c_i\widehat{Y}_{i\cdot})\pm Bs\{\sum_{i=1}^{r}c_i\widehat{Y}_{i\cdot}\}\),其中

\[B=\sqrt{(r-1)F_{1-\alpha,r-1,n_T-r}}\] \[s\{\sum_{i=1}^{r}c_i\widehat{Y}_{i\cdot}\}=\sqrt{\mathrm{MSE}\sum_{i=1}^{r}c_i^2/n_i}\]

双因子模型

模型假设

\[Y_{ijk}=\mu_{ij}+\varepsilon_{ijk}\]

其中\(\varepsilon_{ijk}\sim\mathcal{N}(0,\sigma^2)\)

对于分别有\(a,b\)个水平的因子\(A,B\),我们记\(\mu_{i,j}\)为组内水平,有

\[\mu_{i,j}=\mu+\alpha_i+\beta_j+(\alpha\beta)_{ij}\]

同样的,这些参数中有内部约束 \[\sum_{i=1}^{a}\alpha_i=0\] \[\sum_{j=1}^{b}\beta_j=0\] \[\sum_{i=1}^{a}(\alpha\beta)_{ij}=\sum_{j=1}^{b}(\alpha\beta)_{ij}=0\]

让我们分析一下自由度,模型中一共有了\(1+a+b+ab\)个参数,三条约束分别减少了\(1+1+(a+b-1)\)个自由度,因此剩下\(ab\)个自由度,这是很自然的。

\((\alpha\beta)_{ij}\equiv 0,\forall i,j\),则称为加性模型(additive model). 加性模型最形象化的理解就是一个平行四边形。对于加性模型,其参数减少为\(1+(a-1)+(b-1)=a+b-1\)个。

模型求解

由参数约束不难得到下列信息

\[\mu=\frac{1}{ab}\sum_{i=1}^{a}\sum_{j=1}^{b}\mu_{i,j}=\mu_{\cdot\cdot}\] \[\alpha_i=\frac{1}{b}\sum_{j=1}^{b}\mu_{i,j}-\mu_{\cdot\cdot}=\mu_{i\cdot}-\mu_{\cdot\cdot}\] \[\beta_j=\frac{1}{a}\sum_{i=1}^{a}\mu_{i,j}-\mu_{\cdot\cdot}=\mu_{\cdot j}-\mu_{\cdot\cdot}\] \[(\alpha\beta)_{ij}=\mu_{\cdot\cdot}-\mu_{i\cdot}-\mu_{\cdot j}+\mu_{\cdot\cdot}\]

\(\mu\)代表总体均值,\(\alpha_i,\beta_j\)称为主效应(main effects),\((\alpha\beta)_{ij}\)称为交叉效应(interaction)。

Balanced Design

如果\(n_{ij}\equiv n\),我们称其为平衡的(balanced design),此时我们可以得到极大似然估计

\[\widehat{\mu}_{i,j}=\bar{Y}_{ij\cdot}=\frac{1}{n}\sum_{k=1}^{n}Y_{ijk}\] \[\mu=\frac{1}{ab}\sum_{i=1}^{a}\sum_{j=1}^{b}\widehat{\mu}_{i,j}=\bar{Y}_{\cdot\cdot\cdot}\] \[\widehat{\alpha}_i=\frac{1}{b}\sum_{j=1}^{b}\widehat{\mu}_{i,j}-\widehat{\mu}_{\cdot\cdot}=\bar{Y}_{i\cdot\cdot}-\bar{Y}_{\cdot\cdot\cdot}\] \[\widehat{\beta}_j=\frac{1}{a}\sum_{i=1}^{a}\widehat{\mu}_{i,j}-\widehat{\mu}_{\cdot\cdot}=\bar{Y}_{\cdot j\cdot}-\bar{Y}_{\cdot\cdot\cdot}\] \[(\widehat{\alpha\beta})_{ij}=\bar{Y}_{ij\cdot}-\bar{Y}_{i\cdot\cdot}-\bar{Y}_{\cdot j\cdot}+\bar{Y}_{\cdot\cdot\cdot}\]

此时有\(\mathrm{SSE}=\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{n}(Y_{ijk}-\bar{Y}_{ij\cdot})^2\), 且

\[\mathrm{MSE}=\frac{1}{ab(n-1)}\mathrm{SSE}=\frac{1}{ab}\sum_{i=1}^{a}\sum_{j=1}^{b}s_{ij}^2\]

其中\(s_{ij}^2=\frac{1}{n-1}\sum_{k=1}^{n}(Y_{ijk}-Y_{ij\cdot})^2\)

事实上,此时该模型和单因子模型并无差别。

平衡设计下的加性模型

此时最小化目标为\(Q=\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{n}(Y_{ijk}-\mu-\alpha_i-\beta_j)^2\),得最小二乘解

\[\widehat{\mu}=\bar{Y}_{\cdot\cdot\cdot}\] \[\widehat{\alpha}_i=\bar{Y}_{i\cdot\cdot}-\bar{Y}_{\cdot\cdot\cdot}\] \[\widehat{\beta}_j=\bar{Y}_{\cdot j\cdot}-\bar{Y}_{\cdot\cdot\cdot}\]

此时\(\mathrm{SSE}=\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{n}(Y_{ijk}-\bar{Y}_{i\cdot\cdot}-\bar{Y}_{\cdot j\cdot}+\bar{Y}_{\cdot\cdot\cdot})^2\)\(abn-(a+b+1)\)个自由度,\(\mathrm{MSE}=\mathrm{SSE}/(abn-a-b-1)\)\(\sigma^2\)的无偏估计。

方差分析

\[Y_{ijk}-\bar{Y}_{\cdot\cdot\cdot}=(\bar{Y}_{i\cdot\cdot}-\bar{Y}_{\cdot\cdot\cdot}) + (\bar{Y}_{\cdot j\cdot}-\bar{Y}_{\cdot\cdot\cdot}) + (\bar{Y}_{ij\cdot}-\bar{Y}_{i\cdot\cdot}-\bar{Y}_{\cdot j\cdot}+\bar{Y}_{\cdot\cdot\cdot}) + (\bar{Y}_{ijk}-\bar{Y}_{ij\cdot})=\widehat{\alpha}_i+\widehat{\beta}_j+(\widehat{\alpha\beta})_{ij}+e_{ijk}\] \[\sum_{ijk}(Y_{ijk}-\bar{Y}_{\cdot\cdot\cdot})^2=\sum_{ijk}\widehat{\alpha}_i^2+\sum_{ijk}\widehat{\beta}_j^2+\sum_{ijk}(\widehat{\alpha\beta})_{ij}^2+\sum_{ijk}e_{ijk}^2\] \[\mathrm{SSTO}:(abn-1)=\mathrm{SSA}:(a-1)+\mathrm{SSB}:(b-1)+\mathrm{SSAB}:(a-1)(b-1)+\mathrm{SSA}:ab(n-1)\]

值得一提的是,\(\mathrm{SSA},\mathrm{SSB},\mathrm{SSAB}\)彼此“正交”,因此其不像多元线性回归一样会因为引入变量的先后顺序的不同而导致Type-I SS的不同。

对应地写出均方误差

\[\mathrm{MSA}=\frac{1}{a-1}\mathrm{SSA}=\frac{bn}{a-1}\sum_{i=1}^{a}\widehat{\alpha}_i^2\] \[\mathrm{MSB}=\frac{1}{b-1}\mathrm{SSB}=\frac{an}{b-1}\sum_{j=1}^{b}\widehat{\beta}_j^2\] \[\mathrm{MSA}=\frac{1}{(a-1)(b-1)}\mathrm{SSAB}=\frac{n}{(a-1)(b-1)}\sum_{i=1}^{a}\sum_{j=1}^{b}(\widehat{\alpha\beta})_{ij}^2\]

同时,可以证明

\[\mathbb{E}(\mathrm{MSE})=\sigma^2\] \[\mathbb{E}(\mathrm{MSA})=\sigma^2+\frac{bn}{a-1}\sum_{i=1}^{a}{\alpha_i}^2\] \[\mathbb{E}(\mathrm{MSB})=\sigma^2+\frac{an}{b-1}\sum_{j=1}^{b}{\beta_j}^2\] \[\mathbb{E}(\mathrm{MSAB})=\sigma^2+\frac{n}{(a-1)(b-1)}\sum_{i=1}^{a}\sum_{i=1}^{a}\sum_{j=1}^{b}(\alpha\beta)_{ij}^2\]

因此我们可以做\(F\)检验去检测交叉项以及主效应,例如\(H_0:\sum_{i=1}^{a}\sum_{i=1}^{a}\sum_{j=1}^{b}(\alpha\beta)_{ij}^2=0\),只需使用

\[F^\ast=\frac{\mathrm{MSAB}}{\mathrm{MSE}}\sim F_{(a-1)(b-1),ab(n-1)}\]

如果检测的结果发现交互项的效应确实不明显,说明我们可以用加性模型替代。此时,\(\mathrm{MSAB},\mathrm{MSE}\)合并在一起构成新的\(\mathrm{MSE}\),即

\[\frac{\mathrm{MSAB}+\mathrm{MSE}}{(a-1)(b-1)+ab(n-1)}=\frac{\sum_{i=1}^{a}\sum_{j=1}^{b}\sum_{k=1}^{n}(Y_{ijk}-\bar{Y}_{i\cdot\cdot}-\bar{Y}_{\cdot j\cdot}+\bar{Y}_{\cdot\cdot\cdot})^2}{abn-a-b+1}\]

这和上面谈到的加性模型的结论是一致的。

组间效应的分析

\(\mu_{i\cdot}\)\(\mu_{ij}\)

只需注意到\(\mathbb{E}(\widehat{Y_{i\cdot\cdot}})=\mu_{i\cdot},\,s^2\{\widehat{Y}_{i\cdot\cdot}\}=\frac{1}{bn}\mathrm{MSE}\)即可求解。事实上,对于\(\mu_{i\cdot}\)的线性组合的估计,其方法和单因子模型中提到的完全一样。

而对于\(\mu_{ij}\),只需注意\(s^2\{\widehat{Y}_{ij\cdot}\}=\frac{1}{n}\mathrm{MSE}\)即可,其与方法和单因子模型一致。

两个特殊的变种

单样本情形

\(n=1\)时,也就是说我们只有\(ab\)个样本,此时使用\(ab\)个参数的全模型显然是不合适的,因此其交叉项必须丢掉以减少\((a-1)(b-1)\)个参数

此时,\(\widehat{\mu}=\bar{Y}_{\cdot\cdot},\widehat{\alpha}_i=\bar{Y}_{i\cdot}-\bar{Y}_{\cdot\cdot},\widehat{\beta}_j=\bar{Y}_{\cdot j}-\bar{Y}_{\cdot\cdot},\mathrm{SSE}=\sum_{i=1}^{a}\sum_{j=1}^{b}(Y_{ij}-\bar{Y}_{i\cdot}-\bar{Y}_{\cdot j} +\bar{Y}_{\cdot\cdot})^2,\mathrm{MSE}=\mathrm{SSE}/\left[(a-1)(b-1)\right]\)

然而有时候加性模型这个过强的假设并不能很好地刻画交叉项的效应,一个常用的解决方案是假设\((\widehat{\alpha\beta})_{ij}=D\alpha_i\beta_j\),即引入多一个参数\(D\)。可以求出\(D\)的最小二乘估计为

\[\widehat{D}=\frac{\sum_{i=1}^{a}\sum_{j=1}^{b}\widehat{\alpha}_i\widehat{\beta}_jY_{ij}}{\sum_{i=1}^{a}\widehat{\alpha}_i^2\sum_{j=1}^{b}\widehat{\beta}_j^2}\]

显然我们需要检验这样的一个参数的检验是不是必要的。使用\(F\)检验的思想去检验\(H_0:D=0\),计算

\[\mathrm{SSAB}^\ast=\sum_{i=1}^{a}\sum_{j=1}^{b}\widehat{D}^2\widehat{\alpha}_i^2\widehat{\beta}_j^2\] \[F^\ast=\frac{\mathrm{SSAB}^\ast}{(\mathrm{SSE}-\mathrm{SSAB}^\ast)/(ab-a-b)}\sim F_{1,ab-a-b}\]

嵌套因子模型

我们之前的模型都假设因子\(A\)和因子\(B\)可以随意组合出\(ab\)种可能,但现实中仍有一种情况是因子\(B\)嵌套在因子\(A\)中。一个简单的例子就是当你在研究不同汽车的某个性能指标如油耗时,因子\(A\)为汽车的品牌,因子\(B\)为汽车的型号,显然因子\(B\)依赖于因子\(A\),也就是说此时的模型为

\[Y_{ijk}=\mu+\alpha_i+\beta_{j(i)}+\varepsilon_{ijk}\]

其中\(i=1,\ldots,a,\,j=1,\ldots,b_i,\,k=1,\ldots,n_{ij}\)。约束为\(\sum_{i=1}^{n}\alpha_i=0,\,\sum_{j=1}^{b_i}\beta_{j(i)}=0\)

出于简单我们考虑平衡的设计,可以求得最小二乘估计

\[\widehat{\mu}=\widehat{Y}_{\cdot\cdot\cdot}\] \[\widehat{\alpha}_i=\widehat{Y}_{i\cdot\cdot}-\widehat{Y}_{\cdot\cdot\cdot}\] \[\widehat{\beta}_{j(i)}=\widehat{Y}_{ij\cdot}-\widehat{Y}_{i\cdot\cdot}\] \[\mathrm{SSA}=\sum_{ijk}\widehat{\alpha}_i^2\] \[\mathrm{SSB(A)}=\sum_{ijk}\widehat{\beta}_{j(i)}^2\] \[\mathrm{SSE}=\sum_{ijk}(Y_{ijk}-\bar{Y}_{ij\cdot})^2\]