多元统计之主成分分析

发表于： 更新于： 分类于：
字数： 897 阅读：≈ 5分钟 浏览： 评论：

概述

  • 主成分分析就是设法将原来指标重新组合成一组新的互相无关的几个综合指标来代替原来指标, 同时根据实际需要从中取几个较少的综合指标尽可能多地反映原来指标的信息。
  • 这种将多个指标转化为少数互相无关的综合指标的统计方法叫做主成分分析或称主分量分析。

几何意义

  • 代数观点:
    • p个原始变量的一些特殊的线性组合
  • 几何意义:
    • 这些线性组合通过把由$X_1,X_2,\cdots,X_p$构成的坐标系旋转而产生的新坐标系。这样的新坐标轴使其通过样本变差最大的方向(或者说具有最大的样本方差)。

数学推导

  • 设$X = (X_1,\cdots,X_p)’$为一个p维随机向量, 并假定存在二阶矩, 其均值向量与协差阵分别记为: $$\mu = E(X), \quad \Sigma = D(X)$$ 考虑如下线性变换
    $$
\begin{cases}
    Y_1 = t_{11}X_1 + t_{12}X_2 + \cdots + t_{1p}X_p = T'_1X  \\
    Y_2 = t_{21}X_1 + t_{22}X_2 + \cdots + t_{2p}X_p = T'_2X  \\
    \cdots \cdots \\
    Y_p = t_{p1}X_1 + t_{p2}X_2 + \cdots + t_{pp}X_p = T'_pX
\end{cases}
$$
    用矩阵表示为$Y = (Y_1,Y_2,\cdots,Y_p)’$, 其中$Y = T’X, T=(T_1,T_2,\cdots,T_p)$ 希望找到一组新的变量$Y_1,\cdots,Y_m(m \le p)$, 这组新的变量要求充分地反映原变量$X_1,\cdots,X_p$的信息, 而且相互独立。
    对于$Y_1,\cdots,Y_m$有
    $$
D(Y_i) = D(T'_iX) = T'_iD(X)T''_i = T'_i\Sigma T_i \quad i=1,2,\cdots,m \\
Cov(Y_i,Y_k) = Cov(T'_iX,T'_kX) = T'_iCov(X,X)T''_k = T'_i\Sigma T_k \quad i,k=1,2,\cdots,m
$$
    则问题转化为在新的变量$Y_1,\cdots,Y_m$相互独立的条件下, 求得$T_i$使得$D(Y_i)$达到最大。
    由于$D(Y_i) = T’_i\Sigma T_i$是向量T的增函数, 添加约束条件$T’_iT_i=1$(即将线性组合$Y=T’X$的系数标准化即单位化), 问题变为在约束条件下求最大的问题。
  • 第一主成分满足$T’_1T_1=1$, 使得$D(Y_1) = T’_1\Sigma T_1$达到最大的$Y_1 = T’_1X$。
  • 第二主成分满足$T’_2T_2=1, Cov(Y_2,Y_1) = Cov(T’_2X, T’_1X) = 0$使得$D(Y_2) = T’_2\Sigma T_2$达到最大的$Y_2=T’_2X$。
  • 第k主成分满足$T’_kT_k=1, Cov(Y_k,Y_i) = Cov(T’_kX, T’_iX) = 0 \quad (i < k)$使得$D(Y_k) = T’_k\Sigma T_k$达到最大的$Y_k = T’_kX$。
  • 求第一主成分, 构造目标函数为: $$\quad \varphi_1(T_1, \lambda) = T’_1\Sigma T_1 - \lambda(T’_1T_1 - 1) \tag{1.2-1}$$ 求导有: $$\quad \frac{\partial \varphi_1}{\partial T_1} = 2\Sigma T_1 - 2\lambda T_1 = 0 \tag{1.2-2}$$ 即: $$\quad (\Sigma - \lambda I)T_1 = 0 \tag{1.2-3}$$ 两边左乘$T’_1$得到: $$\quad T’_1\Sigma T_1 = \lambda \tag{1.2-4}$$ 由于X的协方差阵$\Sigma$为非负, 其特征方程(1.2-3)均大于零, 不妨设$\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_p \ge 0$。由(1.2-4)知$Y_1$的方差为$\lambda$。那么, $Y_1$的最大方差值为$\lambda_1$, 其相应的单位化特征向量为$T_1$。
  • 求第二主成分, 由(1.2-2)有$\Sigma T_1 = \lambda T_1,$则$Cov(Y_2,Y_1) = T’_2\Sigma T_1 = \lambda T’_2T_1$。如果$Y_2$与$Y_1$相互独立, 既有$T’_2T_1 = 0$或$T’_1T_2 = 0$。构造第二主成分的目标函数, 即 $$\varphi_2(T_2, \lambda, \rho) = T’_2\Sigma T_2 - \lambda(T’_2T_2 - 1) - 2\rho(T’_1T_2)$$ 对目标函数求导有: $$\frac{\partial\varphi_2}{\partial T_2} = 2\Sigma T_2 - 2\lambda T_2 - 2\rho T_1 = 0$$ 用$T’_1$左乘有$T’_1\Sigma T_2 - \lambda T’_1T_2 - \rho T’_1T_1 = 0$
    由于$T’_1T_2 = 0, T’_1\Sigma T_2 = 0$, 那么$\rho T’_1T_1 = 0$, 即有$\rho = 0$
    从而$(\Sigma - \lambda I)T_2 = 0$, 且$T’_2\Sigma T_2 = \lambda \quad (1.2-5)$
    这样说明, 如果X的协差阵$\Sigma$的特征根为$\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_p \ge 0$, 由(1.2-5)知$Y_2$的最大方差值为第二大特征根$\lambda_2$, 其相应的单位化的特征向量为$T_2$。
  • 针对一般情形, 第k主成分应该是在$T’_kT_k = 1$且$T’_iT_k = 0$或$T’_kT_i = 0(i < k)$的条件下, 使得$D(Y_k) = T’_k\Sigma T_k$达到最大的$Y_k = T’_kX$。构造目标函数为: $$\varphi_k(T_k, \lambda, \rho_i) = T’_k\Sigma T_k - \lambda(T’kT_k - 1) - \sum{i=1}^{k-1}\rho_i(T’iT_k) \tag{1.2-6}$$ 对目标函数求导有: $$\frac{\partial\varphi_k}{\partial T_k} = 2\Sigma T_k - 2\lambda T_k - 2\sum{i=1}^{k-1}\rho_iT_i = 0 \tag{1.2-7}$$ 用$T’_k$左乘(1.2-7)有$T’_i\Sigma T_k - \lambda T’_iT_k - T’i(\sum{i=1}^{k-1}\rho_iT_i) = 0$即$\rho_iT’_iT_i = 0$, 那么$\rho_i = 0 \quad (i=1,2,\cdots,k-1)$。从而$(\Sigma - \lambda I)T_k = 0 \quad (1.2-8)$而且$T’_k\Sigma T_k = \lambda \quad (1.2-8)$
    对于X的协差阵$\Sigma$的特征根为$\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_p \ge 0$, 由(1.2-7)和(1.2-8)知$Y_k$的最大方差值为第k大特征根$\lambda_k$, 其相应的单位化的特征向量为$T_k$。
  • 综上所述, 设$X = (X_1, \cdots, X_p)$的协差矩阵为$\Sigma$, 其特征根为$\lambda_1 \ge \cdots \ge \lambda_k$, 相应的单位化的特征向量为$T_1,T_2,\cdots,T_p$。那么, 由此所确定的主成分为$Y_1 = T’_1X, \cdots, Y_m = T’_mX$, 其方差分别为$\Sigma$的特征根。

性质

  • 设$Y = (Y_1,Y_2,\cdots,Y_p)’$是X的主成分, 由$\Sigma$的所有特征根构成的对角阵为
    $$ \Lambda = 
\begin{bmatrix}
  \lambda_1 &  & 0 \\
            & \ddots  \\
  0         &  & \lambda_p
\end{bmatrix}
$$
    主成分可表示为 $Y = T’X$

性质一

  • 主成分的协方差矩阵是对角阵

性质二

  • 主成分的总方差等于原始变量的总方差

性质三

  • 主成分$Y_k$与原始变量$X_i$的相关系数为: $$ \rho(Y_k,X_i) = \frac{\sqrt{\lambda_k}}{\sqrt{\sigma_{ii}}}t_{ki} \tag{1.3.3-1} $$ 并称之为因子负荷量(或因子载荷量)

性质四

$$ \sum_{i=1}^p\rho^2(Y_k,X_i)\cdot\sigma_{ii} = \lambda_k \quad (k = 1,2,\cdots,p) \tag{1.3.4-1} $$

主成分的偏差贡献率与信息提取率

  • 在选取主成分时, 不仅要考虑累计贡献率, 还应考虑信息提取率。

主成分的偏差贡献率

  • 称 $$\varphi_k = \lambda_k / \sum_{k=1}^p\lambda_k$$ 为第k个主成分$Y_k$的贡献率。第一主成分的贡献率最大, 其他依次减弱。
  • 若只取$m(<p)$个主成分, 则称 $$\psi_m = \sum_{k=1}^m\lambda_k/\sum_{k=1}^p\lambda_k$$ 为主成分$Y_1,\cdots,Y_m$的累计贡献率, 累计贡献率表明$Y_1,\cdots,Y_m$综合$X_1,X_2,\cdots,X_p$的能力。通常取m, 使得累计贡献率达到一个较高的百分数(如85%以上)。

信息提取率

  • 任一原始变量$X_i$可用主成分表示为: $$X_i = t_{1i}Y_1 + t_{2i}Y_2 + \cdots + t_{pi}Y_p$$ 其中各个主成分之间互不相关, 所以原始变量的$X_i$的方差为: $$Var(X_i) = t_{1i}^2\lambda_1 + t_{2i}^2\lambda_2 + \cdots + t_{pi}^2\lambda_p = \sum_{j=1}^pt_{ji}^2\lambda_j = \sigma_{ii}$$
  • 定义比率: $$\Omega_i = \sum_{j=1}^mt_{ji}^2\lambda_j/\sigma_{ii}$$ 为原始变量$X_i$的信息提取率。

实际应用

出发点

  • 为消除量纲的不同可能带来的一些不合理的影响, 常常需要事先对变量标准化。
  • 标准化处理, 即令:
    $$X_i^* = \frac{X_i - E(X_i)}{\sqrt{D(X_i)}} \quad i=1,\cdots,p$$
    $$X^* = (X_1^*, \cdots, X_p^*)'$$的协方差矩阵就是X的相关系数矩阵R.

标准化的实质

  • 标准化后各原始指标的方差均为1, 抹杀了数据本身离散度的信息。对于同度量或相似量级的变量还是使用原始的协方差求解主成分为宜。

注意

  • 若各指标的数量级相差悬殊或者有不同量纲, 较合理做法:
    • 使用原始变量的R代替$\Sigma$做主成分分析;
    • 使用标准化变量的相关矩阵或者协差阵做主成分分析, 二者结果相同。
  • 从原始变量的相关阵求得的主成分与用原始变量的协差阵求得的主成分一般不相同, 有时差异很大。

步骤

  1. 将原始数据标准化
  2. 建立变量的相关系数阵
  3. 求$$R^*$$的特征根为$$\lambda_1^* \ge \cdots \lambda_p^* \ge 0$$, 相应的特征向量为$$T_1^*, T_2^*, \cdots, T_p^*$$
  4. 由累积方差贡献率确定主成分的个数(m), 并写出主成分为: $$Y_i = (T_i^*)'X \quad (i=1,2,\cdots,m)$$

综合评价

  • 一般在进行综合评价时, 需要给各指标赋予权重, 再进行综合评价。
  • 如果是用主成分分析进行综合评价, 则在评价前需要给主成分赋予权重, 此时, 可以以方差贡献率作为参考来确定权重。
  • 设$Y_1,Y_2,\cdots,Y_p$是所求出的p个主成分, 它们的特征分布分别是$\lambda_1,\lambda_2,\cdots,\lambda_p$, 将特征根"归一化", 既有: $$w_i = \lambda_i/\sum_{i=1}^m\lambda_i \quad i=1,2,\cdots,p$$ 记为$W = (w_1,w_2,\cdots,w_p)’$, 由$Y = T’X$, 构造综合评价函数为: $$Z = w_1Y_1 + w_2Y_2 + \cdots + w_pY_p = W’Y = W’T’X = (TW)‘X \tag{3.2-1}$$ 令$TW = w^*_{k \times 1}$, 带入(3.2-1)有
    $$Z = (w^*)'X \tag{3.2-2}$$
    综合评价函数是对原始指标的线性综合, 从计算主成分到对之加权, 经过两次线性运算后得到综合评价函数。

主成分得分

  • 指将第t个样品(标准化后的数据)的值$X_{(t)} = (X_{t1}, X_{t2}, \cdots, X_{tp})’$带入第$i$个主成分$Y_i = T’_iX$, 得到的值称为第$t$个样品在第$i$个主成分的得分。

R应用 

mark =read.table("eg6.1.csv", sep = ",", header = T)
# 计算列与列的相关系数
R = round(cor(mark), 3)
# 使用相关矩阵计算
pca = princomp(mark, cor = T)
summary(pca, loadings = T)
#取前2个主成分，分别为课程差异因子和课程均衡因子
pre = round(predict(pca), 3)
pca$scores
screeplot(pca, types = "lines")
load = loadings(pca)
plot(load[,1:2],xlim=c(-0.6,0.6),ylim=c(-0.6,0.6))
text(load[,1],load[,2],adj=c(0.5,-0.5))
abline(h=0)
abline(v=0)
# 主成分得分
score=as.matrix(pca$scores)
sd=as.vector(pca$sdev)
weight=sd^2/sum(sd^2)
tscore=score%*%weight
outcome=data.frame(pca$scores,tscore)

主成分回归(PCR)

多元线性回归

  • 样本回归模型: $$Y = (I_n \quad X) {\beta_0 \choose \beta} + \mu = \beta_0I_n + X\beta + \mu$$
  • 模型系数向量估计量为: $${\beta_0 \choose \beta} = [(I_nX)’(I_nX)]^{-1}(I_nX)‘Y$$

主成分回归的原理

  • 最小二乘估计的假定前提之一: 自变量不相关, 即自变量之间不存在共线性问题。但在经济管理中, 许多经济变量之间都存在相关关系。
  • 解决办法
    • 先对模型中的自变量进行主成分分析, 并提取一定的主成分;
    • 然后用因变量对所取的主成分进行回归;
    • 最后, 再将因变量对所提取的主成分的回归方程转化为对原自变量的回归方程, 即为主成分回归
  • 对上述回归模型中的自变量$X_1,X_2,\cdots,X_p$的样本观测数据, 若进行主成分分析, 则可得k个主成分分别为:
    $$
  \begin{cases}
    Z_1 = t_{11}X_1 + t_{12}X_2 + \cdots + t_{1p}X_p \\
    Z_2 = t_{21}X_1 + t_{22}X_2 + \cdots + t_{2p}X_p \\
    \cdots  \\
    Z_k = t_{k1}X_1 + t_{k2}X_2 + \cdots + t_{kk}X_k
  \end{cases}
$$
  • 记第i个样品在第j个主成分上的得分为$Z_{ij}$, 并记全部样本的主成分得分矩阵为$Z = (z_{ij})_{n \times k}$, 即有:
    $$Z = 
\begin{bmatrix}
  z_{11} & z_{12} & \cdots & z_{1k} \\
  z_{21} & z_{22} & \cdots & z_{2k} \\
  \vdots & \vdots & \ddots & \vdots \\
  z_{n1} & z_{n2} & \cdots & z_{nk} 
\end{bmatrix} \quad
T = \begin{bmatrix}
  t_{11} & t_{12} & \cdots & t_{1k} \\
  t_{21} & t_{22} & \cdots & t_{2k} \\
  \vdots & \vdots & \ddots & \vdots \\
  t_{n1} & t_{n2} & \cdots & t_{nk} 
\end{bmatrix}
$$
    若主成分分析是根据样本协方差S做出来的, 则主成分的分矩阵为:$Z = X_0T$
    若主成分分析是根据样本相关矩阵R做出来的, 则主成分的分矩阵为:$$Z^* = X^*T^*$$
  • 若将主成分得分矩阵(经过中心化或标准化变换)看作是自变量的样本观测矩阵, 建立因变量Y对k个主成分的线性回归方程:
    $$Y = r_0 + r_1z_1 + \cdot + r_kz_k + \mu \\
或者 \quad Y = r_0^* + r_1^*z_1^* + \cdots + r_k^*z_k^* + \mu$$
    此回归方程的因变量与原回归方程的因变量相同, 所以随机误差也相同, 并且$E(\mu) = 0, E(\mu^2) = \sigma^2$, 即可使用最小二乘法进行系数的估计。 模型系数的估计量为:
    $${\hat{r}_0 \choose \hat{r}} = {n \quad I'_nZ \choose Z'I'_n \quad Z'Z}^{-1}{I_nY \choose Z'Y} \\
或者 {\hat{r}_0^* \choose \hat{r}^*} = {n \quad I'_nZ^* \choose Z^{*'}I'_n \quad Z^{*'}Z}^{-1}{I_nY \choose Z^{*'}Y}
$$
    主成分回归的系数向量估计值与相应的原自变量回归系数向量估计值之间的关系:
    $$\hat{\beta} = T\hat{r} \quad 或者 \quad \hat{\beta}^* = T^*\hat{r}^*$$
    根据上述系数值即可得出因变量与原始变量的回归模型。