线性回归（二）-单因素方差分析

方差分析模型（ANOVA）可以看成是线性模型的一类，但它重点考虑的是离散自变量对连续因变量的影响。

基本模型

我们有 $k$ 组数据，每组有 $n_i$ 个观察值。基本模型为 $$y_{i,j}={\mu }_i+e_{i,j},i=1,\cdots ,k,j=1,\cdots ,n_i,$$ 其中 $e_{i,j}\sim N(0,{\sigma }^2)$.

如果令 $\eta =\sum {\mu }_i/k,{\tau }_i={\mu }_i-\eta$，便可以写成 $$y_{i,j}=\eta +{\tau }_i+e_{i,j}.$$

模型求解

我们直接把他看成多元线性回归问题，自变量 $x$ 使用 one-hot 编码，设计矩阵为 $$X=\left(\begin{array}{c}{X}_1\\ ⋮\\ X_k\end{array}\right),$$ 其中 $$X_i=\left(\begin{array}{ccccc}0& \cdots & 1& \cdots & 0\\ ⋮& \ddots & ⋮& \ddots & ⋮\\ 0& \cdots & 1& \cdots & 0\end{array}\right).$$ 即只有第 $i$ 列是 $1$，其他都是 $0$. 显然 $X$ 是列满秩的。计算可得 $$(X^{T}X{)}^{-1}=\left(\begin{array}{ccc}\frac{1}{n_1}& & \\ & \ddots & \\ & & \frac{1}{n_k}\end{array}\right).$$

于是有最优估计（请自行推导）： $${\mu }_i={\overline{y}}_{i,\cdot }.$$ 同样有方差估计 $${\hat{\sigma }}^2={\displaystyle \frac{1}{N-k}}\sum _i\sum _j(y_{i,j}-{\overline{y}}_{i,\cdot }{)}^2.$$

ANOVA 表

记 $$\begin{array}{rl}{S}^2& =\sum _i\sum _j(y_{i,j}-y_{\cdot ,\cdot }{)}^2,\\ S_1^2& =\sum _i{n}_i(y_{i,\cdot }-y_{\cdot ,\cdot }{)}^2,\\ S_2^2& =\sum _i\sum _j(y_{i,j}-y_{i,\cdot }{)}^2.\end{array}$$ 类似地，有方差分解： $$S^2=S_1^2+S_2^2.$$ 我们可以得到 ANOVA 表：

假设检验

$F$-检验

$H_0:{\mu }_1=\cdots ={\mu }_k$（等价于 ${\tau }_1=\cdots ={\tau }_k=0$）. 构造 $F$-统计量： $$F={\displaystyle \frac{S_1^2/(k-1)}{S_2^2/(N-k)}}\sim F_{k-1,N-k}.$$ $F$ 足够大时拒绝 $H_0$.

$t$-检验

传统的 $t$ 检验是类似的。这里我们更关心两组之间有无明显差异，于是令 $H_0:{\mu }_i={\mu }_j$. 注意到 $\mathrm{V}\mathrm{a}\mathrm{r}({\mu }_i-{\mu }_j)={\displaystyle \frac{1}{n_i}}+{\displaystyle \frac{1}{n_j}}$, 于是构造统计量 $$T_{i,j}={\displaystyle \frac{{\overline{y}}_{i,\cdot }-{\overline{y}}_{j,\cdot }}{\sqrt{{\hat{\sigma }}^2(1/n_i+1/n_j)}}}\sim t_{N-k}.$$ 当 $|t_{i,j}|>t_{N-k,\alpha /2}$ 时拒绝 $H_0$.

独立样本 $t$ 检验和配对样本 $t$ 检验的主要区别是样本之间有无配对关系。如果是把一堆人平均分成两组，就是独立样本，如果是在同一个人身上先后测了两次数据，则应该是配对样本。独立样本 $t$ 检验的模型是 $y_{i,j}\sim N({\mu }_i,{\sigma }^2)$，而配对样本则是 $y_{1,i}-y_{0,i}\sim N(\mu ,{\sigma }^2)$. 这导致了自由度直接差了 2 倍。本文讨论的 ANOVA 可以看成是独立样本 $t$ 检验的拓展，组数是 2 的 ANOVA 和独立样本 $t$ 检验是等价的。

多重检验

如果我们想要检验多组数据之间有无差异，比如 $H_0:{\mu }_1={\mu }_2={\mu }_3={\mu }_4$，就要做 $6$ 组两两之间的 $t$-检验。任意检验犯第一类错误的概率是 $\alpha$, 则整体犯错的概率显然大于 $\alpha$. 这怎么办呢？

Bonferroni Method

简单粗暴的方法：如果一共做 $k^{\prime }$ 组检验，就把单个检验的 significant level 设置为 $\alpha /k^{\prime }$，这样整体的 significant level 就会小于 $\alpha$.

Tukey Method

我们先介绍 studentized range q 统计量。设我们从分布 $N(\mu ,{\sigma }^2)$ 中取出了 $k$ 组样本，每组 $n$ 个 , 则 $$q={\displaystyle \frac{{\overline{y}}_{max}-{\overline{y}}_{min}}{S\sqrt{1/n}}}$$ 服从一个叫做 studentized range ($q$) distribution 的东西。它的分布函数可以在 R 中使用 ptukey 得到。

后来 Karmer 将其推广到每组样本数不同的情况（Karmer，1956）——把 $1/n$ 换成 $1/2(1/n_i+1/n_j)$. 具体来说，构造统计量 $$Q=\sqrt{2}\underset{i,j}{max}{\displaystyle \frac{|{\overline{y}}_{i,\cdot }-{\overline{y}}_{j,\cdot }|}{\hat{\sigma }\sqrt{1/n_i+1/n_j}}}\sim Q_{k,N-k}.$$ 其中 $N$ 为总样本数。然后查表就行了。$|Q|$ 越大，越倾向于拒绝 $H_0$.

对比检验

我们可以推广普通 $t$-检验的概念，令 $$C=c_1{\mu }_1+\cdots +c_k{\mu }_k,$$ 其中 $\sum _i{c}_i=0$. $H_0:C=0$，此时可以构造 $t$-统计量： $$T_C={\displaystyle \frac{\sum _i{c}_i{\hat{\mu }}_i}{\sqrt{{\hat{\sigma }}^2\sum _i{c}_i^2/n_i}}}\sim t_{N-k}.$$

随机效应模型浅谈

基本模型

这个例子来自[3]。假设我们现在要调查一棵萝卜的含钙量。先验知识告诉我们，只要随便取一片叶子测定叶子的含钙量，就可以代表整颗萝卜的含钙量。我们取了 4 片叶子，每片进行了 4 次测试，得到了 16 个数据。

我们关注的并不是对比叶片之间的差异，而是想得到整颗萝卜的含钙量。不同的叶片在这里不视为 treatment，而视为 block。假设一共 $N$ 个观测值，分 $k$ 组，建立模型如下： $$y_{i,j}=\eta +a_i+e_{i,j},$$ 其中 $e_{i,j}\stackrel{i.i.d.}{\sim }N(0,{\sigma }_e^2)$ 是一样的，但 $a_i\stackrel{i.i.d.}{\sim }N(0,{\sigma }_a^2)$ 也变成随机的了。容易得到 $$\mathrm{V}\mathrm{a}\mathrm{r}(y_{i,j})={\sigma }_e^2+{\sigma }_a^2.$$ 同组的不同观察值不再独立，因为 $j\ne k$ 时， $$\mathrm{C}\mathrm{o}\mathrm{v}(y_{i,j},y_{i,k})={\sigma }_a^2.$$ 显然 $$\hat{\eta }={\overline{y}}_{\cdot ,\cdot }.$$

ANOVA 表

这里我们不加证明地给出 ANOVA 表。

证法类似，写矩阵硬算就行，挺麻烦的，还需要计算技巧。这里最后要落到 $N+k$ 维的独立正态分布上，$N$ 个 $e$ 和 $k$ 个 $a$. 所以矩阵就不是正方形的了，不过 $A^{T}A$ 仍然是正方形的。注意只有在 $n_i$ 相等时，$S_1^2$ 才是卡方分布，否则它只是方差不同的几个正态分布的和。

对于 $n_i=n$ 的情况：

如果 $n_i$ 各不相等，平方和将不是卡方分布，但我们仍可以算它的期望，只需把 $n$ 换成 ${\displaystyle \frac{N}{k}}-{\displaystyle \frac{\sum _{i=1}^k{n}_i^2}{kN}}$.

有了这个，我们可以给出方差的无偏估计（$n_i$ 不相等时按照上面的替换）： $$\begin{array}{rl}{\hat{\sigma }}_a^2& ={\displaystyle \frac{1}{n}}({\displaystyle \frac{S_1^2}{k-1}}-{\displaystyle \frac{S_2^2}{N-k}}),\\ {\hat{\sigma }}_e^2& ={\displaystyle \frac{S_2^2}{N-k}}.\end{array}$$

注：${\hat{\sigma }}_a^2$ 有可能是负的，但方差是负的实在天理难容，这时候我们可以用 $0$ 作为估计值，并骂一句“他娘的”。

$F$-检验

此时我们关注的是 $H_0:{\sigma }_a^2=0$. 只要推翻了 $H_0$，就说明单个叶片的观测值可以代表整体了。在 $H_0$ 成立时，请自行证明以下统计量服从 $F$-分布： $$F={\displaystyle \frac{S_1^2/(k-1)}{S_2^2/(N-k)}}\sim F(k-1,N-k).$$ 这和固定效应模型是一样的！当 $F$ 足够大时拒绝 $H_0$.

参考文献

1. https://en.wikipedia.org/wiki/Tukey%27s_range_test
2. Karmer, C. Y. (1956). Extension of multiple range tests to group means with unequal numbers of replication. Biometrics, 12, 307-310.
3. https://faculty.franklin.uga.edu/dhall/sites/faculty.franklin.uga.edu.dhall/files/STAT8200-Fall13-lec2.pdf