「データ解析のための統計モデリング入門」をPythonで写経 Vol.9 ～ ６章「GLMの応用範囲を広げる」①二項分布とロジスティック回帰

６章「GLMの応用範囲を広げる」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」６章「GLMの応用範囲を広げる」の Python写経活動記録 です。 

この記事は 一般化線形モデル（GLM）の一種「ロジスティック回帰」を実践 します。
二項分布とロジットリンク関数を用います。

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

はじめに

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

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

準備

準備

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

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

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

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

# インポート

# 数値計算
import numpy as np
import pandas as pd
from scipy.special import expit, logit, comb

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

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

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

この記事で扱う統計モデリングの概要です。

■ 統計モデル
ロジスティック回帰と呼ばれる統計モデルです。

■ モデリング手続き

ロジスティック回帰（体サイズ＋施肥処理モデル）

データの確認

データを読み込み、データの外観を眺めてから、統計モデリングのためのデータの特徴確認を行います。

■ CSVファイルの読み込み
今回の統計モデルは前回記事と同じデータに当てはめします。
data4a.csv ファイルを pandas データフレームの data に読み込みます。

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

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

【変数の説明】
個体からの観察種子数 N と生存種子数 y、個体の体サイズ x、肥料を与えたかどうかの情報 f です。
今回の統計モデリングではすべての変数を利用します。
目的変数は生存種子数 y です。

生存種子数 y を観察種子数 N で割った $${q = y/N}$$ を生存確率 q と呼びます。

◆ ◆ ◆

■ データの確認
基本的な統計量やチャートでデータを概観します。

① 要約統計量の表示
量的変数の要約統計量を表示します。

# 基本統計量 p.117
data.describe().round(3)

【実行結果】
上からデータの個数、平均、標準偏差、最小値、第１四分位数、中央値、第三四分位数、最大値です。
標準偏差は不偏分散の標準偏差です。

② 標本分散の表示
要約統計量に含まれない分散を確認します。

# 標本分散
data.var(numeric_only=True).rename('var').to_frame().T.round(3)

【実行結果】
生存種子数の標本分散は標本平均 5.08 よりもやや大きな値です。

③ 質的変数の頻度の表示
施肥処理 f の要素とその頻度を確認します

# 変数fの要素の個数 p.117
data.f.value_counts().to_frame()

【実行結果】
施肥なし $${\mathtt{C}}$$、施肥あり $${\mathtt{T}}$$ の両方が 50 ずつ含まれています。

④ 相関係数の表示
体サイズ x と 生存種子数 y の相関係数を施肥処理あり／なしを交えて確認します。

# 相関係数
corrs = [data[['y', 'x']].corr().iloc[0, 1], 
         data[data.f=='C'][['y', 'x']].corr().iloc[0, 1],
         data[data.f=='T'][['y', 'x']].corr().iloc[0, 1]]
pd.DataFrame(
    corrs, index=['全体', '施肥なし', '施肥あり'], columns=['x,yの相関係数']
).round(3)

【実行結果】
施肥処理別の体サイズ x と生存種子数 y の相関は、施肥なしの方がやや相関係数が大きいです。

⑤ 散布図の描画
体サイズ x と生存種子数 y の散布図を施肥処理別に描画します。
書籍 p.117 図 6.2 に相当します。
seaborn の scatterplot() を利用します。

# 体サイズxと生存種子数yの散布図の描画 p.117 図6.2

# 散布図の描画
sns.scatterplot(
    data=data, x='x', y='y', hue='f', s=70, palette=['tab:blue', 'tab:red'],
    alpha=0.7)
# 修飾
plt.xlabel('植物の体サイズ $x_i$', fontsize=14)
plt.ylabel('生存種子数 $y_i$', fontsize=14)
plt.legend(title='施肥処理 f');

【実行結果】
次の２つの傾向を読み取れました。
・体サイズが大きいほど生存種子数が大きい
・体サイズが同じの場合に施肥ありの方が生存種子数が大きい

⑥ 生存種子数のヒストグラムの描画
施肥処理別に生存種子数 y のヒストグラムを確認します。

# 生存種子数yのヒストグラムの描画

# 設定
bins = np.arange(-0.5, 9)
palette=['tab:blue', 'tab:red']

# 施肥処理=Cの場合の生存種子数のヒストグラムの描画
plt.hist(data[data['f']=='C']['y'], bins=bins, edgecolor='white', alpha=0.5,
         label='C：施肥なし')
# 施肥処理=Cの場合の生存種子数のヒストグラムの描画
plt.hist(data[data['f']=='T']['y'], bins=bins, histtype='step', color='tab:red',
         linewidth=2, label='T：施肥あり')
# 修飾
plt.xlabel('生存種子数 $y$', fontsize=14)
plt.ylabel('頻度', fontsize=14)
plt.legend(title='施肥処理');

【実行結果】
施肥なし $${\mathtt{C}}$$ の場合、生存種子数 y は平坦な感じがします。
施肥あり $${\mathtt{T}}$$ の場合、生存種子数 8 がグンと大きいです

⑦ 施肥処理 f と生存種子数 y の箱ひげ図の描画
せっかくなので箱ひげ図も見ます。
seaborn の箱ひげ図 boxplot() に スウォームプロット swarmplot() を重ねて、分布の様子も確認します。

# 施肥処理fと生存種子数yの箱ひげ図の描画

# 設定
palette=['tab:blue', 'tab:red']
# 箱ひげ図の描画
sns.boxplot(data=data, x='f', y='y', hue='f', fill=False, palette=palette)
# スウォームプロットの描画
sns.swarmplot(data=data, x='f', y='y', hue='f', palette=palette)
# 修飾
plt.xlabel('施肥処理 $f_i$', fontsize=14)
plt.ylabel('生存種子数 $y_i$', fontsize=14);

【実行結果】
施肥なし $${\mathtt{C}}$$ の生存種子数の中央値は約５。
施肥あり $${\mathtt{T}}$$ の生存種子数の中央値は約６。
施肥ありのほうが生存種子数が多くなる傾向が見えます。

◆ ◆ ◆

■ データの特徴まとめ
データの特徴を整理します。

データのばらつきは「確率分布」で表現します。
上述のデータの特徴を表現できる確率分布の候補は「二項分布」です。
前回記事までのポアソン分布は上限なしのカウントデータに適していました。
今回のデータは上限ありのカウントデータなのでポアソン分布ではなく、二項分布を選択します。

統計モデルをデータに当てはめ（モデルの理解）

今回は二項分布・ロジットリンク関数を用いるロジスティック回帰の GLM を例題データ＝観測データに当てはめます。
モデル名は「$${\mathtt{x + f}}$$ モデル」です。

◆ ◆ ◆

■ 確率分布、リンク関数、線形予測子
GLMの３要素である「確率分布」、「リンク関数」、「線形予測子」を定義します。

🔷 確率分布と確率質量関数
生存種子数 $${y_i}$$ は観察種子数 $${N_i}$$、生存確率 $${q_i}$$ の二項分布に従います。

$${\displaystyle \binom{N_i}{y_i}}$$ は二項係数です。組み合わせ $${{}_{N_i} \text{C}_{y_i}}$$ と等しいです。

◆

二項分布を深堀りします。

📶 二項分布の概要
二項分布は確率変数がとびとびの値（離散値）をとる離散型確率分布です。
確率変数を$${y}$$、試行回数を $${N}$$、成功確率を $${q}$$ とすると、二項分布の確率質量関数は次の式で表されます。

📶 二項分布のプロフィール

📶 Python で二項分布の確率質量関数を表現
① 確率質量関数の数式を関数化
$${y=3, N=10, q=0.5}$$ で関数の出力をテストします。

# 二項分布の確率質量関数の実装

# 二項分布の確率質量関数の定義 ※comb はscipyの組み合わせ_C_の関数
bimon_pmf = lambda y, N, q: comb(N, y) * q**y * (1-q)**(N-y)

# テスト
y, N, q = 3, 10, 0.5
print(f'y={y}, N={N}, q={q} の二項分布確率 p(y|N,q) = {bimon_pmf(y, N, q):.4f}')

【実行結果】

② scipy.stats の 二項分布クラス利用
binom.pmf() で確率質量関数を算出できます。

# scipyで計算
y, N, q = 3, 10, 0.5
print(f'scipy.statsで計算 p(y|N,q) = {stats.binom.pmf(y, N, q):.4f}')

【実行結果】

📶 二項分布の確率質量関数の可視化
テキスト p.119 図 6.3 に相当します。
確率変数の範囲 $${y=\{0, 1, \cdots, 8\}}$$、パラメータ $${N=8}$$、$${q = \{0.1, 0.3, 0.8\}}$$ のケースで描画します。

# 二項分布の図示 p.119 図6.3

## 設定
N = 8                 # 種子数N
qs = [0.1, 0.3, 0.8]  # 確率qのパターン
ys = range(0, 9)      # x軸の生存種子数yの値

## 描画
# 生起確率qごとに二項分布の確率算出と折れ線グラフ描画を繰り返し処理
for q in qs:
    # 二項分布の確率算出と折れ線グラフの描画 ※scipyのbinomを利用
    plt.plot(ys, stats.binom.pmf(k=ys, n=N, p=q), '-o', label=f'$q=${q}')
# 修飾：x軸ラベル、y軸ラベル、凡例表示
plt.xlabel('生存種子数 $y_i$', fontsize=14)
plt.ylabel(f'確率 $p(y_i \mid N={N}, q)$', fontsize=14)
plt.legend(loc='lower center', ncol=3, bbox_to_anchor=(0.5, 1));

【実行結果】
分布の形状は峰が一つです。
$${q}$$ の値が大きくなるにつれて、峰が上側（右側）に移動し、$${q(1-q)}$$ に対応して尖り具合が変化します。

📶 観測データと二項分布の特徴
観測データの特徴は二項分布と似ています。

二項分布を観測データに当てはめる際の「前提条件」をテキストから引用いたします。

🛸 ちょっと寄り道：確率分布をデータに当てはめ 🛸
scipy の fit 関数で確率分布をデータに当てはめできます。
観測データに二項分布を当てはめして、最尤法で成功確率パラメータ $${p}$$ を推定します。

# 二項分布のパラメータの最尤推定 ※scipyのfitで二項分布に当てはめ
# 結果は一定モデルと同じだと思います。

## 設定
# 種子数N
N = 8
# x軸の生存種子数yの値
y_val = np.arange(0, N + 1)

## 生存種子数yを二項分布に当てはめ（確率パラメータqの推定）
res = stats.fit(stats.binom, data.y, [(N, N), (0, 1)])

## 描画
# 生存種子数yのヒストグラムの描画
plt.hist(data.y, bins=np.arange(-0.5, N+1), density=True, edgecolor='white',
         alpha=0.7)
# 当てはめした二項分布の確率質量関数の折れ線グラフの描画
plt.plot(y_val, stats.binom.pmf(y_val, *res.params), '-o', color='tab:red',
         lw=2, label=f'fitted binom\n$N$=8, $q$={res.params.p:.3f}')
# 修飾
plt.title(f'二項分布の確率パラメータの推定\n最大対数尤度={-res.nllf():.2f}')
plt.xlabel('生存種子数 $y$', fontsize=14)
plt.ylabel('確率', fontsize=14)
plt.legend();

【実行結果】
$${q}$$ の推定値は $${0.635}$$ でした。
この当てはめはロジスティック回帰の一定モデル（切片だけモデル）と同じ結果になります。

◆

🔷 リンク関数と線形予測子
リンク関数は「ロジット」、線形予測子は「$${\beta_1 + \beta_2 x_i + \beta_3 f_i}$$」です。
施肥処理 $${f_i}$$ は説明の便宜上、$${\mathtt{C=0, T=1}}$$ と読み替えます。
線形予測子の説明変数に $${\mathtt{x,\ f}}$$ を用いるので「$${\mathtt{x+f}}$$ モデル」です。

◆

ロジット関数を深堀りします。

📶 ロジット関数の概要
ロジット関数 $${\text{logit}(\cdot)}$$ は次の式で示されます。
$${z_i}$$ は線形予測子を表します。

📶 ロジット関数の可視化
scipy.special の logit 関数を利用してロジット関数を可視化します。

# ロジット関数の曲線

# 設定
qs = np.linspace(0, 1, 101)  # x軸の確率qの値
# ロジット関数の曲線の描画 ※scipyのlogitを利用
plt.plot(qs, logit(qs), label='$z$ = logit($q$) = log $\\frac{q}{1-q}$')
# 修飾
plt.xlabel('確率 $q$', fontsize=14)
plt.ylabel('線形予測子 $z$', fontsize=14)
plt.legend();

【実行結果】
傾きが「急 ⇒ 緩 ⇒ 急」の曲線です。

ロジット関数の逆関数である「ロジスティック関数」を深堀りします。

📶 ロジスティック関数の概要
ロジスティック関数 $${\text{logistic}(\cdot)}$$ は次の式で示されます。
$${z_i}$$ は線形予測子を表します。

線形予測子と成功確率パラメータの関係を表していて、これからお世話になる予感がします。
ロジスティック関数は、機械学習の文脈だと「シグモイド関数」と呼ばれることが多いかもです。

📶 Python でロジスティック関数を表現
① ロジスティック関数の数式を関数化
$${z=2}$$ で関数の出力をテストします。

# ロジスティック関数の実装

# ロジスティック関数の定義
logistic_func = lambda x: 1 / (1 + np.exp(-x))

# テスト
z = 2
print(f'logistic({z}) = {logistic_func(z):.4f}')

【実行結果】

② scipy.special の ロジスティック関数利用
expit 関数を利用します。

# scipyで計算
z = 2
print(f'scipy.statsで計算 logistic({z}) = {expit(z):.4f}')

【実行結果】

📶 ロジスティック関数の可視化
テキスト p.120 図 6.4 に相当します。

# ロジスティック曲線 p.120 図6.4
# ロジスティック関数＝シグモイド関数にはscipyのexpitを利用

# 設定：x軸の線形予測子zの値
zs = np.linspace(-6, 6, 1001)
# ロジスティック曲線の描画
plt.plot(zs, expit(zs), label=r'$q=\frac{1}{1 + exp(-z)}$') # scipy.special.expit
# z=0の垂直線の描画
plt.axvline(0, color='black', ls='--', lw=1)
# 修飾：x軸ラベル、y軸ラベル、凡例
plt.xlabel('線形予測子 $z$', fontsize=14)
plt.ylabel('確率 $q$', fontsize=14)
plt.legend();

【実行結果】
傾きが「緩 ⇒ 急 ⇒ 緩」に変わりました。
$${z = 0}$$ のとき、確率 $${q=0.5}$$ になります。
x 軸と y 軸を交換すると、先ほどのロジット関数の曲線の形状とよく似ています。

📶 ロジスティック関数の可視化：線形予測子との関係
テキスト p.121 図 6.5 に相当します。
線形予測子 $${x = z_i = \beta_1 + \beta_2 x_i}$$ を用いて
・$${\beta_2=2}$$ に固定したときのチャート
・$${\beta_1=0}$$ に固定したときのチャート
を描画します。

# ロジスティック曲線　p.121 図6.5
# ロジスティック関数＝シグモイド関数にはscipyのexpitを利用

## 設定
x = np.linspace(-6, 6, 101)                # xの設定
beta2_fix, beta1_fix = 2, 0                # Aのβ2, Bのβ1の設定
var_params = [[-3, 0, 2], [-1, 2, 4]]      # Aのβ1, Bのβ0の設定
titles = [f'(A) $\\beta_2={beta2_fix}$', f'(B) $\\beta_1={beta1_fix}$']

## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
# AとBごとに描画処理
for i, (var, title, ax) in enumerate(zip(var_params, titles, axes.flat)):

    if i == 0:  # 左の(A)の図の描画
        for beta1 in var:    # β1の値ごとに描画を繰り返し処理
            ax.plot(x, expit(beta1 + beta2_fix * x), label=f'$\\beta_1$={beta1}')

    else:       # 右の(B)の図の描画
        for beta2 in var:    # β2の値ごとに描画を繰り返し処理
            ax.plot(x, expit(beta1_fix + beta2 * x), label=f'$\\beta_2$={beta2}')

    # 修飾：x軸ラベル、y軸ラベル、タイトル、凡例
    ax.set_xlabel('説明変数 $x$', fontsize=14)
    ax.set_ylabel('確率 $q$', fontsize=13)
    ax.set_title(title)
    ax.legend()

plt.show()

【実行結果】
(A) の切片 $${\beta_1}$$ の変化では、曲線の傾きの変化と、曲線の左右の位置の変化がみられます。
(B) の傾き $${\beta_2}$$ の変化では、曲線の傾きの変化がみられます。$${\beta_2}$$ が負の場合、傾きが反転しています。

統計モデルをデータに当てはめ（当てはめと評価）

■ 統計モデルをデータに当てはめ、の準備
GLMの３要素「確率分布」「リンク関数」「線形予測子」を用いる統計モデルをデータに当てはめします。
統計モデルの対数尤度 $${\log L}$$ が最大になるパラメータ $${\beta_1, \beta_2, \beta_3}$$（線形予測子のパラメータ）を推定します。

🔷 今回の統計モデルの尤度関数 $${L}$$ と対数尤度関数 $${\log L}$$

$${n}$$ は標本サイズです（観察種子数 $${N}$$ と被って…）

◆

少々脱線して…二項分布の尤度関数・対数尤度関数を実装します。

📶 Python で尤度関数・対数尤度関数を表現
尤度関数と対数尤度関数の数式を Python の関数で記述します。

# 二項分布の尤度関数、対数尤度関数の実装

# 尤度関数
def binom_likelihood(y, N, q):
    return np.prod(comb(N, y) * q**y * (1 - q)**(N - y))

# 対数尤度関数
def binom_log_likelihood(y, N, q):
    return sum(np.log(comb(N, y)) + y * np.log(q) + (N - y) * np.log(1 - q))

【実行結果】なし

尤度関数を実験します。
観測データと、scipy の fit で得た確率 $${q=0.63 \cdots}$$ で試します。

# 尤度関数の実験 ★とても小さな値になってしまう
# 生存種子数:y, 種子数:N, 確率:さっきフィッティングした二項分布の確率パラメータ
binom_likelihood(data.y, data.N, res.params.p)

【実行結果】
$${3.18 \times 10^{-140}}$$ というものすごく小さな値になりました。

ちなみに尤度関数の対数をとると：

# ちなみに尤度関数の対数をとると...
np.log(binom_likelihood(data.y, data.N, res.params.p))

【実行結果】

今度は対数尤度関数を実験します。

# 対数尤度関数の実験
# 生存種子数:y, 種子数:N, 確率:さっきフィッティングした二項分布の確率パラメータ
binom_log_likelihood(data.y, data.N, res.params.p)

【実行結果】
尤度関数の対数とほぼほぼ同じです。

scipy で fit したときに得た対数尤度を表示します。

# scipyでfitしたときの最大対数尤度
-res.nllf()

【実行結果】
尤度関数の対数と合致しています。

脱線終了です。

◆ ◆ ◆

■ Python ライブラリで統計モデルをデータに当てはめ
statsmodels の glm を利用して統計モデルを実装します。
引数 family には二項分布 sm.families.Binomial() を設定します。
引数 link は未設定とし、デフォルトのリンク関数 logit を使います。
当てはめ結果を変数 result に格納します。

# GLMによるロジスティック回帰 p.122 ※こちらの最大対数尤度のほうが正確

# 設定
family = sm.families.Binomial()   # GLMの引数familyに与える確率分布＝二項分布

# モデルの当てはめ ※目的変数には生存種子数yと死滅種子数N-yを与える
result = smf.glm(formula='y + I(N-y) ~ x + f', data=data, family=family).fit()
result.summary()

【実行結果】
最下３行の Intercept が $${\beta_1}$$、x が $${\beta_2}$$、f[T.T] が $${\beta_3}$$ に対応しており、テキスト p.123  glm() の推定結果 Coefficients に相当します。

【formula 上の目的変数の表現】
二項分布の場合、目的変数に「成功回数」と「失敗回数」を記述します。
このモデルでは $${\mathtt{y + I(N-y)}}$$ のように記述して、生存種子数 $${\mathtt{y}}$$ と死滅種子数 $${\mathtt{N - y}}$$ を与えます。

formula='y + I(N-y) ~ x + f'

◆

ちょっと脱線です２
statsmodels の glm の別の書き方のご紹介です。
ChatGPTが教えてくれました。

特徴は…
・目的変数に生存割合 $${\mathtt{y / N}}$$ を与える
・引数 freq_weights に観察種子数 N を与える

# GLMによるロジスティック回帰 別の書き方 p.122 ※最大対数尤度等が歪む

# 設定
family = sm.families.Binomial()   # GLMの引数familyに与える確率分布＝二項分布

# モデルの当てはめ ※目的変数に成功割合y/Nを与え、freq_weights引数で試行回数Nを与える
result_freq = smf.glm(formula='I(y/N) ~ x + f', data=data, family=family,
                      freq_weights=data['N']).fit()
result_freq.summary()

【実行結果】

注意点は、対数尤度、逸脱度、AIC などの指標があり得ない値になることです。
従いまして、今回の統計モデリングでは不適切なモデリングになります。

ちなみに、成功割合＆freq_weights の設定が適切な例は…
ChatGPTが教えてくれました。

脱線終了です。

◆ ◆ ◆

■ データにロジスティック回帰モデルを当てはめるイメージ
施肥処理 $${f_i = \mathtt{C}}$$ のデータに絞り込んで、(A) 体サイズ x と生存種子数 y の散布図と、(B) ロジスティック回帰モデルの曲線等を描画します。
テキスト p.123 図 6.6 に相当します。

# ロジスティック回帰のモデルの推定の可視化 p.123 図6.6

## 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)

## 左：(A)例題データの一部(f_i=C)の描画
# 散布図の描画
sns.scatterplot(data=data[data.f=='C'], x='x', y='y', s=80, alpha=0.7, ax=ax1)
# 修飾：x軸範囲、y軸範囲
ax1.set(xlim=(7, 12.5), ylim=(-0.5, 8.5), xlabel='$x$', ylabel='y',
        title='(A) 例題データの一部 ($f_i$=C)')

## 右：(B)推定されるモデルの描画
# 設定
N = 8
x_val = np.linspace(7, 12.3, 101)  # x軸の値 
y_val = list(range(0, 9))          # 二項分布のy軸の値
feat_xs = [8.5, 10, 12]            # 二項分布の変数xの値

# 変数xの値ごとに二項分布の棒グラフ・ロジスティック曲線のデータ点の描画を繰り返し処理
for feat_x in feat_xs:

    ## 二項分布の棒グラフの描画
    # 線形予測子の説明変数xとf=Cの値の設定
    feat_data = dict(x=[feat_x], f=['C'])
    # 説明変数の値に応じた確率の予測値を用いて二項分布の確率質量関数の値を算出
    prob = stats.binom.pmf(k=y_val, n=N, p=result.predict(feat_data))
    # 二項分布の垂直線の描画
    ax2.axvline(feat_x, color='tab:orange', lw=1, alpha=0.5)
    # 二項分布の確率分布の棒グラフの描画
    ax2.barh(y_val, prob*1.5, left=feat_x - prob*1.5, color='tab:orange',
             alpha=0.5)

    ## ロジスティック曲線上の点（変数xの確率値）の描画
    ax2.plot([feat_x], [result.predict(feat_data) * N], 'o', ms=7,
             color='tab:blue')

# ロジスティック曲線の描画
ax2.plot(x_val, result.predict(dict(x=x_val, f=['C']*len(x_val))) * N, lw=2)

# 散布図の描画
sns.scatterplot(data=data[data.f=='C'], x='x', y='y', s=50, alpha=0.7, ax=ax2)

# 修飾：
ax2.set(ylim=(-0.5, 8.5), xlabel='$x$', title='(B) 推定されるモデル');

【実行結果】

【 B のチャートのロジスティック回帰モデルの概要】
（B) の曲線は $${\mathtt{x+f}}$$ モデルによる施肥処理 f = $${\mathtt{C}}$$ の場合の生存種子数 y の予測値の平均（期待値）です。
この予測値は、GLM の結果 result に対して predict メソッドを適用して得た生存確率 q の予測値に観察種子数 N = 8 を乗じて算出しています。
薄オレンジの棒グラフは二項分布の確率質量関数であり、$${x = 8.5, 10, 12}$$ のときの生存種子数 y の「ばらつき」を表しています。

⏰️ 当てはめのアディショナルタイム

当てはめ結果の細かな確認と、最適化ライブラリによる最尤推定を行います。

■ 当てはめ結果の分析
当てはめ結果の出力を再掲します。

🔷 係数の推定値
coef に注目します。
パラメータである係数の最尤推定値は 切片（Intercept）$${\beta_1 = -19.5361}$$、x の係数 $${\beta_2 = 1.9525}$$、f の係数 $${\beta_3 = 2.0215}$$ です。
この推定値を線形予測子に当てはめてみます。
テキスト p.123 の glm 関数によるパラメータ推定値に相当します。

🔷 標準誤差の推定値
続いて標準誤差 std err です。
パラメータ $${\beta_1, \beta_2, \beta_3}$$ の推定値の標準偏差であり、それぞれ $${1.414, 0.139, 0.231}$$ です。

🔷 $${z}$$ 値
$${z}$$ 値はパラメータの推定値を標準誤差で割って求めた統計量です。
パラメータ $${\beta_1, \beta_2, \beta_3}$$ の $${z}$$ 値はそれぞれ $${-13.818, 14.059, 8.740}$$ です。

🔷 P>|z|（$${p}$$ 値）
P>|z| は $${z}$$ 値の $${p}$$ 値です。
パラメータ $${\beta_1, \beta_2, \beta_3}$$ の $${p}$$ 値はすべて$${0.000}$$ です。

🔷 パラメータ推定値の 95% 信頼区間
95% 信頼区間は、区間が 95% の確率でパラメータの真値を含む、そんな区間です。
$${z}$$ 値を用いる場合、95% 信頼区間は「パラメータ推定値$${\pm 1.96 \times}$$標準誤差」です。
パラメータ $${\beta_1, \beta_2, beta_3}$$ の 95% 信頼区間はそれぞれ $${[-22.307\ -16.765],\ [1.680,\ 2.225],\ [1.568,\ 2.475]}$$ です。

🔷 パラメータ推定値の評価
$${p}$$ 値が極小であり、統計的には有意だと考えられます。
95% 信頼区間はゼロを含んでおらず、幅は広くなさそうです。
パラメータ推定値は妥当な感じがします。

ChatGPTによると信頼区間の幅の目安には「相対幅」（信頼区間の幅と係数の比率）を用いることがあるそうです。
ChatGPTの回答を貼ります。

こちらは Google 検索「Search Labs」の AI による概要で判明した目安です。

なお、信頼区間の幅の目安には「明確に統一された教科書的な定義・出典が存在するわけではありません」とのことですので、上記の値をご覧になる際にはご留意ください。

◆ ◆ ◆

■ 最適化ライブラリで最尤推定
こちらは趣味のコードです。
ご興味ない方はスルーしてくださって大丈夫です。

statsmodels を使わないで最尤推定を行ってみます。
scipy の最小化問題ソルバー mininize で、負の対数尤度が最小となるパラメータを推定します。

# 最適化ライブラリで最尤推定

## 追加インポート
from scipy.optimize import minimize

## 設定と準備
# 目的関数の定義：負の対数尤度関数 ※負の対数尤度の最小化問題を解く
def llf(params, x, f, y, N):
    # パラメータの分解
    beta1, beta2, beta3 = params
    # qの算出：ロジスティック関数
    q = expit(beta1 + beta2 * x + beta3 * f)
    # 戻り値：負の対数尤度
    return -np.sum(stats.binom.logpmf(k=y, p=q, n=N))

# 施肥処理を C⇒0,T⇒1に変換
data_f = data['f'].apply(lambda x: 0 if x=='C' else 1)

## 最適化の実行
#  β1, β2, β3の最尤推定値を算出
res_optim = minimize(
    fun=llf, x0=[0, 0, 0], args=(data['x'], data_f, data['y'], 8))

# 結果の表示
print(f'β1の最尤推定値 = {res_optim.x[0]:9.4f}')
print(f'β2の最尤推定値 = {res_optim.x[1]:9.4f}')
print(f'β3の最尤推定値 = {res_optim.x[2]:9.4f}')
print(f'最大対数尤度   = {-res_optim.fun:9.4f}')

【実行結果】
statsmodels の当てはめ結果と同じになっています。

生存種子数の予測

テキストにならって、生存種子数 y の予測を行って可視化します。

■ 予測（期待値）を可視化
ロジスティック回帰モデルを当てはめるイメージと同様に、GLM の結果 result に対して predict メソッドを適用して得た生存確率 q の予測値に観察種子数 N = 8 を乗じて算出しています。
p.124 図 6.7 に相当します。

# ロジスティック回帰のモデルの推定の可視化 p.123 図6.7

# 設定
x_val = np.linspace(7, 12.7, 101)  # x軸の値
fs = ['C', 'T']                    # 施肥なし・あり
titles = ['(A) 施肥処理なし ($f_i$=C)', '(B) 施肥処理あり ($f_i$=T)']
N = 8

# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)

# 施肥なり・ありごとに観測値の散布図・ロジスティック関数の描画を繰り返し処理
for f, title, ax in zip(fs, titles, axes.flat):
    # 散布図の描画
    sns.scatterplot(data=data[data.f==f], x='x', y='y', s=80, alpha=0.7,
                    ax=ax)
    # ロジスティック関数の描画
    ax.plot(x_val, result.predict(dict(x=x_val, f=[f]*len(x_val))) * N,
               color='tab:red')
    # 修飾：x軸範囲、y軸範囲
    ax.set(xlim=(7, 12.7), ylim=(-0.5, 8.5), xlabel='植物の体サイズ $x_i$',
           ylabel='生存種子数 $y_i$', title=title)

【実行結果】
（B) の施肥ありの方が、生存種子数 y の予測値の曲線が左側にシフトしています。
施肥処理をすることによって、施肥なしよりも小さな体サイズで生存種子数が多くなる、こんな予測になっています。

◆ ◆ ◆

■ 予測区間も推定してみる
生存種子数の 95% 予測区間の可視化にもチャレンジします。
ChatGPT と相談しながら「正規分布近似」による予測区間の計算を目指します。
上のチャートに予測区間を重ねます。

# ロジスティック回帰のモデルの推定の可視化 p.123 図6.7＋

# 設定
x_val = np.linspace(7, 12.7, 101)  # x軸の値
fs = ['C', 'T']                    # 施肥なし・あり
N = 8                              # 観察種子数
z = stats.norm.isf(q=0.025)        # 正規分布の上側2.5%点 ≈ 1.96
titles = ['(A) 施肥処理なし ($f_i$=C)', '(B) 施肥処理あり ($f_i$=T)']

# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)

# 施肥なり・ありごとに観測値の散布図・ロジスティック関数の描画を繰り返し処理
for f, title, ax in zip(fs, titles, axes.flat):
    # 散布図の描画
    sns.scatterplot(data=data[data.f==f], x='x', y='y', s=80, alpha=0.7,
                    ax=ax)
    # ロジスティック関数の描画
    q_hat = result.predict(dict(x=x_val, f=[f]*len(x_val)))
    ax.plot(x_val, q_hat * N, color='tab:red')
    # 修飾：x軸範囲、y軸範囲
    ax.set(xlim=(7, 12.7), ylim=(-0.5, 8.5), xlabel='植物の体サイズ $x_i$',
           ylabel='生存種子数 $y_i$', title=title)
    # 95%予測区間の描画 ※正規分布近似を利用
    mean = N * q_hat                        # 平均
    std = np.sqrt(N * q_hat * (1 - q_hat))  # 標準偏差
    lower = mean - z * std                  # 95%信頼区間 下端
    upper = mean + z * std                  # 95%信頼区間 上端
    ax.fill_between(x_val, lower, upper, color='lightpink', alpha=0.3)

【実行結果】
薄赤色の塗りつぶし部分が生存種子数 y の 95% 予測区間です。
正規分布近似を用いたので帯が滑らかです。

ロジットリンク関数の解釈

ロジット関数にパラメータ推定値を当てはめて、オッズ $${\cfrac{q_i}{1-q_i}}$$ の形式に変換します。

オッズは 成功確率と失敗確率の比です。
観測データの場合、生存確率と死滅確率の比です。
テキストによると、$${q_i = 0.5}$$ のとき、$${\frac{0.5}{1-0.5}=1}$$ でオッズは１倍、 $${q_i = 0.8}$$ のとき、$${\frac{0.8}{1-0.8}=4}$$ でオッズは４倍、 と言うそうです。

ロジスティック回帰モデルによるオッズは、上のロジット関数の変形を使って次のように変化します。

• 体サイズが１単位大きくなると、オッズは $${\exp(1.952) \approx 7.0}$$ 倍になる

• 施肥を行うと、施肥しない場合に比べてオッズは $${\exp(2.022) \approx 7.5}$$ 倍になる

テキストは「ロジットリンク関数で生存確率を定義することによって、さまざまな要因と応答事象のオッズの解釈が簡単になります」とまとめています。

モデル選択

以下の４つのロジスティック回帰モデルはネストしているモデルになっています。
それぞれのモデルを観測データに当てはめて、AIC でモデル選択を行います。

４つのモデルを一気に当てはめして AIC 等の統計量を算出するコードを実行します！

# 種子の生存確率モデルのAICなど p.127 表6.2
# ちなみに情報：以下のコードでフルモデルの最大対数尤度を算出できる
# np.sum(stats.binom.logpmf(k=data.y, n=data.N, p=data.y / data.N))

## 設定
# モデル名のリスト（最後にフルモデルを表す'フル'を設定）
model_names = ['一定', 'f', 'x', 'x+f', 'フル']
# 各モデルのformulaのリスト（最後にフルモデルで用いる'一定モデル'を設定）
formulas = ['y + I(N-y) ~ 1', 'y + I(N-y) ~ f', 'y + I(N-y) ~ x',
            'y + I(N-y) ~ x + f', 'y + I(N-y) ~ 1']
# GLMの引数familyに与える確率分布＝二項分布
family = sm.families.Binomial()
# 結果を格納するデータフレームの初期化
aic_df = pd.DataFrame()


## 関数定義
# モデルの結果からデータフレーム（１行）を作成する関数
def make_df(model_name, k, llf, deviance, resid_deviance, aic):
    return pd.DataFrame({'モデル': [model_name], '説明変数の数': k,
                         '最大対数尤度': llf, '逸脱度':  deviance,
                         '残差逸脱度': resid_deviance, 'AIC': aic})

# モデルのあてはめ～結果から１行のデータフレームを作成する関数
def model_fitting(data, model_name, formula, family):
    # GLM・ロジスティック回帰のあてはめ
    result = smf.glm(formula=formula, data=data, family=family).fit()
    # フルモデルの場合、一定モデルの結果から各種数値を算出する関数でデータフレーム化
    if model_name == 'フル':
        return full_model_result(result, data, model_name)
    # フルモデル以外の場合、あてはめ結果の各種数値をデータフレーム化
    else:
        return make_df(model_name, len(result.params), result.llf,
                       -2 * result.llf, result.deviance, result.aic)

# フルモデルの各種数値を算出する関数
def full_model_result(result, data, model_name):
    # 各種数値を算出
    k = len(data)                                      # 説明変数の数
    deviance = -2 * result.llf - result.null_deviance  # 逸脱度
    max_llf = deviance / -2                            # 最大対数尤度
    aic = -2 * (max_llf - k)                           # AIC
    # 戻り値：データフレーム化した各種数値
    return make_df(model_name, k, max_llf, deviance, 0, aic)


## モデル比較の実行
# 各モデルのあてはめと各種数値算出を繰り返し処理
for model_name, formula in zip(model_names, formulas):
    # モデルのあてはめと各種数値の算出
    tmp_df = model_fitting(data, model_name, formula, family)
    # 結果を格納するデータフレームに追加
    aic_df = pd.concat([aic_df, tmp_df], axis=0)

# データフレームの最終化と結果表示(最小AICをハイライト)
aic_df = aic_df.reset_index(drop=True)  # インデックスのリセット
aic_df.style.highlight_min(subset='AIC', color='lightpink').format(precision=1)

【実行結果】
AIC 最小の $${\mathtt{x+f}}$$ モデルを選択します！

まとめ

今回はロジスティック回帰と呼ばれる GLM を実践しました。

🔷 確率分布と確率質量関数

🔷 リンク関数と線形予測子

🔷 オッズ $${\cfrac{q_i}{1-q_i}}$$

🔷 statsmodels のロジスティック回帰モデル構築と結果表示

family = sm.families.Binomial()
result = smf.glm(formula='y + I(N-y) ~ x + f', data=data, family=family).fit()
result.summary()

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

次回は、交互作用項を含むロジスティック回帰を実践します。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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