【数理统计基础】 06 - 相关分析和方差分析
1. 相关分析
1.1 相关系数
在一堆变量中，找到并分析它们之间的关系，是复杂环境和模型中的重要任务。由于线性关系的特殊、常见和简单，数学上往往采用线性关系来逼近实际关系。上篇的线性回归以及概率论中的线性回归，更关注的是线性函数的参数估计。如果想单纯地度量随机变量的线性关系，直接讨论相关系数即可，请先复习斜方差的相关概念。
两个变量之间的线性关系，就是之前学过的协方差的概念\(\text{Cov}(X,Y)\)。在得到\(n\)个样本\((X_i,Y_i)\)后，容易得到式（1）的无偏估计，注意其中降低了一个自由度，继而还可以有式（2）的样本相关系数。相关系数是线性关系的直接度量，它可以作为相关假设的检验条件，最常用的就是当\(|r|\leqslant C\)时认为\(X,Y\)是不相关的。
\[\dfrac{1}{n-1}\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y})\approx\text{Cov}(X,Y)\tag{1}\]
\[r=\dfrac{1}{S_XS_Y}\sum_{i=1}^n(X_i-\bar{X})(Y_i-\bar{Y}),\;\;S_X^2=\sum_{i=1}^n(X_i-\bar{X})^2\tag{2}\]
为了能找到关于\(r\)的枢轴变量，这里还是要做一些假设，即\((X,Y)\)是一个二元正态分布。回顾二元正态分布的知识（《初等概率论》第5篇公式（27）），可知\(X,Y\)完全符合一元线性回归的模型。为此这里暂且取定\(X_i\)，而把\(Y_i\)看成随机变量，并对它们进行一元回归分析。比较发现系数估计满足\(\alpha_1=r\cdot\dfrac{S_Y}{S_X}\)，在假设\(\rho=0\)（即系数\(a_1=0\)）的情况下，把这个等式代入上篇公式（12）右的枢轴变量，整理后得到式（3）。由于该结论与\(X_i\)的取值无关，因此它对于变量\(X_i\)也成立，它就是我们要找的枢轴变量。
\[\dfrac{r\sqrt{n-2}}{\sqrt{1-r^2}}\sim t_{n-2}\tag{3}\]
1.2 复相关系数
相关系数度量了两个随机变量之间的线性关系，当系统中的变量很多时，关系也会变得复杂，这时需要引入更多的关系分析。以下记要讨论的\(n\)个变量为\(X_i\)，\(X_i,X_j\)的相关系数为\(\rho_{ij}\)，并记矩阵\(P=[\rho_{ij}]\)，而去除\(i\)行\(j\)列后的子矩阵记作\(P_{ij}\)。在得到样本后，同样可以计算样本相关系数\(r_{ij}\)，并记矩阵\(R=[r_{ij}]\)和子矩阵\(R_{ij}\)。
首先比较容易想到的关系，是一个变量\(X_1\)与多个变量\(X_2,\cdots,X_p\)的整体关系。回顾概率论中的线性回归，假设\(X_1\)对\(X_2,\cdots,X_p\)的线性回归是\(L(X_2,\cdots,X_p)\)，则容易证明\(X_1-L\)与\(X_2,\cdots,X_p\)都不相关。仿照线性空间中的最小二乘法，\(L\)可以看成是\(X_1\)在\(X_2,\cdots,X_p\)空间中“投影”，故用\(X_1\)和\(L\)的关系作为\(X_1\)与\(X_2,\cdots,X_p\)的关系是比较合理的，这个关系被称为\(X_1\)与\(X_2,\cdots,X_p\)的复相关系数（式（4）左）。
\[\rho_{1(23\cdots p)}=\dfrac{\text{Cov}(X_1,L)}{\sqrt{D(X_1)D(L)}}=\sqrt{1-|P|\,/\,|P_{11}|}\tag{4}\]
式（4）右的证明比较繁杂，这里先从一些引论开始。考察随机变量\(Y\)和随机向量\(X=[X_1,\cdots,X_n]\)，为简化讨论，设它们已经中心化。设\(Y\)关于\(X\)的回归函数是\(L(X)=\alpha_1X_1+\cdots+\alpha_nX_n\)，则由最小二乘法可以得到式（5）。求解方程组便得到\(\alpha=[\alpha_1,\cdots,\alpha_n]^T\)的解为\(C_x^{-1}C_y\)，其中\(C_x,C_y\)分别为方程组的系数矩阵和常数列向量。
\[\min\{E[Y-\sum_{i=1}^n\alpha_iX_i]^2\}\;\Rightarrow\;\sum_{i=1}^n\text{Cov}(X_i,X_j)\alpha_i=\text{Cov}(Y,X_j)\tag{5}\]
然后可以计算得\(\text{Cov}(Y,L)=D(L)=C_y^TC_x^{-1}C_y\)，这时再计算复相关系数，并把协方差换算成相关系数，可得式（6）左。其中\(P_y\)是\(Y\)与\(X_i\)的相关系数组成的列向量，而\(P_x\)是\(X_i\)之间的相关系数组成的矩阵。设\(P_x\)的伴随矩阵为\(P_x^*\)，而记\(P\)为\((Y,X_1,\cdots,X_n)\)的相关系数矩阵，则不难发现，\(|P|\)按第\(1\)行、第\(1\)列展开后其实是\(|P_x|-P_y^TP_x^*P_y\)。这样就有了式（6）右成立，同样也有式（4）右成立。
\[\rho_{Y(X)}=\sqrt{P_y^TP_x^{-1}P_y}=\sqrt{1-|P|/|P_x|}\tag{6}\]
在得到样本后，利用\(r_{ij}\)来估计\(\rho_{ij}\)，带入式（4）后算得的估计值称为样本复相关系数\(r_{1(23\cdots p)}\)。当\((X_1,\cdots,X_p)\)是\(p\)维正态分布时，为检验假设\(\rho_{1(23\cdots p)}=0\)，可以证明有式（7）的枢轴变量。
\[\dfrac{n-p}{p-1}\cdot\dfrac{r^2}{1-r^2}\sim F_{(p-1)/2,(n-p)/2}\tag{7}\]
1.3 偏相关系数
有时候两个变量\(X_1,X_2\)的相关性并不是因为它们有直接联系，而是因为它们共同与\(X_3,\cdots,X_p\)相关。所以有必要将\(X_3,\cdots,X_p\)的相关性从\(X_1,X_2\)中去除后再计算\(X_1,X_2\)的相关性，步骤也是比较自然的，先计算出\(X_1,X_2\)对\(X_3,\cdots,X_p\)的线性回归\(L_i(X_3,\cdots,X_p)\)，然后计算\(X'_1=X_1-L_1,X'_2=X_2-L_2\)的相关系数。这样的关系被称为\(X_1,X_2\)对\(X_3,\cdots,X_p\)偏相关系数（式（8）左）。
\[\rho_{12\cdot(3\cdots p)}=\dfrac{\text{Cov}(X'_1,X'_2)}{\sqrt{D(X'_1)D(X'_2)}}=\dfrac{|P_{12}|}{\sqrt{|P_{11}|\cdot|P_{22}|}}\tag{8}\]
上面引理证明过程中的结论，同样可以证明式（8）右，请自行补齐证明过程。另外同样地，可以利用\(r_{ij}\)估计式（8）得到样本偏相关系数\(\rho_{12\cdot(3\cdots p)}\)。当\((X_1,\cdots,X_p)\)是\(p\)维正态分布时，为检验假设\(r_{12\cdot(3\cdots p)}=0\)，可以证明有式（9）的枢轴变量。
\[\dfrac{r\sqrt{n-p}}{\sqrt{1-r^2}}\sim t_{n-p}\tag{9}\]
2. 方差分析
2.1 单因素完全实验
前面的讨论都集中在线性关系上，更一般地还需要讨论一般的关系模型\(Y=f(X)+e\)。确定具体的\(f(x)\)是一个很开放的问题，前面的线性模型算一种，数学中还有很多逼近理论也可以派上用场。这里不深入讨论\(f(x)\)本身，而是只解决最简单的假设检验问题，即\(X\)对\(Y\)是否有显著影响。
以下假设\(X\)有\(k\)个采样值\(X_i\)，任务是检验\(Y_i\)是否受\(X_i\)影响较大。由于\(Y\)还受到随机因素\(e\)的影响，在同一个\(X_i\)下一定要有多个\(Y\)的采样值，才能对\(Y_i\)有个较好的估计。设\(Y_i\)有\(n_i\)个采样值\(Y_{ij}\)，并记\(n=n_1+\cdots+n_k\)，模型可以写成式（10）。把模型中心化会更便于处理，故令\(f(X_i)=\mu+a_i\)，其中\(a_1+\cdots+a_k=0\)。
\[Y_{ij}=f(X_i)+e_{ij}=\mu+a_i+e_{ij},\;\;(e_{ij}\sim e)\tag{10}\]
你可能注意到，\(X_i\)的具体值在这里并不重要，不同的\(X_i\)只是对\(Y_{ij}\)的一个分组，要检验的假设其实是分布并不受分组影响。以下记\(Y_{ij}\)的平均值是\(\bar{Y}\)，而记\(Y_{i1},\cdots,Y_{in_i}\)的平均值是\(\bar{Y}_i\)。想要搞清楚\(Y_{ij}\)是否受分组影响，首先当然要看\(\bar{Y}_i\)的分散程度。然后因为随机值\(e_{ij}\)会影响\(\bar{Y}_i\)的精确性，评估时还要对比\(e_{ij}\)的分散程度。
具体来说，分散程度一般用平方和来度量，这样的统计量一般称为离差平方和。最简单的就是所有样本\(Y_{ij}\)的总离差平方和\(Q_T\)（式（11）左），其次是每个\(f(X_i)\)的组内离差平方和\(Q_E\)（式（11）右）。直观上可以认为总离差平方和\(Q_T\)分为两个部分，一部分是\(f(X_i)\)的组间离差平方和\(Q_X\)，另一部分就是组内离差平方和\(Q_E\)。因此把\(Q_X\)定义为式（12）也是合理的，计算整理后得到的表达式更是有直观的意义。
\[Q_T=\sum_{i=1}^k\sum_{j=1}^n(Y_{ij}-\bar{Y})^2;\;\;Q_E=\sum_{i=1}^k\sum_{j=1}^n(Y_{ij}-\bar{Y}_i)^2\tag{11}\]
\[Q_X=Q_T-Q_E=\sum_{i=1}^kn_i(\bar{Y}_i-\bar{Y})^2\tag{12}\]
然后很容易算到它们的期望值式（13），从中不难发现，\(E[Q_X]\)仍然会含有误差方差的信息，因此必须结合误差信息来度量\(X\)的影响。为度量影响大小，将假设定为\(a_1=\cdots=a_k=0\)，假设成立时称\(X\)对\(Y\)影响显著，否则是影响不显著。当假设成立时，三个离差平方和中都只剩下\(\sigma^2\)项，预感枢轴变量是它们之间相除得到的\(F\)统计量。
\[E[Q_T]=(n-1)\sigma^2+\sum_{i=1}^kn_ka_i^2;\;\;E[Q_E]=(n-r)\sigma^2\tag{13}\]
为寻找枢轴变量，首先假定\(e\)是正态分布，然后将式（10）右带入式（11）（12），由于\(a_i=0\)，得到的结果其实就是把\(Y\)换成\(e\)。考察这些关于\(e_{ij}\)的正定二次型，不难得到\(Q_T,Q_X,Q_E\)的秩分别为\(n-1,k-1,n-k\)，由柯赫伦分解定理可知，\(\dfrac{Q_X}{\sigma^2},\dfrac{Q_E}{\sigma^2}\)分别是自由度为\(n-k,k-1\)的卡方分布，且它们互相独立。
它们正好可以用来生成\(F\)型枢轴变量（式（14）），另外由于假设不成立时，有\(\dfrac{E[Q_X]}{k-1}>\dfrac{E[Q_E]}{n-k}\)，故检验条件选择\(F<C\)。需要强调，检验的结果只是\(X\)相对随机值\(e\)影响\(Y\)大小的一个度量，如果直观上看\(\bar{Y}_i\)的差别十分明显，则说明误差的影响特别大，需要增加实验次数或先提取主要因素。如果假设不成立，还可以继续对\(a_i-a_j\)做区间估计，请自行讨论其枢轴变量。
\[F=\dfrac{(n-k)Q_X}{(k-1)Q_E}\sim F_{k-1,n-k}\tag{14}\]
2.2 两因素完全实验
当\(Y\)有多个影响因素，并且各因素互相独立时，如果针对每个因素进行方差分析，往往需要较多的样本数。这时可以将多个因素合并进一个模型，以两个因素\(A,B\)为例，建立式（15）左的模型。假设\(A\)有\(m\)个采样点\(A_i\)，\(B\)有\(n\)个采样点\(B_j\)，则总共只需要做\(mn\)次试验（式（15）右）。以下记\(\bar{Y}\)为所有\(Y_{ij}\)的平均值，\(\bar{Y}_{*j}\)为\(Y_{1j},\cdots,Y_{mj}\)的平均值，\(\bar{Y}_{i*}\)为\(Y_{i1},\cdots,Y_{in}\)的平均值。
\[Y=A+B+e;\;\;Y_{ij}=a_i+b_j+e_{ij}\tag{15}\]
很快你会发现，想要对\(a_i,b_j\)进行估值，信息量是不够的。上面的几个平均值的期望值如式（16），其中并不能得到具体的\(a_i,b_j\)。但方差分析其实只关注数据的分散性，因此只要有\(a_i,b_j\)的相对关系即可。为此，记\(\mu=\bar{a}+\bar{b}\)，然后把\(a_i,b_j\)中心化，这样就有了式（17）中更有用的结论。
\[E[\bar{Y}]=\bar{a}+\bar{b};\;E[\bar{Y}_{*j}]=\bar{a}+b_j;\;E[\bar{Y}_{i*}]=a_i+\bar{b}\tag{16}\]
\[E[\bar{Y}]=\mu;\;E[\bar{Y}_{*j}]=\mu+\beta_j;\;E[\bar{Y}_{i*}]=\alpha_i+\mu\tag{17}\]
\(\alpha_i,\beta_j\)的方差和与\(a_i,b_j\)的方差和是一样的，类似式（12）可以到式（18）中方差和的估计。然后按照相同的理念，把式（19）作为误差方差和的估计（不用追究其直观意义）。容易知道\(Q_T,Q_A,Q_B\)的自由度分别是\(mn-1,m-1,n-1\)，则\(Q_E\)的自由度是\((m-1)(n-1)\)。接下来可以得到两个类似式（14）的枢轴变量。
\[Q_A=n\sum_{i=1}^m(\bar{Y}_{i*}-\bar{Y})^2;\;\;Q_B=m\sum_{j=1}^n(\bar{Y}_{*j}-\bar{Y})^2\tag{18}\]
\[Q_E=Q_T-Q_A-Q_B=\sum_{i=1}^m\sum_{j=1}^n(Y_{ij}-\bar{Y}_{i*}-\bar{Y}_{*j}+\bar{Y})^2\tag{19}\]
二元方差分析的模型其实可以直接用在单元素的区组设计上，即假定检验的目标是\(A\)，在每个情况\(A_i\)下进行\(n\)次试验。这\(mn\)次试验原本可以随机安排，但如果\(mn\)个试验环境存在可知的差异，在设计试验时就要使得每种情况\(A_i\)尽量出现在不同的环境中。以最理想的场景为例，试验环境正好可以分为\(n\)种，而每种内部的\(m\)个小环境是相同的，这时环境因素就可以看做是因素\(B\)。
区组设计的目的是为了排除随机环境对试验的影响，当环境差距明显时，直接用两因素模型可以得到更准确的检验。但要注意，如果环境差异并不明显，组内离差平方和会被低估，再加上自由度的损失，平均离差平方和更是被严重低估。因此如果检测出环境影响甚微，应当直接采用单因素的方差分析。
【全篇完】
