10 对应分析

《多元统计分析—基于R》课程的R示例, 主要参考书: 费宇《多元统计分析-基于R》，人民大学出版社，2012。

10.1 对应分析介绍

对应分析是列联表分析方法，考察把行类别和列类别分组降维后，行、列之间的关系。

列联表是有两个分类变量交叉分组后形成的频数（个数）表，比如：

性别	吸烟	不吸烟
男	105	225
女	19	300

这是一个\(2 \times 2\)列联表，内容简单，不需要降维，对应分析不适用。更大的列联表则可以用对应分析发现行类与列类之间的关系。

对应分析把行变量的各个不同值进行降维，把列变量的各个不同值进行降维，与主成分分析和因子分析的思想类似。降维后，行、列两个方向的降维结果具有内在联系，可以用关联的投影方向同时表示行和列的降维结果。

在行、列都降维到二维时，可以把降维结果同时画到坐标平面上，每个行画一个点，每个列也画一个点，从中查看每行行与每列接近与否。这样，可以发现行的类与列的类之间的对应关系。

记号:

设列联表保存为矩阵\(M = (x_{ij})_{n \times p}\), \(x_{ij}\)一般是非负整数，表示行变量取第\(i\)种值并且列变量取第\(j\)种值的观测个数（频数）。列联表的每一行对应于行变量的某个类，每一列对应于列变量的某个类。行、列变量可以是分类变量(如性别、行业、所在地区等)，也可以是数值型变量（如家中孩子个数，一般应存在很多重复值）, 或者数值型变量的分段值（比如年龄1—18, 19—60, 60以上）。

令 \[\begin{align*} N =& \sum_{i=1}^n \sum_{j=1}^p x_{ij} \quad\text{(总频数)}, \quad P = (p_{ij})_{n \times p} = \frac{1}{N} M \quad\text{(百分比)}, \\ x_{i\bullet} =& \sum_{j=1}^p x_{ij}, \ i=1,2,\dots,n \ \text{(行和)}, \\ x_{\bullet j} =& \sum_{i=1}^n x_{ij}, \ j=1,2,\dots,p \ \text{(列和)} \\ p_{i\bullet} =& \frac{x_{i\bullet}}{N} \ \text{(行变量边缘百分比),} \quad p_{\bullet j} = \frac{x_{\bullet j}}{N} \ \text{(列变量边缘百分比)} , \\ \boldsymbol a =& (p_{1 \bullet}, p_{2 \bullet}, \dots, p_{n \bullet})^T = P \boldsymbol 1_p, \\ \boldsymbol b =& (p_{\bullet 1}, p_{\bullet 2}, \dots, p_{\bullet p})^T = P^T \boldsymbol 1_n . \end{align*}\]

对应分析与主成分分析的区别:

主成分分析在计算时，把数据矩阵进行某种标准化后，取出能够最好地保留原来的方差的少数指标。对应分析要保留的是行与列的关联指标，通常是列联表独立性卡方统计量值。

10.2 对应分析的原理

10.2.1 线性代数依据

设\(C\)是\(n \times p\)矩阵，则\(C^T C\)是\(p \times p\)矩阵； \(C C^T\)是\(n \times n\)矩阵。 \(C^T C\)和\(C C^T\)都是对称半正定矩阵，所有特征值都非负。 \(C^T C\)和\(C C^T\)的正特征值都是相同的。

如果\(\boldsymbol\delta\)是\(C^T C\)的正特征值\(\lambda\)的特征向量，则\(C \boldsymbol\delta\)是\(C C^T\)的特征值\(\lambda\)的特征向量。如果\(\boldsymbol\gamma\)是\(C C^T\)的正特征值\(\lambda\)的特征向量，则\(C^T \boldsymbol\gamma\)是\(C^T C\)的特征向量。

证明:

设\(n \geq p\), \(C^T C\)的正特征值为 \(\lambda_1 \geq \lambda_2 \geq \dots \lambda_r > 0\), \(r \leq p \leq n\)。设\(\boldsymbol\delta_k\)是\(C^T C\)的对应于\(\lambda_k\)的特征向量，即 \[\begin{align*} & C^T C \boldsymbol\delta_k = \lambda_k \boldsymbol\delta_k, \\ & \text{则} \\ & C C^T C \boldsymbol\delta_k = \lambda_k C \boldsymbol\delta_k, \\ & \text{即} \\ & (C C^T) (C \boldsymbol\delta_k) = \lambda_k (C \boldsymbol\delta_k) \end{align*}\] 所以\(\lambda_k\)是\(C C^T\)的特征值，\(C \boldsymbol\delta_k\)是对应的特征向量。

反之类似。若\(\lambda_k>0\)和\(\boldsymbol\gamma_k\)是\(C C^T\)的正特征值和对应的特征向量，则 \[\begin{align*} C C^T \boldsymbol\gamma_k =& \lambda_k \boldsymbol\gamma_k \\ \text{于是} & \\ C^T C C^T \boldsymbol\gamma_k =& \lambda_k C^T \boldsymbol\gamma_k \\ \text{即} & \\ (C^T C) (C^T \boldsymbol\gamma_k) =& \lambda_k (C^T \boldsymbol\gamma_k) \end{align*}\] 所以\(\lambda_k\)也是\(C^T C\)的特征值，\(C^T \boldsymbol\gamma_k\)是对应的特征向量。

所以，\(C^T C\)与\(C C^T\)的正特征值完全相同，对应的特征向量可以互相决定。

10.2.1.1 奇异值分解

对\(n \times p\)矩阵\(C\)，设\(C^T C\)和\(C C^T\)共同的正特征值为\(\lambda_1, \dots, \lambda_r\), \(1 \leq r \leq \min(n,p)\)， \(C^T C\)对应于\(\lambda_k\)的特征向量为\(\boldsymbol\delta_k\), \(C C^T\)对应于\(\lambda_k\)的特征向量为\(\boldsymbol\gamma_k\), 令\(\Delta = (\boldsymbol\delta_1, \dots, \boldsymbol\delta_r)_{p \times r}\), \(\Gamma = (\boldsymbol\gamma_1, \dots, \boldsymbol\gamma_r)_{n \times r}\)。则\(C\)有如下的奇异值分解： \[\begin{align} C = \Gamma \Lambda \Delta^T = \sum_{k=1}^r \lambda_k^{1/2} \boldsymbol\gamma_k \boldsymbol\delta_k^T, \end{align}\] 其中\(\Lambda = \text{diag}(\lambda_1^{1/2}, \dots, \lambda_r^{1/2})\)。

注意在主成分分析中，特征向量用来作为各个行向量的线性组合系数（R型主成分），或各个列向量的线性组合系数（Q型主成分）。这样，如果同时得到了每行的前两个因子得分(R型主成分)，以及每列的前两个因子得分(R型主成分), 可以将两种第一主成分放在同一个横轴，将两种第二主成分放在同一个纵轴，同时把代表各行的散点与代表各列的散点画在同一坐标系中。

10.2.2 卡方分解

设\(M=(x_{ij})_{n \times p}\)是列联表，通常用如下的卡方统计量进行列联表独立性检验： \[\begin{align} V = \sum_{i=1}^n \sum_{j=1}^p \frac{(x_{ij} - N p_{i\bullet} p_{\bullet j} )^2}{N p_{i\bullet} p_{\bullet j}} = N \sum_{i=1}^n \sum_{j=1}^p \frac{(p_{ij} - p_{i\bullet} p_{\bullet j})^2}{p_{i\bullet} p_{\bullet j}} . \tag{10.1} \end{align}\]

记\(E_{ij}= p_{i\bullet} p_{\bullet j}\), 这是是在行、列变量独立情况下的联合概率， \(E = (E_{ij})_{n\times p} = \boldsymbol a \boldsymbol b^T\)。当行、列变量独立时，\(V\)近似服从\(\chi^2((n-1)(p-1))\)分布。

行、列独立时，对应分析没有意义。当行、列不独立时，对应分析可以发现哪些行与哪些列有正向或负向的关联，即代表行的变量与代表列的变量在哪些组合中比独立情形出现更多或更少。

令矩阵\(C = (c_{ij})_{n \times p}\), 其中 \[\begin{align*} c_{ij} = \frac{x_{ij} - N p_{i\bullet} p_{\bullet j}}{(N p_{i\bullet} p_{\bullet j})^{1/2}} = \sqrt{N} \frac{p_{ij} - p_{i\bullet} p_{\bullet j}}{(p_{i\bullet} p_{\bullet j})^{1/2}}, \ i=1,2,\dots,n, \, j=1,2,\dots, p . \end{align*}\] 这是(10.1)中通项的平方根，易见 \[\begin{align} V = \sum_{i=1}^n \sum_{j=1}^p c_{ij}^2 = \text{trace}(C^T C) . \end{align}\] 令\(A = \text{diag}(\boldsymbol a)\), \(B = \text{diag}(\boldsymbol b)\), 则\(A \boldsymbol 1_n = \boldsymbol a\), \(B \boldsymbol 1_p = \boldsymbol b\)。易见 \[\begin{align} C = \sqrt{N} A^{-1/2} (P - \boldsymbol a \boldsymbol b^T) B^{-1/2} . \tag{10.2} \end{align}\]

记 \[\begin{align*} \boldsymbol a^{1/2} =& (p_{1 \bullet}^{1/2}, p_{2 \bullet}^{1/2}, \dots, p_{n \bullet}^{1/2})^T, \\ \boldsymbol b^{1/2} =& (p_{\bullet 1}^{1/2}, p_{\bullet 2}^{1/2}, \dots, p_{\bullet p}^{1/2})^T, \\ \end{align*}\] 容易验证 \[\begin{align} C \boldsymbol b^{1/2} = \boldsymbol 0, \quad C^T \boldsymbol a^{1/2} = \boldsymbol 0 . \tag{10.3} \end{align}\]

事实上， \[\begin{align*} C \boldsymbol b^{1/2} =& \sqrt{N} A^{-1/2} (P - \boldsymbol a \boldsymbol b^T) B^{-1/2} B^{1/2} \boldsymbol 1_p \\ =& \sqrt{N} A^{-1/2} (P \boldsymbol 1_p - \boldsymbol a \boldsymbol b^T \boldsymbol 1_p) \\ =& \sqrt{N} A^{-1/2} ( \boldsymbol a - \boldsymbol a ) = \boldsymbol 0, \\ C^T \boldsymbol a^{1/2} =& \sqrt{N} B^{-1/2} (P^T - \boldsymbol b \boldsymbol a^T) A^{-1/2} A^{1/2} \boldsymbol 1_n \\ =& \sqrt{N} B^{-1/2} (P^T \boldsymbol 1_n - \boldsymbol b \boldsymbol a^T \boldsymbol 1_n) \\ =& \sqrt{N} B^{-1/2} ( \boldsymbol b - \boldsymbol b ) = \boldsymbol 0 . \end{align*}\]

设\(r = \text{rank}(C) \leq \min(n-1, p-1)\), 矩阵\(C\)有奇异值分解 \[\begin{align} C = \Gamma \Lambda \Delta^T = \sum_{k=1}^r \lambda_k^{1/2} \boldsymbol\gamma_k \boldsymbol\delta_k^T, \tag{10.4} \end{align}\] 其中\(\Lambda = \text{diag}(\lambda_1^{1/2}, \lambda_2^{1/2}, \dots, \lambda_r^{1/2})\), \(\lambda_1, \lambda_2, \dots, \lambda_r\)是\(C^T C\)和\(C C^T\)共同的正特征值， \(\Gamma^T \Gamma = I_r\), \(\Delta^T \Delta = I_r\), \(\Gamma_{n \times r}\)的第\(k\)列\(\boldsymbol\gamma_k\)是\(C C^T\)对应于\(\lambda_k\)的特征向量， \(\Delta_{p \times r}\)的第\(k\)列\(\boldsymbol\delta_k\)是\(C^T C\)对应于\(\lambda_k\)的特征向量。奇异值分解(10.4)说明 \[\begin{align} c_{ij} = \sum_{k=1}^r \lambda_k^{1/2} \gamma_{ik} \delta_{jk}, \ i=1,2,\dots, n, \, j=1,2,\dots, p . \tag{10.5} \end{align}\]

这里有 \[\begin{align} V =& \text{trace}(C^T C) = \text{trace}(\Delta \Lambda \Gamma^T \Gamma \Lambda \Delta^T) = \text{trace}(\Delta \Lambda^2\Delta^T) \nonumber\\ =& \text{trace}(\Lambda^2 \Delta^T \Delta) = \text{trace}(\Lambda^2) \nonumber\\ =& \sum_{k=1}^r \lambda_k . \end{align}\] 可见列联表独立性卡方统计量可以分解为\(C^T C\)的特征值之和，可以类似主成分分析解释方差那样解释列联表统计量的分解。称\(V = \sum_{k=1}^r \lambda_k\)为数据的总惯性（total inertia）。

当\(\sum_{k=1}^K \lambda_k / \sum_{k=1}^r \lambda_k\)较大时，在(10.4)中可以用前\(K\)项逼近\(C\)。由于\(\lambda_1 \geq \lambda_2 \geq \dots\), 对应于\(\lambda_1\)和\(\lambda_2\)的特征向量\(\boldsymbol\gamma_1, \boldsymbol\gamma_2\) 和\(\boldsymbol\delta_1, \boldsymbol\delta_2\)是更重要的特征向量, 常取\(K=2\)。

在奇异值分解式(10.4)中，\(\boldsymbol\gamma_k\)使得 \[\begin{align*} (C C^T) \boldsymbol\gamma_k =& \lambda_k \boldsymbol\gamma_k, \\ (C^T C) (C^T \boldsymbol\gamma_k) =& \lambda_k C^T \boldsymbol\gamma_k, \\ \text{而} & \\ (C^T C) \boldsymbol\delta_k =& \lambda_k \boldsymbol\delta_k, \\ \text{所以} & \\ \boldsymbol\delta_k =& d C^T \boldsymbol\gamma_k \quad(d\text{是待定系数}) \end{align*}\]

由 \[\begin{align*} 1 =& \boldsymbol\delta_k^T \boldsymbol\delta_k = d^2 \boldsymbol\gamma_k^T C C^T \boldsymbol\gamma_k \\ =& d^2 \gamma_k^T \lambda_k \boldsymbol\gamma_k = d^2 \lambda_k \end{align*}\] 得\(d = \pm \lambda_k^{-1/2}\)，所以 \[\begin{align} \boldsymbol\delta_k =& \pm \lambda_k^{-1/2} C^T \boldsymbol\gamma_k, \quad \boldsymbol\gamma_k = \pm \lambda_k^{-1/2} C \boldsymbol\delta_k , k=1,2,\dots, r . \tag{10.6} \end{align}\]

由(10.3)可知 \[\begin{align} \boldsymbol\delta_k^T \boldsymbol b^{1/2} = 0, \quad \boldsymbol\gamma_k^T \boldsymbol a^{1/2} = 0, \ k=1,2,\dots,r . \tag{10.7} \end{align}\] 这一点是因为\(\boldsymbol\delta_k\)是\(C^T C\)的正特征值\(\lambda_k\)的特征向量，而(10.3)说明\(\boldsymbol b^{1/2}\)是\(C^T C\)的零特征值的特征向量，两者必须正交。

类似于主成分分析中把\(\boldsymbol\delta_k\)作为计算主成分得分的线性组合系数，把\(C\)的每个行向量分别用\(\boldsymbol\delta_1, \boldsymbol\delta_2, \dots\)作线性组合，得到降维表示 \[\begin{align*} C \boldsymbol\delta_j = \lambda_j^{1/2} \boldsymbol\gamma_j, \ j=1,2,\dots \end{align*}\] 其中的前面几个可以作为\(C\)的各列的降维, 把原来\(C\)中每行的\(p\)个元素压缩到较低维数。

类似地，把\(C\)的每个列向量用\(\boldsymbol\gamma_1, \boldsymbol\gamma_2, \dots\)作线性组合，得到降维表示 \[\begin{align*} C^T \boldsymbol\gamma_i = \lambda_i^{1/2} \boldsymbol\delta_i, \ i=1,2,\dots \end{align*}\] 其中前面几个可以作为\(C\)的各行的降维表示，把原来\(C\)中每列的\(n\)个元素降维到较低维数。

例如，当\(\lambda_1\)远大于其它特征值时， \(n\)维列向量\(\boldsymbol\gamma_1\)可以作为\(C\)的\(p\)列的降维表示，可以用\(\boldsymbol\gamma_1\)的\(n\)个元素作为\(C\)的行变量的\(n\)个类的代表。

同时，\(p\)维行向量\(\boldsymbol\delta_1^T\)可以作为\(C\)的各行的降维表示，可以用\(\boldsymbol\delta_1\)的\(p\)个元素作为\(C\)的列变量的\(p\)个类的代表。且这时 \[\begin{align*} c_{ij} \approx \lambda_1^{1/2} \gamma_{i1} \delta_{j1} . \end{align*}\]

当\(\gamma_{i1}\)和\(\delta_{j1}\)同号且都比较大时， \(c_{ij}\)也比较大，这时行变量的第\(i\)类与列变量的第\(j\)类存在正向的关联；当\(\gamma_{i1}\)和\(\delta_{j1}\)异号且绝对值都比较大时， \(c_{ij}\)为负且绝对值较大，这时行变量的第\(i\)类与列变量的第\(j\)类存在负向的关联。这里正向的关联是指这两个类同时出现概率比独立情形明显地大。当\(\gamma_{i1}\)和\(\delta_{j1}\)绝对值较小时，行变量的第\(i\)类与列变量的第\(j\)类没有什么关联。

在更一般的情形中，\(\lambda_1\)和\(\lambda_2\) 在卡方统计量值\(V=\sum_{k=1}^r \lambda_k\)中占据较大比重。这时，用 \[\begin{align*} (\boldsymbol\gamma_1, \boldsymbol\gamma_2)_{n\times 2} =& \left(\begin{array}{cc} \gamma_{11} & \gamma_{12} \\ \gamma_{21} & \gamma_{22} \\ \vdots & \vdots \\ \gamma_{n1} & \gamma_{n2} \end{array}\right) \end{align*}\] 作为\(C\)的每个行向量的降维表示，用这\(n\)个二维行向量近似表示\(C\)的\(n\)行（原来每行是\(p\)维，现在每行用二维近似表示）。

用 \[\begin{align*} \left(\begin{array}{c} \boldsymbol\delta_1^T \\ \boldsymbol\delta_2^T \end{array}\right)_{2 \times p} =& \left(\begin{array}{ccccc} \delta_{11} & \delta_{21} & \cdots & \delta_{p1} \\ \delta_{12} & \delta_{22} & \cdots & \delta_{p2} \end{array}\right) \end{align*}\] 作为\(C\)的每个列向量的降维表示，用这\(p\)个二维列向量近似表示\(C\)的\(p\)个列（原来每列是\(n\)维，现在每列用二维近似表示）。

把\((\boldsymbol\gamma_1, \boldsymbol\gamma_2)_{n\times 2}\)和 \((\boldsymbol\delta_1, \boldsymbol\delta_2)_{p\times 2}\)适当标准化后，可以分别作散点图并画在同一坐标平面中。这样的散点图可以用来考察行变量类\(i\)与列变量\(j\)的关联性。

10.2.3 标准化的降维坐标

为了将\(C\)的每个行向量降维，计算\(C \boldsymbol\delta_1, C \boldsymbol\delta_2\), 即用\(\boldsymbol\delta_1\)和\(\boldsymbol\delta_2\)对\(C\)的每行作线性组合。因为\(\| \delta_k \|=1\)，所以这样降维得到的第\(i\)行两个新坐标受到\(P\)的第\(i\)行的行边缘百分比\(p_{i \bullet}\)影响。

类似地，用\(C^T \boldsymbol\gamma_1, C^T \boldsymbol\gamma_2\)对\(C\)的每个列向量降维，得到的第\(j\)列的两个新坐标受到\(P\)的第\(j\)列的列边缘百分比\(p_{\bullet j}\)影响。

我们希望降维后的行坐标仅代表该行的“行轮廓”(row profile)即行百分比 \(\left( \frac{p_{i1}}{p_{i\bullet}}, \frac{p_{i2}}{p_{i\bullet}}, \dots, \frac{p_{ip}}{p_{i\bullet}} \right)\)，希望降维后的列坐标仅代表该列的“列轮廓”(column profile)即列百分比 \(\left( \frac{p_{1j}}{p_{\bullet j}}, \frac{p_{2j}}{p_{\bullet j}}, \dots, \frac{p_{nj}}{p_{\bullet j}} \right)\)。

为了把\(C\)的各个行向量进行标准化降维，需要改造\(C\)使得各行不再受行和的影响。

原来 \[\begin{align*} c_{ij} = \sqrt{N} \frac{p_{ij} - p_{i\bullet} p_{\bullet j}} { p_{i\bullet}^{1/2} p_{\bullet j}^{1/2}}, \end{align*}\] 令\(\tilde C = N^{-1/2} A^{-1/2} C\), 则 \[\begin{align} \tilde c_{ij} = \frac{1}{\sqrt{N} p_{i\bullet}^{1/2}} c_{ij} = \frac{ \frac{p_{ij}}{p_{i\bullet}} - p_{\bullet j}}{p_{\bullet j}^{1/2}} . \tag{10.8} \end{align}\] 这里\((\frac{p_{ij}}{p_{i\bullet}}, j=1,2,\dots,p)\)是第\(i\)行的行百分比，代表了第\(i\)行的特殊分布，而\((p_{\bullet j}, j=1,2,\dots,p)\)代表列变量的边缘分布。 \((\tilde c_{ij}, j=1,2,\dots, p)\) 代表了第\(i\)行的\(p\)个比例与列变量的边缘分布的差异。

于是，为了把\(C\)的各个行向量进行标准化降维，令 \[\begin{align} \boldsymbol r_k = \tilde C \boldsymbol\delta_k = N^{-1/2} A^{-1/2} C \boldsymbol\delta_k = N^{-1/2} \lambda_k^{1/2} A^{-1/2} \boldsymbol\gamma_k, \ k=1,2,\dots \tag{10.9} \end{align}\] 这样，仅取前两个特征值， \[\begin{align*} (\boldsymbol r_1, \boldsymbol r_2)_{n\times 2} =& \left(\begin{array}{cc} r_{11} & r_{12} \\ r_{21} & r_{22} \\ \vdots & \vdots \\ r_{n1} & r_{n2} \\ \end{array}\right) \end{align*}\] 可以用来代表原来的\(n\)行，把每行从\(p\)维降低到二维。

以\((r_{i1}, r_{i2})\), \(i=1,2,\dots,n\)作为行变量的\(n\)个类的降维坐标画散点图。这里，\(r_{i1} = N^{-1/2} \lambda_1^{1/2} p_{i\bullet}^{-1/2} \gamma_{i1}\), \(r_{i2} = N^{-1/2} \lambda_2^{1/2} p_{i\bullet}^{-1/2} \gamma_{i2}\)。

类似地，为了把\(C\)的各个列向量进行标准化降维，令 \[\begin{align} \boldsymbol s_k =& N^{-1/2} B^{-1/2} C^T \boldsymbol\gamma_k = N^{-1/2} \lambda_k^{1/2} B^{-1/2} \boldsymbol\delta_k . \tag{10.10} \end{align}\] 这样，仅取前两个特征值， \[\begin{align*} \left(\begin{array}{c} \boldsymbol s_1^T \\ \boldsymbol s_2^T \end{array}\right)_{2 \times p} =& \left(\begin{array}{ccccc} s_{11} & s_{21} & \cdots & s_{p1} \\ s_{12} & s_{22} & \cdots & s_{p2} \end{array}\right) \end{align*}\] 可以用来代表原来的\(p\)列，每列从\(n\)维降低到二维。

以\((s_{j1}, s_{j2})\), \(j=1,2,\dots,p\)作为列变量的\(p\)个类的降维坐标画散点图。与行变量的降维散点图画在同一坐标平面中。这里，\(s_{j1} = N^{-1/2} \lambda_1^{1/2} p_{\bullet j}^{-1/2} \delta_{j1}\), \(s_{j2} = N^{-1/2} \lambda_2^{1/2} p_{\bullet j}^{-1/2} \delta_{j2}\)。

由(10.7)可知 \[\begin{align} \boldsymbol r_k^T \boldsymbol a = 0, \quad \boldsymbol s_k^T \boldsymbol b = 0, k=1,2,\dots,r . \tag{10.11} \end{align}\] 所以，上述散点图中的四种横、纵坐标都是满足加权平均为零的，其取值都在原点左右（或上下）相对平衡地分布，加权系数为\(\boldsymbol a\)或\(\boldsymbol b\)。

由(10.6)易见 \[\begin{align*} \begin{cases} \boldsymbol r_k = \pm \lambda_k^{-1/2} A^{-1/2} C B^{1/2} \boldsymbol s_k, \\ \boldsymbol s_k = \pm \lambda_k^{-1/2} B^{-1/2} C^T A^{1/2} \boldsymbol r_k . \\ \end{cases} \end{align*}\] 这又可以化简为 \[\begin{align} \begin{cases} \boldsymbol r_k = \pm N^{1/2} \lambda_k^{-1/2} A^{-1} P \boldsymbol s_k, \\ \boldsymbol s_k = \pm N^{1/2} \lambda_k^{-1/2} B^{-1} P^T \boldsymbol r_k . \end{cases} \tag{10.12} \end{align}\] 这说明压缩后的\(n\)行的\(n\)个第\(1\)坐标\(\boldsymbol r_1\) 与压缩后的\(p\)列的\(p\)个第\(1\)坐标\(\boldsymbol s_1\)存在一一对应关系, 其它坐标类似。

证明:

来证明(10.12)。事实上， \[\begin{align*} \boldsymbol r_k =& \pm \lambda_k^{-1/2} A^{-1/2} C B^{1/2} \boldsymbol s_k \\ =& \pm N^{1/2} \lambda_k^{-1/2} A^{-1/2} \; A^{-1/2} (P - \boldsymbol a \boldsymbol b^T) B^{-1/2} \; B^{1/2} \boldsymbol s_k \\ =& \pm N^{1/2} \lambda_k^{-1/2} A^{-1} (P - \boldsymbol a \boldsymbol b^T) \boldsymbol s_k \\ =& \pm N^{1/2} \lambda_k^{-1/2} A^{-1} \left( P \boldsymbol s_k - \boldsymbol a \boldsymbol b^T \, \boldsymbol s_k \right) \\ =& \pm N^{1/2} \lambda_k^{-1/2} A^{-1} P \boldsymbol s_k \end{align*}\] (注意\(\boldsymbol b^T \boldsymbol s_k = 0\))。而 \[\begin{align*} \boldsymbol s_k =& \pm \lambda_k^{-1/2} B^{-1/2} C^T A^{1/2} \boldsymbol r_k \\ =& \pm N^{1/2} \lambda_k^{-1/2} B^{-1/2} \, B^{-1/2} (P^T - \boldsymbol b \boldsymbol a^T) A^{-1/2} \, A^{1/2} \boldsymbol r_k \\ =& \pm N^{1/2} \lambda_k^{-1/2} B^{-1} ( P^T \boldsymbol r_k - \boldsymbol b \boldsymbol a^T \boldsymbol r_k) \\ =& \pm N^{1/2} \lambda_k^{-1/2} B^{-1} P^T \boldsymbol r_k \end{align*}\] (注意\(\boldsymbol a^T \boldsymbol r_k = 0\))。

10.3 对应分析做法

从\(M=(x_{ij})_{n\times p}\)得\(P = \frac{1}{N} M\), 计算得到\(C=(c_{ij})_{n \times p}\)。
求\(C^T C\)的前两个特征值\(\lambda_1, \lambda_2\)(比例不足时可以取更多个), 得到\(C^T C\)的前两个特征向量\(\boldsymbol\delta_1, \boldsymbol\delta_2\), 然后利用(10.6)计算\(C C^T\)的前两个特征向量 \(\boldsymbol\gamma_1, \boldsymbol\gamma_2\)。
用(10.9)和(10.10)得到标准化的行得分 \(\boldsymbol r_1, \boldsymbol r_2\)和列得分\(\boldsymbol s_1, \boldsymbol s_2\)。
以\((r_{i1}, r_{i2})\), \(i=1,2,\dots, n\)为坐标画行类的散点图，以\((s_{j1}, s_{j2})\), \(j=1,2,\dots,p\)为坐标画列类的散点图。

10.3.1 对应分析结果图形的解释—两个行类之间或两个列类之间

对于图上的两个行类，其散点相近，表示这两行的行百分比很相近。
对于图上两个列类，其散点相近，表示这两列的列百分比很相近。
两个行类的散点距离远则代表这两行的行百分比取值差别大。两个列类的散点距离远则代表这两列的列百分比取值差别大。
以上解释是来自于标准化的行类、列类的降维方法。
因为横坐标\(\boldsymbol r_1, \boldsymbol s_1\)的加权平均为零，纵坐标\(\boldsymbol r_2, \boldsymbol s_2\)的加权平均为零，所以落在原点附近的行类的行百分比取值与列变量的边缘百分比 \(\frac{1}{N} \boldsymbol b^T\)相近，落在原点附近的列类的列百分比取值与行变量的边缘百分比\(\frac{1}{N} \boldsymbol a\)相近。

10.3.2 对应分析结果图形的解释—一个行类与一个列类之间

当行类\(i\)与列类\(j\)的散点距离近而且都远离原点时，这个行类与这个列类有较强的正向关联；
当行类\(i\)与列类\(j\)的散点到原点的距离较远且位于相反方向时，这个行类与这个列类有较强的负向关联；
当行类\(i\)与列类\(j\)至少有一个靠近原点时，这两个类的关联不大。
当\(\lambda_1 + \lambda_2\)占据\(V = \sum_{k=1}^r \lambda_k\)的很大比例时，对应分析的散点图是有意义的，否则散点图意义不明显。
事实上，\(c_{ij}\)绝对值大小代表了第\(i\)行类与第\(j\)行类的关联性强弱，当\(\lambda_1 + \lambda_2\)占据\(V = \sum_{k=1}^r \lambda_k\)的很大比例时， \[\begin{align} c_{ij} =& \sum_{k=1}^r \lambda_k^{1/2} \gamma_{ik} \delta_{jk} \approx \sum_{k=1}^2 \lambda_k^{1/2} \gamma_{ik} \delta_{jk} \tag{10.13} \\ =& \sum_{k=1}^2 \lambda_k^{1/2} \cdot N^{1/2} \lambda_k^{-1/2} p_{i\bullet}^{1/2} r_{ik} \cdot N^{1/2} \lambda_k^{-1/2} p_{\bullet j}^{1/2} s_{jk} \nonumber\\ =& N (p_{i \bullet} p_{\bullet j})^{1/2} \sum_{k=1}^2 \lambda_k^{-1/2} r_{ik} s_{jk} . \tag{10.14} \end{align}\]
\(\sum_{k=1}^2 \lambda_k^{-1/2} r_{ik} s_{jk}\)是向量\((r_{i1}, r_{i2})\) 和向量\((s_{ji}, s_{j2})\)的一种内积。注意(10.13)中的\(\gamma_{ik}\)和\(\delta_{jk}\) 都来自单位特征向量，所以(10.13)中的加权 \(\lambda_1^{1/2}\)和\(\lambda_2^{1/2}\)使得内积更侧重于第一特征向量； (10.14)中的\(r_{ik}\)和\(s_{ik}\)并不是来自单位向量。
当\((r_{i1}, r_{i2})\)与\((s_{j1}, s_{j2})\)两个散点接近而且离原点较远时，这两个向量同方向，内积较大，使得\(c_{ij}\)较大，这时第\(i\)行类与第\(j\)列类有较强的正向关联。

10.3.3 推广用法

因为对应分析计算中只要求列联表\(M = (x_{ij})_{n \times p}\)所有元素非负，所以对应分析也可以对一般的数据集\(M\)计算，但是要求元素都为正而且各变量之间是可比的。
这样的推广用法可以作为主成分分析的补充。主成分分析关心的是解释更多的方差，而对应分析则考虑更一般的关联性。

10.4 案例8.1：现金支出定位问题

P. 112案例8.1。2012年我国城镇居民按收入等级分组，每种收入等级内人均年消费的不同种类分配情况表。行变量为收入分组，共7个级别: ‘最低’, ‘较低’, ‘中下’, ‘中等’, ‘中上’, ‘较高’, ‘最高’ 列变量为不同的消费项目（肉、蛋、鱼、奶、衣、耐用、文娱）。

library(MASS)
incomeUse <- readr::read_csv(
  "data/不同收入组支出分配.csv",
  locale=locale(encoding = "GB18030"),
  show_col_types = FALSE)
d0 <- as.data.frame(incomeUse)
rownames(d0) <- c('最低', '较低', '中下', '中等', '中上', '较高', '最高')
d <- d0[,-1]
colnames(d) <- c('肉', '蛋', '鱼', '奶', '衣', '耐', '文')

knitr::kable(d)

	肉	蛋	鱼	奶	衣	耐	文
最低	767.51	84.30	173.39	125.75	494.80	119.37	144.68
较低	946.71	96.64	235.60	169.02	746.36	185.89	209.60
中下	1088.30	112.02	308.75	208.34	1018.94	272.85	284.93
中等	1249.37	125.68	412.72	260.07	1288.06	395.16	404.79
中上	1341.13	133.33	522.70	308.80	1637.68	532.17	565.38
较高	1480.40	142.42	630.61	365.39	2067.91	738.16	793.42
最高	1555.67	147.05	768.17	423.34	3019.31	1195.76	1187.24

作列联表独立性卡方检验。虽然这里不是列联表但还是可以用的。

chisq.test(d)

## 
##  Pearson's Chi-squared test
## 
## data:  d
## X-squared = 1030.1, df = 36, p-value < 2.2e-16

用corresp()作对应分析, 用nf指定因子个数。为作散点图一般取2:

co <- corresp(d, nf=2)

## Warning in corresp.matrix(as.matrix(x), ...): negative or non-integer entries in
## table

co

## First canonical correlation(s): 0.17933228 0.01896212 
## 
##  Row scores:
##             [,1]       [,2]
## 最低 -1.96830644 -1.3285223
## 较低 -1.44856190 -0.9070070
## 中下 -0.97068022 -0.4603836
## 中等 -0.54845375  0.5670742
## 中上 -0.03232808  1.3639620
## 较高  0.43507458  0.9414503
## 最高  1.25763257 -1.0453509
## 
##  Column scores:
##           [,1]       [,2]
## 肉 -1.36707828 -0.4637594
## 蛋 -1.59707661 -1.1112356
## 鱼  0.02725855  2.6175325
## 奶 -0.46697640  1.5261695
## 衣  0.48491880 -0.2497690
## 耐  1.33266208 -0.8879205
## 文  1.13861727 -0.1015123

输出中典型相关系数的解释：把行变量的每个类当作两点分布结果，共有\(n\)个二值行变量；把列变量的每个类当作两点分布结果，共有\(p\)个二值列变量； \(c_{ij}\)近似等于第\(i\)个二值行变量和第\(j\)个二值列表量的相关系数，把所有\(n\)个二值行变量当成一组，所有\(p\)个二值列变量当成一组，可以计算两组的典型相关系数。

作对应分析的biplot():

biplot(co)
abline(h=0, v=0, lty=3, col='gray')

结果分析：

纵坐标都很小，说明第二因子基本不起作用。
在第一坐标上，比较均匀地分布了从最低到最高的七个收入组，次序完全正确，反映了七个行类在对应分析图中的接近关系是完全正确的。
关于行类（收入组）与列类（消费类型）之间的关系：
最高收入组，其支出结构更偏向耐用消费品和文娱。
较高收入组，其支出结构更偏向服装。
较低收入组，其支出结构在肉类、蛋类占得较多。注意，并不是说这个收入组花费在肉类、蛋类的钱比高收入组还多，这里比的是相对比例。
中等收入组，其支出结构在奶和奶制品占得相对较多。
中上收入组靠近原点，其支出结构和平均的结构相近。

10.4.1 多元图形辅助分析

收入组与支出分配的对应分析的多元图形:

library(graphics)
library(aplpack)

用星图验证分类:

stars(log(d), scale=TRUE,
      len=0.6, cex=1, key.loc=c(5, 2),
      draw.segments=TRUE)

因为归一化，所以不如对应分析结果可解释。

脸谱图：

faces(d)

## effect of variables:
##  modified item       Var 
##  "height of face   " "肉"
##  "width of face    " "蛋"
##  "structure of face" "鱼"
##  "height of mouth  " "奶"
##  "width of mouth   " "衣"
##  "smiling          " "耐"
##  "height of eyes   " "文"
##  "width of eyes    " "肉"
##  "height of hair   " "蛋"
##  "width of hair   "  "鱼"
##  "style of hair   "  "奶"
##  "height of nose  "  "衣"
##  "width of nose   "  "耐"
##  "width of ear    "  "文"
##  "height of ear   "  "肉"

这些图没有给出很好的解释。

10.5 苏格兰Caithness地区的人的眼睛颜色和头发颜色组合情况

眼睛颜色和头发颜色列联表:

data(caith)
d <- caith

knitr::kable(d)

	fair	red	medium	dark	black
blue	326	38	241	110	3
light	688	116	584	188	4
medium	343	84	909	412	26
dark	98	48	403	681	85

作列联表独立性卡方检验:

chisq.test(d)

## 
##  Pearson's Chi-squared test
## 
## data:  d
## X-squared = 1240, df = 12, p-value < 2.2e-16

10.5.1 作对应分析，取两个特征值

co <- corresp(d, nf=2)
co

## First canonical correlation(s): 0.4463684 0.1734554 
## 
##  Row scores:
##               [,1]       [,2]
## blue   -0.89679252  0.9536227
## light  -0.98731818  0.5100045
## medium  0.07530627 -1.4124778
## dark    1.57434710  0.7720361
## 
##  Column scores:
##               [,1]       [,2]
## fair   -1.21871379  1.0022432
## red    -0.52257500  0.2783364
## medium -0.09414671 -1.2009094
## dark    1.31888486  0.5992920
## black   2.45176017  1.6513565

10.5.2 作对应分析的biplot

biplot(co)
abline(h=0, v=0, lty=3, col='gray')

从眼睛颜色与头发颜色关系看：

发色fair则眼睛颜色多blue, light;
发色medium, 眼睛颜色多medium;
发色black, 眼睛颜色多dark。

10.6 奶牛场规模与患病

奶牛场规模小、中、大与某种牛传染病发病无、轻、重的列联表。

d1 <- data.frame(
    herdsize=factor(
      rep(c('小', '中', '大'), each=3),
      levels=c('小', '中', '大')),
    disease=factor(
      rep(c('无', '轻', '重'), 3),
      levels=c('无', '轻', '重')),
    count=c(9, 5, 9,  18, 4, 19,  11, 88, 136))
## 下面把d1转置成列联表格式, 用xtabs()函数
d <- xtabs(count ~ herdsize + disease, data=d1)

knitr::kable(d)

	无	轻	重
小	9	5	9
中	18	4	19
大	11	88	136

作对应分析，取两个特征值:

co <- corresp(d, nf=2)
co

## First canonical correlation(s): 0.47152500 0.04338242 
## 
##  herdsize scores:
##          [,1]        [,2]
## 小 -1.6594198  3.04077718
## 中 -2.0490587 -1.44708027
## 大  0.5199067 -0.04513866
## 
##  disease scores:
##          [,1]       [,2]
## 无 -2.5727754  0.4992470
## 轻  0.6396978  1.2935459
## 重  0.2177730 -0.8807642

作对应分析的biplot:

biplot(co)
abline(h=0, v=0, lty=3, col='gray')

图形解释:

第二因子作用很小。
病情轻、重很接近，与无方向相隔很远，是合理的。
大的牛场倾向于有病；中、小牛场倾向于无病。

10.7 美国50个州数据

state.x77数据，美国50个州在1970年代的基本情况数据。变量:

Population: 人口
Illiteracy: 文盲比例
Life Exp: 期望寿命
Murder: 每十万人杀人犯罪数
HS Grad: 高中毕业生比例
Frost: 日最低温度在结霜点以下平均天数
Area: 面积，平方英里

data(state.x77)

## Warning in data(state.x77): data set 'state.x77' not found

d <- state.x77
rownames(d) <- state.abb
colnames(d) <- c('人口', '收入', '文盲',
                 '寿命', '谋杀', '高中生',
                 '无霜期', '面积')

knitr::kable(d)

	人口	收入	文盲	寿命	谋杀	高中生	无霜期	面积
AL	3615	3624	2.1	69.05	15.1	41.3	20	50708
AK	365	6315	1.5	69.31	11.3	66.7	152	566432
AZ	2212	4530	1.8	70.55	7.8	58.1	15	113417
AR	2110	3378	1.9	70.66	10.1	39.9	65	51945
CA	21198	5114	1.1	71.71	10.3	62.6	20	156361
CO	2541	4884	0.7	72.06	6.8	63.9	166	103766
CT	3100	5348	1.1	72.48	3.1	56.0	139	4862
DE	579	4809	0.9	70.06	6.2	54.6	103	1982
FL	8277	4815	1.3	70.66	10.7	52.6	11	54090
GA	4931	4091	2.0	68.54	13.9	40.6	60	58073
HI	868	4963	1.9	73.60	6.2	61.9	0	6425
ID	813	4119	0.6	71.87	5.3	59.5	126	82677
IL	11197	5107	0.9	70.14	10.3	52.6	127	55748
IN	5313	4458	0.7	70.88	7.1	52.9	122	36097
IA	2861	4628	0.5	72.56	2.3	59.0	140	55941
KS	2280	4669	0.6	72.58	4.5	59.9	114	81787
KY	3387	3712	1.6	70.10	10.6	38.5	95	39650
LA	3806	3545	2.8	68.76	13.2	42.2	12	44930
ME	1058	3694	0.7	70.39	2.7	54.7	161	30920
MD	4122	5299	0.9	70.22	8.5	52.3	101	9891
MA	5814	4755	1.1	71.83	3.3	58.5	103	7826
MI	9111	4751	0.9	70.63	11.1	52.8	125	56817
MN	3921	4675	0.6	72.96	2.3	57.6	160	79289
MS	2341	3098	2.4	68.09	12.5	41.0	50	47296
MO	4767	4254	0.8	70.69	9.3	48.8	108	68995
MT	746	4347	0.6	70.56	5.0	59.2	155	145587
NE	1544	4508	0.6	72.60	2.9	59.3	139	76483
NV	590	5149	0.5	69.03	11.5	65.2	188	109889
NH	812	4281	0.7	71.23	3.3	57.6	174	9027
NJ	7333	5237	1.1	70.93	5.2	52.5	115	7521
NM	1144	3601	2.2	70.32	9.7	55.2	120	121412
NY	18076	4903	1.4	70.55	10.9	52.7	82	47831
NC	5441	3875	1.8	69.21	11.1	38.5	80	48798
ND	637	5087	0.8	72.78	1.4	50.3	186	69273
OH	10735	4561	0.8	70.82	7.4	53.2	124	40975
OK	2715	3983	1.1	71.42	6.4	51.6	82	68782
OR	2284	4660	0.6	72.13	4.2	60.0	44	96184
PA	11860	4449	1.0	70.43	6.1	50.2	126	44966
RI	931	4558	1.3	71.90	2.4	46.4	127	1049
SC	2816	3635	2.3	67.96	11.6	37.8	65	30225
SD	681	4167	0.5	72.08	1.7	53.3	172	75955
TN	4173	3821	1.7	70.11	11.0	41.8	70	41328
TX	12237	4188	2.2	70.90	12.2	47.4	35	262134
UT	1203	4022	0.6	72.90	4.5	67.3	137	82096
VT	472	3907	0.6	71.64	5.5	57.1	168	9267
VA	4981	4701	1.4	70.08	9.5	47.8	85	39780
WA	3559	4864	0.6	71.72	4.3	63.5	32	66570
WV	1799	3617	1.4	69.48	6.7	41.6	100	24070
WI	4589	4468	0.7	72.48	3.0	54.5	149	54464
WY	376	4566	0.6	70.29	6.9	62.9	173	97203

把各列都归一化到[1,100]内:

cmin <- apply(d, 2, min)
cmax <- apply(d, 2, max)
cra <- cmax - cmin
d <- 1 + scale(d, center=cmin, scale=cra)*99

作对应分析，取两个特征值:

co <- corresp(d, nf=2)

## Warning in corresp.matrix(d, nf = 2): negative or non-integer entries in table

co

## First canonical correlation(s): 0.4409611 0.2444077 
## 
##  Row scores:
##           [,1]         [,2]
## AL -2.20149269 -0.542965161
## AK -0.01576789 -0.789803960
## AZ -0.55117605 -0.326692416
## AR -1.37991754 -1.145878975
## CA -0.34790488  2.917753726
## CO  0.73194052 -0.287074654
## CT  0.69466978 -0.118727593
## DE  0.41142720 -0.564367149
## FL -0.65670163  1.351785206
## GA -1.82576130 -0.263388506
## HI -0.07737277 -0.374989147
## ID  0.88627116 -0.586078435
## IL -0.14211633  1.380675923
## IN  0.35846989  0.466472638
## IA  1.21710818  0.008725745
## KS  0.92379266 -0.062928170
## KY -1.15502945 -0.683091362
## LA -2.55092046 -0.901072520
## ME  1.03160918 -1.013970323
## MD  0.10018053  0.359506400
## MA  0.51371677  0.530696459
## MI -0.13877369  0.961650956
## MN  1.08950670  0.047468793
## MS -2.55064321 -1.553753389
## MO -0.02780331  0.273019338
## MT  0.85220769 -0.747805794
## NE  1.12446618 -0.370803801
## NV  0.48584817 -0.639999952
## NH  1.06242729 -0.847482114
## NJ  0.21419679  0.833466286
## NM -0.71656957 -1.431216293
## NY -0.62005972  2.332760102
## NC -1.52184972 -0.151884256
## ND  1.11033958 -0.805523565
## OH  0.12284257  1.444811688
## OK  0.08235278 -0.273159323
## OR  0.84322394  0.362572895
## PA -0.03360095  1.605533210
## RI  0.46497447 -0.942590729
## SC -2.37399259 -1.281210099
## SD  1.35270900 -0.841481787
## TN -1.21730665 -0.335524928
## TX -1.21222663  0.966800790
## UT  1.06582655 -0.496582396
## VT  0.94318448 -0.930617172
## VA -0.53953828  0.182976274
## WA  0.80914377  0.800396474
## WV -0.72738284 -1.152727330
## WI  0.88958550  0.164651407
## WY  0.78056003 -0.798408387
## 
##  Column scores:
##              [,1]        [,2]
## 人口   -0.7387762  3.42345189
## 收入    0.4110482  0.47045368
## 文盲   -2.0307400 -1.05538283
## 寿命    0.7094660 -0.02244498
## 谋杀   -1.3476948 -0.12249617
## 高中生  0.8287947  0.02609934
## 无霜期  0.6810850 -0.81204484
## 面积   -0.1071694 -0.31923610

作对应分析的biplot:

biplot(co)
abline(h=0, v=0, lty=3, col='gray')

print(cbind(state.abb, state.name))

##       state.abb state.name      
##  [1,] "AL"      "Alabama"       
##  [2,] "AK"      "Alaska"        
##  [3,] "AZ"      "Arizona"       
##  [4,] "AR"      "Arkansas"      
##  [5,] "CA"      "California"    
##  [6,] "CO"      "Colorado"      
##  [7,] "CT"      "Connecticut"   
##  [8,] "DE"      "Delaware"      
##  [9,] "FL"      "Florida"       
## [10,] "GA"      "Georgia"       
## [11,] "HI"      "Hawaii"        
## [12,] "ID"      "Idaho"         
## [13,] "IL"      "Illinois"      
## [14,] "IN"      "Indiana"       
## [15,] "IA"      "Iowa"          
## [16,] "KS"      "Kansas"        
## [17,] "KY"      "Kentucky"      
## [18,] "LA"      "Louisiana"     
## [19,] "ME"      "Maine"         
## [20,] "MD"      "Maryland"      
## [21,] "MA"      "Massachusetts" 
## [22,] "MI"      "Michigan"      
## [23,] "MN"      "Minnesota"     
## [24,] "MS"      "Mississippi"   
## [25,] "MO"      "Missouri"      
## [26,] "MT"      "Montana"       
## [27,] "NE"      "Nebraska"      
## [28,] "NV"      "Nevada"        
## [29,] "NH"      "New Hampshire" 
## [30,] "NJ"      "New Jersey"    
## [31,] "NM"      "New Mexico"    
## [32,] "NY"      "New York"      
## [33,] "NC"      "North Carolina"
## [34,] "ND"      "North Dakota"  
## [35,] "OH"      "Ohio"          
## [36,] "OK"      "Oklahoma"      
## [37,] "OR"      "Oregon"        
## [38,] "PA"      "Pennsylvania"  
## [39,] "RI"      "Rhode Island"  
## [40,] "SC"      "South Carolina"
## [41,] "SD"      "South Dakota"  
## [42,] "TN"      "Tennessee"     
## [43,] "TX"      "Texas"         
## [44,] "UT"      "Utah"          
## [45,] "VT"      "Vermont"       
## [46,] "VA"      "Virginia"      
## [47,] "WA"      "Washington"    
## [48,] "WV"      "West Virginia" 
## [49,] "WI"      "Wisconsin"     
## [50,] "WY"      "Wyoming"

人口变量：CA和NY都是人口多的。

寿命变量与高中毕业生比例相近，说明这两个变量相关性强。

图中离谋杀变量近的并不是谋杀犯罪率最高的。

10.8 附录：用ggplot2改写biplot

require(ggplot2)
require(ggrepel)
require(tibble)
ggbiplot.corresp <- function(obj) {
  rscore <- tibble(
    label = rownames(obj$rscore),
    x = obj$rscore[,1],
    y = obj$rscore[,2]
  )
  cscore <- tibble(
    label = rownames(obj$cscore),
    x = obj$cscore[,1],
    y = obj$cscore[,2]
  )
  p <- ggplot(mapping = aes(x = x, y = y, label=label)) +
    geom_text_repel(data = rscore, color="black") +
    geom_text_repel(data = cscore, color="red") +
    geom_hline(yintercept=0, linetype=3, col="gray") +
    geom_vline(xintercept=0, linetype=3, col="gray") +
    labs(x=NULL, y=NULL)
  p
}

现金支出问题的测试：

library(MASS)
d0 <- as.data.frame(incomeUse)
rownames(d0) <- d0[[1]]
d <- d0[,-1]
co <- corresp(d, nf=2)

## Warning in corresp.matrix(as.matrix(x), ...): negative or non-integer entries in
## table

ggbiplot.corresp(co)

10.9 使用FactoMineR包

例子来自FactoMineR包文档。

library(FactoMineR)
data(children)
res.ca <- CA (children, row.sup = 15:18, col.sup = 6:8)

summary(res.ca)

## 
## Call:
## CA(X = children, row.sup = 15:18, col.sup = 6:8) 
## 
## The chi square of independence between the two variables is equal to 98.80159 (p-value =  9.748064e-05 ).
## 
## Eigenvalues
##                        Dim.1   Dim.2   Dim.3   Dim.4
## Variance               0.035   0.013   0.007   0.006
## % of var.             57.043  21.132  11.764  10.061
## Cumulative % of var.  57.043  78.175  89.939 100.000
## 
## Rows (the 10 first)
##                       Iner*1000    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## money               |     3.759 | -0.115  4.550  0.428 |  0.020  0.371  0.013 |
## future              |     8.690 |  0.176 17.567  0.716 | -0.098 14.587  0.220 |
## unemployment        |     9.151 | -0.212 22.616  0.875 | -0.071  6.779  0.097 |
## circumstances       |     3.804 |  0.401  6.274  0.584 |  0.331 11.544  0.398 |
## hard                |     1.199 | -0.250  2.994  0.884 |  0.068  0.592  0.065 |
## economic            |     8.787 |  0.354 12.005  0.484 |  0.321 26.604  0.397 |
## egoism              |     3.287 |  0.060  0.681  0.073 | -0.026  0.338  0.013 |
## employment          |     5.648 | -0.137  2.621  0.164 |  0.215 17.555  0.408 |
## finances            |     3.576 | -0.237  2.790  0.276 | -0.206  5.690  0.209 |
## war                 |     1.025 |  0.217  2.169  0.749 | -0.075  0.694  0.089 |
##                      Dim.3    ctr   cos2  
## money                0.101 16.884  0.328 |
## future              -0.053  7.568  0.064 |
## unemployment        -0.004  0.046  0.000 |
## circumstances       -0.016  0.046  0.001 |
## hard                 0.060  0.845  0.051 |
## economic             0.084  3.280  0.027 |
## egoism               0.179 29.496  0.655 |
## employment          -0.213 30.815  0.398 |
## finances            -0.044  0.469  0.010 |
## war                 -0.098  2.139  0.152 |
## 
## Columns
##                       Iner*1000    Dim.1    ctr   cos2    Dim.2    ctr   cos2  
## unqualified         |    13.146 | -0.209 25.110  0.676 | -0.081 10.082  0.101 |
## cep                 |    10.044 | -0.139 18.297  0.645 |  0.056  8.079  0.105 |
## bepc                |     7.670 |  0.109  6.758  0.312 | -0.028  1.251  0.021 |
## high_school_diploma |    17.732 |  0.274 37.976  0.758 | -0.121 20.099  0.149 |
## university          |    13.468 |  0.231 11.859  0.312 |  0.318 60.488  0.589 |
##                      Dim.3    ctr   cos2  
## unqualified          0.073 14.659  0.081 |
## cep                 -0.018  1.520  0.011 |
## bepc                -0.147 59.874  0.570 |
## high_school_diploma  0.077 14.407  0.059 |
## university           0.094  9.540  0.052 |
## 
## Supplementary rows
##                       Dim.1  cos2   Dim.2  cos2   Dim.3  cos2  
## comfort             | 0.210 0.069 | 0.703 0.775 | 0.071 0.008 |
## disagreement        | 0.146 0.131 | 0.119 0.087 | 0.171 0.180 |
## world               | 0.523 0.876 | 0.143 0.065 | 0.084 0.023 |
## to_live             | 0.308 0.139 | 0.502 0.369 | 0.521 0.397 |
## 
## Supplementary columns
##                        Dim.1   cos2    Dim.2   cos2    Dim.3   cos2  
## thirty              |  0.105  0.138 | -0.060  0.044 | -0.103  0.132 |
## fifty               | -0.017  0.011 |  0.049  0.090 | -0.016  0.009 |
## more_fifty          | -0.177  0.286 | -0.048  0.021 |  0.101  0.093 |

## Ellipses for all the active elements
ellipseCA(res.ca)

## Warning: ggrepel: 5 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## Ellipses around some columns only
ellipseCA(res.ca,ellipse="col",col.col.ell=c(rep("blue",2),rep("transparent",3)),
     invisible=c("row.sup","col.sup"))

10.10 使用FactoMineR包和factoextra包

例子来自factoextra包的使用文档。

library(factoextra)

## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBa

library("FactoMineR")
data(housetasks)
head(housetasks)

##            Wife Alternating Husband Jointly
## Laundry     156          14       2       4
## Main_meal   124          20       5       4
## Dinner       77          11       7      13
## Breakfeast   82          36      15       7
## Tidying      53          11       1      57
## Dishes       32          24       4      53

res.ca <- CA(housetasks, graph=FALSE)

# Biplot of rows and columns
# ++++++++++++++++++++++++++
# Symetric Biplot of rows and columns
fviz_ca_biplot(res.ca)

# Asymetric biplot, use arrows for columns
fviz_ca_biplot(res.ca, map ="rowprincipal",
 arrow = c(FALSE, TRUE))

# Keep only the labels for row points
fviz_ca_biplot(res.ca, label ="row")

# Keep only labels for column points
fviz_ca_biplot(res.ca, label ="col")

# Select the top 7 contributing rows
# And the top 3 columns
fviz_ca_biplot(res.ca,  
               select.row = list(contrib = 7),
               select.col = list(contrib = 3))

# Graph of row variables
# +++++++++++++++++++++
   
# Control automatically the color of row points
   # using the "cos2" or the contributions "contrib"
   # cos2 = the quality of the rows on the factor map
   # Change gradient color
   # Use repel = TRUE to avoid overplotting (slow if many points)
fviz_ca_row(res.ca, col.row = "cos2",
   gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
   repel = TRUE)

# You can also control the transparency 
# of the color by the "cos2" or "contrib"
fviz_ca_row(res.ca, alpha.row="contrib")

# Select and visualize some rows with select.row argument.
 # - Rows with cos2 >= 0.5: select.row = list(cos2 = 0.5)
 # - Top 7 rows according to the cos2: select.row = list(cos2 = 7)
 # - Top 7 contributing rows: select.row = list(contrib = 7)
 # - Select rows by names: select.row = list(name = c("Breakfeast", "Repairs", "Holidays"))
 
 # Example: Select the top 7 contributing rows
fviz_ca_row(res.ca, select.row = list(contrib = 7))

# Graph of column points
# ++++++++++++++++++++++++++++

 
# Control colors using their contributions
fviz_ca_col(res.ca, col.col = "contrib",
   gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"))

# Select columns with select.col argument
   # You can select by contrib, cos2 and name 
   # as previously described for ind
# Select the top 3 contributing columns
fviz_ca_col(res.ca, select.col = list(contrib = 3))