列联表¶
statsmodels 支持多种方法来分析列联表，包括用于评估独立性、对称性、同质性以及用于处理来自分层总体的一组表格的方法。
此处描述的方法主要用于双向表。可以使用对数线性模型分析多维表。statsmodels 目前还没有专门用于对数线性建模的 API，但是 statsmodels.genmod.GLM 中的泊松回归可用于此目的。
列联表是一个多维表，它描述了一个数据集，其中每个观察值属于几个变量中的每一个的某个类别。例如，如果有两个变量，一个有 \(r\) 个水平，另一个有 \(c\) 个水平，那么我们就有 \(r \times c\) 个列联表。该表可以用落在表中某个单元格的观察次数来描述，例如 \(T_{ij}\) 是具有第一个变量的水平 \(i\) 和第二个变量的水平 \(j\) 的观察次数。请注意，每个变量都必须具有有限数量的水平（或类别），这些水平可以是有序的也可以是无序的。在不同的上下文中，定义列联表轴的变量可能被称为 **分类变量** 或 **因子变量**。它们可以是 **名义的**（如果它们的水平是无序的）或 **序数的**（如果它们的水平是有序的）。
列联表的底层总体用 **分布表** \(P_{i, j}\) 描述。\(P\) 的元素是概率，\(P\) 中所有元素的总和为 1。用于分析列联表的方法使用 \(T\) 中的数据来了解 \(P\) 的属性。
代码 statsmodels.stats.Table 是用于处理列联表的两个最基本的类之一。我们可以直接从任何包含列联表单元格计数的矩形类数组对象创建 Table 对象。
In [1]: import numpy as np
In [2]: import pandas as pd
In [3]: import statsmodels.api as sm
In [4]: df = sm.datasets.get_rdataset("Arthritis", "vcd").data
In [5]: df.fillna({"Improved":"None"}, inplace=True)
In [6]: tab = pd.crosstab(df['Treatment'], df['Improved'])
In [7]: tab = tab.loc[:, ["None", "Some", "Marked"]]
In [8]: table = sm.stats.Table(tab)
或者，我们可以传入原始数据，让 Table 类为我们构建单元格计数数组。
In [9]: data = df[["Treatment", "Improved"]]
In [10]: table = sm.stats.Table.from_data(data)
独立性¶
**独立性** 是行因子和列因子独立发生的属性。**关联** 是缺乏独立性。如果联合分布是独立的，它可以写成行和列边缘分布的外积
我们可以获得观察数据的最佳拟合独立分布，然后查看残差，这些残差识别出最强烈地违反独立性的特定单元格。
In [11]: print(table.table_orig)
Improved Marked None Some
Treatment
Placebo 7 29 7
Treated 21 13 7
In [12]: print(table.fittedvalues)
Improved Marked None Some
Treatment
Placebo 14.333333 21.5 7.166667
Treated 13.666667 20.5 6.833333
In [13]: print(table.resid_pearson)
Improved Marked None Some
Treatment
Placebo -1.936992 1.617492 -0.062257
Treated 1.983673 -1.656473 0.063758
在这个例子中，与来自行和列独立的总体的样本相比，我们在安慰剂/无改善和治疗/明显改善单元格中观察到太多观察值，而在安慰剂/明显改善和治疗/无改善单元格中观察到太少观察值。这反映了治疗的明显益处。
如果表的行和列是无序的（即为名义因子），那么正式评估独立性的最常见方法是使用皮尔逊的 \(\chi^2\) 统计量。查看单元格对 \(\chi^2\) 统计量的贡献通常很有用，以了解依赖性的证据来自何处。
In [14]: rslt = table.test_nominal_association()
In [15]: print(rslt.pvalue)
0.0014626434089526352
In [16]: print(table.chi2_contribs)
Improved Marked None Some
Treatment
Placebo 3.751938 2.616279 0.003876
Treated 3.934959 2.743902 0.004065
对于具有有序行和列因子的表，我们可以使用 **线性对线性** 关联检验来获得更多关于尊重排序的备择假设的效力。线性对线性关联检验的检验统计量为
其中 \(r_i\) 和 \(c_j\) 是行和列得分。通常，这些得分被设置为序列 0、1、…。这将给出“科克伦-阿米塔奇趋势检验”。
In [17]: rslt = table.test_ordinal_association()
In [18]: print(rslt.pvalue)
0.023644578093923983
我们可以通过构建一系列 \(2\times 2\) 表格并计算它们的优势比来评估 \(r\times x\) 表格中的关联性。有两种方法可以做到这一点。**局部优势比** 从相邻的行和列类别构建 \(2\times 2\) 表格。
In [19]: print(table.local_oddsratios)
Improved Marked None Some
Treatment
Placebo 0.149425 2.230769 NaN
Treated NaN NaN NaN
In [20]: taloc = sm.stats.Table2x2(np.asarray([[7, 29], [21, 13]]))
In [21]: print(taloc.oddsratio)
0.14942528735632185
In [22]: taloc = sm.stats.Table2x2(np.asarray([[29, 7], [13, 7]]))
In [23]: print(taloc.oddsratio)
2.230769230769231
**累积优势比** 通过在每个可能的点对行和列因子进行二分法来构建 \(2\times 2\) 表格。
In [24]: print(table.cumulative_oddsratios)
Improved Marked None Some
Treatment
Placebo 0.185185 1.058824 NaN
Treated NaN NaN NaN
In [25]: tab1 = np.asarray([[7, 29 + 7], [21, 13 + 7]])
In [26]: tacum = sm.stats.Table2x2(tab1)
In [27]: print(tacum.oddsratio)
0.18518518518518517
In [28]: tab1 = np.asarray([[7 + 29, 7], [21 + 13, 7]])
In [29]: tacum = sm.stats.Table2x2(tab1)
In [30]: print(tacum.oddsratio)
1.0588235294117647
马赛克图是一种图形方法，用于非正式地评估双向表中的依赖性。
In [31]: from statsmodels.graphics.mosaicplot import mosaic
In [32]: fig, _ = mosaic(data, index=["Treatment", "Improved"])
对称性和同质性¶
**对称性** 是 \(P_{i, j} = P_{j, i}\) 对于每个 \(i\) 和 \(j\) 的属性。**同质性** 是行因子和列因子的边缘分布相同的属性，这意味着
请注意，为了使这些属性适用，表 \(P\)（和 \(T\)）必须是正方形的，并且行和列类别必须相同，并且必须以相同的顺序出现。
为了说明这一点，我们加载一个数据集，创建一个列联表，并计算行和列边缘。代码 Table 类包含用于分析 \(r \times c\) 列联表的方法。下面加载的数据集包含对人们左右眼视力评估的结果。我们首先加载数据并创建一个列联表。
In [33]: df = sm.datasets.get_rdataset("VisualAcuity", "vcd").data
In [34]: df = df.loc[df.gender == "female", :]
In [35]: tab = df.set_index(['left', 'right'])
In [36]: del tab["gender"]
In [37]: tab = tab.unstack()
In [38]: tab.columns = tab.columns.get_level_values(1)
In [39]: print(tab)
right 1 2 3 4
left
1 1520 234 117 36
2 266 1512 362 82
3 124 432 1772 179
4 66 78 205 492
接下来，我们从列联表创建一个 SquareTable 对象。
In [40]: sqtab = sm.stats.SquareTable(tab)
In [41]: row, col = sqtab.marginal_probabilities
In [42]: print(row)
right
1 0.255049
2 0.297178
3 0.335295
4 0.112478
dtype: float64
In [43]: print(col)
right
1 0.264277
2 0.301725
3 0.328474
4 0.105524
dtype: float64
代码 summary 方法打印对称性和同质性检验过程的结果。
In [44]: print(sqtab.summary())
Statistic P-value DF
--------------------------------
Symmetry 19.107 0.004 6
Homogeneity 11.957 0.008 3
--------------------------------
如果我们在一个名为 data 的数据框中具有单个案例记录，我们也可以通过使用 SquareTable.from_data 类方法传入原始数据来执行相同的分析。
sqtab = sm.stats.SquareTable.from_data(data[['left', 'right']])
print(sqtab.summary())
单个 2x2 表格¶
代码 sm.stats.Table2x2 类中提供了几种用于处理单个 2x2 表格的方法。代码 summary 方法显示了表格行和列之间关联性的几种度量。
In [45]: table = np.asarray([[35, 21], [25, 58]])
In [46]: t22 = sm.stats.Table2x2(table)
In [47]: print(t22.summary())
Estimate SE LCB UCB p-value
-------------------------------------------------
Odds ratio 3.867 1.890 7.912 0.000
Log odds ratio 1.352 0.365 0.636 2.068 0.000
Risk ratio 2.075 1.411 3.051 0.000
Log risk ratio 0.730 0.197 0.345 1.115 0.000
-------------------------------------------------
请注意，风险比是不对称的，因此如果分析转置后的表格，将会获得不同的结果。
In [48]: table = np.asarray([[35, 21], [25, 58]])
In [49]: t22 = sm.stats.Table2x2(table.T)
In [50]: print(t22.summary())
Estimate SE LCB UCB p-value
-------------------------------------------------
Odds ratio 3.867 1.890 7.912 0.000
Log odds ratio 1.352 0.365 0.636 2.068 0.000
Risk ratio 2.194 1.436 3.354 0.000
Log risk ratio 0.786 0.216 0.362 1.210 0.000
-------------------------------------------------
分层 2x2 表格¶
分层发生在当我们有一组由相同的行和列因子定义的列联表时。在下面的示例中，我们有一组 2x2 表格，反映了中国几个地区的吸烟和肺癌的联合分布。有可能所有表格都具有共同的优势比，即使边缘概率在不同层之间变化。‘Breslow-Day’ 过程测试数据是否与共同的优势比一致。它在下面显示为 常数 OR 检验。Mantel-Haenszel 过程测试此共同的优势比是否等于 1。它在下面显示为 OR=1 检验。还可以估计共同的优势比和风险比，并获得它们的置信区间。代码 summary 方法显示所有这些结果。可以从类方法和属性中获得单个结果。
In [51]: data = sm.datasets.china_smoking.load_pandas()
In [52]: mat = np.asarray(data.data)
In [53]: tables = [np.reshape(x.tolist(), (2, 2)) for x in mat]
In [54]: st = sm.stats.StratifiedTable(tables)
In [55]: print(st.summary())
Estimate LCB UCB
-----------------------------------------
Pooled odds 2.174 1.984 2.383
Pooled log odds 0.777 0.685 0.868
Pooled risk ratio 1.519
Statistic P-value
-----------------------------------
Test of OR=1 280.138 0.000
Test constant OR 5.200 0.636
-----------------------
Number of tables 8
Min n 213
Max n 2900
Avg n 1052
Total n 8419
-----------------------
模块参考¶
|
一个双向列联表。 |
|
可以在 2x2 列联表上执行的分析。 |
|
用于分析平方列联表的方法。 |
|
对一系列 2x2 列联表的分析。 |
|
McNemar 同质性检验。 |
|
Cochran Q 检验，用于检验二项式比例是否相同。 |
参见¶
Scipy 提供了几个用于分析列联表的函数，包括 Fisher 精确检验，目前 statsmodels 中尚未提供。