「データ解析のための統計モデリング入門」をPythonで写経 Vol.8 ～ ５章「GLMの尤度比検定と検定の非対称性」 GLMの尤度比検定

５章「GLMの尤度比検定と検定の非対称性」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」５章「GLMの尤度比検定と検定の非対称性」の Python写経活動記録 です。 

５章はちょいとクセの強い内容になっています。

この記事は GLMの尤度比検定の計算手続き に絞って学びます。
統計学的な検定に関する枠組みと批評については触れないようにします。

では書籍を開いて統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは、書籍「データ解析のための統計モデリング入門 一般化線形モデル・階層ベイズモデル・MCMC」（岩波書店、「テキスト」と呼びます）の Python 写経を通じて得た「統計モデリングの楽しさ」をご紹介します。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備・サマリー

準備

■ 記事の範囲
この記事はテキスト５章の以下の節を取り扱います。

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
▶️ サポートサイト

Jupyter Notebook ファイルと同一フォルダ内に「data」フォルダを用意して、data フォルダ配下の章別フォルダにデータファイルを格納しています。

■ ライブラリのインポート
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計計算
import scipy.stats as stats
import statsmodels.api as sm
import statsmodels.formula.api as smf

# Rデータセットの読み込み
import rdata

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib

統計モデリング・サマリー

統計モデリングの概要です。
この記事は２つの統計モデルを取り扱います。

■ 統計モデル
一定モデルと $${\mathtt{x}}$$ モデルの２つのポアソン回帰を用います。
（この２モデルが続いていますね…）

■ モデリング手続き
今回は モデル選定「尤度比検定」にフォーカスするので、モデリングの詳細説明は割愛します。

モデリング

データの確認

３章の例題データを利用します。

■ データの読み込み
data3a.csv ファイルを pandas データフレームの data に読み込みます。

# データの読み込み
data = pd.read_csv('./data/ch03/data3a.csv')
print('data.shape: ', data.shape)
data.head()

【実行結果】
データの個数（標本サイズ）は 100 です。
100 個体の植物に関する仮想実験の観測データです。

【変数の説明】
個体から採れた種子数 y、個体の体サイズ x、肥料を与えたかどうかの情報 f です。
今回の統計モデリングで利用する変数は、種子数 y と体サイズ x です。

◆ ◆ ◆

■ データの概要の確認
散布図で種子数 y と体サイズ x の関係を確認します。

# 散布図でデータの確認
fig, ax = plt.subplots(figsize=(5, 5))
sns.scatterplot(data=data, x='x', y='y', s=70, alpha=0.7);

【実行結果】
相関関係があるのでしょうか…

相関係数を確認します。

# 相関係数
data.corr(numeric_only=True).round(3)

【実行結果】
体サイズ x と種子数 y には弱い正の相関関係があります。

統計モデルをデータに当てはめ

２つのポアソン回帰モデルをフィッティングします。
違いは線形予測子です。

◆ ◆ ◆

１つ目は「一定モデル」です。
線形予測子が切片のみのポアソン回帰です。

🔷 確率分布と確率質量関数
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$（平均種子数）のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1}$$」です。
推定するパラメータは $${\beta_1}$$ の１つなので、パラメータ数 $${k=1}$$ です。
もう一つのモデル「$${\mathtt{x}}$$ モデル」との違いは、パラメータ $${\beta_2 = 0}$$ になっている点です。

🔷 GLMの当てはめ
statsmodels の glm クラスを利用してポアソン回帰を実装します。

# ２つのポアソン回帰モデルの当てはめ

# 共通設定
family=sm.families.Poisson()   # GLMの引数familyに与える確率分布＝ポアソン分布

#  一定モデル：パラメータ数 k=1（傾きβ2=0）
result1_k1 = smf.glm(formula='y ~ 1', data=data, family=family).fit()
result1_k1.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-237.64}$$ です。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

◆ ◆ ◆

２つ目は「$${\mathtt{x}}$$ モデル」です。
線形予測子が切片と体サイズのポアソン回帰です。

🔷 確率分布と確率質量関数
個体 $${i}$$ の種子数 $${y_i}$$ はパラメータ $${\lambda_i}$$ のポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は「対数」、線形予測子は「$${\beta_1 + \beta_2 x_i}$$」です。
変数 $${x_i}$$ は体サイズです。
推定するパラメータは $${\beta_1, \beta_2}$$ の２つなので、パラメータ数 $${k=2}$$ です。
一定モデルとの違いは、パラメータ $${\beta_2 \neq 0}$$ になっている点です。

🔷 GLMの当てはめ

#  xモデル：パラメータ数 k=2（傾きβ2≠0）
result1_k2 = smf.glm(formula='y ~ x', data=data, family=family).fit()
result1_k2.summary()

【実行結果】
当てはまりの良さの指標「最大対数尤度」は $${-235.39}$$ です。
こちらのモデルの方が最大対数尤度が大きいので、「当てはまりの良いモデル」になっています。

平均種子数 $${\lambda_i}$$ の予測式は次のようになります。

予測

２つのモデルの平均種子数 $${\lambda_i}$$ の予測値を可視化します。
テキスト p.97 図 5.2(B) に相当します。

# 100個体分の観測データ、一定モデルとxモデル p.97 図5.2(B)

# x軸の値の設定設定
x_val = np.linspace(data.x.min(), data.x.max(), len(data))
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 5))
# 観測データの散布図の描画
sns.scatterplot(data=data, x='x', y='y', s=70, alpha=0.7, ax=ax)
# xモデルの予測値の赤い曲線の描画
ax.plot(x_val, result1_k2.predict(dict(x=x_val)), color='tab:red',
         label=f'xモデル: 逸脱度={-2 * result1_k2.llf:.1f}')
# 一定モデルの予測値の黒点線の描画
ax.plot(x_val, result1_k1.predict(), color='black', ls='--',
         label=f'一定モデル: 逸脱度={-2 * result1_k1.llf:.1f}')
# 修飾：x軸ラベル、y軸ラベル、凡例（枠の外に出す）
ax.set_xlabel('体サイズ $x_i$', fontsize=12)
ax.set_ylabel('種子数 $y_i$', fontsize=12)
ax.legend(loc='lower center', bbox_to_anchor=(0.5, 1));

【実行結果】
黒点線が一定モデルの予測、赤実線が $${\mathtt{x}}$$ モデルの予測です。

尤度比検定

尤度比検定とは

ChatGPT に相談して作成した「尤度比検定の概要」です。

🔷 尤度比検定

尤度比検定（Likelihood Ratio Test, LRT） は、ネストしている２つの統計モデルを比較して、「より複雑なモデルが、統計的に有意に良い説明力を持つかどうか」を調べる検定方法です。

「ネストしているモデル」とは、単純なモデルが、複雑なモデルの特別な場合（部分モデル）として含まれている関係を意味します。

🔷 統計モデルの構成

• 単純モデル：一部のパラメータに制約（たとえばゼロ）が課されたモデル

• 複雑モデル：制約がなく、すべてのパラメータを自由に推定するモデル

🔷 仮説

• 帰無仮説 $${H_0}$$：単純モデルが適切である
  （＝複雑モデルの追加的なパラメータはすべて 0 である）

• 対立仮説 $${H_1}$$：複雑モデルが適切である
  （＝追加されたパラメータのうち、少なくとも１つは 0 でない）

🔷 検定統計量

単純モデルと複雑モデルの逸脱度の差（または最大対数尤度の差）に基づく検定統計量 $${\Lambda}$$ を用います。

検定統計量 $${\Lambda}$$ は、追加されたパラメータの数を自由度とするカイ二乗分布に近似的に従うとされます（大標本の場合）。
なおテキストは、① パラメトリックブートストラップ法（PB 法）と ② カイ二乗分布を使った近似計算法で、棄却限界値や $${p}$$ 値を求めます。

🔷 解釈

• $${\Lambda}$$ の値が大きく、対応する $${p}$$ 値が有意水準（たとえば 0.05）以下であれば、単純モデルでは不十分であり、複雑モデルの方が統計的に優れていると判断します。

• 一方、$${p}$$ 値が有意水準を上回る場合は、単純モデル・複雑モデルのどちらが優れているとも不十分とも言えません。
  ２つのモデルのどちらが良いかの判断を保留することになります。

💡 ２つの統計モデル・バージョン 💡
今回取り扱う２つの統計モデルを上述の文章に当てはめてみます！

🔷 モデルの構成

• 一定モデル：$${\beta_2=0}$$ の制約が課されたパラメータ数 $${k=1}$$ のモデル

• $${\mathtt{x}}$$ モデル：$${\beta_2 \neq 0}$$ とするパラメータ数 $${k=2}$$ のモデル

🔷 仮説

• 帰無仮説 $${H_0}$$：一定モデルが適切である（$${\beta_2=0}$$）

• 対立仮説 $${H_1}$$：$${\mathtt{x}}$$ モデルが適切である（$${\beta_2 \neq 0}$$）

🔷 検定統計量

一定モデルと $${\mathtt{x}}$$ モデルの逸脱度の差（または最大対数尤度の差）に基づく検定統計量 $${\Lambda}$$ を用います。

検定統計量 $${\Lambda}$$ は、追加されたパラメータの数 $${1}$$を自由度とするカイ二乗分布に近似的に従うとされます（大標本の場合）。
なおテキストは、① パラメトリックブートストラップ法（PB 法）と ② カイ二乗分布を使った近似計算法で、棄却限界値や $${p}$$ 値を求めます。

🔷 解釈

• $${\Lambda}$$ の値が大きく、対応する $${p}$$ 値が有意水準（たとえば 0.05）以下であれば、一定モデルでは不十分であり、$${\mathtt{x}}$$ モデルの方が統計的に優れていると判断します。

• 一方、$${p}$$ 値が有意水準を上回る場合は、一定モデル・$${\mathtt{x}}$$ モデルのどちらが優れているとも不十分とも言えません。
  ２つのモデルのどちらが良いかの判断を保留することになります。

２つのモデルの検定統計量

■ 逸脱度
検定統計量と関連する逸脱度を振り返ります。
一定モデルと $${\mathtt{x}}$$ モデルの逸脱度等を確認します。
テキスト p.98 表 5.1（３章p.77 表 4.3）に相当します。

# 一定モデルとxモデルの対数尤度・逸脱度・AIC p.98 表5.1

# フルモデルの最大対数尤度の算出
max_logL_full = sum(stats.poisson.logpmf(k=data.y, mu=data.y))

# 3モデルの各列の値の設定・算出 
model_names = ['一定', 'x', 'フル']          # モデル名
ks = np.array([1, 2, 100])                  # k:変数の数
max_logLs = np.array(
    [result1_k1.llf, result1_k2.llf, max_logL_full]) # 最大対数尤度
deviances = -2 * max_logLs                  # 逸脱度
resid_deviances = deviances - deviances[2]  # 残差逸脱度
AICs = -2 * (max_logLs - ks)                # AIC

# データフレーム化
stats_df = pd.DataFrame({
    'k': ks, '最大対数尤度': max_logLs, '逸脱度': deviances, 
    '残差逸脱度': resid_deviances, 'AIC': AICs},
    index=pd.Series(model_names, name='モデル'))
stats_df.round(1)

【実行結果】
逸脱度は当てはまりの悪さを表すので、一定モデルの方が当てはまりが悪いです。
一定モデルと $${\mathtt{x}}$$ モデルの逸脱度の差はおよそ $${4.5}$$ です。

■ 尤度比検定の検定統計量 $${\Lambda}$$
尤度比検定の検定統計量を求めます。

# 検定統計量（逸脱度の差）
stats_val = stats_df.loc['一定', '逸脱度'] - stats_df.loc['x', '逸脱度']
print(f'検定統計量 Λ: {stats_val:.2f}')

【実行結果】
検定統計量の実現値は 4.5 です。
先ほど計算した逸脱度の差 4.5 と一致しています。

続いて、「帰無仮説が正しいと仮定したときに検定統計量が従う確率分布」の学びに進みます。
テキストは ①パラメトリックブートストラップ法と ②カイ二乗分布を使った近似計算法 の２つの分布による検定方法を説明しています。
検定統計量が従う確率分布に基づいて $${p}$$ 値や棄却限界値を求め、「検定統計量の実現値 4.5」で帰無仮説を棄却できるかどうか、判断します。

なお有意水準は 5% とします。

パラメトリックブートストラップ法

■ パラメトリックブートストラップ法とは
パラメトリックブートストラップ法（略して PB 法）は、仮定された統計モデルに従って乱数を生成し、検定統計量の分布をシミュレーションによって推定する方法です。

■ テキストのパラメトリックブートストラップ法
ここでは、帰無仮説である一定モデルを真のモデルと仮定し、真のモデルから乱数生成して、検定統計量＝逸脱度の差の分布を推定します。
この分布を用いて $${p}$$ 値と棄却限界値を算出します。

テキストのシミュレーション方法は次のとおりです。

• 推定したパラメータ $${\hat{\beta_0}}$$ を使って、平均パラメータ $${\hat{\lambda} = \exp(\hat{\beta_0})}$$ のポアソン分布乱数を 100 個生成して逸脱度の差を計算します。

• この処理を 1000 回繰り返し実施して、取得した 1000 個の逸脱度の差を検定統計量の分布にします。

◆ ◆ ◆

■ 一定モデルの内容
一定モデルの $${\hat{\beta_0}}$$ と $${\hat{\lambda}}$$ を確認します。

# 帰無仮説である一定モデルで推定された平均種子数の計算 p.103
k1_beta1 = result1_k1.params['Intercept'].round(2)  # ≒ 2.06
print(f'一定モデルのβ0　　　 : {k1_beta1:.2f}')
print(f'一定モデルの平均種子数: {np.exp(k1_beta1):.2f}')

【実行結果】
平均種子数 $${\hat{\lambda} = 7.85}$$ です。

ところで、種子数 y が従うポアソン分布の平均パラメータ＝平均種子数の最尤推定値は 種子数 y の標本平均 $${\overline{\bm Y}}$$ です。

# 平均種子数の最尤推定量は種子数 y の標本平均
print(f"種子数 y の標本平均: {data['y'].mean()}")

【実行結果】

テキストは、一定モデルの平均種子数の推定値 $${\hat{\lambda}}$$ に種子数 y の標本平均 $${\overline{\bm Y} = 7.83}$$ を使っています。
この記事もテキストに同調して、$${\hat{\lambda}=7.83}$$ を使います。

◆ ◆ ◆

■ １回のお試しシミュレーション
平均パラメータ $${\lambda = 7.83}$$ のポアソン分布で乱数を 100 個生成して、一定モデルと $${\mathtt{x}}$$ モデルの逸脱度の差を算出する練習をします。

ポアソン分布乱数 100 個を生成します。
テキスト p.103 の１つ目の R 実行コード（rpois）に相当します。

# ポアソン乱数生成関数を使って、新のモデルから100個体分のデータを新しく生成 p.103
# テキストと乱数が異なるため、結果はテキストと異なります

# 乱数シードの設定
np.random.seed(3124)
# 観測データの種子数yの標本平均★をパラメータにするポアソン分布乱数を100個生成
data['y_rnd'] = stats.poisson.rvs(mu=data['y'].mean(), size=100)
# 結果の確認
data.head()

【実行結果】
y_rnd 列が生成した新しいデータです。

新しいデータに一定モデルと $${\mathtt{x}}$$ モデルを当てはめて、逸脱度の差を求めます。
テキスト p.103 の２つ目の R 実行コード（glm）に相当します。
乱数がテキストと異なるので、結果もテキストと異なります。

# 新データに一定モデル・xモデルをあてはめる p.103
# テキストと乱数が異なるため、結果はテキストと異なります

# 共通設定
family=sm.families.Poisson()   # GLMの引数familyに与える確率分布＝ポアソン分布

# モデルの当てはめ
fit1 = smf.glm(formula='y_rnd ~ 1', data=data, family=family).fit() # 一定モデル
fit2 = smf.glm(formula='y_rnd ~ x', data=data, family=family).fit() # xモデル

# 逸脱度の差の算出
print(f'逸脱度の差: {fit1.deviance - fit2.deviance:.6f}')

【実行結果】
新データの逸脱度の差は $${1.92}$$ です。

◆ ◆ ◆

■ PB 法による逸脱度の差の分布データの作成
尤度比検定のためのデータ生成を行います。
まずテキスト p.104 の PB 法実行関数に似せた pb 関数を定義します。

# PB法実行関数pb()の定義 p.104

def pb(df, n_bootstrap=1000, seed=123):

    ## 設定
    n_sample = len(df)              # データ数
    y_mean = df.y.mean()            # 標本平均
    np.random.seed(seed)            # 乱数シード
    family = sm.families.Poisson()  # GLMの引数familyの確率分布＝ポアソン分布
    stats_list = []                 # 逸脱度の差を格納するリストの初期化

    ## パラメトリックブートストラップ法で逸脱度の差を取得（bootstrap回繰り返し処理）
    for _ in range(n_bootstrap):
        # λ＝標本平均のポアソン乱数の取得
        df['y_rnd'] = stats.poisson.rvs(mu=y_mean, size=n_sample)
        # 一定モデルのあてはめ
        fit1 = smf.glm(formula='y_rnd ~ 1', data=df, family=family).fit()
        # xモデルのあてはめ
        fit2 = smf.glm(formula='y_rnd ~ x', data=df, family=family).fit()
        # 逸脱度の差をリストに格納
        stats_list.append(fit1.deviance - fit2.deviance)
    
    ## 戻り値：逸脱度の差（bootstrap個）
    return pd.DataFrame(stats_list, columns=['dd12'])

【実行結果】なし

繰り返し 1000 回のパラメトリックブートストラップ法を実行します。
テキスト p.104 の R 実行コードに相当します。

# パラメトリックブートストラップ法の実行 p.104
# テキストと乱数が異なるため、結果はテキストと異なります

dd12 = pb(data, n_bootstrap=1000, seed=10)

【実行結果】なし

得られた逸脱度の差の要約統計量を確認しましょう。
テキスト p.105 の R 実行コード（summary）に相当します。
乱数がテキストと異なるので、結果もテキストと異なります。

# 得られた逸脱度の差の要約統計量の表示 p.105
# テキストと乱数が異なるため、結果はテキストと異なります

dd12.describe().T.round(3)

【実行結果】
検定統計量の実現値 4.5 は第３四分位数から最大値の間に位置することが分かりました。

ヒストグラムで逸脱度の差の分布を確認しましょう。
検定統計量の実現値 4.5 もプロットします。
テキスト p.105 図 5.4 に相当します。

# 得られた逸脱度の差の要約統計量の表示 p.105 図5.4
# テキストと乱数が異なるため、結果はテキストと異なります

# 逸脱度の差のヒストグラムの描画
dd12['dd12'].hist(bins=45, edgecolor='white', alpha=0.7, grid=False,
                  figsize=(7, 4))
# 観察された逸脱度の差4.5の垂直線の描画
plt.axvline(4.5, color='tab:red', ls='--',
            label='観察された逸脱度の差 $\Delta D_{1,2}=4.5$')
# 修飾：x軸ラベル、凡例表示
plt.xlabel('一定モデルとxモデルの逸脱度の差 $\Delta D_{1,2}$', fontsize=12)
plt.legend();

【実行結果】
果たして検定統計量の実現値 4.5 で帰無仮説を棄却できるのでしょうか…

◆ ◆ ◆

■ $${p}$$ 値と棄却限界値の算出
検定統計量の実現値 4.5 以上のデータの個数を数えて、全データ数 1000 に対する割合＝ $${p}$$ 値を計算します。
テキスト p.105 の３つ目のR 実行コードに相当します。
乱数がテキストと異なるので、結果もテキストと異なります。

# 観察された逸脱度の差4.5以上になるデータの個数 p.105
# テキストと乱数が異なるため、結果はテキストと異なります

num_rare = sum(dd12.dd12 >= 4.5)
print(f'4.5以上となるデータの個数  　: {num_rare}')
print(f'4.5以上となるデータの確率 p値: {num_rare / len(dd12):.3f}')

【実行結果】
4.5 以上のデータ個数は 37 であり、その割合＝ $${p}$$ 値は 0.037 です。

$${p}$$ 値 $${=0.05}$$ となる逸脱度の差、つまり、逸脱度の差の分布の 95% 点を調べます。
この値が棄却限界値です。
テキスト p.106 の R 実行コードに相当します。
乱数がテキストと異なるので、結果もテキストと異なります。

# P=0.05となる逸脱度の差の値を調べる p.106
# テキストと乱数が異なるため、結果はテキストと異なります

p005 = dd12.dd12.quantile(q=0.95)
print(f'P=0.05となる逸脱度の差: {p005:.6f}')

【実行結果】
棄却限界値は 3.59 です。

◆ ◆ ◆

■ 尤度比検定の結論
検定統計量の実現値である逸脱度の差 4.5 の $${p}$$ 値は 0.038 は有意水準 5% で有意です。
あるいは、検定統計量の実現値である逸脱度の差 4.5 は棄却限界値 3.59 以上なので、棄却域に含まれます。
したがって、一定モデルは棄却され、$${\mathtt{x}}$$ モデルを採択します。

先ほどの逸脱度の差のヒストグラムに棄却限界値を重ねたチャートを描画します。

# 得られた逸脱度の差の要約統計量ｔ棄却限界値の表示
# テキストと乱数が異なるため、結果はテキストと異なります

# 逸脱度の差のヒストグラムの描画
dd12['dd12'].hist(bins=45, edgecolor='white', alpha=0.7, grid=False,
                  figsize=(7, 4))
# 観察された逸脱度の差4.5の垂直線の描画
plt.axvline(4.5, color='tab:red', ls='--',
            label='観察された逸脱度の差 $\Delta D_{1,2}=4.5$')
# 有意水準 5% のときの棄却限界値 3.59 の垂直線の描画
plt.axvline(p005, color='gray', ls='--', label=f'棄却限界値 ${p005:.3}$')
# 修飾：x軸ラベル、凡例表示
plt.xlabel('一定モデルとxモデルの逸脱度の差 $\Delta D_{1,2}$', fontsize=12)
plt.legend();

【実行結果】
グレイの点線が棄却限界値です。
検定統計量の実現値 4.5 が棄却限界値を超えている様子が分かります。

カイ二乗分布を使った近似計算法

■ 概要
検定統計量である逸脱度の差は、一定モデルと $${\mathtt{x}}$$ モデルのパラメータ数の差１を自由度とするカイ二乗分布に近似的に従います。

テキストによると、このカイ二乗分布近似は、標本サイズが大きい場合に有効であり、例題の観測データの標本サイズ 100 程度では $${p}$$ 値が正確でない可能性がある、とのことです。

◆ ◆ ◆

■ $${p}$$ 値と棄却限界値の算出
カイ二乗分布の近似計算法を実践しましょう。
自由度１のカイ二乗分布において、検定統計量の実現値＝逸脱度の差 4.5 になる $${p}$$ 値と棄却限界値を算出します。
テキスト p.107 の２つ目の R 実行コードに相当します。

# カイ二乗分布を用いた尤度比検定（近似計算）の実行 p.107
# anova関数に相当するPythonライブラリを見つけられませんでしたので

# 自由度の算出 ※大きなモデルのパラメータ数－小さなモデルのパラメータ数
df_diff = len(result1_k2.params) - len(result1_k1.params)

# 逸脱度の差の算出 ※小さなモデルの残差逸脱度－大きなモデルの残差逸脱度
deviance_diff = result1_k1.deviance - result1_k2.deviance 

# カイ二乗分布よりp値を算出 ※scipy.stats.chi2利用
# deviance_diffのカイ二乗分布の上側確率を取得
p_value = stats.chi2.sf(x=deviance_diff, df=df_diff)

# カイ二乗分布よりp値を算出 ※scipy.stats.chi2利用
# deviance_diffのカイ二乗分布の上側5%点
c_value = stats.chi2.isf(q=0.05, df=df_diff)

# 結果の表示
print(f'Df: 自由度\t\t= {df_diff}')
print(f'Deviance: 逸脱度の差\t= {deviance_diff:.3f}')
print(f'P(>|Chi|): p値\t\t= {p_value:.3f}') 
print(f'棄却限界値:\t\t= {c_value:.3f}')

【実行結果】
$${p}$$ 値は 0.034、棄却限界値は 3.84 です。

◆ ◆ ◆

■ 尤度比検定の結論
検定統計量の実現値である逸脱度の差 4.5 の $${p}$$ 値は 0.034 は有意水準 5% で有意です。
あるいは、検定統計量の実現値である逸脱度の差 4.5 は棄却限界値 3.84 以上なので、棄却域に含まれます。
したがって、一定モデルは棄却され、$${\mathtt{x}}$$ モデルを採択します。

自由度１のカイ二乗分布と逸脱度の差、棄却限界値を可視化します。

# カイ二乗分布近似による逸脱度の差と棄却限界値の可視化

# 自由度１のカイ二乗分布の確率密度関数の描画
x_val = np.linspace(0, 5, 101)
plt.plot(x_val, stats.chi2.pdf(x_val, df=df_diff),
         label='自由度1のカイ二乗分布')
# 検定統計量：逸脱度の差の垂直点線の描画
plt.axvline(deviance_diff, color='tab:red', ls='--',
            label='検定統計量：逸脱度の差')
# 棄却限界値の垂直点線の描画
plt.axvline(c_value, color='gray', ls='--', label='棄却限界値')
# 修飾
plt.xlabel('逸脱度の差', fontsize=12)
plt.ylabel('確率密度', fontsize=12)
plt.legend(loc='lower center', bbox_to_anchor=(0.5, 1));

【実行結果】
検定統計量：逸脱度の差が棄却限界値を超えている様子が分かります。

尤度比検定の結果の解釈

尤度比検定は、帰無仮説の「単純モデルが適切である」を棄却できるか・棄却できないかのみを取り扱います。
帰無仮説を棄却して対立仮説「複雑モデルが適切である」と結論づけたとして、しかし、複雑モデルの適切さ・良さの程度は分かりません。
テキストの言葉をお借りすると「$${p}$$ 値は効果の大きさを表さない」のです。

まとめ

今回は尤度比検定を学びました。

■ 尤度比検定
尤度比検定（Likelihood Ratio Test, LRT） は、ネストしている２つの統計モデルを比較して、「より複雑なモデルが、統計的に有意に良い説明力を持つかどうか」を調べる検定方法です。

■ 検定統計量
尤度比検定の検定統計量は２つのモデルの逸脱度の差です。

■ 逸脱度の差の分布
今回は次の２つの方法を学びました。
① パラメトリックブートストラップ法
　帰無仮説の単純モデルを真のモデルと仮定して乱数を生成して、逸脱度の差の分布データを取得します。
② カイ二乗分布を使った近似計算法
　標本サイズが大きい場合、逸脱度の差が、２つのモデルのパラメータ数の差を自由度とするカイ二乗分布に近似的に従うことを利用します。

■ 検定結果の解釈
尤度比検定の結果に基づくモデル選択は、単純モデルを棄却して複雑モデルを採択するか、単純モデルを棄却できないかの「二択」で判断します。
選択したモデルの良さの程度や効果の大きさは分かりません。

今回のブログは以上です。

次回は、ロジスティック回帰です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。