4.1　列联分析

列联分析是利用列联表（Contingency Table）来研究两个分类变量之间关系的统计方法。以皮尔逊提出的卡方检验为代表列联分析技术在数据分析领域具有非常广泛的应用。

4.1.1　类别数据与列联表

1912年4月15日，豪华巨轮泰坦尼克号在她的处女航中因不幸与冰山相撞而沉入冰冷的大西洋。这场空前的海难令世界震动。根据维基百科所提供的数据，当时船上有908名船员和1316名乘客，共计2224人。事故发生后，共有710名乘客获救，约占总人数的32%。具体伤亡统计数据请见表4-1。

数据是冰冷和残酷的，但它背后的故事却是值得我们深入思考的。比如死亡是否与性别或年龄有关？面对死亡，人是否有高低贵贱之分？面对突如其来的灾难，船员是否忠于职守？这些看似干瘪的数据是否能够反映当时人们的价值取向以及面对死亡的态度？基于本章的学习，相信读者应该有能力解开这些疑问。

我们都知道，统计数据有类别型数据和数值型数据之分。对于类别型数据来说，尽管最终结果也是以数值来表示的，但不同数值所描述的对象特征却是彼此区别的。例如在泰坦尼克号伤亡情况分析的例子中，若想讨论死亡是否与性别有关，就可以把成年人群体分成男性和女性两类，并用1来表示男性，而用0来表示女性。如果想研究死亡是否与舱位有关，又可将乘客分位头等舱乘客（用1表示）、二等舱乘客（用2表示）和三等舱乘客（用3表示）。显然，对于上述问题的分析，是在以统计数据的汇总分类为基础的。而且为了便于后续分析工作的开展，选择一种有效的方式来对数据进行组织也很有必要，这种有效的方式就是所谓的列联表。

列联表是由两个以上的变量进行交叉分类的频数分布表。例如在研究泰坦尼克号中乘客的舱位与其最终是否获救之间关系的问题上，可以建立如表4-2所示的列联表。表中的行（row）是生存变量，它被划分成两类：获救或者罹难；表中的列（column）是舱位变量，它被划分为三类，即头等舱、二等舱和三等舱。因此表4-2是一个2×3的列联表。表中的每个数据，都反映着来自生成变量和舱位变量两个方面的信息。列联表是进行列联分析的基础，下一小节我们就将通过实际的用例来介绍基于χ²检验的列联分析方法。

4.1.2　皮尔逊（Pearson）的卡方检验

皮尔逊的卡方统计量（Pearson's Chi-Square Statistic），或者简写成χ²统计量，是由皮尔逊于1899年提出的用于检验实际分布与理论分布配合程度，即配合度检验的统计量。皮尔逊认为，不管理论分布构造得如何好，它与实际分布之间总存在着或多或少的差异。这些差异可能是由于观察次数不充分、随机误差太大而引起的，也可能是因所选配的理论分布本身与实际分布有实质性差异而导致的。要甄别导致差异的原因，还需要用一种方法来检验。为此，皮尔逊提出了著名的χ²统计量，用来检验实际值的分布数列与理论值数列是否在合理范围内相符合，换句话说，χ²统计量可被用以测定观察值与期望值之间的差异显著性。卡方检验提出后得到了广泛的应用，在现代统计理论中占有重要地位。

卡方统计量由各项实际观测数值与理论分布数值之差的平方除以理论数值，然后再求和而得出的。若用f_o表示观察值频数，用f_e表示期望值频数，χ²统计量可以写为

作为一种统计方法，χ²检验主要用于对两个定类变量之间关系的分析。对定类变量进行分析，一般是把检验问题进行转化，通过考察频数与其期望频数之间的吻合程度，达到检验目的。χ²检验还依赖于χ²分布的自由度，自由度定义为类别数量与限制数量之差，具体的计算我们在后续结合例子加以说明。

假设有一枚骰子，投掷120次并记录其结果如表4-3所示，请问该枚骰子是否是无偏的？

首先提出原假设和备择假设。

H₀：骰子是无偏的，即所有投掷结果出现的可能性大致是均等的。

H₁：原假设是错误的，即某些投掷结果的可能性较其他结果更大。

在原假设基础上，可以得到期望的投掷结果如表4-4所示。

据此可以计算χ²统计量如下

因此P值应该是Pr（X²≥7.90），其中。显然，骰子掷出后可能的结果有6种，而在我们的例子中限制条件的数量只有1个，即所有观察频数之和等于120，所以自由度为5。查询统计表可知，P值将介于0.1和0.25之间，因此在5%的水平下我们无法拒绝原假设。

在R中，自然无须这样烦琐的计算，chisq.test（）函数将执行χ²检验。该函数以记录观测频数的向量为输入，而且默认情况下以均等占比为原假设。对于掷骰子的例子，即假设掷出每个面的比例都是1/6。下面给出示例代码

因此得到P值为0.1616。我们不能得出骰子是有偏的这个结论。注意如果要修改默认为均等占比的原假设，可以通过调整函数中的参数p来实现，具体细节请参阅R的帮助文档，这里不再赘述。

掷骰子的例子其实只是让读者体验了一下χ²检验的基本思想，引入列联表之后将有能力处理更为复杂的例子。注意掷骰子的例子中所给出的表格并不是列联表，因为它并未涉及多分类变量之间的交叉。

下面就来看一个基于列联表的例子。众所周知，妇女怀孕期间饮酒或抽烟将会对胎儿造成不良影响。有人认为饮酒和吸烟之间存在某种联系，例如通常酗酒的人都有抽烟的嗜好。理解两者之间的关系对于研究孕妇相关行为给胎儿可能带来的影响十分重要。1984年，研究人员对452名母亲进行了调查，根据她们在得知自己怀孕前的酒精和烟草摄入量得出了如表4-5所示的列联表。请问饮酒和吸烟之间是否有关联？

列联表是两个因素（变量）从横向和纵向交叉而成的，因此以列联表为基础的假设检验中，原假设H₀通常为两个因素之间是没有联系的，即彼此独立的。相应地，备择假设H₁为原假设为假。在掷骰子的例子中，H₀确定了每个可能输出的概率，彼时H₀仅仅指定了概率之间的关系。对于饮酒和吸烟关系的例子，我们可以提出下列原假设和备择假设。

H₀：吸烟和饮酒之间没有关系，即两者彼此独立。

H₁：原假设是错误的。

回想一下概率论中的有关结论，即如果事件A和B彼此独立，则当且仅当Pr（A∩B）=Pr（A）×Pr（B）。所以如果原假设为真，那么对于本例而言必然有“酒精日均摄入超过1oz并且尼古丁日均摄入超过16mg的概率”就等于Pr（酒精日均摄入≥1oz）×Pr（尼古丁日均摄入≥16mg）。

这个原理也为我们计算期望值列联表提供了依据，相应的期望值就等于行和与列和之积再除以表中数据总和。例如，从表4-5中可知，有90名母亲日均饮酒量超过1oz，在原假设基础上，则其中应该有83/452的人日均尼古丁摄入量超过16mg。因此相应的期望值就应该是90×83/452=16.53。按照此方法，最终可以得出期望值的列联数据如表4-6所示。由此可以算得χ²统计量如下

再来考虑一下用于检验的χ²分布的自由度，对于列联表而言，一个通常的计算公式为

df=（r-1）×（c-1）=rc-（r+c-1）

其中，r表示行数，c表示列数，所以rc就是表中所给出的类别总数。r同时给出了r个限制条件，列数c同时给出了c个限制条件。但总行和＝总列和＝表中数值总和，所以在计算由行限制与列限制给出的限制条件数目时，有一个重复计算，我们应该将其减去。最终限制数量为r+c-1。针对当前所讨论的问题，自由度为

df=（4-1）×（3-1）=6

查表可知的临界值，由此得到的P值小于0.001。计算P值的代码如下。

整个χ²检验执行过程的R代码如下。

由此，便可以果断地拒绝原假设，并推得结论：饮酒与吸烟之间确实有联系。

4.1.3　列联分析应用条件

在卡方检验中使用χ²分布来获取P值，其实隐含地使用了一个条件，即用正态分布来近似二项分布。为了保证卡方检验的有效性，下列执行条件应当予以满足

（1）每个单元格中的数据都是确切的频数（而非占比）。

（2）类别不可相互交织。

（3）所有的期望频数应当都不小于1。

（4）至少80%的期望频数都应该不小于5。

如果上述条件无法都满足，我们就不得不通过合并单元格的方法来满足这些条件。但合并单元格的做法也会令自由度下降，进而削弱检验的效力。

来看一个研究铝元素摄入与阿尔兹海默病之间关系的例子。研究人员选择了一组阿尔兹海默病患者。作为对照实验，又选择了一组没有患阿尔兹海默病的人，但其他方面与实验组中的病患非常形似。参与实验的对象，他们的含铝抗酸剂使用情况如表4-7所示。

下面的代码对上述数据执行了以皮尔逊卡方检验为基础的列联分析，其中将chisq.test（）函数的输出赋给了一个对象，即a.by.a.test。而且我们用一个括号来把赋值语句括了起来，如果不这样做，程序将仅会对函数的结果进行储存但并不会将其输出。

下面想检查一下每个期望频数的大小，于是在R中输入下面的代码。

易见8个期望频数中有两个都小于5，因此得到的P值可能不是十分可靠。尽管它也比较小，但是在5%的显著水平下，并不显著。

为了说明这个问题，可以通过多种方法来陈述这个问题。其中一种方法就是使用模拟的方法来获得一个更加精确的P值，例如

将参数simulate.p.value的值置为TRUE时，就表示要通过蒙特卡洛（Monte Carlo）模拟的方法来计算P值。具体来说，这个过程会产生2000个随机表的独立数据，并以此来评估观察值表在原假设前提下的极端性。通过调整参数B的值（缺省情况下为2000），可以改变蒙特卡洛模拟的重复量。

另外一种可以把问题陈述清楚的方法是对数据进行重新分类。可能更有问题的观察值位于那些表示使用了中等剂量抗酸剂的单元格，所以将中等和高等两类数据进行合并是比较合理的做法。于是便得到了如表4-8所示的结果。

然后再以此为基础执行卡方检验，于是可得下面的结果。

可见所有的期望频数都已经不再小于5了，此时给出的P值为0.037 38，因此可以得出抗酸剂的使用和阿尔茨海默病之间存在某种联系。当然，将低等和中等两类数据进行合并也比较合理，读者不妨尝试这种做法，然后再观察一下其对最终结果的影响。

4.1.4　费希尔（Fisher）的确切检验

如果在2×2的列联表中观察值太小，χ²检验因近似程度较差，易导致分析的偏性（尤其是当所得概率接近检验水准时）。1934年，统计学家费希尔提出了一种新的检验方法，即费希尔确切检验（Fisher's exact test），这是一种专门用来对2×2的列联表进行检验的方法。该方法不属于χ²检验的范畴，但可作为2×2表格的χ²检验的补充。

假设在一项有47名学生参与的调查中，研究人员试图检验性别与惯用左手还是右手之间是否存在联系。调查数据见表4-9，从中易见女生中惯用左手的比例为1/31=0.032，该值小于男生中惯用左手的比例4/16=0.25。为了检验这两个比例是否有显著的不同，自然会想到使用之前介绍的χ²检验，但是四个单元格中有两个所包含的期望值都小于5。根据之前的讨论，我们知道此时应该χ²检验并不明智。

费希尔确切检验假设仅知道2×2列联表中边界上的加和值，但对表中的详细数据一无所知。此时可以得到如表4-10所示的一张残缺表。

在仅知道上面这些边缘加和值的情况下，可以推得总共的可能情况有六种，如表4-11所示。这是因为调查数据中左撇子的数量一共只有5个，那么表格中左上角位置的取值就仅可能是0，1，2，3，4，5这几种情况。据此我们就可以推出全部可能的结果。

注意在表4-11中还计算出了每一种可能情况的概率。费希尔确切检验的基础是超几何分布，超几何分布是统计学上的一种离散型概率分布。假设N件产品中有M件次品，不放回的抽检中，抽取n件时得到X=k件次品的概率分布就是超几何分布，它的概率质量函数PMF为

对应到表4-12中，现在有a+b+c+d=N件产品，其中次品有a+b件，进行不放回的抽检，共抽取了a+c件产品，其中得到X=a件次品的概率即为

于是表4-11中的各个概率值可以使用R中的内嵌分布函数来计算，代码如下。

由于超几何概率分布的非对称性，一个双尾的P值并没有被唯一并确切的定义。但在统计分析中，双尾的P值更为常用。一种计算方法是将两个方向上的单尾P值都算出来，然后将其中的较小者乘以2作为双尾P值使用。另外一种方法，也是R中所使用的，就是将输出结果中小于等于观察值概率的所有概率进行加总。在这种方法中，比观察值概率更小的概率值被看成是比远离原假设的观察值更加极端，这也与我们对P值的定义相吻合。例如在当前所讨论的问题中，P值就应该为

在5%的显著水平下，可以拒绝（惯用哪只手与性别无关的）原假设，并认为女生中左撇子的比例要低于男生中左撇子的比例。

在R中，可以用下面的代码来执行费希尔确切检验，易见其中得出的P值与之前算得的结果一致。