方差分析模型(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$.

参考文献

  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