「データ解析のための統計モデリング入門」をPythonで写経 Vol.21 ～ 10章「階層ベイズモデル」②ベイズ統計モデリング～階層ベイズモデル（個体差＋場所差）

10章「階層ベイズモデル」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」10章「階層ベイズモデル」の Python写経活動記録 です。 

この記事は前回記事に引き続き GLMM の 階層ベイズモデル 化に取り組みます。
前回からの変化点は、ランダム切片に「場所差」（植木鉢差）が加わることです。
二項分布・ロジットリンク関数・ランダム切片（個体差と場所差）の GLMM がベイズ統計モデルに！

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

はじめに

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

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

準備

準備

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

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

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

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

# インポート

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

# 統計計算
import scipy.stats as stats

# PyMC
import pymc as pm
import pytensor.tensor as pt
import arviz as az

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

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

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

■ 統計モデル
GLMM のベイズ統計モデルです。
ランダム切片は「個体差」と「植木鉢差」です。

■ モデリング手続き

データの確認

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

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

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

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

【変数の説明】
植物の個体 id ごとの生存種子数 y です。

【施肥有無と植木鉢の関係】
植木鉢 A ～ E は施肥なし、植木鉢 F ～ J は施肥ありです。

🍀🍀🍀

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

① 要約統計量の表示

# 要約統計量
data.describe().round(3)

【実行結果】
生存種子数 y の値は 0 ～ 37 個です。

② 標本分散の表示

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

【実行結果】

y の分散は 54.394 です。
y の分布にポアソン分布を仮定すると標本分散は標本平均と近い値になるはずですが、平均 5.520 よりも分散はとても大きな値であり、過分散が発生することになります。

③ y の標本平均と標本分散の深堀り
施肥有無別の y の標本平均と標本分散を確認します。

# 施肥有無別の標本平均と標本分散
data.groupby(['f'])['y'].agg(['mean', 'var']).round(3)

【実行結果】
施肥あり $${\mathtt{T}}$$ の方が種子数 y が小さいとは…（波乱の予感）
施肥有無別に見ても、標本平均に比べて標本分散はかなり大きな値になっています。
個体ごとのばらつきが大きいのかもしれません。
そうならば、個体差が過分散の要因になるでしょう。

続いて植木鉢別の y の標本平均と標本分散を確認します。

# 植木鉢別の標本平均と標本分散
data.groupby(['pot'])['y'].agg(['mean', 'var']).round(3)

【実行結果】
施肥なしの植木鉢 A ～ E、施肥ありの植木鉢 F ～ J の標本平均は一体どんな傾向があるというのでしょう（無いですね無いですね）。
標本平均と標本分散が近い鉢があれば、遠い鉢もあり、植木鉢差が過分散の要因になっている感じがします。

④ y のヒストグラムの描画
種子数 y のヒストグラムで分布の形状を確認します。

# yのヒストグラムの描画
sns.histplot(data=data, x='y', edgecolor='white');

【実行結果】
ゼロ付近に多くの観測値が集中し、右に裾の長い強い正の歪みをもつ単峰の分布です。

⑤ 個体別・植木鉢別のばらつきの可視化
テキスト p.135 図 10.7 の散布図・箱ひげ図に相当するチャートを描画します。

# 観測値の可視化 p.235 図10.7

## 設定と準備
# 色（無処理：青、施肥処理：赤）
colors = ['tab:blue', 'tab:red']

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

## (A)個体ごとの描画
# 施肥なしの平均値の水平点線（青）の描画
ax1.plot(range(50), [data[data.f=='C']['y'].mean()] * 50,
         color=colors[0], ls='--', lw=2)
# 施肥ありの平均値の水平点線（赤）の描画
ax1.plot(range(50, 100), [data[data.f=='T']['y'].mean()] * 50,
         color=colors[1], ls='--', lw=2)
# 個体iと種子数yの散布図の描画（テキスト表示のためのダミー描画）
ax1.scatter(data.id, data.y, color='none')
# 個体iと種子数yの植木鉢の文字による散布図の描画
for i in range(100):
    ax1.text(x=data.id[i], y=data.y[i], s=f'{data.pot[i]}',
             ha='center', va='center', color=colors[i//50])
# 施肥なし・ありの文字の表示
ax1.text(x=25, y=40, s='無処理', ha='center', color=colors[0], fontsize=12)
ax1.text(x=75, y=40, s='施肥処理', ha='center', color=colors[1], fontsize=12)
# 修飾
ax1.set_xlabel('個体 $i$', fontsize=12)
ax1.set_ylabel('種子数 $y_i$', fontsize=12)
ax1.set(ylim=(None, 45), title='(A) 個体ごと')

## (B)植木鉢ごとの描画
# 箱ひげ図の描画
sns.boxplot(data=data, x='pot', y='y', fill=False,
            hue='f', palette=colors, legend=False, ax=ax2)
# 施肥なし・ありの文字の表示
ax2.text(x=2, y=40, s='無処理', ha='center', color=colors[0], fontsize=12)
ax2.text(x=7, y=40, s='施肥処理', ha='center', color=colors[1], fontsize=12)
# 修飾
ax2.set_xlabel('植木鉢(pot) $j$', fontsize=12)
ax2.set(ylabel='', title='(B) 植木鉢ごと')
ax2.set_ylim(top=45);

【実行結果】
（A) 個体ごとのチャートの水平点線は施肥有無別の種子数の平均です。
個体のばらつきも植木鉢のばらつきもさまざまです。
モデルに個体差・植木鉢差を組み込みたくなるデータです。

🍀🍀🍀

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

テキストは p.225 ～ で「個体差だけでなく場所差なども考慮する GLMM の場合は、パラメータの最尤推定が困難になります」と解説しています。
ということでベイズ統計モデルの登場です！

特徴①②より「ポアソン分布」「対数リンク関数」とし、特徴③より「個体差と植木鉢差（ランダム切片）」を考慮する GLMM をベイズ統計モデル化します。
今回取り組むモデルは 階層ベイズモデル です。

階層ベイズモデル（個体差＋植木鉢差）

ベイズ統計モデルの理解

◼️ パラメータの事後分布
ベイズ統計モデリングでは「パラメータの事後分布の推定」を中心に置いて動きます。
ベイズ統計モデルの事後分布は尤度と事前分布の積に比例します。

◼️ 今回モデルの概観
今回のベイズ統計モデリングで推定するパラメータは $${\beta_1, \beta_2, \{r_i\}, \{r_{pj(i)}\}, s, sp}$$ です（詳細は後ほど！）。
データを $${\bm Y}$$ とし、事後分布を数式化します。
$${n}$$ は標本サイズであり、例題データの場合は 100 です。

◼️ GLMM の構成要素
個体 $${i}$$ の種子数 $${y_i}$$ はポアソン分布 $${p(y_i \mid \lambda_i)}$$ に従うとします。
平均パラメータ $${\lambda_i}$$ は線形予測子と対数リンク関数を用いて $${\log \lambda_i = \beta_1 + \beta_2 f_i + r_i + rp_{j(i)}}$$ とします。
$${\beta_1, \beta_2}$$ は全個体共通のパラメータです。
個体差 $${r_i}$$ は平均パラメータ 0、標準偏差パラメータ $${s}$$ の正規分布 $${\text{Normal}(0, s^2)}$$ に従うとします。
植木鉢差 $${r_{pj(i)}}$$ は平均パラメータ 0、標準偏差パラメータ $${s_p}$$ の正規分布 $${\text{Normal}(0, s_p^2)}$$ に従うとします。

🍀🍀🍀

◼️ 尤度
尤度は次の数式で表されます。

PyMC のコードに寄せて、次のように表現します。

🍀🍀🍀

◼️ 事前分布
事後分布の数式化に現れた６つの事前分布を具体化します。

1️⃣ $${\beta_1}$$ と $${\beta_2}$$ の事前分布
無情報事前分布を指定します。
平均パラメータ 0、標準偏差パラメータ 100 の「押しつぶされた正規分布」です。

2️⃣ $${r_i}$$ と $${r_{pj(i)}}$$ の事前分布
平均パラメータ 0、標準偏差パラメータ $${s, s_p}$$ の正規分布を指定します。
階層事前分布に該当します。

3️⃣ $${s}$$ と $${s_p}$$ の事前分布
無情報事前分布を設定します。
幅が十分に広い連続一様分布です。

🍀🍀🍀

今回のベイズ統計モデルの数式をまとめます。

$${\sim}$$ は左側の確率変数が右側の確率分布に従うことを意味します。
$${=}$$ で表現された変数は、等号で結ばれた数式どおりに特定の値が決定される「決定論的変数」です。

パラメータの事後分布の推定

PyMC ライブラリでベイズ統計モデリングを実装します。
テキストの WinBUGS の設定を解釈しつつ、PyMC コードに書き換えていきます。

◼️ 設定
カテゴリデータを数値に変換します。
・施肥処理を $${\texttt{C:0, T:1}}$$ のダミー変数化
・植木鉢番号（A ～ J）を ０から始まる整数連番に変換

# 設定と準備

# 施肥処理：文字列のカテゴリデータを数値に変換 C:0, T:1
F_val = (data['f'] == 'T').astype(int).values

# 植木鉢番号：文字列のカテゴリデータを数値に変換 A:0～J:9
Pot_val = data['pot'].astype('category').cat.codes        # 数値化した植木鉢番号
Pot_cat = data['pot'].astype('category').cat.categories   # 元の文字列

◼️ モデルの定義

# モデルの定義

# coordsの設定
coords = {
    'id': data.id.values,   # 個体番号
    'pot': Pot_cat,         # 植木鉢番号（文字列）
}

# モデリング
with pm.Model(coords=coords) as model:
    
    # dataの定義
    # 目的変数: 種子数Y
    Y = pm.Data('Y', value=data['y'].values, dims='id')
    # 説明変数: id別の施肥処理の有無: 0:施肥なし, 1:施肥あり
    F = pm.Data('F', value=F_val, dims='id')
    # 説明変数: id別の植木鉢番号: 0:A～9:J
    Pot = pm.Data('Pot', value=Pot_val, dims='id')

    # 事前分布
    beta1 = pm.Normal('beta1', mu=0, sigma=100)       # β1: 無情報 N(0,100)
    beta2 = pm.Normal('beta2', mu=0, sigma=100)       # β2: 無情報 N(0,100)
    s = pm.Uniform('s', lower=0, upper=10**4)         # s:  無情報 U(0,10^4)
    sp = pm.Uniform('sp', lower=0, upper=10**4)       # sp: 無情報 U(0,10^4)
    r = pm.Normal('r', mu=0, sigma=s, dims='id')      # r:  階層   N(0,s)
    rp = pm.Normal('rp', mu=0, sigma=sp, dims='pot')  # rp: 階層   N(0,sp)

    # 線形予測子: 指数関数でlog(λ)を平均λに変換
    lam = pm.Deterministic(
        'lam', pt.exp(beta1 + beta2 * F + r + rp[Pot]), dims='id')

    # 尤度関数: 平均λのポアソン分布
    obs = pm.Poisson('obs', mu=lam, observed=Y, dims='id')

【実行結果】なし

【コードの補足：線形予測子中の植木鉢差 rp の書き方】
植木鉢差 rp は植木鉢番号ごと（10個）の配列的な変数です。
平均値 lam は個体単位（100個）です。
どの個体がどの植木鉢に植えられているかを変数 Pot（データの pot ）に持っているので、$${\texttt{rp[Pot]}}$$ と指定することで、個体ごとに植木鉢差 rp を識別できます。

pt.exp(beta1 + beta2 * F + r + rp[Pot])

🍀🍀🍀

◼️ モデルの確認
モデルの数式と有向グラフを可視化します。
数式を表示します。

# モデルの表示
model

【実行結果】

モデルの有向グラフを描画します。

# モデルの可視化
pm.model_to_graphviz(model)

【実行結果】
個体差 $${\texttt{r}}$$ は個体 id 単位で推定されます。
植木鉢差 $${\texttt{rp}}$$ は植木鉢 pot 単位で推定されます。

🍀🍀🍀

◼️ MCMC の実行
MCMC サンプリングを行います。
NUTS サンプラーに nutpie を使います。
合計 20000 個の MCMC サンプルを得ます。

%%time
# MCMCサンプリング
# chain=4, draws=5000（テキストは50000）, tune=1000, thinなし, パラメータ初期値なし,
# NUTSサンプラー利用

with model:
    idata = pm.sample(
        draws=5000, tune=1000, chains=4, random_seed=123,
        nuts_sampler='nutpie',  # nutpieを使わない場合はこの行を削除
    )

【実行結果】
Divergences は０個です。

パラメータ推定値の確認

ざっくり収束の確認などを行います。

◼️ $${\widehat{R}}$$ の確認
$${\widehat{R}}$$ が 1.01 以下になっていることを確認します。
全パラメータの「$${\widehat{R}}$$ ＞1.01」の個数が０になればOKです。

# r_hat>1.01の確認
# 設定
idata_in = idata         # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
すべてのパラメータの $${\widehat{R}}$$ は1.01 以下です。

◼️ 有効サンプルサイズ（ESS）
「実質的に独立なサンプル数」であるESS が 400 以上になっていることを確認します。
全パラメータの「ESS ＜ 400」の個数が０になればOKです。

# 有効サンプルサイズ しきい値を400とした
print((az.ess(idata) < 400).sum())

【実行結果】
すべてのパラメータの ESS は400 以上です。

◼️ トレースプロットの確認

# トレースプロットの表示
var_names = ['beta1', 'beta2', 's', 'sp', 'rp', 'r']
pm.plot_trace(idata, var_names=var_names, backend_kwargs={'tight_layout': True});

【実行結果】
右のトレースプロットは、各チェーンが「毛糸玉」のようにゲジゲジと混ざり合い、ドリフトがないのでＯＫ（収束ＯＫ）としましょう。

🍀🍀🍀

◼️ 事後分布の要約統計量
パラメータの事後分布の推定値を確認します。
推定値は対数スケールですので、$${\exp(\cdot)}$$ することで平均種子数 $${\lambda_i}$$ の単位に変換できます。

# 推論データの要約統計情報の表示
pm.summary(idata, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】
$${\beta_2}$$ 推定値の $${95\%}$$ HDI 区間が０を含んでいる点が気になります。
後ほど詳細を確認します。

🍀🍀🍀

◼️ 事後分布の確認
推定したパラメータの事後分布を可視化します。
こちらは $${\beta_1, \beta_2}$$ です。

# 事後分布プロット β1, β2
pm.plot_posterior(
    idata, var_names=['beta1', 'beta2'], hdi_prob=0.95, round_to=3,
    ref_val=0, ref_val_color='tab:red',
    backend_kwargs={'tight_layout': True, 'sharex': True, 'sharey': True},
    figsize=(10, 3));

【実行結果】

続いて $${s, s_p}$$ です。

# 事後分布プロット s, sp
pm.plot_posterior(
    idata, var_names=['s', 'sp'], hdi_prob=0.95, round_to=3,
    backend_kwargs={'tight_layout': True, 'sharex': True, 'sharey': True},
    figsize=(10, 3));

【実行結果】
植木鉢差のばらつきパラメータ $${s_p}$$ は右に裾が長い偏った分布になっています。

最後は植木鉢差 $${r_{pi(j)}}$$ です。

# 事後分布プロット rp
pm.plot_posterior(
    idata, var_names=['rp'], hdi_prob=0.95, round_to=3,
    ref_val=0, ref_val_color='tab:red',
    textsize=10, grid=(2, 5), figsize=(10, 5),
    backend_kwargs={'tight_layout': True, 'sharex': True, 'sharey': True}
);

【実行結果】
分布の形状はよく似ています。
０の赤い補助線を手がかりにすると、植木鉢差はばらついてることが分かります。

フォレストプロットも見てみましょう。

# フォレストプロットの描画 rp: 植木鉢差
pm.plot_forest(
    idata, var_names=['rp'], hdi_prob=0.95, combined=True, figsize=(5, 4)
)
plt.axvline(0, color='tab:red', ls='--');

【実行結果】
植木鉢差がバラバラとしている様子がわかります。
また、植木鉢 I を除いて 95% HDI は０を含んでいるようです。
植木鉢差が無い（植木鉢差が０である）ことを否めないかもです。

🍀🍀🍀

◼️ 事後予測チェック
y の事後予測値と観測値を可視化して比べます。

# 事後予測チェック

# 事後予測
with model:
    idata_pp = pm.sample_posterior_predictive(idata)

# 事後予測チェックプロット
pm.plot_ppc(idata_pp, num_pp_samples=100, random_seed=123);

【実行結果】
オレンジ点線は事後予測値の平均、薄青色線は MCMCサンプルのうち 100 個分の事後予測値、黒実線は観測値です。
事後予測値のオレンジ点線と観測値の黒実線はまあまあ似ていて、ベイズ統計モデルに大きな問題は無さそうです。

（注意）
アディショナルは趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

アディショナル：事後分布・事後予測の考察

事後分布と事後予測のサンプルを用いて、簡単な考察を行います。
ChatGPT にサポートしてもらいました。

idata 内の MCMC サンプルを取り出すヘルパー関数を定義します。

# ---- idataからMCMCサンプルを取り出すヘルパー関数 ----

# xarray を shape=(samples,) の numpy 配列に変換
def stack_1d(xa):  # 例：beta2, s, sp
    return xa.stack(sample=('chain','draw')).values  # shape (S,)

# xarray を shape=(samples, id) の numpy 配列に変換
def stack_by_id(xa):  # 例：posterior_predictive['obs'], r
    return xa.stack(sample=('chain','draw')).transpose('sample','id').values

【実行結果】なし

🍀🍀🍀

◼️ 切片 $${\beta_1}$$

# =========================================
# 1) 切片 β1 （対数スケール、観測値スケール）
# =========================================

# beta1のMCMCサンプルの平坦化
beta1 = stack_1d(idata.posterior['beta1'])  # shape (S,)

# beta1が0以上の確率、beta1の中央値・95%HDIの算出
pr_beta1_pos = float((beta1 > 0).mean())
beta1_med = np.quantile(beta1, 0.5)
beta1_lo, beta1_hi = az.hdi(beta1, hdi_prob=0.95)

# 観測値スケールでの中央値・95%HDIの算出
exp_beta1 = np.exp(beta1)
exp_beta1_med = np.quantile(exp_beta1, 0.5)
exp_beta1_lo, exp_beta1_hi = az.hdi(exp_beta1, hdi_prob=0.95)

# 結果の表示
print('【切片 β1（対数スケール）】')
print(f'  Pr(β1 > 0) = {pr_beta1_pos:.3f}')
print(f'  β1: median={beta1_med:.3f}, 95%HDI=[{beta1_lo:.3f}, {beta1_hi:.3f}]')
print('\n【観測値スケール】')
print(f'  rate ratio: median={exp_beta1_med:.3f}, '
      f'95%HDI=[{exp_beta1_lo:.3f}, {exp_beta1_hi:.3f}]')

# （お好み）観測値スケールのヒストグラム
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(exp_beta1, bins=40, alpha=0.85, edgecolor='white')
ax.set_xlabel('exp(β1)', fontsize=12)
ax.set_ylabel('頻度', fontsize=12)
ax.set_title('$\\beta_1$ の事後分布（観測値スケール）')
plt.show()

【実行結果】

【考察】
観測値スケールの $${\exp(\beta_1)}$$ は「施肥なし・個体差なし・植木鉢差なし」の場合の平均種子数に相当します。
中央値は 3.9 個、95% HDI は 0.9 ～ 8.8 個と幅が広い感じです。
この「個数」がベースとなり、その他のパラメータの指数値$${\exp(\beta ほにゃらら)}$$ で「倍」計算されます。

🍀🍀🍀

◼️ 傾き $${\beta_2}$$

# =========================================
# 2) 施肥効果 β2 と 発生率比 exp(β2)
# =========================================

# beta2のMCMCサンプルの平坦化
beta2 = stack_1d(idata.posterior['beta2'])  # shape (S,)

# beta2が0以上の確率、beta2の中央値・95%HDIの算出
pr_beta2_pos = float((beta2 > 0).mean())
beta2_med = np.quantile(beta2, 0.5)
beta2_lo, beta2_hi = az.hdi(beta2, hdi_prob=0.95)

# 施肥あり／なしの発生率比とその中央値・95%HDIの算出
rr = np.exp(beta2)
rr_med = np.quantile(rr, 0.5)
rr_lo, rr_hi = az.hdi(rr, hdi_prob=0.95)

# 結果の表示
print('【施肥効果 β2（対数スケール）】')
print(f'  Pr(β2 > 0) = {pr_beta2_pos:.3f}')
print(f'  β2: median={beta2_med:.3f}, 95%HDI=[{beta2_lo:.3f}, {beta2_hi:.3f}]')
print('\n【発生率比 exp(β2)（観測値スケールでの倍率）】')
print(f'  rate ratio: median={rr_med:.3f}, 95%HDI=[{rr_lo:.3f}, {rr_hi:.3f}]  '
      f'# 1なら無効果')

# （お好み）発生率比のヒストグラム
fig, ax = plt.subplots(figsize=(5, 3))
ax.hist(rr, bins=40, alpha=0.85, edgecolor='white')
ax.axvline(1.0, color='tab:red', ls='--')
ax.set_xlabel('exp(β2) 発生率比 [倍]', fontsize=12)
ax.set_ylabel('頻度', fontsize=12)
ax.set_title('施肥効果 $\\beta_2$ の発生率比（事後分布）')
ax.set_xticks(range(int(rr.max())+1))
plt.show()

【実行結果】

【考察】
対数スケールの $${\beta_2}$$ は施肥ありのときの固定効果です。
０超となる確率は 0.12、０以下となる確率は 0.88 です。
施肥効果はきっと無いのでしょう。

観測値スケールの $${\exp(\beta_2)}$$ は「施肥あり」の場合の平均種子数の「倍率」です。
中央値で 0.45 倍、95% HDI は 0.04 ～ 1.53 倍です。
肥料をやることで種子数は減るかも・増えるかも、のどっちつかずな状態です。
下のチャートでは $${\exp(\beta_2)=1}$$（１倍、つまり影響なし）に赤点線を引いています。
１倍以下の確率が高そうです。

ChatGPT の考察は以下のとおりです。

ちなみにテキスト p.239 に「肥料の効果をゼロと設定して架空データを生成した」と書かれています。

$${\beta_2=0}$$ ですって！

対数スケールの 95% HDI が０を含むのは「妥当な推定」（テキストより）なのです。

🍀🍀🍀

◼️ 個体差 $${r}$$ と植木鉢差 $${r_{p}}$$

# =========================================
# 3) ばらつきの強さ：個体差の s と 植木鉢差の sp
#  ・s, sp はいずれも「log(λ) 尺度の標準偏差」
#  ・観測値スケール（λ）の“1SDシフト倍率”は exp(s), exp(sp)
# =========================================

# sとspのMCMCサンプルを平坦化してnumpy配列化
s  = stack_1d(idata.posterior['s'])
sp = stack_1d(idata.posterior['sp'])

# s,spの中央値、95%HDIを算出
s_med = np.quantile(s, 0.5)
s_lo,  s_hi  = az.hdi(s, hdi_prob=0.95)
sp_med = np.quantile(sp, 0.5)
sp_lo, sp_hi = az.hdi(sp, hdi_prob=0.95)

# “1SD シフト時の倍率”の直感用（観測値スケール）：exp(sd)
s_factor_med  = float(np.exp(s_med))
sp_factor_med = float(np.exp(sp_med))

print('【個体差 s（対数スケールのSD）】')
print(f'  s: median={s_med:.2f}, 95%HDI=[{s_lo:.2f}, {s_hi:.2f}]')
print(f'  1SDシフト倍率（λの直感）: exp(median s) ≈ ×{s_factor_med:.2f}')

print('\n【植木鉢差 sp（対数スケールのSD）】')
print(f'  sp: median={sp_med:.2f}, 95%HDI=[{sp_lo:.2f}, {sp_hi:.2f}]')
print(f'  1SDシフト倍率（λの直感）: exp(median sp) ≈ ×{sp_factor_med:.2f}')

# （お好み）s と sp のMCMCサンプルのヒストグラム
fig, axes = plt.subplots(1,2, figsize=(8.5,3))
axes[0].hist(s,  bins=40, alpha=0.85, edgecolor='white')
axes[0].set_title('s の事後分布')
axes[0].set_xlabel('$s$', fontsize=12)
axes[1].hist(sp, bins=40, alpha=0.85, edgecolor='white')
axes[1].set_title('sp の事後分布')
axes[1].set_xlabel('$s_p$', fontsize=12)
for ax in axes:
    ax.set_ylabel('頻度', fontsize=12)
plt.tight_layout()
plt.show()

【実行結果】

【考察】
1SDシフト倍率は 1 標準偏差の倍率です（中央値利用）。
個体差 $${r}$$ は 2.75 倍、植木鉢差 $${r_p}$$ は 2.65 倍であり、1 標準偏差で 平均種子数が 2.6 ～ 2.8 倍となるようです。
下のチャートは２つのパラメータの事後分布です。
$${s_p}$$ は右に裾の長い分布ですね（デジャヴ）。

ChatGPT の考察は以下のとおりです。
SD は標準偏差です。

🍀🍀🍀

◼️ y の事後予測の施肥あり・なしの差
施肥ありの事後予測と施肥なしの事後予測の差を確認します。
種子数の予測値による施肥ありの効果を見ます。

# =========================================
# 4) 施肥の群差 Δ = E[y~|F=1] − E[y~|F=0] （観測の単位で）
#    事後予測 y~ を使うので「観測の世界」での差がそのまま出ます
# =========================================

# 施肥有無をBool値に変換（施肥なし：False, 施肥あり：True）
F_bool = (data['f'].values == 'T')

# 事後予測サンプルを取得：y~ （shape (Sample, N)）
y_pp = stack_by_id(idata_pp.posterior_predictive['obs'])

# サンプルごとに群平均を取り、その差を作る shape=(Sample, )
mu_T = y_pp[:, F_bool].mean(axis=1)                  # F=1 の平均（各サンプル）
mu_C = y_pp[:, np.logical_not(F_bool)].mean(axis=1)  # F=0 の平均（各サンプル）
delta = mu_T - mu_C                                  # 差 Δ（各サンプル）

# 要約（確率・中央値・95%HDI）の算出
pr_delta_pos = float((delta > 0).mean())
d_med = np.quantile(delta, 0.5)
d_lo, d_hi = az.hdi(delta, hdi_prob=0.95)

# 結果の表示
print('【yの事後予測：施肥あり・なしの差 Δ = E[y~|T] − E[y~|C]（観測値スケール）】')
print(f'  Pr(Δ > 0) = {pr_delta_pos:.3f}')
print(f'  Δ: median={d_med:.3f}, 95%HDI=[{d_lo:.3f}, {d_hi:.3f}]  '
      '# 0 を跨ぐなら差は不確か')

# （お好み）Δ のヒスト（0 の縦線付き）
fig, ax = plt.subplots(figsize=(5.2,3))
ax.hist(delta, bins=40, alpha=0.85, edgecolor='white')
ax.axvline(0.0, color='tab:red', ls='--')
ax.set_xlabel('Δ = E[y~|T] − E[y~|C] [個]', fontsize=12)
ax.set_ylabel('頻度', fontsize=12)
ax.set_title('$y$ の事後予測：施肥あり・なしの差 Δ')
plt.show()

【実行結果】

【考察】
Pr(Δ > 0 ) は、差が０超である、つまり、施肥ありの y の事後予測の方が大きい確率であり、なんと 0.000 ！
この「差の事後予測サンプル」のヒストグラムは０未満に正規分布に似たきれいな分布を描いています。
ちなみに、差の事後予測サンプル 20000 個のうち ０超は 7 個だけです。

delta[delta > 0]

【実行結果】

ChatGPT の考察は以下のとおりです。

（注意）
こちらの予測は趣味的な深堀りとコードであり、テキストに掲載はありません。
ご興味ない方はスルーしてくださって大丈夫です。

種子数の予測

平均種子数 $${\lambda}$$ のMCMC サンプルを用いて、予測を可視化します。

描画およびヘルパー関数群を定義します。

# チャート描画ヘルパー関数 by ChatGPT

def _summary_band(arr, q_low=0.025, q_high=0.975, axis=0):
    '''
    arr に対し、中央値・下側・上側のCI分位点を返す（軸は axis）。
    返り値: (median, q_low, q_high)
    '''
    med = np.quantile(arr, 0.5, axis=axis)
    lo  = np.quantile(arr, q_low, axis=axis)
    hi  = np.quantile(arr, q_high, axis=axis)
    return med, lo, hi

def _histo_counts(values, kmax):
    '''
    離散値（0..kmax）のヒストグラムカウントを返す。
    values: 1次元整数配列
    '''
    c = np.bincount(values, minlength=kmax + 1)
    return c[:kmax + 1]

def plot_ppc_frequency_band(
        y_obs, y_pp_samples, ax, title='生存種子数ごとの個体数'):
    '''
    各カウント k=0..kmax の度数について、事後予測の中央値・95%区間のバンドと
    観測の度数を重ねて表示。
    '''
    # x軸の上限：観測とPPCの 99.5% 分位の大きい方まで
    kmax = int(max(np.max(y_obs), np.quantile(y_pp_samples, 0.995)))
    S, N = y_pp_samples.shape

    # 各サンプルごとに、0..kmax のカウントベクトルを作る
    counts = np.vstack([_histo_counts(y_pp_samples[s], kmax) for s in range(S)])

    # 中央値・95%CI
    med, lo, hi = _summary_band(counts, axis=0)
    
    # 観測の度数
    obs_counts = _histo_counts(y_obs, kmax).astype('float')
    obs_counts[obs_counts == 0] = np.nan

    # 描画
    xs = np.arange(kmax+1)
    # 観測値の散布図
    ax.scatter(xs, obs_counts, s=40, label='観測値', zorder=1)
    # 平均予測の散布図
    ms = np.where(med==0, 8, 30)  # 予測値==0の場合、小さいマーカーを描画
    ax.scatter(xs, med, s=ms, color='tab:red', alpha=0.5,
            label='$\\lambda$:平均予測（中央値）', zorder=2)
    # 95%CIの塗りつぶし
    ax.fill_between(
        xs, lo, hi, color='tomato', alpha=0.1, 
        label='$\\lambda$:平均予測（95%CI）', zorder=0)
    # 修飾
    ax.set_xlabel('種子数')
    ax.set_ylabel('度数（個体数）')
    ax.legend(loc='best')
    ax.set_title(title)
    # 戻り値：axes
    return ax

【実行結果】なし

描画用のデータを作成します。

# 観測値と平均生存種子数λの取得
y_obs = data.y.values
print('y_obs.shape:', y_obs.shape)

lam_samples = stack_by_id(idata.posterior.lam).astype('int')
print('lam_samples.shape:', lam_samples.shape)

【実行結果】なし

🍀🍀🍀

◼️ 施肥有無別の平均予測の描画
横軸が平均種子数、縦軸が各平均種子数の個体数（頻度）を描画します。

# 施肥有無別の平均種子数ごとの個体数（度数）の描画

# 設定
title_head = '平均種子数'
# 描画領域の設定
fig, axes = plt.subplots(1, 2, sharey=True, figsize=(10, 4), tight_layout=True)
# 施肥有無ごとにチャート描画を繰り返し処理
for treat, ax in zip(['C', 'T'], axes.flat):
    # 施肥有無を絞り込むインデックスの取得
    idx = data[data['f']==treat].index.values
    # グラフタイトルの設定
    title = f'{title_head} : 施肥 {treat}'
    # 生存種子数ごとの個体数（度数の描画）
    plot_ppc_frequency_band(y_obs[idx], lam_samples[:, idx], ax=ax, title=title);

【実行結果】
青い点が y の観測値、赤い点が平均予測の中央値、薄赤色が平均予測の 95% 信用区間です。
赤い小さい点は「その種子数に該当する個体数（中央値）が０」を示します。

【考察】
施肥あり T の方が平均種子数 $${0, 1}$$の個体数が多い印象です。
それ以外の種子数では施肥なし C と施肥あり T の違いが分かりません。

🍀🍀🍀

◼️ 植木鉢別の平均予測の描画

# 植木鉢別の平均種子数ごとの個体数（度数）の描画

# 設定
title_head = '平均種子数'

# 描画領域の設定
fig, axes = plt.subplots(5, 2, sharex=True, sharey=True, figsize=(10, 12),
                         tight_layout=True)
# 植木鉢ごとにチャート描画を繰り返し処理
for pot, ax in zip(Pot_cat, axes.T.flat):
    # 植木鉢を絞り込むインデックスの取得
    idx = data[data['pot']==pot].index.values
    # グラフタイトルの設定
    title = f'{title_head} : 植木鉢 {pot}'
    # 生存種子数ごとの個体数（度数の描画）
    plot_ppc_frequency_band(y_obs[idx], lam_samples[:, idx], ax=ax, title=title);

【実行結果】
左列が施肥なしの植木鉢、右列が施肥ありの植木鉢です。

【考察】
植木鉢ごとに種子数の大きさが異なっていることが分かります。
一方で、施肥なしの植木鉢と施肥ありの植木鉢の違いを見つけることはできません。

まとめ

今回のベイズ統計モデルをまとめます。

🔷 尤度
観測データはポアソン分布に従います。

🔷 リンク関数と線形予測子
リンク関数は対数です。
線形予測子は、切片、傾き×施肥あり、個体差（ランダム切片）、植木鉢差（ランダム切片）の和です。

🔷 事前分布
線形予測子の切片 $${\beta_1}$$ と傾き $${\beta_2}$$ の事前分布は、平均０、標準誤差 100 の「押しつぶされた正規分布」です。

個体差 $${r_i}$$ と植木鉢差 $${r_{pj(i)}}$$ の事前分布は、平均０、標準誤差 $${s , s_p}$$ 正規分布であり、階層事前分布です。

個体差のばらつきパラメータ $${s}$$ と植木鉢差のばらつきパラメータ $${s_p}$$ は階層事前分布に関するハイパーパラメータであり、事前分布は $${0}$$ から $${10^4}$$ までの一様分布です。

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

次回は 空間構造のある階層ベイズモデル を学びます。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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