方差分析模型（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= \begin{pmatrix} X_1\\ \vdots\\ X_k \end{pmatrix}, $$ 其中 $$ X_i= \begin{pmatrix} 0 & \cdots & 1 & \cdots & 0\\ \vdots&\ddots&\vdots&\ddots&\vdots\\ 0 & \cdots & 1 & \cdots & 0\\ \end{pmatrix}. $$ 即只有第 $i$ 列是 $1$，其他都是 $0$. 显然 $X$ 是列满秩的。计算可得 $$ (X^TX)^{-1}= \begin{pmatrix} \frac{1}{n_1} &&\\ &\ddots&\\ &&\frac{1}{n_k} \end{pmatrix}. $$
于是有最优估计（请自行推导）： $$ \mu_i=\bar y_{i,\cdot}. $$ 同样有方差估计 $$ \hat\sigma^2=\dfrac{1}{N-k}\sum_i\sum_j (y_{i,j}-\bar y_{i,\cdot})^2. $$
ANOVA 表
记 $$ \begin{align} 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{align} $$ 类似地，有方差分解： $$ S^2=S_1^2+S_2^2. $$ 我们可以得到 ANOVA 表：
| 来源 | 自由度 | 平方和 | 平方和均值 |
|---|---|---|---|
| 处理组（Treatments） | $k-1$ | $S_1^2$ | $S_1^2/(k-1)$ |
| 误差（Error） | $N-k$ | $S_2^2$ | $S_2^2/(N-k)$ |
| 合计 | $N-1$ | $S^2$ |
假设检验
$F$-检验
$H_0:\mu_1=\cdots=\mu_k$（等价于 $\tau_1=\cdots=\tau_k=0$）. 构造 $F$-统计量： $$ F=\dfrac{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$. 注意到 ${\rm Var}(\mu_i-\mu_j)=\dfrac{1}{n_i}+\dfrac{1}{n_j}$, 于是构造统计量 $$ T_{i,j}=\dfrac{\bar y_{i,\cdot}-\bar y_{j,\cdot}}{\sqrt{\hat\sigma^2\left(1/n_i+1/n_j\right)}}\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’$ 组检验，就把单个检验的 significant level 设置为 $\alpha/k’$，这样整体的 significant level 就会小于 $\alpha$.
Tukey Method
我们先介绍 studentized range q 统计量。设我们从分布 $N(\mu,\sigma^2)$ 中取出了 $k$ 组样本，每组 $n$ 个 , 则
$$
q=\dfrac{\bar y_\max-\bar 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}\max_{i,j}\dfrac{|\bar y_{i,\cdot}-\bar 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=\dfrac{\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}\overset{i.i.d.}\sim N(0,\sigma_e^2)$ 是一样的，但 $a_i\overset{i.i.d.}\sim N(0,\sigma_a^2)$ 也变成随机的了。容易得到 $$ {\rm Var}(y_{i,j})=\sigma_e^2+\sigma_a^2. $$ 同组的不同观察值不再独立，因为 $j\ne k$ 时， $$ {\rm Cov}(y_{i,j},y_{i,k})=\sigma_a^2. $$ 显然 $$ \hat\eta=\bar y_{\cdot,\cdot}. $$
ANOVA 表
这里我们不加证明地给出 ANOVA 表。
证法类似，写矩阵硬算就行，挺麻烦的，还需要计算技巧。这里最后要落到 $N+k$ 维的独立正态分布上，$N$ 个 $e$ 和 $k$ 个 $a$. 所以矩阵就不是正方形的了，不过 $A^TA$ 仍然是正方形的。注意只有在 $n_i$ 相等时，$S_1^2$ 才是卡方分布，否则它只是方差不同的几个正态分布的和。
对于 $n_i=n$ 的情况：
| 来源 | 自由度 | 平方和 | 平方和均值（MS） | E(MS) |
|---|---|---|---|---|
| 区组（Block） | $k-1$ | $S_1^2$ | $S_1^2/(k-1)$ | $\sigma_e^2+n\sigma_a^2$ |
| 误差（Error） | $N-k$ | $S_2^2$ | $S_2^2/(N-k)$ | $\sigma_e^2$ |
| 合计 | $N-1$ | $S^2$ |
如果 $n_i$ 各不相等，平方和将不是卡方分布，但我们仍可以算它的期望，只需把 $n$ 换成 $\dfrac{N}{k}-\dfrac{\sum_{i=1}^k n_i^2}{kN}$.
有了这个，我们可以给出方差的无偏估计（$n_i$ 不相等时按照上面的替换）： $$ \begin{align} \hat\sigma_a^2&=\dfrac{1}{n}\left(\dfrac{S_1^2}{k-1}-\dfrac{S_2^2}{N-k}\right),\\ \hat\sigma_e^2&=\dfrac{S_2^2}{N-k}. \end{align} $$
注：$\hat\sigma_a^2$ 有可能是负的，但方差是负的实在天理难容，这时候我们可以用 $0$ 作为估计值，并骂一句“他娘的”。
$F$-检验
此时我们关注的是 $H_0:\sigma_a^2=0$. 只要推翻了 $H_0$，就说明单个叶片的观测值可以代表整体了。在 $H_0$ 成立时，请自行证明以下统计量服从 $F$-分布： $$ F=\dfrac{S_1^2/(k-1)}{S_2^2/(N-k)}\sim F(k-1,N-k). $$ 这和固定效应模型是一样的！当 $F$ 足够大时拒绝 $H_0$.
参考文献
- https://en.wikipedia.org/wiki/Tukey%27s_range_test
- Karmer, C. Y. (1956). Extension of multiple range tests to group means with unequal numbers of replication. Biometrics, 12, 307-310.
- https://faculty.franklin.uga.edu/dhall/sites/faculty.franklin.uga.edu.dhall/files/STAT8200-Fall13-lec2.pdf