首页 > 代码库 > 【数理统计基础】 06 - 相关分析和方差分析

【数理统计基础】 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\)。

  区组设计的目的是为了排除随机环境对试验的影响,当环境差距明显时,直接用两因素模型可以得到更准确的检验。但要注意,如果环境差异并不明显,组内离差平方和会被低估,再加上自由度的损失,平均离差平方和更是被严重低估。因此如果检测出环境影响甚微,应当直接采用单因素的方差分析。


全篇完

【数理统计基础】 06 - 相关分析和方差分析