这部分我写的比较抽象,估计只有学过的才看得懂,没学过的可能云里雾里,不过我功力实在有限,难以用自然语言表达清楚这件事……

模型

Multiple way ANOVA

这是一种特殊的 multiple way ANOVA. 假设我们有 $k$ 个因子,每个只有 2 个 level,记作 $-$ 和 $+$.

以 $2^3$ 为例,我们有三个因子 $A,B,C$,模型为: $$ \begin{align} y_i=&\mu +\alpha x_{i,A}+\beta x_{i,B} +\gamma x_{i,C}\ &+(\alpha\beta)x_{i,AB}+(\alpha\gamma)x_{i,AC}+(\beta\gamma)x_{i,BC}\ &+(\alpha\beta\gamma)x_{i,ABC}+\varepsilon_i, \end{align} $$ 其中

  • $x_{i,A}$ 在 $A+$ 时取值 $1$,在 $A-$ 时取 $-1$,
  • $x_{i,AB}=x_{i,A}x_{i,B}$,也只有 $1$ 或 $-1$ 两个值,
  • $x_{i,ABC}=x_{i,A}x_{i,B}x_{i,C}$,
  • $\varepsilon_i\overset{i.i.d.}\sim N(0,\sigma^2)$.

注:这里面每个交叉项只有一个,例如 $\alpha\beta$,而不是像一般的模型那样有 $(\alpha\beta)_{ij}$. 这样正好 $2^k$ 个自由参数,就不必额外加约束条件了。

一般地,用 $A_1,\cdots,A_k$ 表示 $k$ 个 factor,模型可以写为 $$ y_i=\sum_{\mathcal{I}\in 2^{{1,2,\cdots,k}}} \theta_{\mathcal{I}}x_{i,\mathcal{I}}. $$ 这里 $\mathcal{I}$ 是指标集。例如 $\mathcal{I}={1,2,3}$,则 $x_{i,\mathcal{I}}$ 指的就是前面的 $x_{i,A_1A_2A_3}$. 特别地,$x_{i,\varnothing}=1,\theta_{\varnothing}=\mu$.

有了这种记号,我们可以方便的写出参数估计(就不推导了): $$ \theta_{\mathcal{I}}=\frac{1}{2}\left(\bar y(x_{\mathcal{I}}=1)-\bar y(x_{\mathcal{I}}=-1) \right). $$ 这么看太抽象,还是回到前面 $2^3$ 的例子,有:

  • $\alpha=\bar y(A+)-\bar y(A-)$,
  • $(\alpha\beta)=\dfrac{1}{2}\left(\bar y(A+,B+)+\bar y(A-,B-)\right)-\dfrac{1}{2}\left(\bar y(A+,B-)+\bar y(A-,B+)\right)$.

我们定义 effect 是系数的二倍,这样就没有 $1/2$ 了。即 $$ E(A_{i_1}\cdots A_{i_m})=2\theta_{i_1,\cdots,i_m}. $$ 例如 $A$ 组和 $B$ 组的交互作用项 $E(AB)=\left(\bar y(A+,B+)+\bar y(A-,B-)\right)-\left(\bar y(A+,B-)+\bar y(A-,B+)\right)$.

注:这里的 $E(AB)$ 就是书上的 $AB$,表示 effect,$E$ 只是个记号,和期望没关系。我们一会要用 $AB$ 表示群中的元素,所以只好这样了。

群论视角

我们想知道,对于一般情况,如何得知某个 effect $E(A_{i_1}\cdots A_{i_m})$ 是怎么算的 —— 即哪些是加,哪些是减?为此,我们需要新的武器。

表一:

  • 第一列是观测值,ab 表示 $A+,B+$ 组的观测结果(每组只有一个),(1) 表示全 $-$ 的观测结果;
  • 右边 8 列是 8 个参数,也就是各种 effect。例如 $AB$ 列表示 $AB$ 的联合 effect 是 $E(AB)=((1)-a-b+ab+c-ac-bc+abc)/4$.

可以看出,每个 effect 都是各个观测值的组合,要么是加要么是减。所以研究这个表格的行和列的性质至关重要。

令 $\Omega={(a_1,\cdots,a_{2^k}):a_i\in{+,-}}$,即由 $+,-$ 组成的 $k$ 维向量。定义 $\Omega$ 上的乘法: $$ (a_1,\cdots,a_{2^k})(b_1,\cdots,b_{2^k})=(a_1b_1,\cdots,a_{2^k}b_{2^k}), $$ 其中 $++=+, +-=-, -+=-,–=+$. 记向量 $X=(2^k-1,2^k-2,\cdots,0)^T$. 记 $B:\mathbb{N}\to{-,+}$. 若 $m$ 的二进制表示中 1 的个数为奇数,则 $B(m)=-$,否则 $B(m)=+$. 我们令 $$ A_i=B(X\ &\ 2^{i-1}). $$ 其中 $&$ 是二进制按位与。$A_1,\cdots, A_k$ 就是表一中的前几列(蓝色的那几列)。利用位运算的性质,有 $$ A_{i_1}\cdots A_{i_m}=\prod_{s=1}^m B(X\ &\ 2^{i_s-1})=B\left( X\ &\ \sum_{s=1}^m 2^{i_s-1}\right).\qquad (1) $$ 设 $$ G=\left{\prod_{i\in\mathcal{I}} A_{i}:\mathcal{I}\subset {1,\cdots,k}\right}. $$

$G$ 上的乘法继承自 $\Omega$. 特别地,定义 $\prod_{i\in\varnothing} A_{i}=\mu=(+,\cdots,+)$. 我们可以迅速得到 $G$ 的一些基本性质(自行证明):

  1. $G$ 在乘法下构成一个 Abel 群,$\mu$ 是单位元;

  2. $G$ 中除单位元外,所有元素的阶都是 2;

  3. $|G|=2^k$;($G$ 中的元素就是表一中的列,它定义了各个 effect 的计算方式)

  4. 更进一步地,有 $G\cong (\mathbb{Z}_2)^k$. ($k$ 个 $\mathbb{Z}_2$ 的直积)

对 $g\in G$,定义重量 $w(g)=g\text{中'-‘的个数}$。其中 $i_k$ 各不相等。根据式(1)不难看出 $w(A_{i_1}\cdots A_{i_m})=1/2^{k-1}$,即 $G$ 中除了 $\mu$ 外的元素都是一半 $+$ 一半 $-$. 这意味着每个 effect 都是一半观测值的均值减掉另一半观测值的均值

Yates 算法

我们有一个快速计算各个 effect 的算法:

  1. 把观测值按标准顺序排成一列
  2. 在右边新建 $k$ 列:
    • 每一列的前半部分:第 $i$ 个元素是前一列的第 $2i-1$ 个和第 $2i$ 个之和;
    • 每一列的后半部分:第 $i$ 个元素是前一列的第 $2i$ 个和第 $2i-1$ 个之差;
  3. 在右边在写一列($k+1$ 列),第一个元素是 $2^k$,其它全是 $2^{k-1}$,
  4. 此时用新建第 $k$ 列除以 $k+1$ 列,就是对应的 effect.

如下表:

标准顺序就是最右边的顺序,$ab$ 表示 $A+,B+,C-$,其它类似。

这个过程是不是感到莫名熟悉?其实它就是快速傅里叶算法的早期版本。原始论文:

Yates, Frank (1937). “The design and analysis of factorial experiments”. Technical Communication No. 35 of the Commonwealth Bureau of Soils. 142 (3585): 90–92. Bibcode:1938Natur.142…90F. doi:10.1038/142090a0. S2CID 23501205.

显著性分析

方案 A

每个 effect 都是一半观测值的均值减掉另一半观测值的均值,所以 $\mathrm{Var}(\dfrac{4\sigma^2}{N})$,$N$ 为总观测样本数。只要有了 $\sigma^2$ 的估计,置信区间就有了。

可设某些高阶 interaction 为 0. 例如对 $2^4$,设三阶和四阶 interaction 为 0,此时我们算出来的 $ABC,ABD,ACD,BCD, ABCD$ 就服从 $N(0,\dfrac{4\sigma^2}{N})$,且独立(这是因为设计矩阵的列正交)。于是可以进行方差估计: $$ \hat\sigma^2=\frac{N}{4}\frac{ABC^2+ABD^2+ACD^2+BCD^2+ ABCD^2}{5}. $$

自由度是 5。注意这里均值是已知的 0,所以自由度就是 $n$。之前是因为均值未知,用样本均值替代,自由度就是 $n-1$ 了。

某个 effect 的置信区间是 $$ \hat\theta\pm t_{5,\frac{\alpha}{2}}\sqrt\frac{4\hat\sigma^2}{N}. $$

方案 B (Daniel)

我们假设所有 effect 都是 0,此时只有均值和误差项。那么共有 $2^k-1$ 个 effect,它们服从 $N(0,\dfrac{4\sigma^2}{N})$. 我们对这 $2^k-1$ 个点进行正态性检验,哪些不合群哪些就显著。

具体地,把 $\hat\theta_i,1\le i\le 2^k-1$ 加上绝对值,然后从小到大排序,然后画 half normal plot。它类似 Q-Q 图,只是有两点不同:

  • 只考虑 $x>0$,相当于截尾正态分布
  • 横坐标是标准 half normal quantile,但纵坐标是 $N(0,\dfrac{4\sigma^2}{N})$ 的 quantile

令 $I=2^k-1$,图上包含的点是 $$ (z_{0.5+0.5[i-0.5]/I},|\hat\theta|_{(i)}). $$ 然后目测一条直线,拟合 0 附近的点,看谁不合群。比如这个例子:

B 不合群,说明 B 显著。

方案 C (Lenth)

这里吧,感觉写的巨详细,我就不抄一遍了。

Ofat Approach

暂略,之后补上

混淆(confounding)

2^k Design in 2 Blocks

Confounding, 直译为“混淆”。在因果推断中表示干扰因素:当我们研究自变量 X 与因变量 Y 的相关性时,存在一个共同的原因 Z 同时影响了 X 和 Y。最典型的就是我们常说的“相关不等于因果”:冰淇淋的销量和溺水率正相关,但其原因是天气炎热导致冰淇淋销量上升,同时也导致人们游泳变多,溺水事故自然也变多了,而非“冰淇淋导致溺水”。

这也是一个经典的逻辑谬误,详见我的另一篇文章中的因果谬误

我们通常考虑 $2^k$ 因子设计。以 $2^2$ 设计为例,两个因子为 $A,B$,此外还有一个干扰变量 $C$,它有两个取值 $C1,C2$. 设计时应该考虑到 $C$ 的分配,我们给出三种方案:

Runs A B AB 方案1 方案2 方案3
$(1)$ $-$ $-$ $+$ $C1$ $C1$ $C1$
$a$ $+$ $-$ $-$ $C1$ $C2$ $C1$
$b$ $-$ $+$ $-$ $C2$ $C1$ $C2$
$ab$ $+$ $+$ $+$ $C2$ $C2$ $C1$

直观上看,方案 1 中 $C$ 的分配和 $B$ 的分配一样,这会导致因子 $B$ 出现 confounding;方案 2 会导致因子 $A$ confounding;而方案 3 则没有 confounding。

你能提出一个使得 $AB$ confounding 而 $A$ 和 $B$ 都没有 confounding 的方案吗?

数学解释

模型: $$ y=\mu+\frac{A}{2}x_A+\frac{B}{2}x_B+\frac{AB}{2}x_Ax_B. $$ 四个未知数刚好解出来。而当我们引入干扰因子 $C$ 时,便乘了 $$ y=\mu_1I(C1)+\mu_2I(C2)+\frac{A}{2}x_A+\frac{B}{2}x_B+\frac{AB}{2}x_Ax_B. $$ 五个未知数,而我们只有四次实验,所以无法全部解出。以方案 1 为例,$C$ 的行为和 $B$ 一样,我们可以设 $$ \begin{gather} \tau_1=\frac{B}{2}-\mu_1,\ \tau_2=\frac{B}{2}+\mu_2. \end{gather} $$ 此时有 $$ y=\frac{\tau_2-\tau_1}{2}+\frac{A}{2}x_A+(\tau_1+\tau_2)x_B+\tau_1\frac{AB}{2}x_Ax_B. $$ 可见我们解不出 $B,\mu_1,\mu_2$,但其它的都能解。所以 $B$ 被“混淆”了。

以上讨论给我们的启发:混淆是无法避免的。有干扰因子时,未知数个数多于方程个数,肯定无法全部解出,所以实验设计时应当直接牺牲一个 effect,把它和干扰因子进行混淆,让其它的顺利解出。如果一个都不牺牲,结局就是一个也解不出。

一般理论:$2^n$ Blocks ($n\ge2$)

当混淆因子有 4 个 level 时,我们就得牺牲两列,如下:

Effect 1 ($D1$) Effect 2 ($D2$) Blocks
$-$ $-$ 1
$+$ $-$ 2
$-$ $+$ 3
$+$ $+$ 4

正好对应。类似地,干扰因子有 $2^n$ 个 level,就得牺牲 $n$ 列。

但事实没有那么美好,对于上例,如果 $D_1=AB$, $D_2=A$,我们可以证明 $B$ 也惨遭混淆。如果 $D_1=A,D_2=B$ 呢?那 $AB$ 会被混淆。更一般地,我们如果选取 $G$ 中的非单位元 ${h_1,\cdots,h_n}$ 作为牺牲契约,可以证明 ${h_1,\cdots,h_n}$ 生成的子群 $H$ 中的 $2^n$ 个元素都会被混淆,而不在 $H$ 中的元素都不会被混淆。

具体证明比较繁琐,和群没啥关系,需要研究每一列生成的规律(二进制按位与),就略去了。有人有好的证法可以评论一手。

分式析因设计

分式析因设计(Fractional Factorial Experiments, F.F.)是 $2^k$ 因子设计的变种。对 block 来说,“混淆”相当于集合的减法,而分式析因设计则相当于除法,也就是商群。

先来说说动机。根据层次原则,高阶的交叉项显著的概率很低。比如一个 10 因子设计,需要做 1024 次试验,而如果我们假定大于等于 4 阶的 effect 都不显著,我们可以忽略掉整整 968 个。这不禁让我们动起歪脑筋:既然忽略了 968 个 effect,有没有办法少做 968 次试验呢?这不是活活省钱吗?

做法

我们可以从 $G$ 中选取若干非单位元 $h_1,\cdots, h_p$,记它们生成的子群为 $H$,显然 $|H|=2^p$. 我们认为 $H$ 中的元素所代表的 effect 可以忽略,此时我们只需做 $2^{k-p}$ 次实验。具体来说,只要考虑 $h_1,\cdots, h_p$ 中全为正的行即可。这样做相当于考虑商群 $G/H$. 对于 $gH\in G/H$,我们认为 $gH$ 中的元素所代表的 effect 都相等。$H$ 称为 defining contrast subgroup.

为了方便,人们规定了 $H$ 的几种表达方式。

$H$ 有 $p$ 个生成元,如果依次记作 $h_1,\cdots,h_p$,则可以写 $I=h_1,I=h_2,\cdots,I=h_p$. 这组 $p$ 个式子称为 defining relation. ($\mu$ 也写作 $I$)

注意到 $A_i$ 是 $G$ 的生成元,我们总可以将其划分为不交的两部分,$\mathcal{L}={A_{l_1},\cdots,A_{l_p}}$,$\mathcal{R}={A_{r_1},\cdots,A_{r_{k-p}}}$,使得对于 $i=1,\cdots,p$ 有 $$ A_{l_i}H=A_{r_{i_1}}\cdots A_{r_{i_{s_i}}}H $$ 这种形式。也就是左边的每个元素都能用右边的表示出来(可以自己试着证明下)。这种形式叫 generators。

F.F. 的评价

如何评价 F.F. 的好坏?分辨率是一个指标。我们商掉了 $H$,自然希望 $H$ 中尽可能是高阶项。$H$ 中元素对应 effect 阶数称为字长(word-length),例如 $ABC$ 字长为 3. 分辨率 定义为 $H$ 中非单位元的最短字长。

Clear effect: 对于字长为 1 或 2 的元素 $g$,称它是 clear 的,如果陪集 $gH$ 中没有其他的字长为 1 或 2 的元素;若 $gH$ 中没有其他的字长为 1,2,3 的元素,则称为 strongly clear.

更进一步,我们可以考虑 $H$ 中元素的字长,记 $n_i$ 为 $H$ 中字长为 $i$ 的元素个数,则 $W=(n_3,n_4,\cdots)$ 称为 wordlength pattern. (通常我们要求 $n_2=0$,否则糟糕透了,两个字长为 1 的 effect 在同一个陪集中)

选择 $W$ 的字典序最小的设计。