列联表¶
statsmodels 支持多种方法来分析列联表，包括用于评估独立性、对称性、同质性的方法，以及用于处理来自分层人口的表格集合的方法。
这里描述的方法主要用于双向表。多向表可以使用对数线性模型进行分析。statsmodels目前没有专门的对数线性建模API，但可以使用statsmodels.genmod.GLM中的泊松回归来实现这一目的。
列联表是一种多维表格，用于描述一个数据集，其中每个观测值属于多个变量中的每个变量的一个类别。例如，如果有两个变量，一个有\(r\)个水平，另一个有\(c\)个水平，那么我们有一个\(r \times c\)的列联表。该表可以用落在表中给定单元格的观测值数量来描述，例如\(T_{ij}\)是具有第一个变量的水平\(i\)和第二个变量的水平\(j\)的观测值数量。请注意，每个变量必须有有限数量的水平（或类别），这些水平可以是顺序的或无序的。在不同的上下文中，定义列联表轴的变量可能被称为分类变量或因子变量。它们可以是名义的（如果它们的水平是无序的）或有序的（如果它们的水平是有序的）。
列联表的底层总体由一个分布表 \(P_{i, j}\) 描述。\(P\) 的元素是概率，\(P\) 中所有元素的总和为 1。分析列联表的方法使用 \(T\) 中的数据来了解 \(P\) 的属性。
The 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, …。这给出了“Cochran-Armitage趋势检验”。
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"])
对称性和均匀性¶
对称性是指对于每一个\(i\)和\(j\)，\(P_{i, j} = P_{j, i}\)。同质性是指行因子和列因子的边际分布相同，这意味着
请注意，为了使这些属性适用，表格 \(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
The 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’程序测试数据是否与共同的比值比一致。它在下面显示为共同比值比检验。Mantel-Haenszel程序测试这个共同的比值比是否等于一。它在下面显示为比值比=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's Q检验。 |
另请参阅¶
Scipy 提供了几个用于分析列联表的函数， 包括目前不在 statsmodels 中的 Fisher 精确检验。