UA MATH571B 试验设计III 单因素试验设计
单因素ANOVA模型设定与假设ANOVA F检验单因素ANOVA
方差分析（Analysis of Variance，ANOVA）是两样本均值的检验的扩展，其作用在于同时比较多组样本的均值。最基础最简单的ANOVA是单因素ANOVA，这一讲会介绍这个模型是什么，有什么用，以及使用这个模型需要的假设、如何验证这些假设是否成立、不成立的话要怎么处理，之后所有的试验设计模型都按这个套路来介绍。回到第一篇博客的例子，要研究二氧化碳浓度对天竺葵光合作用强度的影响，设置对照组（Control Group）二氧化碳浓度是410ppm，实验组（Treatment Group）的二氧化碳浓度是400ppm、420ppm。假设对照组叶片累积的淀粉量均值为μ1\mu_1μ1，实验组累积的淀粉量为μ2\mu_2μ2和μ3\mu_3μ3，关于二氧化碳浓度是否影响天竺葵光合作用强度的假设检验可以写成
H0:μ1=μ2=μ3H_0:\mu_1=\mu_2=\mu_3 H0:μ1=μ2=μ3
要两两比较的话需要C32C_3^2C32次比较，设置n−1n-1n−1个实验组就需要做Cn2C_n^2Cn2次比较，因此两两比较的效率非常低，我们需要一种能同时处理多组样本均值的假设检验方法。
模型设定与假设
对于有aaa组样本均值的假设检验问题，可以将两样本均值假设检验的DGP推广为
yij=μi+ϵij,ϵij∼iidN(0,σ2)i=1,2,⋯,a;j=1,2,⋯,ny_{ij} = \mu_i+ \epsilon_{ij},\epsilon_{ij}\sim_{iid}N(0,\sigma^2)\\ i = 1,2,\cdots,a; j=1,2,\cdots,n yij=μi+ϵij,ϵij∼iidN(0,σ2)i=1,2,⋯,a;j=1,2,⋯,n
其中μi\mu_iμi是第iii个level的组内平均或者treatment mean，ϵij\epsilon_{ij}ϵij是试验误差，这个模型叫做均值模型（mean model）。另外一种等价的写法是
yij=μ+τi+ϵij,ϵij∼iidN(0,σ2)i=1,2,⋯,a;j=1,2,⋯,ny_{ij} = \mu + \tau_i + \epsilon_{ij},\epsilon_{ij}\sim_{iid}N(0,\sigma^2)\\ i = 1,2,\cdots,a; j=1,2,\cdots,n yij=μ+τi+ϵij,ϵij∼iidN(0,σ2)i=1,2,⋯,a;j=1,2,⋯,n
其中μ\muμ是总体均值（grand mean），τi\tau_iτi是第iii种factor level的treatment effect，μi=μ+τi\mu_i=\mu+\tau_iμi=μ+τi，这个模型叫效应模型（effect model）。这里想研究的是单个treatment factor的不同level对response的影响，所以这个模型叫单因素ANOVA（one-way ANOVA）。假设只研究给定的这些factor level，那么这个模型叫固定效应模型（fixed effect model）；如果这些factor level视为是factor level总体的一组随机样本，则treatment effect就是随机的，这种模型叫随机效应模型（random effect model）。本讲只研究固定效应模型。ANOVA模型最终想做的检验是
H0:μ1=μ2=⋯=μaH_0:\mu_1=\mu_2=\cdots=\mu_a H0:μ1=μ2=⋯=μa
也可以写成
H0:τ1=τ2=⋯=τaH_0:\tau_1=\tau_2=\cdots=\tau_a H0:τ1=τ2=⋯=τa
先定义几个符号
yi.=∑j=1nyij,yˉi.=yi.ny..=∑i=1ayi.,yˉ..=y..ay_{i.} = \sum_{j=1}^n y_{ij}, \bar{y}_{i.} = \frac{y_{i.}}{n} \\ y_{..} = \sum_{i=1}^a y_{i.}, \bar{y}_{..} = \frac{y_{..}}{a} yi.=j=1∑nyij,yˉi.=nyi.y..=i=1∑ayi.,yˉ..=ay..
要构造这个检验需要先估计这个模型，固定效应模型可以用最小二乘法来估计
L=∑i=1a∑j=1n(yij−μ^−τ^i)2L = \sum_{i=1}^a \sum_{j=1}^n (y_{ij}-\hat{\mu}-\hat{\tau}_i)^2 L=i=1∑aj=1∑n(yij−μ^−τ^i)2
分别对μ^\hat{\mu}μ^与τ^i\hat{\tau}_iτ^i求偏导，并令其为零可得：
−2∑i=1a∑j=1n(yij−μ^−τ^i)=0⇒anμ^+∑i=1aτ^i=y..−2∑j=1n(yij+μ^−τ^i)=0⇒nμ^+nτ^i=yi.,i=1,⋯,a-2\sum_{i=1}^a \sum_{j=1}^n (y_{ij}-\hat{\mu}-\hat{\tau}_i)=0 \Rightarrow an\hat{\mu}+\sum_{i=1}^a \hat{\tau}_i=y_{..}\\ -2 \sum_{j=1}^n (y_{ij}+\hat{\mu}-\hat{\tau}_i)=0 \Rightarrow n\hat{\mu}+ n\hat{\tau}_i=y_{i.},i=1,\cdots,a −2i=1∑aj=1∑n(yij−μ^−τ^i)=0⇒anμ^+i=1∑aτ^i=y..−2j=1∑n(yij+μ^−τ^i)=0⇒nμ^+nτ^i=yi.,i=1,⋯,a
注意到后aaa个方程加起来就是第一个方程，因此这其实是一个超静定的线性系统。这是由于从两样本到多样本的自然推广其实是均值模型，均值模型恰好能让方程数目与未知参数数目一致，而效应模型与均值模型相比多了一个参数，所以方程数就少了一个。对均值模型用最小二乘估计，简单计算一下就知道
nμ^i=yi.⇒nμ^+nτ^i=yi.,i=1,⋯,an \hat{\mu}_i = y_{i.} \Rightarrow n\hat{\mu}+ n\hat{\tau}_i=y_{i.},i=1,\cdots,a nμ^i=yi.⇒nμ^+nτ^i=yi.,i=1,⋯,a
所以明显均值模型正好存在唯一的估计。为了估计效应模型，最常见的做法是增加一个约束
∑i=1aτ^i=0\sum_{i=1}^a \hat{\tau}_i = 0 i=1∑aτ^i=0
增加这个约束后，μ^=y../an\hat{\mu}=y_{..}/anμ^=y../an正好是总体样本均值，而τi\tau_iτi的含义也可以解读为因为第iii种treatment factor level，导致的treatment mean相对grand mean的偏离程度，正好可以用来衡量treatment的作用，所以增加这个约束是合理的。那么在有了这个约束的前提下，效应模型的参数最小二乘估计为
μ^=y..an=yˉ..,τ^i=yi.n−y..an=yˉi.−yˉ..\hat{\mu}=\frac{y_{..}}{an}=\bar{y}_{..} , \hat{\tau}_i = \frac{y_{i.}}{n}-\frac{y_{..}}{an}=\bar{y}_{i.} -\bar{y}_{..} μ^=any..=yˉ..,τ^i=nyi.−any..=yˉi.−yˉ..
残差的估计量为
eij=yij−yˉi.e_{ij}=y_{ij}-\bar{y}_{i.} eij=yij−yˉi.
在最小二乘估计的基础上，可以用广义线性检验方法完成ANOVA检验，具体参考回归那个系列。
ANOVA F检验
Response数据中的信息含量可以用总平方和表示
SST=∑i=1a∑j=1n(yij−yˉ..)2SST = \sum_{i=1}^a \sum_{j=1}^n (y_{ij}-\bar{y}_{..})^2 SST=i=1∑aj=1∑n(yij−yˉ..)2
现在考虑对总平方和做分解
SST=∑i=1a∑j=1n(yij−yˉi.+yˉi.−yˉ..)2=∑i=1a∑j=1n(yij−yˉi.)2+∑i=1a∑j=1n(yˉi.−yˉ..)2+∑i=1a∑j=1n(yˉi.−yˉ..)(yij−yˉi.)=∑i=1a∑j=1neij2+n∑i=1aτ^i2+0SST = \sum_{i=1}^a \sum_{j=1}^n (y_{ij}-\bar{y}_{i.}+\bar{y}_{i.}-\bar{y}_{..})^2 \\ = \sum_{i=1}^a \sum_{j=1}^n (y_{ij}-\bar{y}_{i.})^2 + \sum_{i=1}^a \sum_{j=1}^n (\bar{y}_{i.}-\bar{y}_{..})^2 + \sum_{i=1}^a \sum_{j=1}^n (\bar{y}_{i.}-\bar{y}_{..})(y_{ij}-\bar{y}_{i.}) \\ = \sum_{i=1}^a \sum_{j=1}^n e_{ij}^2 + n\sum_{i=1}^a \hat{\tau}_{i}^2 + 0 SST=i=1∑aj=1∑n(yij−yˉi.+yˉi.−yˉ..)2=i=1∑aj=1∑n(yij−yˉi.)2+i=1∑aj=1∑n(yˉi.−yˉ..)2+i=1∑aj=1∑n(yˉi.−yˉ..)(yij−yˉi.)=i=1∑aj=1∑neij2+ni=1∑aτ^i2+0
其中第一项是残差平方和（SSE），第二项是试验平方和（SS of Treatment）记为SSM，第三项为零，因为
∑i=1a∑j=1n(yˉi.−yˉ..)(yij−yˉi.)=∑i=1a(yˉi.−yˉ..)(∑j=1n(yij−yˉi.))∑j=1n(yij−yˉi.)=nyˉi.−nyˉi.=0\sum_{i=1}^a \sum_{j=1}^n (\bar{y}_{i.}-\bar{y}_{..})(y_{ij}-\bar{y}_{i.}) = \sum_{i=1}^a (\bar{y}_{i.}-\bar{y}_{..}) (\sum_{j=1}^n (y_{ij}-\bar{y}_{i.}) ) \\ \sum_{j=1}^n (y_{ij}-\bar{y}_{i.}) = n\bar{y}_{i.}- n\bar{y}_{i.}=0 i=1∑aj=1∑n(yˉi.−yˉ..)(yij−yˉi.)=i=1∑a(yˉi.−yˉ..)(j=1∑n(yij−yˉi.))j=1∑n(yij−yˉi.)=nyˉi.−nyˉi.=0
因此总平方和可以分解为
SST=SSM+SSESST = SSM + SSE SST=SSM+SSE
现在研究这个分解的分布特征。定义N=anN=anN=an为总样本数。首先考虑，每一组的样本方差为
Si2=∑j=1n(yij−yˉi.)2n−1S_i^2 = \frac{\sum_{j=1}^n (y_{ij}-\bar{y}_{i.})^2}{n-1} Si2=n−1∑j=1n(yij−yˉi.)2
考虑每一组样本方差的平均
(n−1)S12+⋯+(n−1)Sa2(n−1)+⋯+(n−1)=∑i=1a∑j=1n(yij−yˉi.)2N−a=SSEN−a\frac{(n-1)S_1^2 + \cdots + (n-1)S_a^2}{(n-1)+\cdots+(n-1)}=\frac{\sum_{i=1}^a \sum_{j=1}^n (y_{ij}-\bar{y}_{i.})^2}{N-a}=\frac{SSE}{N-a} (n−1)+⋯+(n−1)(n−1)S12+⋯+(n−1)Sa2=N−a∑i=1a∑j=1n(yij−yˉi.)2=N−aSSE
因为残差独立同方差，因此上式是方差的无偏估计，也就是说
E[SSEN−a]=σ2E \left[ \frac{SSE}{N-a}\right] = \sigma^2 E[N−aSSE]=σ2
这个证明比较直接，就是把SSESSESSE仔细展开求期望就好。定义MSE=SSE/dfEMSE=SSE/df_EMSE=SSE/dfE，dfEdf_{E}dfE是其对应的自由度。类似地可以定义MSM=SSM/dfMMSM=SSM/df_MMSM=SSM/dfM，dfM=a−1df_{M}=a-1dfM=a−1是其对应的自由度，
E[MSM]=σ2+n∑i=1aτi2a−1E[MSM]=\sigma^2 + \frac{n\sum_{i=1}^a \tau_i^2}{a-1} E[MSM]=σ2+a−1n∑i=1aτi2
当ANOVA检验的原假设成立时，显然MSMMSMMSM也是方差的无偏估计。根据Cochran定理（参考概率论那个系列的文章），可以构造统计量
F0=SSM/a−1SSE/N−a=MSMSSE∼F(a−1,N−a)F_0 = \frac{SSM/a-1}{SSE/N-a} = \frac{MSM}{SSE} \sim F(a-1,N-a) F0=SSE/N−aSSM/a−1=SSEMSM∼F(a−1,N−a)
由此可以对原假设做ANOVA F检验。ANOVA Table可以表示为
如果觉得《UA MATH571B 试验设计III 单因素试验设计1》对你有帮助，请点赞、收藏，并留下你的观点哦！
