分割表¶

statsmodelsは、独立性、対称性、均一性を評価するための方法、および層化母集団からのテーブルの集合を操作するための方法など、分割表を分析するためのさまざまなアプローチをサポートしています。

ここで説明する方法は、主に2元表用です。多元表は、対数線形モデルを使用して分析できます。 statsmodelsには現在、対数線形モデリング専用のAPIはありませんが、statsmodels.genmod.GLMのポアソン回帰をこの目的で使用できます。

分割表は、各観測値が複数の変数のそれぞれについて1つのカテゴリに属するデータセットを記述する多元表です。たとえば、2つの変数があり、一方が\(r\)レベルで、もう一方が\(c\)レベルである場合、\(r \times c\)分割表になります。テーブルは、テーブルの特定のセルに分類される観測値の数で表すことができます。たとえば、\(T_{ij}\)は、最初の変数のレベル\(i\)と2番目の変数のレベル\(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)

独立性¶

**独立性**とは、行と列の因子が独立して発生するという性質です。**関連性**とは、独立性の欠如です。結合分布が独立している場合、行と列の周辺分布の外積として記述できます。

\[P_{ij} = \sum_k P_{ij} \cdot \sum_k P_{kj} \quad \text{for all} \quad i, j\]

観測データに最適な独立分布を取得し、独立性に最も強く違反している特定のセルを識別する残差を表示できます。

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

順序付けられた行と列の因子を持つテーブルの場合、**線形対線形**関連検定を使用して、順序付けを尊重する代替仮説に対してより大きな検出力を持つことができます。線形対線形関連検定の検定統計量は、

\[\sum_k r_i c_j T_{ij}\]

ここで、\(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つの方法があります。**ローカルオッズ比**は、隣接する行と列のカテゴリから\(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

モザイクプロットは、2元表の依存性を非公式に評価するためのグラフィカルなアプローチです。

In [31]: from statsmodels.graphics.mosaicplot import mosaic

In [32]: fig, _ = mosaic(data, index=["Treatment", "Improved"])

対称性と均一性¶

**対称性**とは、すべての\(i\)と\(j\)について\(P_{i, j} = P_{j, i}\)であるという性質です。**均一性**とは、行因子の周辺分布と列因子の周辺分布が同一であるという性質であり、

\[\sum_j P_{ij} = \sum_j P_{ji} \forall 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テーブル¶

個々の2x2テーブルを操作するためのいくつかのメソッドは、sm.stats.Table2x2クラスで提供されています。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」手順は、データが共通のオッズ比と一致するかどうかを検定します。これは、以下にTest of constant ORとして表示されます。Mantel-Haenszel手順は、この共通のオッズ比が1に等しいかどうかを検定します。これは、以下にTest of 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  
-----------------------

モジュールリファレンス¶

`Table`(table[, shift_zeros])	2元分割表。
`Table2x2`(table[, shift_zeros])	2x2分割表に対して実行できる分析。
`SquareTable`(table[, shift_zeros])	正方分割表を分析するためのメソッド。
`StratifiedTable`(tables[, shift_zeros])	2x2分割表の集合に対する分析。
`mcnemar`(table[, exact, correction])	同質性のマクネマー検定。
`cochrans_q`(x[, return_object])	同一の二項比率に対するコクランのQ検定。