カイ二乗検定が有意だった場合におこなう残差分析をPythonで実装する

『Rで学ぶ確率統計学 多変量統計編』を読んでいて、分割表の検定の章で残差分析という手法が出てきました。

カイ二乗検定の説明は市販の統計学の参考書で頻出しているものの、有意な検定結果を得られた後に実施する下位検定（post hoc test）の方法までを述べている本は少ないと思います。

今回はそんな下位検定の一つである残差分析を理解するために、その簡単な理論とともにPythonでの実装コードを紹介します。

そもそも残差分析とは

残差分析はカイ二乗検定で有意な \( \chi^2 \) 値が得られた後に行われる分析方法（下位検定）です。
そもそもカイ二乗検定だけでは、分割表全体に対しての検定はできますが、仮に独立であることがわかってもどのセルがその結果に起因しているのかまではわかりません。
この問題を解決する方法の1つが残差分析です。
残差分析によって、どのセルがカイ二乗検定の有意な検定結果に寄与しているのかを把握できます。

残差分析の理論

残差分析では、カイ二乗検定で使用する観測値 \( O_{ij} \)、期待度数 \( E_{ij} \) に加えて、調整済み標準化残差を用います。
残差分析ではセルごとに調整済み標準化残差を計算し、それを大標本と考えて正規近似して \(p\) 値を計算します。
※カイ二乗検定については市販の統計学書籍やWebの記事に詳しい解説があるためここでは省略します。

調整済み標準化残差

今回の残差分析では、以下で定義される調整済み残差標準化残差をもとに検定を実施します。
$$
\mbox{Adj Residual} = \frac{O_{ij} – E_{ij}}{\sqrt{E_{ij} \cdot (1 – n_{i.}/N)(1 – n_{.j}/N)}}
$$
ここで、各記号は以下とします。

• \( O_{ij} \)： \(i\) 行 \(j\) 列の観測度数
• \( E_{ij} \)： \(i\) 行 \(j\) 列の期待度数
• \(n_{i.}\): \(i\) 行の観測値の合計値
• \(n_{.j}\): \(j\) 列の観測値の合計値
• \(N\): 観測値の合計値

ちなみに、標準残差は以下で定義されます。
$$
\mbox{Std Residual} = \frac{O_{ij} – E_{ij}}{\sqrt{E_{ij}}}
$$
期待度数 \( E_{ij} \) が十分に大きい場合は標準残差でも標準正規分布によく近似しますが、期待度数 \( E_{ij} \) が十分に大きい場合は少ないために調整済み標準残差を使うことが一般的です。

標準正規分布を用いた検定

調整済み残差標準化残差は標準正規分布に近似するため、標準正規分布を用いて検定します。
したがって、調整済み標準化残差が0、すなわち、観測度数\( O_{ij} \) と期待度数 \( E_{ij} \) が等しいことを帰無仮説とした統計的仮説検定をおこないます。

調整済み標準化残差が上記の標準正規分布の薄い青の部分にある場合、帰無仮説が棄却されることになります。

ちなにみ、上記のグラフは以下のコードで描画しています。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# https://pythonforundergradengineers.com/plotting-normal-curve-with-python.html
reject_inf, reject_sup = norm.ppf([0.025, 0.975],0,1)
x = np.arange(reject_inf, reject_sup, 0.001)
x_all = np.arange(-10, 10, 0.001)
y = norm.pdf(x,0,1)
y2 = norm.pdf(x_all,0,1)
# build the plot
fig, ax = plt.subplots(figsize=(9,6))
plt.style.use('fivethirtyeight')
ax.plot(x_all,y2)

ax.fill_between(x,y,0, alpha=0.3, color='b')
ax.fill_between(x_all,y2,0, alpha=0.1)
ax.set_xlim([-4,4])
ax.set_xlabel('# of Standard Deviations Outside the Mean')
ax.set_yticklabels([])
ax.set_title('Normal Gaussian Curve')
plt.show()

残差分析の実装

上記で紹介した調整済み標準残差をセルごとに計算し、標準正規分布をもとに検定を実施します。

import pandas as pd
import numpy as np
from scipy.stats import norm, chi2, chi2_contingency


def residual_analysis(table: pd.DataFrame, p_value: int=0.05):
    """
    クロス集計結果に対して残差分析を実施し、指定したp値以下の組み合わせを取得するメソッド。
    
    Parameters
    -------
    table : pd.DataFrame
        クロス集計結果。インデックス、カラム名を指定すること。
    p_value : int
        p値。
    
    Returns
    -------
    pair list : list
        インデックス、カラム名の組み合わせtupleのlist
    
    """
    
    # numpy.arrayに変換
    np_data = np.array(data)
    
    # カイ二乗検定
    chi_sqared, chi_p_value, df, exp = chi2_contingency(np_data)
    if chi_p_value < p_value:
        print(f'カイ二乗検定：有意水準{p_value}で有意差あり。')
    else:
        print(f'カイ二乗検定：有意水準{p_value}で有意差なし。')
    
    # インデックスとカラム名
    index = data.index
    column = data.columns
    
    # 行数と列数を取得
    row_num, col_num = np_data.shape
    # 合計
    total = np_data.sum()
    # 行と列ごとの合計
    total_by_row = [np_data[i, :].sum() for i in range(row_num)]
    total_by_col = [np_data[:, i].sum() for i in range(col_num)]
    
    # 期待値
    exp = np.array(exp)
    
    pairs = list()
    # 期待値と残差分散を算出
    for i in range(row_num):
        for j in range(col_num):
            # 残差分散
            res_var = (1 - total_by_row[i]/total)*(1 - total_by_col[j]/total)
            # 調整済み標準化残差
            std_res = (np_data[i, j] - exp[i, j])/np.sqrt(exp[i, j] * res_var)
            # 両側検定
            p = norm.sf(x=abs(std_res), loc=0, scale=1)*2
            # p値を下回るペア
            if p <= p_value:
                pairs.append((index[i], column[j]))
    return pairs

上記で定義したメソッドに、引数にクロス集計のpandas.DataFrameと \( p\) 値を指定して実行します。
このメソッドでは、カイ二乗検定を実施後に、残差分析によってどのセルが検定結果に寄与しているのかをtupleで返却します。

# sampleデータ作成
# 参考：https://www.jstage.jst.go.jp/article/jraps/28/2/28_KJ00004980771/_pdf
data = pd.DataFrame(
    [
        [15, 25, 13],
        [8, 18, 20],
        [5, 8, 27]
    ],
    columns=['A', 'B', 'C'],
    index=['年少', '年中', '年長']
)

residual_analysis(table=data, p_value=0.05)

カイ二乗検定：有意水準0.05で有意差あり。
[('年少', 'B'), ('年少', 'C'), ('年長', 'B'), ('年長', 'C')]

最後に

今回は下位検定の一つである残差分析について紹介しました。
今後はそもそものカイ二乗検定だったり、カイ二乗検定が不適切な場合の検定法やその他の下位検定などについてもご紹介できればと思います。

参考

• Chi-Square Test is Statistically Significant: Now What?
• クロス集計表に対する統計分析の手法
• 分散分析の下位に多重検定を置くな