前言
方差分析方法可以用于比较两组或多组均数，其处理因素一般是可以控制的。方差分析要求各比较组除了所施加的处理因素不同外，其它对观察指标有影响的因素应齐同或均衡。
例如，我想要研究某个药物的不同剂量对小鼠的体重影响时，我们要保证每组的小鼠体重在给药前都相同或接近。
但是在实际工作中，有时候有些因素根本无法控制，或者由于实验设计的疏忽、实验条件的限制等原因造成对观察指标有影响的个别因素未加控制或难以它抑制。
例如，在降压药物疗效考核的临床试验中，患者的初始血压水平对服药一段时间后的血压下降量有影响，但患者的初始血压水平是难以控制的。如果不考虑各组患者的初始血压水平的差异，直接用方差分析的方法比较不同处理组患者的平均血压下降量，以评价药物的降压效果是不恰当的。如何在比较两组或多组均数的同时扣除或均衡这些不可控制因素的影响，可以考虑协方差分析的方法。
协方差矩阵
举个例子，如下所示：
|
|
上图是x与y对应的散点图，我们把这个图形分为4个象限，即1、2、3、4，现在我们求一组值，分别为：
$y_{i}-\bar{y}$：每一个y值与y均值的偏离；
$x_{i}-\bar{x}$：每一个x值与x均值的偏离；
$(y_{i}-\bar{y})\times (x_{i}-\bar{x})$。
这3组值分别如下所示：
|
|
计算结果如下所示：
|
|
这里面有正有负，可以得出以下结果：
1、2象限中，$y_{i}-\bar{y}$为正，3、4象限中，$y_{i}-\bar{y}$为负，
1、4象限中，$x_{i}-\bar{x}$为正，2、3象限中，$x_{i}-\bar{x}$为负。
画成表格就是下面的样子：
| 象限 | $y_{i}-\bar{y}$ | $x_{i}-\bar{x}$ | $(y_{i}-\bar{y})\times (x_{i}-\bar{x})$ |
|---|---|---|---|
| 1 | + | + | + |
| 2 | + | - | - |
| 3 | - | - | + |
| 4 | - | + | - |
通过上面的结果可以得出以下结论：
如果x和y变动一致，此时
根据这3组值，可以得到以下结论：
其中红色方框的部分是每列的方差，蓝色表示的是两两组合后（也就是xy，xz，yz）的协方差，其中我们可以发现，左下角用蓝色圈起来的部分与右上的未圈部分是一样的。用红色圈起来的部分好理解，在R中直接使用var(x)或sd(x)^2就可以直接计算，现在讲一下协方差如何计算，它的公式如下所示：
$cov(X,Y)=\frac{\sum_{i=1}^n (X_{i}-\bar{X})(Y_{i}-\bar{Y})}{n-1}$
这个公式如果不好理解的话，可以看一个案例，如下所示：
|
|
以上面的这个矩阵为例说明一下cov(X,Y)是如何计算的，计算公式表示，x的值减去x的均值得到一个向量，y的值减去y的均值，得到一个向量，这两个向量相乘，并除以n-1（这里就是9了），就得到了cov(X,Y)，现在计算一下：
|
|
那么此时，cov_data协方差矩阵的第1列，第2行或第2列第1行的值就是-96.55556，也就是下面的黄色部分，如下所示：
接着再计算cov(z,x)和cov(z,y)，如下所示：
|
|
把它们填到这个矩阵中，就是下面的这个样子：
剩下的绿色部分就是每列的方差的，很好计算，如下所示：
|
|
再补充好矩阵，如下所示：
上面只是过计算过程，在R中直接使用cov()函数即可，如下所示：
|
|
在R中协方差矩阵的计算原理就搞清楚了。
协方差分析基本思想
协方差分析(analysis of covariance)是将线性回归分析与方差分析结合起来的一种统计方法。在协方差分析中，影响观察指标Y（应变量）的往往是一些定性变量；而在线性回归分析中，影响Y的都是定量变量。协方差分析的基本思想就是将那些定量产变量X（指未加或难以控制的因素）对Y的影响看作协变量(covariate)，建立应变量Y随协变量X变化的线性回归关系，并利用这种回归关系把X值化为相等后再进行各级Y的修正均数(adjusted means)的比较，其实质是从Y的总平方和中扣除协变量X对Y的回归平方和，对残差平方和作进一步分解后再进行方差分析，以更好地评价处理因素的效应。
按方差分析不同实验设计类型，相应有不同的协方差分析，如有完全随机设计、随机区组（配伍组）设计、拉丁方设计和析因设计等类型的协方差分析，而且协变量可以有一个、两个或多个，分析方法略有不同，但其解决问题的基本思想相同。
应用条件
协方差分析有两个重要的应用条件：一是与方差分析的应用条件相同，即理论上要求观察变量服从正态分布，各观察变量相互独立，各样本的总体方差齐性；二是各总体客观存在应变量对协变量的线性回归关系且斜率相同（回归线平行），即要求各样本回归系数b本身有统计学意义而各样本回归系数b间的差别无统计学意义。因此，进行协方差分析时，必须先对样本资料进行方差齐性检验和回归系数的假设检验，若满足这两个条件或经变量变换后满足这两个条件，才可作协方差分析。
协方差分析要求协变量是连续变量，而且不能是影响处理的变量，其取值应在研究之前被观察；或者虽在研究中被观察，但不能受到处理的影响。
协方差分析中比较的修正均数，而修正均数实际上是假设协变量取值固定在其总均数时的观察变量的均数。修正均数间的差别与实际均数间的差别并不是一回事，只是用作合理比较。当各比较组的协变量相关悬殊时，协变量的总均数可能不落在各处理组协变量的实测范围内，这时的修正均数可能是对回归线的一种外推，但是这种外推是否仍满足线性和平等的条件无人知晓。因此，协方差分林的时最好先对协变量间的差别作假设检验，协方差分析应用于协变量均数间差别不大的资料较好，利于对分析结果作出较恰当的解释。
如果在多因素研究中有一个或多个因素（协变量）不能或难以控制，而这个或这些协变量对观察变量可能有影响，解决这一类问题可以用多元协方差分析方法(analysis of multiple covariance)。
分析步骤
现在我们以完全随机资料为例，设有G组回归样本，每组有n对观察值，总观察值对数为N=nG（一般情况为各组有$n_{j}$对观察值，总观察值对数为$N=n_{1}+n_{2}+\cdots+n_{G})$，其基础数据与协方差分析计算模式如表13-1和表13-2。
计算各组$\sum X_{j},\sum Y_{j},\sum X_{j}^2,\sum Y_{j}^2,\sum X_{j},\sum Y_{j},\bar{X_{j}},\bar{Y_{j}}$
利用合计项各数据计算较正数$C_{1},C_{2},C_{3}$，以及总变异的离均差平方和$l_{XX},l_{YY}$积和$l_{XY}$及自由度$\nu$，如下所示：
$C_{1}=\frac{(\sum X)^2}{N},C_{2}=\frac{\sum Y)^2}{N},C_{3}=\frac{\sum X \cdot \sum Y}{N}$
$l_{XX}=\sum X^2=C_{1},l_{YY}=\sum Y^2=C_{2}$
$l_{XY}=\sum XY=C_{3}$
$\nu=N-1$
计算各处理组间的离均差平方和$l_{XX},l_{YY}$，积和$l_{XY}$及自由度$\nu$，如下所示：
$l_{XX}=\sum \frac{(\sum X_{j})^2}{n_{j}}-C_{1}$
$l_{XX}=\sum \frac{(\sum Y_{j})^2}{n_{j}}-C_{2}$
$l_{XX}=\sum \frac{(\sum X_{j})\sum y_{j})}{n_{j}}-C_{3}$
$\nu=G-1$
列出协方差分析计算表13-2，将上述第2步，第3步的计算结果列入表的左侧部分，再由总变量的$l_{XX},l_{XY}$及$\nu$减去与处理组（组间）相应的各值，得到组内的$l_{XX},l_{XY}$及$\nu$。
计算回归估计估计误差平方和$\sum(Y-\bar{Y})^2$及徙度，如表13-2右侧部分，其中总的平方和、组内平方和分别由下面的公式计算：
$\sum(Y-\bar{Y})^2=l_{YY}-\frac{l_{XY}^2}{l_{XX}}$
总的平方和减去组内平方和即为“修正零数”的平方和，总的及组内的自由度分别为表13-2左侧相应部分的自由度减去1，再由总的自由度减去组内的自由度得修正均数的自由度。
以修正均数和组内的估计误差平方和分别除以相应的自由度得到修正均数和组内的估计误差均方，简称为修正均数均方和组内均方，再按下面公式计算F值：
$F=\frac{修正均数均方}{组内均方}$
以上述公式中的公子和分母两均方的自由度查表F表，得到P值，按所设立的检验水准α作出统计推断，若$P\leq \alpha$，则可认为各总体修正均数问总的来说有差别，但空间哪两组间有差别，哪两组间没有差别，必要的时候还要做进一步的多重比较。
若需要计算各修正均数，可先按下面的公式计算出公共回归系数(common regression coefficient)$b_{c}$，再由下面的公式计算出各修正均数$\bar{y_{j}}^$
$b_{c}=\frac{组内l_{XY}}{组内l_{XX}},\quad \bar{Y_{j}}^=\bar{Y_{j}}-b_{c}(\bar{X_{j}}-\bar{X})$修正均数间的多重比较，使用q检验，公式如下所示：
上式中$\bar{Y_{A}}^,\bar{Y_{B}}^$为多重比较中任两个对比组的修正均数，分母为其差值的标准误，$S_{Y \cdot X}^2$为组内估计误差的均方，$n_{0}$为每组平均例数，a为将修正均数按大小顺序排列后两对比修正均数间所包含的组数，求得q之后，先按协方差分析表“估计误差”栏中组内自由度和a，查表得到p值，然后再按所设立的检验水准a作推断结论。
完全随机设计资料的协方差分析
现在看一个案例A，如下所示：
案例A
为研究某降血糖药物的有效性及其合用盐酸二甲双胍片的有效性，选择收治90名2型糖尿病患者，并采用随机对照试验，分为三个治疗组，第一组为该降糖药组，第二组为盐酸二甲双胍片组，第三组为该降糖药+盐酸二甲双胍片组，每组30名患者，治疗3个月，主要有效性指标为糖化血红蛋白。测得每个患者入组前（X）和3个月后（Y）的糖化血红蛋白含量（％）见表13-3的上部，试分析三种治疗降糖化血红蛋白的效果是否不同（《医学统计学》（第四版，孙振球，例13-1，P204），数据如下所示：
各种计算结果如下所示：
在此例中，给予何种药物治疗是人为可以控制的定性因素，称为定性变量；患者初始（入组前）的糖化血红蛋白含量是难以控制的定量因素，称为协变量X；试验的观察指标是患者治疗3个月后的糖化血红蛋白含量，称为应变量Y，如果不考虑3组患者初始糖化血红蛋白含量不同的影响，那么此例就是一个典型的完全随机设计类的方差分析问题。
现在我们使用常规的方差分析来计算一下数据（这个案例不能用方差分析，这里使用方差分析只是为了说明一下为什么不能用方差分析），计算过程如下所示：
|
结果如下所示：
|
|
由上面的结果我们可以认为这3种治疗组降糖化血红蛋白的效果不同，但是这样得出的结论是不恰当的。因为从前面的原始数据我们可以得知，3种治疗组患者的初始糖化血红蛋白含量均数不同，第一组最高（10），第三组最低（9.947）。
仅对给药后的糖化血红蛋白进行方差分析，再经SNK-q的两两比较我们可以知道，除了第一组与第二组间的差异无统计学意义外（P>0.05），第一组和第二组分别与第三组间的差异均有统计学意义(P<0.01)，计算过程如下所示：
|
结果如下所示（只截取了各组均分与两两比较的结果）：
|
|
但是，这是在没有扣除X对Y影响的情况下的计算的结果，也就是说第二组降糖均数最低，第三组最高（看Q50），提示患者初始的糖化血红蛋白含量与治疗药物的效应混杂，即患者治疗3个月后的糖化血红蛋白含量的多少除了与治疗的药物有关外，还受其初始糖化血红蛋白含量的影响。对此，可以采用完全随机设计且有一个协变量的协方差分析方法，将三组患者初始的糖化血红蛋白含量化为相等，以扣除其影响，再比较三种治疗组的降糖效果，即检验三组修正均数间的差别有无统计学意义。具体的计算步骤如下：
- $H_{0}:$ 各组降糖的总体修正数相等，$H_{1}$：各组降糖总体修正均数不全相等，$\alpha=0.05$
- 列表，并计算初步结果，也就是下面的这一部分：
- 根据前面的公式，计算相应的校正数，总的，组间和组内的离均差平方和、积和自由度，计算结果如下所示：
其中，校正系数为：
总的离均差平方和、积和及自由度为：
组间离均差平方和、积和及自由度为：
组内离均差平方和、积和及自由度：
- 根据前面的公式计算出总的、组内和修正均数的估计误差平方和、自由度、均方及F值，见前面表13-5右侧部分，如下所示：
$\nu_{1}=2,\nu_{2}=86$查表可得，Ｐ<0.01,按α=0.05水准,拒绝H0，接受H1，可以认为在扣除初始糖化血红蛋白含量因素的影响后，3组患者的总体降糖均数有差别。
根据前面的表13-3，表13-5及公式，求得公共回归系数$b_{c}$及各组修正均数$\bar{Y}^*$为
$b_{c}=13.26/53.40=0.25$
第一组$\bar{Y_{1}}^*=1.61-0.25(10.09-9.970)=1.61$
第一组$\bar{Y_{2}}^*=1.07-0.25(9.88-9.970)=1.09$
第一组$\bar{Y_{3}}^*=2.24-0.25(9.95-9.970)=2.24$
根据前面的q值公式，对这三个修正均数间差别进行两两比较的结果显示，第一组与第二组、第一组与第三组、第二组与第三组总体修正均数间均有差别(P<0.01)，可以认为在扣除患者初始糖化血红蛋白含量因素的影响后，第三种治疗患者的平均降糖量最多，第一种治疗次之，第二种治疗最少。
使用R进行协方差分析时，需要安装HH包，用这个包中的ancova()函数进行协方差分析。
计算过程如下所示：
|
|
结果如下所示：
|
|
从这个协方差的最终结果图上我们可以知道，第三种疗法对血红蛋白的降低程度最大，第一种次之，第二种降低程度最弱。
参考资料
- 孙振球. 医学统计学.第3版[M]. 人民卫生出版社, 2010.
- 浅谈协方差矩阵
- 协方差与协方差矩阵
- 从协方差到相关系数.冯国双.微信公众号文章：小白学统计