列聯表¶
statsmodels 支援多種分析列聯表的方法，包括評估獨立性、對稱性、同質性的方法，以及處理來自分層母體的表格集合的方法。
這裡描述的方法主要用於雙向表格。多向表格可以使用對數線性模型來分析。statsmodels 目前沒有專門用於對數線性建模的 API，但可以使用 statsmodels.genmod.GLM 中的 Poisson 迴歸來達到此目的。
列聯表是一種多向表格，用於描述資料集，其中每個觀測值都屬於多個變數的每個類別之一。例如，如果有兩個變數，一個有 \(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
在這個範例中，與來自行列獨立的母體的樣本相比，我們在安慰劑/沒有改善和治療/顯著改善單元格中有過多的觀測值，而在安慰劑/顯著改善和治療/沒有改善單元格中觀測值過少。這反映了治療的明顯益處。
如果表格的行和列是未排序的（即，是名義因子），那麼正式評估獨立性的最常見方法是使用 Pearson 的 \(\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
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's Q 檢定，用於檢定相同的二項式比例。 |
另請參閱¶
Scipy 有幾個用於分析列聯表的函數，包括目前 statsmodels 中沒有的費雪精確檢定。