「データ解析のための統計モデリング入門」をPythonで写経 Vol.16 ～ ７章「一般化線形混合モデル（GLMM）」③一般化線形混合モデルの挑戦～Pythonライブラリで実装

７章「一般化線形混合モデル（GLMM）」

書籍の著者　久保拓弥 先生

書籍「データ解析のための統計モデリング入門」７章「一般化線形混合モデル（GLMM）」の Python写経活動記録 です。 

書籍は第７章で 一般化線形混合モデル（GLMM） をやり切ります！

前回記事では ChatGPT との共同作業で GLMM コードを作り込みました。
処理速度が爆速になって、とても満足したのですが…

やはり、Python のライブラリを活用して、手軽に分析したいです。
けれども Python には GLMM を扱うメジャーなライブラリが無さそうで…。

そこで本記事では、Python にとって難関の GLMM について、なんとか Python の既存ライブラリを代替的に利用 することを試みます。
今回は、statsmodels、GPBoost、PyMC を使い比べてみます。

ChatGPT に頼りっきりの記事ですが、どうぞお読みください。
では書籍を開いて統計モデリングの旅に出かけましょう🚀

（注意）
今回の記事は趣味的な深堀りとコードです。
ご興味ない方はスルーしてくださって大丈夫です。

はじめに

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

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

準備

準備

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

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

# インポート

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

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

# gpboost
import gpboost as gpb
from patsy import dmatrix

# PyMC
import pymc as pm
import arviz as az

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

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

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

■ 統計モデル
GLMM を３つの Python ライブラリで実装します。

３つのライブラリは「二項分布・ロジットリンク関数・最尤推定を直接扱えない」ので、この記事で紹介するのは「代用的な利用方法」になります。

ちなみにChatGPTによると、３ライブラリの特徴は…

statsmodelsは軽量、GPBoostは高性能、PyMCは正確さを重視。

データの準備

前回記事と同じ例題データを利用します。
データの確認などの様子はぜひ前回記事をご覧ください。

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

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

【実行結果】
データの個数（標本サイズ）は 100 です。
植物の個体（個体数 100）に関する仮想の観測データです。
葉数 x の値ごとに 20 の個体を調査しています。

【変数の説明】
植物の個体ごとに採取した調査種子数 N、生存種子数 y（目的変数）、葉数 x です。

■ データの特徴

statsmodels

概要

statsmodels ライブラリの BinomialBayesMixedGLM クラスでモデリングします。
変分ベイズ（略：VB）でパラメータを推定します。

２点の制限事項があります。

このクラスは「二項分布・ロジットリンク関数」をモデリングできません。
一方で「ベルヌーイ分布・ロジットリンク関数」はモデリングできます。

そこで、代用的に「ベルヌーイ分布・ロジットリンク関数」でモデリングします。

データの変換

例題データの目的変数「生存種子数」が二値変数になるように、データを変換します。
１行のデータを 調査種子数 $${N=8}$$ の８行に増幅します。
８行のうち、生存種子数 $${y}$$ 個の行を $${y=1}$$、残りの $${N-y}$$ 個の行を $${y=0}$$ にします。

# BinomialBayesMixedGLM に乗せてみる

## データの整形

# 結果を格納するデータフレームの初期化
data2 = pd.DataFrame()

# １行毎にN行増やす処理を繰り返す
for idx, row in data.iterrows():
    # yは生存種子数分の1と残りを0に設定
    y = [1] * row.y + [0] * (row.N - row.y)
    # xとidはN行分の同じ値をリピート
    x = [row.x] * row.N
    id = [row.id] * row.N
    # 一時データフレームにy,x,idを統合
    temp_df = pd.DataFrame({'ex_index':idx , 'y': y, 'x': x, 'id': id})
    # 結果を格納するデータフレームに一時データフレームを行結合
    data2 = pd.concat([data2, temp_df], axis=0, ignore_index=True)      

# 結果の表示
data2

【実行結果】
変換前に 100 行あったデータは $${\times 8}$$ の 800 行になっています。

モデリング

ひとまずデータにモデルを当てはめて、パラメータを推定しましょう。

# statsmodels で GLMM 体験

# モデリング

# GLMMモデルの定義
model = sm.BinomialBayesMixedGLM.from_formula(
    formula='y ~ x',                  # 固定効果の式
    vc_formulas={'id': '0 + C(id)'},  # ランダム効果の式 C()で囲って質的変数表現
    data=data2,                       # データフレーム
    vcp_p=100,                        # ランダム効果の事前分布のスケールパラメータ
    fe_p=100,                         # 固定効果の事前分布のスケールパラメータ
)

# モデルの当てはめ ※変分ベイズ fit_vb() で推定
result = model.fit_vb()
result.summary()

【実行結果】
モデリングの結果が表示されました。

モデリング結果の分析

ChatGPT に結果を分析してもらいました。

📊 表の列の意味

statsmodels（$${\mathtt{fit\_vb}}$$） の結果と、テキストの R：$${\mathtt{glmmML}}$$ の結果を比較します。

📊 推定値の比較表（固定効果）

📊 推定値の比較表（ランダム効果）

※真値 3 にやや小さめ（変分ベイズ特有の収縮）

🔍 ポイント

• 固定効果（Intercept, x）は真値に非常に近い
  fe_p, vcp_p を大きくして事前をほぼフラットにした効果

• ランダム効果SDはRの結果と近いが、真値より小さい
  これは変分ベイズ（VB）の近似特性でよくある過小推定

• R の SE と VBの 事後 SD は意味が異なる

  • R：最尤推定における標準誤差

  • VB：事後分布の標準偏差（ただしVBは過小評価しやすい）

予測

■ 生存種子数 y の予測（個体差＝０）
生存種子数 y の予測を行って、テキスト p.161 図 7.10 の可視化を行います。
まず共通設定と描画関数の定義を行います。

# GLMM化したロジスティック回帰の推定にもとづく予測 p.161 図7.10

## ３モデル共通の設定
# 真の係数パラメータ
beta1_true, beta2_true = -4, 1
# 調査種子数
n_survey = 8
# x軸の値
x_val_a = np.linspace(1, 7, 100)

## このモデルの左のチャート用の共通設定
# 真の生存確率の生存種子数y
y_val_a = expit(beta1_true + beta2_true * x_val_a) * n_survey

## このモデルの右のチャート用の共通設定
# xiの値=4
xi = 4
# x=4のデータ個数
num_xi = len(data[data.x==xi])
# 生存種子数yごとの個体数
xi_data = data[data.x==xi]['y'].value_counts().to_frame()
# xi=4をリピートしたデータ
x_val_b = np.repeat(xi, n_survey + 1)
# 二項分布の横軸のyの値
y_val_b = list(range(n_survey + 1))


## ２チャートを描画する関数の定義
def plot_pred_twinchart(beta1_hat, beta2_hat, y_pred, r_sigma, y_mix=None):
    
	## 描画領域の設定
	fig, ax = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)

	## 左の(A)の描画
	# 観測値の散布図の描画（jitter処理込み）
	sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15,
				scatter_kws={'alpha': 0.5}, label='観測値', ax=ax[0])
	# 真の生存確率の赤い点線の描画
	ax[0].plot(x_val_a, y_val_a, color='tab:red', ls='--', label='真の生存確率')
	# GLMMの予測値(ランダム切片を除く平均予測値)の緑実線の描画
	ax[0].plot(x_val_a, y_pred, color='green', lw=2)
	# 修飾
	ax[0].set(xlim=(1.6, 6.4), xlabel='葉数 $x_i$', ylabel='生存種子数 $y_i$',
			  title='(A) 葉数と生存種子数の関係')

	## 右の(B)の描画
	# 計算
	# 正規分布に従うrの200個のパーセンタイル点の取得
	rs2 = stats.norm.ppf(
			q=np.linspace(0.001, 0.999, 200), loc=0, scale=r_sigma
	)
	# 生存確率logistic(β1+β2*x+r)の二項分布に従う生存種子数yの取得 shape=(200, 9)
	if y_mix is None:
		y_mix = np.array(
		[stats.binom.pmf(
				k=y_val_b, n=n_survey, 
				p=expit(beta1_hat + beta2_hat * x_val_b + r))
		for r in rs2]
		) * num_xi

	# 描画
	# x=4のときの生存種子数yの個体数の描画
	sns.scatterplot(data=xi_data, x='y', y='count', s=60, alpha=0.5, ax=ax[1])
	# x=4のときの種子数分布の描画
	ax[1].plot(y_val_b, y_mix.mean(axis=0), '-o', ms=7, color='green', mec='white')
	# 修飾
	ax[1].set(xlabel='生存種子数 $y_i$', ylabel='個体数',
			title=f'(B) 葉数 $x_i=${xi} での種子数分布')
	ax[1].set_ylim(0, 6.2)

	return fig, ax

【実行結果】なし

予測値を描画します。

# 描画処理の実行

## 設定
# パラメータ推定値の取得
beta1_hat, beta2_hat = result.fe_mean
# yの平均予測値の算出
y_pred = expit(beta1_hat + beta2_hat * x_val_a) * n_survey
# 個体差rの標準偏差sの推定値
r_sigma = np.exp(result.vcp_mean[0])

## 描画
fig, ax = plot_pred_twinchart(beta1_hat, beta2_hat, y_pred, r_sigma)

【実行結果】

【チャートの読み取り】
左のチャートは個体差 $${r=0}$$ として、葉数 x を用いた生存種子数 y の予測値（平均）を描画しています。
パラメータ真値による予測値の赤点線に対して、変分ベイズで推定したパラメータ値による予測値の緑実線は近い感じになっています。
右のチャートは生存種子数 y の分布です。
変分ベイズで推定したパラメータ値による予測値の緑実線は、観測値の分布（青い点）に近い感じになっています。

🍀🍀🍀

■ 生存種子数 y の予測（個体差ごと）
個体差 $${r_i}$$ を個別に反映した予測値（平均）を描画します。

変分ベイズの結果 result に対して $${\texttt{random\_effects()}}$$ メソッドで個体差（ランダム切片）を取得できます。

# ランダム切片rの表示
result.random_effects()

【実行結果】
平均 Mean と標準偏差 SD を取得できます。

描画関数を定義します。

# １チャートを描画する関数の定義
def plot_pred_chart(beta1_hat, beta2_hat, y_pred, r_mean):

	# 描画領域の設定
	fig, ax = plt.subplots(figsize=(6.4, 4.8))
	# 観測値の散布図の描画（jitter処理）
	sns.regplot(data=data, x='x', y='y', fit_reg=False, x_jitter=0.15,
				scatter_kws={'alpha': 0.5}, label='観測値', ax=ax)
	# 真の生存確率による予測値の赤い点線の描画
	ax.plot(x_val_a, y_val_a, color='tab:red', ls='--', label='真の生存確率')
	# GLMMの予測値(ランダム切片を除く平均予測値)の緑実線の描画
	ax.plot(x_val_a, y_pred, color='tab:green', lw=2,
		    label='予測値（パラメータ推定値）')
	# GLMMの予測値のオレンジ実線の描画
	for r in r_mean:
		ax.plot(x_val_a, expit(beta1_hat + beta2_hat * x_val_a + r) * n_survey,
				color='tab:orange', lw=0.5, alpha=0.3)
	ax.plot(0, 0, color='tab:orange', lw=1, alpha=0.5,
		    label='予測値（個体差含む）')

	# 修飾
	ax.set_title('(A) 葉数と生存種子数の関係')
	ax.set_xlabel('葉数 $x_i$', fontsize=14)
	ax.set_ylabel('生存種子数 $y_i$', fontsize=14)
	ax.set_xlim(1.6, 6.4)
	ax.legend(bbox_to_anchor=(1, 1), loc='upper left')

	return fig, ax

【実行結果】なし

描画します。

# 描画処理の実行

## 設定
# 個体差rの推定値（平均）
r_mean = result.random_effects().Mean

## 描画
fig, ax = plot_pred_chart(beta1_hat, beta2_hat, y_pred, r_mean)

【実行結果】
薄オレンジ線が 100 個体の個体差（ランダム切片）を反映した個体ごとの平均予測線です。

モデルとコードの補足説明

ChatGPT による補足説明です。難解でしたら飛ばし読みしてください。

1️⃣ 前提となる GLMM のモデル（なにを学習している？）

• 目的変数 $${y}$$ は０または１をとる二値変数です。

• 期待確率 $${p}$$ はロジットリンクで説明変数＋ランダム効果から決まります。

• ここで 引数 $${\mathtt{vc\_formulas={'id': '0 + C(id)'}}}$$ は、「id の各水準ごとに独立なランダム切片を1つずつ入れる」指定です。
  （0 は固定効果の切片を含めない、C(id) は id をカテゴリとして扱う、の意味）

• statsmodels の Bayes 混合モデルは 二項分布（実質的にはベルヌーイ分布）とポアソン分布をサポートし、推定法として Laplace近似（fit_map） と 変分ベイズ（fit_vb） を備えます。(statsmodels.org)

2️⃣ 尤度（ざっくり数式）

ランダム効果 $${u_g}$$ を積分消去した周辺尤度で学習します：

ベイズ枠なので、固定効果 $${\beta}$$ や「分散パラメータ（$${\log\sigma }$$ など）」にガウス系の事前分布を置き、事後分布（likelihood×prior）を近似します。
fe_p と vcp_p がその事前分布のスケール引数です（下で解説）。(statsmodels.org)

3️⃣ コードのここがポイント

model = sm.BinomialBayesMixedGLM.from_formula(
    formula='y ~ x',                  # 固定効果：切片＋x
    vc_formulas={'id': '0 + C(id)'},  # idごとのランダム切片
    data=data2,
    vcp_p=100,                        # ランダム効果“標準偏差の対数”への事前のSD
    fe_p=100,                         # 固定効果 β への事前のSD
)

• formula='y ~ x'
  y を x で説明（切片は自動で入る）。固定効果は $${\beta_0,\beta_1}$$。

• vc_formulas={'id': '0 + C(id)'}
  id の各レベルごとに独立なランダム切片 $${u_g}$$ を持つ指定。

• fe_p=100
  固定効果 $${\beta}$$ の事前分布（正規分布）の標準偏差。
  大きいほど弱情報＝ゆるい事前。(statsmodels.org)

• vcp_p=100
  “ランダム効果の標準偏差の対数”に対する事前分布（正規分布）の標準偏差。
  "標準偏差の対数" $${\log \sigma}$$ の事前分布は正規分布。
  → $${\sigma}$$ の事前分布は対数正規分布。
  これも値が大きいほど弱情報事前分布です。(statsmodels.org)

4️⃣ パラメータの推定方法（変分ベイズ fit_vb）

fit_vb() は平均場（mean-field）ガウス近似で事後分布を近似します。
ざっくり言うと、真の事後 $${p(\theta\mid y)}$$ を

という独立ガウスの積で近似し、そのELBO（Evidence Lower BOund）を最大化（≒ 事後分布への KL ダイバージェンスを最小化）します：

• こうして得られるのは、$${q}$$ の平均（≒推定値）と分散（≒不確実性）です。

• ドキュメントでも「因子化ガウス近似を探し、ELBO を最大化する」と説明されています。(statsmodels.org)

5️⃣ result.summary() の見方（なにが出る？）

• 固定効果の推定値（近似事後平均）と標準誤差（近似事後SD）。

• 分散成分（$${\sigma}$$ やそのログパラメータ）に対する近似事後の要約。

• 結果オブジェクト BayesMixedGLMResults には fe_mean（固定効果の事後平均）や、対応する SD などが入ります。(statsmodels.org)

6️⃣ 使うときのコツと注意点

• VB は分散をやや小さく見積もりがち（不確実性がタイトになりやすい）。重要な結論は fit_map() と照合するのがおすすめ。(statsmodels.org)

• 事前スケール（fe_p, vcp_p）は結果に影響します。よくわからなければ 大きめ（弱情報） から始め、診断しながら調整。(statsmodels.org)

• データが大きい・グループが多いとき、VB は速くて安定。一方、厳密なベイズ（MCMC）の代替ではない点は理解しておくと◎。

• vc_formulas はどのランダム効果を入れるかの核心。最初は「ランダム切片」から始め、必要に応じて拡張。

7️⃣ まとめ（このコードで得られるもの）

• モデル：ロジスティック混合モデル（固定効果 xx、id ごとのランダム切片）。

• 推定：変分ベイズにより 固定効果と分散パラメータの近似事後（平均・SD）。

• ハンドル：fe_p と vcp_p が事前の強さ、vc_formulas がランダム効果の構造。

• 実務：fit_vb()（速い・スケールしやすい）で当たりをつけ、重要な分析は fit_map() でも確認。

GPBoost

概要

GPboost ライブラリでモデリングします。
最尤推定でパラメータを推定します。
最尤推定中の積分にラプラス近似を用います。

こちらも２点の制限事項があります。

このライブラリはstatsmodelsと同様に、「二項分布・ロジットリンク関数」をモデリングできません。
一方で「ベルヌーイ分布・ロジットリンク関数」はモデリングできます。

そこで、代用的に「ベルヌーイ分布・ロジットリンク関数」でモデリングします。
さきほど変換したデータを利用します。

モデリング

変換データにモデルを当てはめて、パラメータを推定しましょう。

# gpboost で GLMM 体験

# モデリング

# 設定
group = data2.id                   # ランダム切片（個体差）
X = dmatrix('x', data=data2)       # 線形予測子（β1 + β2 x）

# GLMMモデルの定義 ※group_data：ランダム切片、likelihood：尤度＆リンク関数
gp_model = gpb.GPModel(group_data=group, likelihood='bernoulli_logit')

# フィッティング
gp_model.fit(y=data2['y'], X=X);

# 結果の表示
gp_model.summary();

【実行結果】
モデリングの結果が表示されました。

モデリング結果の分析

ChatGPT に結果を分析してもらいました。

📊 推定結果の整理

📝 解釈のポイント

1. 固定効果（Intercept, x）

   • 真値（-4, 1）に非常に近い。

   • 負の切片で、x=0では成功確率が非常に低い設定。

   • x が大きくなると成功確率は急速に上昇する。

2. ランダム効果（id）

   • 分散 ≈ 5.799 → 標準偏差 ≈ 2.408

   • この標準偏差は、群（id）ごとの平均的なロジット値のばらつきの大きさを表す。

   • 真値（3）よりやや小さいが、statsmodels（VB）の推定値とほぼ同じ。

   • このばらつきにより、群によって成功確率の基準水準が大きく異なることを示す。

3. モデル全体

   • 二項（0/1）データに対する一般化線形混合モデル（GLMM）のロジット版。

   • 推定はラプラス近似ベースで行われ、VBよりやや信頼できる不確実性推定が可能。

   • 固定効果はほぼ真値通り、ランダム効果はやや収縮傾向。

予測

■ 生存種子数 y の予測（個体差＝０）
生存種子数 y の予測を行って、テキスト p.161 図 7.10 の可視化を行います。

# 描画設定の実行

## 設定
# パラメータ推定値の取得
beta1_hat, beta2_hat = gp_model.get_coef(False)
# yの平均予測値の算出
y_pred = expit(beta1_hat + beta2_hat * x_val_a) * n_survey
# 個体差rの標準偏差sの推定値
r_sigma = np.sqrt(gp_model.get_cov_pars(format_pandas=False)[0])

## 描画
fig, ax = plot_pred_twinchart(beta1_hat, beta2_hat, y_pred, r_sigma)

【実行結果】

（参考：statsmodels の予測）

【チャートの考察】
statsmodels の予測値と近い結果になっています。

🍀🍀🍀

■ 生存種子数 y の予測（個体差ごと）
個体差 $${r_i}$$ を個別に反映した予測値（平均）を描画します。

GPBoost のモデル gp_model に対して 学習データのランダム効果を算出する $${\texttt{predict\_training\_data\_random\_effects()}}$$ メソッドで個体差（ランダム切片）を取得できます。

# ランダム切片rの推定
gp_model.predict_training_data_random_effects(predict_var=True).iloc[::8]

【実行結果】
平均と分散 var を取得できます。

描画します。

# 描画処理の実行

## 設定
# 個体差rの推定値
r_mean = gp_model.predict_training_data_random_effects(
    predict_var=False).iloc[::8].values

## 描画
fig, ax = plot_pred_chart(beta1_hat, beta2_hat, y_pred, r_mean)

【実行結果】
薄オレンジ線が 100 個体の個体差（ランダム切片）を反映した個体ごとの平均予測線です。

モデルとコードの補足説明

ChatGPT による補足説明です。難解でしたら飛ばし読みしてください。

1️⃣ 前提となる GLMM（モデルと数式）

• 目的変数は 0/1（二値変数）。リンク関数はロジット。

• 個体（id）ごとにランダム切片 $${u_{g(i)}}$$（平均0・分散 $${\tau^2}$$）を入れます。

GPBoost の GPModel は、こうしたグループ化ランダム効果（必要に応じて GP：ガウス過程も）を扱え、二項ロジットは likelihood="bernoulli_logit" で指定します。(gpboost.readthedocs.io)

2️⃣ コードの引数（なにを渡している？）

group = data2.id                 # ランダム切片のグループID
X = dmatrix('1 + x', data=data2) # 固定効果：[切片, x]

gp_model = gpb.GPModel(
    group_data=group,            # グループ化ランダム効果
    likelihood='bernoulli_logit' # 二項ロジット
)

gp_model.fit(y=data2['y'], X=X)
gp_model.summary()

• group_data=group
  グループIDのベクトル（または列を増やせば複数のランダム切片）。
  要するに「どの観測がどの個体（群）に属するか」を渡します。
  型は数値でも文字列でもOK。
  列数＝ランダム切片の数です。(gpboost.readthedocs.io)

• likelihood='bernoulli_logit'
  ベルヌーイ×ロジットリンクを指定。
  ほかに bernoulli_probit や poisson 等もあります。(gpboost.readthedocs.io)

• X
  固定効果のデザイン行列（ここでは切片＋説明変数 x）。
  fit(y, X=X) で一緒に推定します。

3️⃣ パラメータの推定方法（中身はどう動く？）

◆ なにを最適化？

• ざっくりいうと、周辺対数尤度

　を最大化します（= 負の周辺対数尤度を最小化）。

◆ 非ガウス（二項）への近似

• $${\int}$$ が解析的に解けないので、ラプラス近似で周辺化します。
  内部ではラプラス近似のモード探索にニュートン法が使われ、行列表現のトレースの近似にランダムベクトルを再利用するオプションなどが実装されています（高速化の工夫）。(gpboost.readthedocs.io)

◆ 具体的な最適化器

• 共分散（＝ランダム効果の分散など）や補助パラメータの推定には lbfgs / fisher_scoring / newton / gradient_descent / nelder_mead などが利用可能（デフォルトは多くのケースで L-BFGS）。
  固定効果の係数も WLS や L-BFGS 等で推定されます。収束判定や学習率のパラメータも用意されています。(gpboost.readthedocs.io)

4️⃣ 推定後にできること（便利ワザ）

• 要約の確認
  gp_model.summary() で、推定された係数や分散パラメータのサマリを確認。
  get_cov_pars() で共分散パラメータ一覧を DataFrame でも取れます。(gpboost.readthedocs.io)

• 確率 or 潜在スコアの予測
  gp_model.predict(...) は既定で応答（確率）を返し、predict_response=True でリンク前の線形予測子（固定＋ランダム）を返せます。

• 学習データのランダム効果をまとめて取得
  gp_model.predict_training_data_random_effects(predict_var=False) で、学習に使った各観測のランダム効果の推定値（事後平均／EBLUP に相当）が取れます。
  グループごとのユニーク値が欲しければ、同一グループの先頭行を抽出すればOK。
  必要なら分散も predict_var=True で。(gpboost.readthedocs.io)

• 未観測グループの扱い
  予測時に学習で見ていない group_data_pred を渡すと、そのグループのランダム効果は0（事前平均）に縮む扱いになります。(Cross Validated)

5️⃣ 実務でのコツ・注意

• X の設計：dmatrix('x') は切片を自動で入れます。

• グループIDの型：数値でも文字列でもOKですが、欠損がないか確認（形が揃っていれば内部でハンドリングされます）。(gpboost.readthedocs.io)

• 速度と安定性：非ガウスの大規模データや GP を入れた時には、内部でラプラス近似＋反復解法（必要に応じて）などが使われ、パラメータ optimizer_cov・maxit・delta_rel_conv などで調整できます。(gpboost.readthedocs.io)

• 目的の出力を選ぶ：

  • 「確率が欲しい」→ predict(predict_response=False)

  • 「線形予測子が欲しい」→ predict_response=True

  • 「学習データのランダム効果を見たい」→ predict_training_data_random_effects()。

6️⃣ ひとことで総括

このコードは、「二項ロジット＋グループ別ランダム切片」の標準的なGLMMを、GPBoost の GPModel で（近似）周辺尤度最大化として解いています。非ガウスの周辺化はラプラス近似、最適化は L-BFGS 等。
学習後は summary() で要約、predict で確率/潜在値、predict_training_data_random_effects でランダム効果推定値が得られます。(gpboost.readthedocs.io)

PyMC

概要

第８章以降で取り組む「ベイズ統計モデル」を先取りします。
PyMC ライブラリでモデリングします。

二項分布・ロジットリンク関数・ランダム切片をモデリングできます。
そこで、「二項分布・ロジットリンク関数」でモデリングします。
変換前の例題データを利用します。

なお、パラメータ推定方法は MCMC 法です。
GLMM の最尤推定とは異なります。

モデリング

さきほどの２つのモデルと比べて、モデリングのステップが長くなります。

■ 数式
ベイズ統計モデルの数式イメージです。

■ モデルの定義
PyMC のベイズモデルを定義します。

# モデルの定義

# coordsの設定
coords = {'data': data.index.values,     # データの識別子
          'id_index': data.id.unique(),  # 個体番号（１始まり）
}

# モデリング
with pm.Model(coords=coords) as model:
    
    ## dataの定義
    # 目的変数 Y
    Y = pm.Data('Y', value=data['y'].values, dims='data')
    # 説明変数 X
    X = pm.Data('X', value=data['x'].values, dims='data')
    # 調査種子数 N
    N = pm.Data('N', value=data['N'].values, dims='data')
    # 個体番号 id：配列の添字になるので０始まりの値に変換
    id = pm.Data('id', value=data['id'].values - 1, dims='data')

    ## 事前分布
    # パラメータβ1, β2：平均0, 標準偏差100000の正規分布
    beta1 = pm.Normal('beta1', mu=0, sigma=1e5) 
    beta2 = pm.Normal('beta2', mu=0, sigma=1e5)
    # ランダム切片効果 r：平均0, 標準偏差sigma_rの正規分布
    sigma_r = pm.HalfNormal('sigma_r', sigma=1e5)
    r = pm.Normal('r', mu=0, sigma=sigma_r, dims='id_index')

    ## 線形予測子: q=invlogit(線形予測子)
    q = pm.Deterministic('q', pm.invlogit(beta1 + beta2 * X + r[id]),
                         dims='data')

    ## 尤度関数: 二項分布
    obs = pm.Binomial('obs', n=N, p=q, observed=Y, dims='data')

【実行結果】なし

モデルの数式を可視化します。

# モデルの表示
model

【実行結果】

モデルを図示（グラフィカルモデル）します。

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

【実行結果】

■ MCMCの実行
パラメータ推定を実行します。

%%time
# MCMCサンプリング
with model:
    idata = pm.sample(
        draws=2000, tune=2000, chains=4, random_seed=42,
        nuts_sampler='nutpie'   # NUTSサンプラーを使用
    )

【実行結果】
処理時間はおよそ 43 秒です。
ダイバージェンス（Divergences）は検出されず、事後分布の探索は安定していると考えられます。

■ 収束確認
$${\widehat{R}}$$ 値がしきい値 $${1.01}$$ 以下になっていることを確認します。
パラメータごとに $${\widehat{R} > 1.01}$$ の個数をカウントし、０になっていればＯＫとします。

# r_hat>1.01の確認

# 設定
idata_in = idata         # idata名
threshold = 1.01         # しきい値

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

【実行結果】
すべてのパラメータがしきい値以下になっています。
収束の証拠になると思われます。

続いてトレースプロットを描画します。

# トレースプロットの表示
var_names = ['beta1', 'beta2', 'sigma_r']
pm.plot_trace(idata, var_names=var_names, compact=False, figsize=(7, 6))
plt.tight_layout();

【実行結果】
推定したパラメータ $${\beta_1, \beta_2, \sigma_r}$$ を描画しました。
左のチャートの分布図では４つの Chain の線が単峰でほぼ同じになっています（OK）。
右のチャートのサンプリングの軌跡は同じ範囲を乱雑に行き来しています（OK）。
収束の証拠になると思われます。

■ パラメータ推定値の表示

# 要約統計量の表示
az.summary(idata, var_names=var_names, hdi_prob=0.95, round_to=3)

【実行結果】
mean 列がパラメータ推定値の平均、sd がパラメータ推定値の標準誤差です。

モデリング結果の分析

ChatGPT に結果を分析してもらいました。

1️⃣ 数値表（`az.summary`の結果）

パラメータの解釈例です。

2️⃣ 収束診断

• r_hat 全て ≈ 1.00 → チェーン間の収束良好

• ess_bulk と ess_tail が十分大きい → サンプルの独立性も高い

3️⃣ 事後分布の形（左図）

• beta1, beta2, sigma_r いずれもなめらかな単峰性（混合モードなし）

• beta1 はやや広がり大きめ（sd ≈1）、sigma_r も幅広め

• beta2 は比較的狭く、真値周辺に集中

4️⃣ トレースプロット（右図）

• 4チェーンとも同じレンジを行き来 → 収束問題なし

• 飛び跳ねやドリフトはほぼ無し

5️⃣ モデル設定の良さ

• 広い事前（sigma=1e5）を設定しているので、ほぼ非情報的事前

• Binomial(N, p) で集計データを直接モデリングしており、展開不要

• ランダム効果も dims 指定で個体ごとに割り当てられ、拡張性あり

予測

■ 生存種子数 y の予測（個体差＝０）
生存種子数 y の予測を行って、テキスト p.161 図 7.10 の可視化を行います。
MCMCで個体差 $${r}$$ のサンプリングデータを取得しています。
生存種子数 y の平均予測値の 95% 信用区間を塗りつぶし描画しましょう。

# 描画設定の実行

## パラメータβ1,β2,sのMCMCサンプルの取得
beta1_samples = az.extract(idata.posterior)['beta1'].data
beta2_samples = az.extract(idata.posterior)['beta2'].data
sigma_r_samples = az.extract(idata.posterior)['sigma_r'].data

## 設定
# パラメータ推定値の取得
beta1_hat = beta1_samples.mean()
beta2_hat = beta2_samples.mean()
# yの平均予測値の算出
y_val_a_samples = (                        # 生存種子数の平均サンプル, ※r=0である
        expit(beta1_samples + np.outer(x_val_a, beta2_samples)) * n_survey)
y_pred = y_val_a_samples.mean(axis=1)
# 個体差rの標準偏差sの推定値
r_sigma = np.sqrt(gp_model.get_cov_pars(format_pandas=False)[0])

# 生存確率logistic(β1+β2*x+r)の二項分布に従う生存種子数yの取得 shape=(200, 9)
# 正規分布に従うrの200個のパーセンタイル点の取得
rs2 = stats.norm.ppf(
        q=np.linspace(0.001, 0.999, 200), loc=0, scale=sigma_r_samples.mean()
)
y_mix = np.array(
    [stats.binom.pmf(
        k=y_val_b, n=n_survey, 
        p=expit(beta1_samples.mean() + beta2_samples.mean() * x_val_b + r))
     for r in rs2]
    ) * num_xi


## 描画
# メインチャートの描画
fig, ax = plot_pred_twinchart(beta1_hat, beta2_hat, y_pred, r_sigma, y_mix)
# GLMMの予測値(ランダム切片を除く平均予測値)の95%区間緑塗りつぶしの描画
ax[0].fill_between(
        x_val_a, *np.quantile(y_val_a_samples, q=[0.025, 0.975], axis=1),
        color='green', alpha=0.1, zorder=0);

【実行結果】
左チャートの薄緑塗りつぶしが平均予測値の 95%ベイズ信用区間です。

（参考：statsmodels の予測）

（参考：GPBoost の予測）

【チャートの考察】
statsmodels、GPBoost の予測値と近い結果になっています。

🍀🍀🍀

■ 生存種子数 y の予測（個体差ごと）
個体差 $${r_i}$$ を個別に反映した予測値（平均）を描画します。
引き続き 95% ベイズ信用区間を重ね描きします。

# 描画処理の実行

## 設定
# 個体差rの推定値
r_mean = gp_model.predict_training_data_random_effects(
    predict_var=False).iloc[::8].values

## 描画
fig, ax = plot_pred_chart(beta1_hat, beta2_hat, y_pred, r_mean)

【実行結果】
薄オレンジ線が 100 個体の個体差（ランダム切片）を反映した個体ごとの平均予測線です。

モデルとコードの補足説明

ChatGPT による補足説明です。難解でしたら飛ばし読みしてください。

1️⃣ 前提となる GLMM（二項分布×ランダム切片）

各観測 $${i=1,\dots,n}$$・グループ $${g(i)}$$

2️⃣ コードの各行は何を表してる？

coords = {'data': data.index.values, 'id_index': data.id.unique()}

• ラベルを付けて配列次元を管理（結果が読みやすく）。

Y = pm.Data('Y', data['y'].values, dims='data')
X = pm.Data('X', data['x'].values, dims='data')
N = pm.Data('N', data['N'].values, dims='data')
id = pm.Data('id', data['id'].values - 1, dims='data')  # 0始まりに

• pm.Data は後で差し替え可能なデータ容器。
  dims で行方向に data ラベル。

• id は 0始まりの整数で、r[id] の添字に使うため。

beta1 = pm.Normal('beta1', mu=0, sigma=1e5)
beta2 = pm.Normal('beta2', mu=0, sigma=1e5)
sigma_r = pm.HalfNormal('sigma_r', sigma=1e5)   # σ_r > 0 の事前
r = pm.Normal('r', mu=0, sigma=sigma_r, dims='id_index')  # 各グループの u_g

• 事前分布：固定効果は正規、ランダム切片は $${u_g\sim \text{Normal}(0,\sigma_r^2)}$$、その $${\sigma_r}$$ は $${\text{HalfNormal}}$$。

• r はベクトル（次元 id_index）。

q = pm.Deterministic('q', pm.invlogit(beta1 + beta2 * X + r[id]), dims='data')

• $${\eta=\beta_1+\beta_2X+u_{g(i)}}$$ を逆ロジット（invlogit）で 成功確率 $${p_i=q_i}$$ に。

obs = pm.Binomial('obs', n=N, p=q, observed=Y, dims='data')

• 尤度：$${y_i\sim \text{Binomial}(N_i, q_i)}$$。

3️⃣ 推定はどうやる？（PyMC の中身）

• with model: idata = pm.sample(...) と呼ぶと、PyMC がNUTS（MCMC）で $${\{\beta_1,\beta_2,\sigma_r, u_g\}}$$ の事後分布をサンプリングします。

• 返るのはフルのサンプル分布（平均・SD だけでなく、相関や非対称性も持ったまま）。
  近似（VB/Laplace）と違い、不確実性が過小になりにくいのが利点です。

4️⃣ 実務メモ（より良く・より安全に）

(A) 事前は弱情報でも“現実的なスケール”に
1e5 は極端に広く、発散やダイバージェンスの原因になりがち。
まずは下のような弱情報が安全です（標準化前提ならもっと小さく）：

beta1 = pm.Normal('beta1', 0, 5)
beta2 = pm.Normal('beta2', 0, 5)
sigma_r = pm.HalfNormal('sigma_r', 1)  # 例

(B) 非中心化パラメータ化（divergence 減らし）
階層モデルの定番テク。

z = pm.Normal('z', 0, 1, dims='id_index')
r = pm.Deterministic('r', sigma_r * z, dims='id_index')

$${u_g = \sigma_r \times z_g}$$ にしてから r[id] を使うとサンプルが安定しやすいです。

(C) 数値安定：logit_p 引数の活用
PyMC の Binomial は p の代わりに logit_p=eta を直接渡せます：

eta = beta1 + beta2 * X + r[id]
obs = pm.Binomial('obs', n=N, logit_p=eta, observed=Y, dims='data')

invlogit を明示しなくてよく、極端な確率でも安定。

(D) サンプリング実行と要約

with model:
    idata = pm.sample(target_accept=0.9, chains=4, cores=4, random_seed=0)
    ppc = pm.sample_posterior_predictive(idata)
# 結果閲覧（ArviZ）
import arviz as az
az.summary(idata, var_names=['beta1','beta2','sigma_r'])

(E) 予測・新しい id への外挿

• 既存 id の予測 → pm.set_data で X/N/id を差替 → sample_posterior_predictive。

• 新規 id を予測したいなら、z_new ~ Normal(0,1) を生成し σ_r*z_new を線形予測子に足す、という“階層的予測”が必要です（既知 id の r[id] を使い回すと流用になってしまう）。

5️⃣ ひとことで

この PyMC モデルは、

• 二項GLMM（ロジットリンク、ランダム切片）を

• Bayesian MCMC（NUTS）でフル事後として推定し、

• coords と dims でグループ別の $${u_g}$$ をきれいに扱える設計です。

まとめ

３つのモデルの比較でまとめとします。
こちらも ChatGPT にまとめてもらいました。

3つの推定結果（statsmodels: 変分ベイズ、GPBoost: ラプラス近似、PyMC: MCMC）を、真値との比較という観点で整理します。

1️⃣ 固定効果（Intercept, x）

$$
\begin{array}{lrrrr}
モデル & \text{Intercept} & 真値との差 & x & 真値との差 \\
\hline
\text{statsmodels (VB)} & -4.061 & +0.061 & 0.975 & -0.025 \\
\text{GPBoost (Laplace)} & -4.190 & -0.190 & 1.005 & +0.005 \\
\text{PyMC (MCMC)} & -4.336 & -0.336 & 1.039 & +0.039 \\
\end{array}
$$

【コメント】

$$
\begin{array}{ll}
モデル & コメント \\
\hline
\text{statsmodels (VB)} & 切片はほぼ真値、x も良好 \\
\text{GPBoost (Laplace)} & x はほぼ真値、切片はやや過小 \\
\text{PyMC (MCMC)} & 傾きは良好、切片は3つの中で最も過小寄り \\
\end{array}
$$

【評価】

• 固定効果は3モデルとも高精度。

• 切片は statsmodels が最も近く、傾きは GPBoost と PyMC がほぼ真値。

🍀🍀🍀

2️⃣ ランダム切片の標準偏差（σ）

$$
\begin{array}{lrrr}
モデル & 推定値 & 真値 & 差 \\
\hline
\text{statsmodels (VB)} & 2.410 & 3.0 & -0.590 \\
\text{GPBoost (Laplace)} & 2.408 & 3.0 & -0.592 \\
\text{PyMC (MCMC)} & 2.653 & 3.0 & -0.347 \\
\end{array}
$$

【コメント】

$$
\begin{array}{ll}
モデル & コメント \\
\hline
\text{statsmodels (VB)} & 明確な過小推定 \\
\text{GPBoost (Laplace)} & \text{statsmodels} とほぼ同じ過小 \\
\text{PyMC (MCMC)} & 過小だが最も真値に近い \\
\end{array}
$$

【評価】

• VB と Laplace はほぼ同じ結果で、標準偏差を約20%小さく推定。

• MCMC（PyMC）は過小傾向がやや弱い。これは近似法よりも分散の過小評価が少ないため。

🍀🍀🍀

3️⃣ 推定方法の影響

• statsmodels (VB)：

  • 計算は速いが、変分ベイズ特有の「分散の過小評価」がσ推定に影響。

  • 固定効果は安定して良好。

• GPBoost (Laplace)：

  • 固定効果は非常に良好。

  • σはVBと同じく過小推定傾向。

  • 大規模データや勾配ブースティングとの併用が強み。

• PyMC (MCMC)：

  • 計算負荷は最も高いが、事後分布を忠実にサンプリングできる。

  • σの推定が最も真値に近く、不確実性評価も信頼度高め。

🍀🍀🍀

4️⃣ 総合コメント

• 固定効果の推定精度は3モデルとも非常に高く、差はごくわずか。

• ランダム効果の分散推定では、近似法（VB, Laplace）がMCMCより過小になりやすい傾向が確認できる。

• 不確実性の定量化（Post. SDやHDIなど）を重視するならPyMC、計算効率重視ならstatsmodelsやGPBoostが有力。

5️⃣ ３ツールのざっくりまとめ

• 計算速さ重視 → statsmodels

• 中規模で実務的に安定 → GPBoost

• 厳密さ・不確実性重視 → PyMC

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

次回からベイズ統計モデルの章に進みます。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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