「ベイズ統計モデリングによるデータ分析入門」をPythonとStanで写経 ～ Vol.17 階層ベイズモデルと一般化線形混合モデルの基本
書籍の著者 馬場真哉 先生
この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第４部第１章「階層ベイズモデルと一般化線形混合モデルの基本」の Python 写経活動記録です。
書籍の第４部から、いよいよ「応用編」に進みます！
第４部は一般化線形混合モデル（通称：GLMM）の階層ベイズモデルです。
今回は過分散に対処するための ポアソンGLMM を題材にして階層ベイズモデルの第一歩を踏み出します👣
では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀
はじめに
このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。
テキストの紹介と引用表記はリンク先の記事に掲載しています。
準備
■ 記事の範囲
この記事はテキスト第４部第１章の以下の節を取り扱います。
1.3 分析の準備
1.4 通常のポアソン回帰モデルを適用した結果
1.8 GLMMのためのStanファイルの実装
1.9 MCMCの実行
1.10 brmsによるGLMMの推定
ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部および本章の理論面の記載をご覧いただき、本記事との繋がりをご確認下さいませ。
■ コード記述法
Jupyter Notebook 形式でコードを記述します。
■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。
■ ライブラリのインポート
この記事で用いるライブラリをインポートします。
# インポート
# 数値計算
import numpy as np
import pandas as pd
# ベイズ統計モデリング
from cmdstanpy import CmdStanModel # stan
import arviz as az # 分析・可視化
# デザイン行列
from patsy import dmatrices
# 統計処理
import scipy.stats as stats
# ユーティリティ
import os
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme() # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'第１章 階層ベイズモデルと一般化線形混合モデルの基本
一般化線形混合モデルと階層ベイズモデルとは？
① 一般化線形混合モデル
一般化線形混合モデル（Generalized Linear Mixed Model：GLMM）を一言で表してみます。
GLMMは、GLMの線形予測子にランダム切片やランダム傾きを追加し、それらを多くの場合、平均0・分散 $${\sigma^2}$$（または共分散行列 $${\Sigma}$$）の正規分布に従うと仮定したモデルです。
前回までの記事は GLM：一般化線形モデルに焦点を当てました。
今回は「混合」が加わった GLMM を学びます。
GLM に混合するのは「ランダム効果」です。
上の文章のランダム切片、ランダム傾きがランダム効果を表現しています。
ちなみにランダム効果を除くと、残りは「固定効果」です。
固定効果とランダム効果が「混合している」モデルなので一般化線形混合モデルなのです。
本記事で学ぶ「ポアソン回帰を拡張したGLMM」の線形予測子を題材にしてランダム効果を見てみましょう。
◆ ポアソン回帰の線形予測子
$$
\text{log}(\lambda_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}
$$
◆ 上の式にランダム切片 $${r_i}$$ を追加した GLMM の線形予測子
$$
\text{log}(\lambda_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_i
$$
◆ 固定効果とランダム効果を区別すると…
$$
\text{log}(\lambda_i) =
\underbrace{\underbrace{\beta_0} + \underbrace{\beta_1} x_{i1} + \underbrace{\beta_2} x_{i2}}_{固定効果}
+ \underbrace{r_i}_{\substack{ランダム効果\\（ランダム切片）}}
$$
ここではパラメータを固定効果と呼びましたが、文脈によっては、固定効果のパラメータにかかる線形予測子の部分 $${\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2}}$$ を固定効果と呼ぶこともあるそうです。
ところで、一般化線形混合モデルはベイズでなくてもモデリングできます。今回のポアソンGLMM は「GPBoost ライブラリ」で分析できます。
最後尾の「アディショナルタイム」で GPBoost によるポアソンGLMM を実験しましょう。
🔵
② 階層ベイズモデル
テキストの階層ベイズモデルの説明をお借りします。
階層ベイズモデルは、その名の通り階層構造を持つモデルです。上位の層の確率変数の実現値が、下位の層の確率分布の母数（確率分布のパラメータ）となります。
パラメータが階層構造になっているモデルですね！
今回取り組む階層ベイズモデルに照らして確認します。
$$
\begin{align*}
\sigma_r &\sim \text{HalfNormal}((1e5)^2) \\
r_i &\sim \text{Normal}(0, \sigma_r^2) \\
\beta_0, \beta_1, \beta_2 &\sim \text{Normal}(0, (1e5)^2) \\
\log(\lambda_i) &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_i \\
y_i &\sim \text{Poisson}(\lambda_i) \\
\end{align*}
$$
上位の層の確率変数 $${\sigma_r}$$ の実現値が、下位の層の確率変数 $${r_i}$$ が従う正規分布のパラメータ $${\sigma_r^2}$$ になっています！
📝 ちなみに情報
$${\sigma_r}$$ の事前分布と $${r_i}$$ の事前分布の関係も階層になっています。
この階層構造の事前分布に関して、$${r_i}$$ の事前分布を「階層事前分布」、$${\sigma_r}$$ を「ハイパーパラメータ」、$${\sigma_r}$$ の事前分布を「ハイパー事前分布」と呼ぶこともあります。
ではテキストの実装に進みましょう。
🔵🔵🔵
データの読み込みと外観の確認
テキスト 1.3 節に相当します。
テキストの仮想の魚の釣獲尾数データを引用いたします。
csv ファイルを pandas のデータフレーム形式で変数 fish_num_climate_2 に読み込みます。
# p.246 分析対象データの読み込み
# ファイルの読み込み
fish_num_climate_2 = pd.read_csv('./data/4-1-1-fish-num-2.csv')
# id列を数値ではなくcategoryとして扱う（後続処理で活かせるか！？）
fish_num_climate_2 = fish_num_climate_2.astype({'id': 'category'})
# 結果の表示
print('fish_num_climate_2.shape: ', fish_num_climate_2.shape)
fish_num_climate_2.head(3)【実行結果】
標本サイズ 100のデータです。
fish_num：釣った魚の数（釣獲尾数）
weather：天気の種類（晴れ：sunny、曇り：cloudy）
temperature：気温
id：データの連番（1～100）
天気ごとのデータ件数をカウントします。
pandas の value_counts メソッドを利用します。
# 天気ごとのデータ件数
fish_num_climate_2['weather'].value_counts().to_frame()【実行結果】
曇り、晴れそれぞれ 50 件のデータがあります。
データの要約統計量を確認します。
まずは全体（量的変数のみ）です。
# データの要約統計量
fish_num_climate_2.describe().T.round(2)【実行結果】
釣獲尾数は平均 2.7、最小値 0、最大値 15 です。
気温は平均 14.5、最小値 0.3、最大値 29.8 です。
次は天気別釣獲尾数の要約統計量です。
# 天気別の釣獲尾数の要約統計量
fish_num_climate_2.groupby(['weather'])['fish_num'].describe().round(2)【実行結果】
曇りの方が大きい感じがします。
データを可視化しましょう。
seaborn の scatterplot を利用して散布図を描画します。
# 釣獲尾数と気温・天気の散布図
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 散布図の描画
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}
sns.scatterplot(data=fish_num_climate_2, x='temperature', y='fish_num',
hue='weather', palette=colors)
# 修飾
plt.title('気温・天気と釣獲尾数の散布図', fontsize=14)
plt.xlabel('気温', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12)
plt.xticks(range(0, 31, 2))
plt.legend(loc='upper left', title='天気');【実行結果】
気温の上昇により釣獲尾数の増加する関係が見られます。
ただし、個々のデータ点の「ばらつき」が激しいので、増加傾向の説明は一筋縄ではいかないのです。
🔵🔵🔵
ポアソン分布と過分散
過分散について、少々掘り下げてみましょう。
今回モデルでは、釣獲尾数 $${y_i}$$ がポアソン分布に従うと仮定します。
ポアソン分布に従う確率変数の期待値と分散は、ポアソン分布の平均パラメータ $${\lambda}$$ そのものです。
$$
E[y] = \lambda, \quad V[y] = \lambda \\
$$
GLMやGLMMではデータ１件ごとに平均パラメータ $${\lambda_i}$$ を推定しますが、ここでは、データ全体の平均と天気別の平均でざっくり過分散を体感したいと思います。
釣獲尾数全体の標本平均と標本不偏分散を算出します。
# 釣獲尾数の標本平均と標本不偏分散
fish_num_climate_2[['fish_num']].agg(['mean', 'var']).T.round(2)【実行結果】
標本平均 2.66、標本不偏分散 11.46 です。
ポアソン分布で予測される分散は平均と同じになって欲しいところ…
データでは分散が平均の 4.3 倍になっています！
予測される分散よりもデータの分散が大きい状態を過分散といいます。
このデータはまさに過分散の予感がします。
天気別の釣獲尾数の標本平均と標本不偏分散を算出します。
# 天気別の釣獲尾数の標本平均と標本不偏分散
fish_num_climate_2.groupby('weather')[['fish_num']].agg(['mean', 'var']).round(2)【実行結果】
曇り・晴れの両方で、分散が平均よりも大きい過分散を感じられます。
天気別の過分散の状況を可視化してみましょう。
「データのヒストグラム」（観測ベース）と「データの平均をパラメータにしたポアソン分布の確率質量関数」（理論ベース）を重ね描きします。
# 天気別釣獲尾数のヒストグラムとポアソン分布の確率質量関数（理論値）の描画
# 設定
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'} # 天気の色
x_line = np.arange(0, fish_num_climate_2['fish_num'].max()+2) # bins, x軸の値
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(10, 4), tight_layout=True)
# 天気ごとにヒストグラム描画とポアソン分布確率質量関数の描画を繰り返し処理
for weather, ax in zip(['cloudy', 'sunny'], axes.flat):
# 天気に合致するデータを抽出
tmp_df = fish_num_climate_2[fish_num_climate_2['weather'] == weather]
# 天気の釣獲尾数の平均値を算出
fish_mean = tmp_df['fish_num'].mean()
# 天気別釣獲尾数のヒストグラムの描画
ax.hist(
tmp_df['fish_num'], bins=x_line-0.5, density=True,
color=colors[weather], alpha=0.5,
label=f'ヒストグラム'
)
# 天気別ポアソン分布の確率質量関数の描画
pois_pmf = stats.poisson.pmf(k=x_line, mu=fish_mean)
ax.plot(
x_line[:-1], pois_pmf[:-1], '-o', color=colors[weather],
label=f'Poisson(λ={fish_mean:.2f})'
)
# 修飾
ax.set_title(f'天気 {weather} のヒストグラムとポアソン分布')
ax.set_xlabel('釣獲尾数', fontsize=12)
ax.set_ylabel('確率', fontsize=12)
ax.legend()【実行結果】
棒グラフがデータのヒストグラム、折れ線グラフがポアソン分布の確率質量関数です。
データのヒストグラムとポアソン分布の折れ線グラフは以下の点でかけ離れている印象です。
・ポアソン分布と比べて０匹データが非常に多い
・ポアソン分布の右裾が確率０の部分で釣れたデータが点在する
生成AIにこのグラフを過分散の文脈で読み解いてもらいましょう。
可視化結果には典型的な 「過分散」 の特徴が出ています。
理論分布（ポアソン分布）と経験分布（実際のヒストグラム）を比較した際、以下の3つの乖離が見られる場合、それは過分散であると言えるでしょう。
1. 「0」の頻度が理論値より明らかに高い
ポアソン分布は、平均値（$${\lambda}$$）が決まると「0」が出る確率も自動的に決まってしまいます。
今回のケース：晴れ（$${1.94}$$）も曇り（$${3.38}$$）も、理論線が示す「0」の確率よりも、実際の青や赤のバーの方がずっと高くなっています。
解釈：「全然釣れない（0匹）」という結果が、理論上の確率を超えて頻発している状態です。
2. 「ピーク（山）」が理論値ほど尖っていない
ポアソン分布は平均付近にギュッと確率が集中する性質があります。
今回のケース：理論線（点つき折れ線）は平均付近でこんもりと高い山を作っていますが、実際のヒストグラムはそこまで高くありません。
解釈：平均付近のデータ密度が薄く、その分が他の数値へ「散らばって」しまっています。これが「分散が大きい」状態の正体です。
3. 「長い裾」が理論値を超えて伸びている
ポアソン分布は、平均から離れるほど確率が急激に（指数関数的に）ゼロへ収束します。
今回のケース：理論線は 7〜8 匹あたりで「確率0」に張り付いていますが、ヒストグラムは 10 匹、12 匹、15 匹…と、遥か右側まで「島」のようにデータが存在しています。
解釈：理論上は「まず起こりえない」はずの大漁が、現実には時々起きている。この「裾の長さ」が分散を跳ね上げています。
結論としてどう認識するか
一言でいうと、「真ん中（平均付近）が凹んで、両端（0と大漁）が膨らんでいる」 状態です。
ポアソン分布という型に、無理やりデータを押し込もうとしたけれど、両端からはみ出してしまった
可視化された図で、折れ線グラフがヒストグラムの「0のバー」と「右側の遠くのバー」を拾いきれていないことこそが、過分散の視覚的な定義そのものです。
この「はみ出した部分」を、データ1件単位のランダム効果が「個別のゆらぎ」として拾い上げてくれることで、モデルはようやくこのヒストグラムの形に寄り添うことができるようになります。
この先、テキストは過分散対策の目的で一般化線形混合モデルのベイズ統計モデリングに進みます。
🔵🔵🔵
ベイズモデリング by Stan
テキスト 1.4、1.8、1.9、1.10 節に相当します。
デザイン行列を作成して CmdStanPy でモデリングします。
🚀
「通常のポアソン回帰モデル」編
テキスト 1.4 節では brms でポアソン回帰モデルを実装しています。
この記事は CmdStanPy でコツコツ実装します。
① モデルの概要
デザイン行列と係数ベクトルを用いて、次のモデルを実装します。
$$
\begin{align*}
\bm Y &\sim \text{Poisson}(\bm \lambda) \\
\log(\bm \lambda) &= \bm{X \beta}
\end{align*}
$$
目的変数 $${\bm Y}$$（fish_num）はポアソン分布に従うと仮定しています。
ポアソン分布の平均パラメータベクトルの対数 $${\log(\bm \lambda)}$$ は、デザイン行列 $${\bm X}$$ と係数ベクトル $${\bm \beta}$$ の積で示される線形予測子と等しいです。
係数ベクトル $${\bm \beta}$$ には事前分布を明示的に設定しません。
② Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
事後予測サンプル（サンプル内予測）を取得する目的で、generated quantities ブロックで 予測値 Y_pred を設定します。
📑ファイル名：4-1-1-glm-pois-design-matrix.stan
data {
int N; // 標本サイズ
int K; // デザイン行列の列数（説明変数の数＋１）
array[N] int Y; // 目的変数
matrix[N, K] X; // 説明変数
}
parameters {
vector[K] b; // 切片を含む係数ベクトル
}
transformed parameters {
vector[N] lam = X * b;
}
model {
Y ~ poisson_log(lam);
}
generated quantities {
array[N] int Y_pred; // Posterior Predictive Check用
for (n in 1:N) {
Y_pred[n] = poisson_log_rng(lam[n]);
}
}【実行結果】なし
③ データセットの作成
patsy ライブラリを用いてデザイン行列等を作成します。
# デザイン行列の作成
# formula構文の設定
formula_pois = 'fish_num ~ weather + temperature'
# 目的変数Y, デザイン行列（説明変数）Xの作成
Y, X = formula_lm = dmatrices(formula_pois, fish_num_climate_2,
return_type='dataframe')
# int型の設定
X = X.astype({'Intercept': int, 'weather[T.sunny]': int})
Y = Y.astype(int)
# デザイン行列の先頭5行の表示
X.head()【実行結果】
デザイン行列は定数項、天気のダミー変数（晴れ）、気温で構成されます。
# 目的変数の先頭5行の表示
Y.head()【実行結果】
目的変数も作ってくれました。
標本サイズ・説明変数の数を算出して、Stan に渡すデータセットを辞書にまとめます。
# データセットの準備
# サンプルサイズ、デザイン行列の列数（説明変数の数＋１）
N, K = X.shape
# 辞書にまとめる ※Y:pd.Series, X:pd.DataFrame
data_dict_design = dict(N=N, K=K, Y=Y['fish_num'], X=X)【実行結果】なし
④ モデルのコンパイル
モデルのコンパイルを実行します。
%%time
# モデルのコンパイル
# stanプログラムファイルのパス指定
stan_file = '4-1-1-glm-pois-design-matrix.stan'
current_dir = os.path.abspath(os.getcwd())
stan_path = os.path.join(current_dir, 'stan', stan_file)
# モデルオブジェクトの作成（exeの作成）
model_glm = CmdStanModel(stan_file=stan_path) # stanのpathを設定【実行結果】
MCMC の準備が整いました！
⑤ MCMC の実行
MCMCを実行しましょう。
%%time
# MCMCの実行
fit_glm = model_glm.sample(
data=data_dict_design, # 対象データ
seed=1, # 乱数の種
sig_figs=18, # 出力CSV等に適用する数値精度
)【実行結果】
⑥ 収束確認
収束の確認をします。
診断メソッド diagnose を利用します。
# 事後分布の診断
print(fit_glm.diagnose())【実行結果】（１行目のファイルパスは記載省略）
問題は検出されませんでした（no problems detected.）。
CmdStanPy の 出力結果 fit から MCMCサンプルの要約統計量を把握しましょう。
# 結果の表示
(
fit_glm.summary(percentiles=[2.5, 50, 97.5])
.loc[['b[1]', 'b[2]', 'b[3]'], :]
.round(2)
)【実行結果】
$${\widehat{R}}$$（R_hat）、有効サンプル数（N_Eff）に問題はなさそうです。
b[1] は切片、b[2] は晴れ、b[3] は気温の係数に対応しています。
fit を arviz の idata に変換します。
事後予測サンプルデータを idata に適切に変換するために、呪文のような引数を３つ設定しています（コードのコメント参照）。
# arvizのidataに変換
idata_glm = az.from_cmdstanpy(
posterior=fit_glm,
posterior_predictive='Y_pred', # 事後予測 posterior_predictive グループを指定
observed_data={'Y': Y}, # 観測値 observed_data グループを指定
dtypes={'Y_pred': int}, # Y_predがfloat型で変換されないように int 指定
log_likelihood='lp__',
)
idata_glm【実行結果】
事後予測 Y_pred を格納した posterior_predictive グループを表示します。
トレースプロットを描画します。
# トレースプロットの描画
var_names = ['b']
az.plot_trace(idata_glm, var_names=var_names, compact=False,
backend_kwargs={'tight_layout': True});【実行結果】
左側のチャートの４本の Chain はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。
以上のチェックに基づいて、収束していると考えましょう。
⑦ 可視化で過分散の様子を確認
テキストの図 4.1.1 に相当します。
99% 予測区間を描画して、予測区間からはみ出る観測値のデータ点を眺めに行きましょう。
テキストの図に合わせるため、先ほど生成した事後予測サンプルデータを使わず、このコード内で事後予測サンプルを作成します。
# p.247 図4.1.1 通常のポアソン回帰の結果
## 予測釣獲尾数データの作成
# MCMCサンプルからβ0, β1, β2を取り出し
intercept, w_sunny, temper = az.extract(idata_glm.posterior).b.values
# x軸の値
x_val = np.linspace(X['temperature'].min(), X['temperature'].max(), 121)
# MCMCサンプルから曇り、晴れの平均釣獲尾数λを算出
cloudy_exp_lam = np.exp(intercept + w_sunny*0 + np.outer(x_val, temper)).T
sunny_exp_lam = np.exp(intercept + w_sunny*1 + np.outer(x_val, temper)).T
# 曇り、晴れの平均釣獲尾数λから予測値を算出：ポアソン分布乱数を生成
cloudy_pred = stats.poisson.rvs(mu=cloudy_exp_lam, random_state=123)
sunny_pred = stats.poisson.rvs(mu=sunny_exp_lam, random_state=123)
# 平均釣獲尾数データと描画色を天気別の辞書に格納
fishes = {'cloudy': cloudy_pred, 'sunny': sunny_pred}
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}
## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=fish_num_climate_2, x='temperature', y='fish_num',
hue='weather', palette=colors)
# 天気ごとに平均釣獲尾数の平均値と99%予測区間の描画を繰り返し処理
for weather in fishes.keys():
# 平均釣獲尾数の中央値の描画
plt.plot(x_val, np.median(fishes[weather], axis=0), color=colors[weather])
# 平均釣獲尾数の99%予測区間の塗りつぶし描画
ci99 = np.quantile(fishes[weather], q=[0.005, 0.995], axis=0).T
plt.fill_between(x_val, ci99[:, 0], ci99[:, 1], color=colors[weather],
alpha=0.2)
# 修飾
plt.title('釣獲尾数の99%予測区間')
plt.xlabel('気温', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12)
plt.legend(loc='upper left', title='天気');【実行結果】
薄赤塗りは曇り、薄青塗りは晴れの 99% 予測区間です。
いくつかのデータ点は 99% 予測区間の範囲外にあります。
99% 予測区間に収まらない状況をテキストは、データのばらつきが想定された分散よりも大きい「過分散」の状態だと説明しています。
そして GLMM 降臨✨️
個々のデータ点が持つ「特別なばらつき」をランダム切片で表現します！
🔵🔵🔵
「ポアソン GLMM モデル」編
テキスト 1.8、1.9 節の Stan 実装を CmdStanPy でコツコツ実装します。
① 過分散対策の立役者「ランダム切片」
冒頭の階層ベイズモデルの説明で示したモデル数式を再掲いたします。
$$
\begin{align*}
\sigma_r &\sim \text{HalfNormal}((1e5)^2) \\
r_i &\sim \text{Normal}(0, \sigma_r^2) \\
\beta_0, \beta_1, \beta_2 &\sim \text{Normal}(0, (1e5)^2) \\
\log(\lambda_i) &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_i \\
y_i &\sim \text{Poisson}(\lambda_i) \\
\end{align*}
$$
線形予測子 $${\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_i}$$ の $${r_i}$$ がランダム切片です。
添字 $${i}$$ が付くとおり、データ点１件ごとに推定されます。
このランダム切片 $${r_i}$$ は、ポアソン分布が持つ理論上のばらつき（誤差構造）を超えた個別のばらつきの表現を担当します。
過分散をランダム切片が捉えるのです！
そして、ランダム切片の事前分布に使われる $${\sigma}$$ パラメータは１つだけであり、全データ点で共通で使います。
唯一の $${\sigma}$$ パラメータでランダム切片のばらつきをコントロールしているかのようです。
② モデルの概要
デザイン行列、係数ベクトル、ランダム切片を用いて、次のモデルを実装します。
$$
\begin{align*}
\bm Y &\sim \text{Poisson}(\bm \lambda) \\
\log(\bm \lambda) &= \bm{X \beta} + \bm r \\
\bm r &\sim \text{Normal}(0, \sigma_r^2) \\
\end{align*}
$$
目的変数 $${\bm Y}$$（fish_num）はポアソン分布に従うと仮定しています。
ポアソン分布の平均パラメータベクトルの対数 $${\log(\bm \lambda)}$$ は、デザイン行列 $${\bm X}$$ と係数ベクトル $${\bm \beta}$$ の積にランダム切片 $${\bm r}$$ を加えた線形予測子と等しいです。
ランダム切片 $${\bm r}$$ は正規分布に従うと仮定しています。
係数ベクトル $${\bm \beta}$$ と ランダム切片の標準偏差パラメータ $${\sigma_r^2}$$ には事前分布を明示的に設定しません。
③ Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
デザイン行列を用いる書き方を採用しています。
事後予測サンプル（サンプル内予測）を取得する目的で、generated quantities ブロックで 予測値 Y_pred を設定します。
📑ファイル名：4-1-1-glmm-pois-design-matrix.stan
data {
int N; // 標本サイズ
int K; // デザイン行列の列数（説明変数の数＋１）
array[N] int Y; // 目的変数
matrix[N, K] X; // 説明変数
}
parameters {
vector[K] b; // 切片を含む係数ベクトル
vector[N] r; // ランダム効果
real<lower=0> sigma_r; // ランダム効果の標準偏差
}
transformed parameters {
vector[N] lam = X * b + r; // 線形予測子
}
model {
r ~ normal(0, sigma_r);
Y ~ poisson_log(lam);
}
generated quantities {
array[N] int Y_pred; // Posterior Predictive Check用
for (n in 1:N) {
Y_pred[n] = poisson_log_rng(lam[n]);
}
}
【実行結果】なし
④ データセットの作成
GLM モデル構築時に作成したデザイン行列や設定辞書を再利用します。
⑤ モデルのコンパイル
モデルのコンパイルを実行します。
%%time
# モデルのコンパイル
# stanプログラムファイルのパス指定
stan_file = '4-1-1-glmm-pois-design-matrix.stan'
current_dir = os.path.abspath(os.getcwd())
stan_path = os.path.join(current_dir, 'stan', stan_file)
# モデルオブジェクトの作成（exeの作成）
model_glmm = CmdStanModel(stan_file=stan_path) # stanのpathを設定【実行結果】
MCMC の準備が整いました！
⑥ MCMC の実行
MCMCを実行しましょう。
%%time
# p.251 MCMCの実行
fit_glmm = model_glmm.sample(
data=data_dict_design, # 対象データ
seed=1, # 乱数の種
sig_figs=18, # 出力CSV等に適用する数値精度
)【実行結果】
⑦ 収束確認
収束の確認をします。
診断メソッド diagnose を利用します。
# 事後分布の診断
print(fit_glmm.diagnose())【実行結果】（１行目のファイルパスは記載省略）
問題は検出されませんでした（no problems detected.）。
CmdStanPy の 出力結果 fit から MCMCサンプルの要約統計量を把握しましょう。
# p.252 結果の表示
(
fit_glmm.summary(percentiles=[2.5, 50, 97.5])
.loc[['b[1]', 'b[2]', 'b[3]', 'sigma_r'], :]
.round(2)
)【実行結果】
$${\widehat{R}}$$（R_hat）、有効サンプル数（N_Eff）に問題はなさそうです。
b[1] は切片、b[2] は晴れ、b[3] は気温の係数に対応しています。
fit を arviz の idata に変換します。
事後予測サンプルデータを idata に適切に変換するために、呪文のような引数を３つ設定しています（コードのコメント参照）。
# arvizのidataに変換
idata_glmm = az.from_cmdstanpy(
posterior=fit_glmm,
posterior_predictive='Y_pred', # 事後予測 posterior_predictive グループを指定
observed_data={'Y': Y}, # 観測値 observed_data グループを指定
dtypes={'Y_pred': int}, # Y_predがfloat型で変換されないように int 指定
log_likelihood='lp__',
)
idata_glmm【実行結果】
事後予測 Y_pred を格納した posterior_predictive グループを表示します。
収束の追加確認をします。
テキストは 図 4.1.3 の可視化で大量パラメータの $${\widehat{R}}$$ を確認しています。
ここでは arviz の rhat 関数を用いて、$${\widehat{R} > 1.01}$$ のパラメータがないことを数値で確かめます。
# p.251 収束の確認 r_hat>1.01の確認
# 設定
idata_in = idata_glmm # idata名
threshold = 1.01 # しきい値
# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())【実行結果】
4B の右側の数値が $${\widehat{R} > 1.01}$$ のパラメータの個数です。
すべて $${0}$$ ですので、$${\widehat{R} \leq 1.01}$$ と言えます。
トレースプロットを描画します。
# トレースプロットの描画
var_names = ['b', 'sigma_r']
az.plot_trace(idata_glmm, var_names=var_names, compact=False,
backend_kwargs={'tight_layout': True});【実行結果】
左側のチャートの４本の Chain はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。
以上のチェックに基づいて、収束していると考えましょう。
⑧ 推定されたモデルの解釈
線形予測子の数式表現を再掲します。
$$
\log(\lambda_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_i \\
$$
両辺の指数 $${\exp}$$ をとると次のようになります。
$$
\begin{align*}
\lambda_i &= \exp(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_i) \\
&= \exp(\beta_0) \times \exp(\beta_1 x_{i1}) \times \exp(\beta_2 x_{i2}) \times \exp(r_i)
\end{align*}
$$
係数（および変数の値）の変動によって、平均パラメータ $${\lambda_i}$$ は指数 $${\exp(\cdot)}$$ 倍変動します。
これを踏まえて、GLMM モデルの事後分布要約統計量を確認しましょう。
ここはひとつ、Gemini に解説をお願いしましょう！
1. 気温の影響は「確実」にある
気温 $${\texttt{b[3]}}$$ の 95%信用区間が 0.04 〜 0.11 となっており、0 を含まずにプラス側に振り切っています。
気温が上がると釣果が伸びるという傾向は、この過分散なデータの中でも統計的にしっかりと検出できています。
2. ランダム効果の標準偏差 $${\texttt{sigma\_r}}$$ の存在感
ここが今回の肝です。ランダム切片の標準偏差の平均が 1.10 です。
対数リンク関数を使っているので、この数値は「個体ごとのゆらぎによって、釣果が平均的に（1 標準偏差的に）$${\exp(1.10) \approx 3}$$ 倍程度、あるいは $${\exp(-1.10) \approx 1/3}$$ 程度まで変動する」という、ノイズの激しさを物語っています。
3. 天気の効果 $${\texttt{b[2]}}$$
天気（晴れ） $${\texttt{b[2]}}$$ が -0.73 です。指数をとると $${\exp(-0.73) \approx 0.48}$$。
つまり、「晴れ」は「曇り」に比べて、気温などの条件が同じでも釣果が半分くらいに落ち込む傾向がある、と読み解けます。
ランダム切片で個々のデータ特有のばらつきを表現できているようですね。
それにしてもランダム切片を解釈して、個々のデータが $${1/3}$$ ～ $${3}$$ 倍の変動するというのは驚きです！
🔵🔵🔵
GLMとGLMMを比べる
ポアソン回帰モデル（GLM）とポアソンGLMMモデルを比べましょう。
ランダム効果を取り入れた GLMM が過分散をうまくモデルに取り入れていることを追確認します。
① 事後予測チェックプロット（plot_ppc）
２モデルの事後予測チェックプロットで比べます。
観測データとモデルの予測データの分布が一致しているかを確認します。
# 事後予測チェック
# 描画領域の設定
fig, ax = plt.subplots(1, 2, figsize=(12, 5), sharey=True)
# ポアソンGLM (過分散を捉えきれず、実測 y とズレやすい)
az.plot_ppc(
idata_glm, kind='kde', data_pairs={'Y': 'Y_pred'}, num_pp_samples=100,
ax=ax[0]
)
ax[0].set_title('Poisson GLM')
# ポアソンGLMM (OLREにより過分散を吸収し、実測 y にフィットしやすい)
az.plot_ppc(
idata_glmm, kind='kde', data_pairs={'Y': 'Y_pred'}, num_pp_samples=100,
ax=ax[1]
)
ax[1].set_title('Poisson GLMM');【実行結果】
観測データの黒実線と事後予測データのオレンジ破線に注目します。
左側の GLM は両者のズレが明確です。
一方で、右側の GLMM は観測データと予測データが一致している感じがします。GLMM が過分散にうまく対応できていると思います。
ただ、予測データは 20 以上の右裾の予測も含んでますね…
② 99% 予測 HDI 区間のプロット
GLM と GLMM の 99% 予測HDI 区間が観測データを含んでいるかどうかを確認します。
ここでは天気別ではなく、全体の事後予測を描画します。
# GLMとGLMMの事後予測分布のチェック
# 事後予測分布の描画ヘルパー関数
def plot_post_pred(
idata, x_data, y_data, y_obs_col, color, label, title, y_label, ax,
hdi_prob=0.95
):
# 1. 事後予測分布の取得
pred_y = idata.posterior_predictive[y_obs_col]
pred_y_stack = pred_y.stack(sample=('chain', 'draw'))
# 2. 中央値と95%高密度区間(HDI)の算出
y_med = pred_y_stack.median('sample')
y_hdi = az.hdi(pred_y, hdi_prob=hdi_prob)[y_obs_col]
# 3. x軸データでソートするためのインデックスの取得
idx = np.lexsort((np.arange(len(x_data)), x_data))
# 4. 予測分布の可視化
# 観測値の散布図
ax.scatter(
x_data, y_data, s=40, color=color, edgecolor='white', label='観測値'
)
# 事後予測のHDI区間 (ソート順に描画)
ax.fill_between(
x_data[idx], y_hdi[idx, 0], y_hdi[idx, 1], color=color, alpha=0.2,
label=f'{label} {hdi_prob:.0%} HDI区間'
)
# 事後予測の中央値 (ソート順に描画)
ax.plot(x_data[idx], y_med[idx], color=color, label=f'{label} 中央値')
# 修飾
ax.set_title(title, fontsize=12)
ax.set_ylabel(y_label, fontsize=12)
ax.legend(loc='upper left')
## 設定と準備
# x軸のデータ
x_data = fish_num_climate_2['temperature'].to_numpy()
# y軸のデータ
y_data = fish_num_climate_2['fish_num'].to_numpy()
## 描画
# 描画領域の設定
fig, ax = plt.subplots(2, 1, figsize=(8, 6), sharey=True, tight_layout=True)
# GLMの事後予測の描画
plot_post_pred(idata_glm, x_data, y_data, 'Y_pred', 'tab:green', 'GLM',
'GLMの事後予測', '釣獲尾数', ax[0], 0.99)
# GLM<の事後予測の描画
plot_post_pred(idata_glmm, x_data, y_data, 'Y_pred', 'tab:green', 'GLMM',
'GLMMの事後予測', '釣獲尾数', ax[1], 0.99)
# 修飾
ax[1].set_xlabel('気温', fontsize=12)
plt.show()【実行結果】
GLM は先ほど見たのと同様に、観測データの幾つかは 99% 予測 HDI 区間の外側にあります。モデルがデータの過分散に対応できていません。
一方で GLMM は観測データが 99% 予測 HDI 区間の範囲内に含まれています。過分散対応ができていると思います。
③ GLMM の天気別 99% 予測 HDI 区間のプロット
GLMM を天気別で見てみましょう。
# GLMMの天気別の事後予測分布のチェック
# GLMMの天気別の事後予測分布の描画ヘルパー関数
def plot_post_pred_glmm(
idata, df, x_col, y_col, y_obs, cat_col, cat_name, color, title, y_label,
ax, hdi_prob=0.99
):
# 0. 設定と準備
x_data = df[x_col].to_numpy() # x軸のデータ
y_data = df[y_col].to_numpy() # y軸のデータ
# 1. 事後予測分布の取得
y_pred = idata.posterior_predictive[y_obs]
# 2. 該当する要素のデータのインデックスを取得
cat_idx = (df[cat_col] == cat_name).to_numpy()
# 2. 該当する要素の観測値 x, y を抽出
cat_x = df.loc[cat_idx, x_col].values
cat_y = df.loc[cat_idx, y_col].values
# 3. 事後予測分布から該当データのみを切り出す
y_pred_sub = y_pred[:, :, cat_idx]
# 4. xでソートするためのインデックスを取得（同値の場合は元のインデックス順を維持）
# lexsortは後ろに指定したキーが優先されるため、(index, temp) の順で指定
sort_idx = np.lexsort((np.arange(len(cat_x)), cat_x))
cat_x_sorted = cat_x[sort_idx]
# 5. HDIと中央値を計算（ソート済みデータに対して）
y_pred_stacked = y_pred_sub.stack(sample=('chain', 'draw'))
hdi = az.hdi(y_pred_sub, hdi_prob=hdi_prob)[y_obs]
median = y_pred_stacked.median(dim='sample')
# 6. 可視化
# 観測値の散布図
ax.scatter(cat_x, cat_y, s=20, color=color, alpha=0.8,
label=f'観測値 ({cat_name})')
# 事後予測のHDI区間 (ソート順に描画)
ax.fill_between(cat_x_sorted, hdi[sort_idx, 0], hdi[sort_idx, 1],
color=color, alpha=0.2, label=f'{hdi_prob:.1%} 事後予測 HDI')
# 事後予測の中央値 (ソート順に描画)
ax.plot(cat_x_sorted, median[sort_idx], color=color, label='事後予測中央値')
# 修飾
ax.set_title(f'{title}（{cat_col} = {cat_name}）', fontsize=12)
ax.set_ylabel(y_label, fontsize=12)
ax.legend(loc='upper left')
## 描画
# 描画領域の設定
fig, ax = plt.subplots(2, 1, figsize=(8, 6), sharey=True, tight_layout=True)
# GLMの事後予測の描画
plot_post_pred_glmm(
idata_glmm, fish_num_climate_2, 'temperature', 'fish_num', 'Y_pred',
'weather', 'cloudy', 'tomato', 'GLMMの事後予測チェック', '釣獲尾数', ax[0]
)
# GLM<の事後予測の描画
plot_post_pred_glmm(
idata_glmm, fish_num_climate_2, 'temperature', 'fish_num', 'Y_pred',
'weather', 'sunny', 'tab:blue', 'GLMMの事後予測チェック', '釣獲尾数', ax[1]
)
# 修飾
ax[1].set_xlabel('気温', fontsize=12)
plt.show()【実行結果】
GLMM は天気別においても、観測データが 99% 予測 HDI 区間の範囲内に含まれています。過分散対応ができていると思います。
GLMM モデルはデータの過分散にうまく対処できています！
🔵🔵🔵
アディショナルタイム：GPBoost 実装例
話は GLMM の最尤推定（つまり非ベイズ）に移ります。
GPBoost ライブラリのパラメータ推定を Gemini に訊いてみました。
さっそく GPBoost でポアソン GLMM モデルを構築しましょう。
① 追加インポート
# ライブラリの追加インポート
import gpboost as gpb
gpb.__version__【実行結果】
今回利用する GPBoost のバージョンは 1.6.1 です（最新版では無いです）。
② ポアソン GLMM のモデリング
さっそくクライマックスです！
ランダム切片の変数は 引数 group_data に設定します。
ランダム切片と説明変数に pandas の DataFrame、目的変数に pandas の Series を渡すと、分析結果に変数名が表示されるのでおすすめです。
このモデリングでは先ほど作成したデザイン行列を再利用しています。
# ポアソンGLMMモデルによる最尤推定：ラプラス近似を用いた（制限付き）最尤法
# 1. データセットの作成
# ランダム効果（グループ）を作成 ※pandas DataFrame
group_data = fish_num_climate_2[['id']]
# 2. GLMMの実行
# モデルの定義：ランダム効果はグループ（ID）を指定
gp_model = gpb.GPModel(group_data=group_data, likelihood='poisson')
# 学習 ※説明変数には先ほど作成したデザイン行列を利用
gp_model.fit(y=Y_dm.iloc[:, 0], X=X_dm, params={'std_dev': True})
# 結果サマリーの表示 ※ランダム効果のCovarianceは分散（標準偏差でない）
print(gp_model.summary())【実行結果】
このモデルサマリーの見方を Gemini に教えてもらいました。
【GPBoostによるパラメータ推定結果の解釈】
1. 固定効果（Fixed effects）
平均的な傾向を決定する「集団共通のルール」です。
weather[T.sunny] (-0.7193)：晴れの日（sunny）は、基準となる天気よりも対数期待値が約0.72減少します（P値 < 0.01で有意）。
temperature (0.0737)：気温が1度上がるごとに、対数期待値が約0.07増加します（P値 < 0.001で非常に有意）。
2. ランダム効果（Random effects）
データの個別のゆらぎを担当する部分です。
id (1.0263)：データ点単位のランダム切片 $${r_i}$$ の分散です。
ランダム切片 $${r_i}$$ は、ポアソン分布が持つ理論上のばらつき（誤差構造）を超えた個別のばらつきの表現を担当します。今回のモデルでは、この 1.0263 という数値が「過分散」をランダム切片が捉えた結果となります。
3. 学習手法の背景
この結果は、周辺尤度（ランダム効果を積分消去した尤度）を最大化する「最尤法」によって導かれました。 ポアソン分布（非ガウス尤度）の周辺化には、計算負荷を抑えつつ精度を確保する「ラプラス近似」が用いられています。
4. ブログ用まとめ
ランダム切片のばらつきは「分散」で表示されています。
ベイズモデリングでは「標準偏差」を求めています。
分散の平方根 $${\sqrt{1.0263}=1.013}$$ が標準偏差です。
③ モデルの学習で得られた情報を表示
学習済みモデル gp_model から情報を取り出してみましょう。
◆ パラメータ推定値と標準誤差
# パラメータ推定値の表示 推定値と標準誤差
gp_model.get_coef().T【実行結果】
◆ ランダム切片の分散（の平方根を取って標準偏差を算出）
# ランダム切片の分散推定値を標準偏差に変換して表示
gp_model.get_cov_pars()**(1/2)【実行結果】
④ 予測分布の可視化
パラメータ推定値とランダム切片 $${r}$$ の推定値を使って予測分布を算出し、可視化します。
◆ 描画用データの算出
ランダム切片の推定値は学習済みモデル gp_model に対して predict_training_data_random_effects メソッドを適用して取得します。
# 描画用データの算出
# 設定
x_line = np.linspace(0, 30, 100)
# パラメータ推定値の取得
beta0_hat, beta1_hat, beta2_hat = gp_model.get_coef(False)[0]
# yの平均予測値の算出
y_pred_cloudy = np.exp(beta0_hat + beta1_hat * 0 + beta2_hat * x_line)
y_pred_sunny = np.exp(beta0_hat + beta1_hat * 1 + beta2_hat * x_line)
# ランダム切片rの推定値の取得
r_mean = gp_model.predict_training_data_random_effects(predict_var=True)
r_mean_df = pd.concat([fish_num_climate_2['weather'], r_mean], axis=1)
r_mean_df.head()【実行結果】
こちらはランダム切片 $${r}$$ の推定値（先頭 5 行）です。
では描画しましょう！
# 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(8, 4), tight_layout=True)
# 観測値の散布図の描画
sns.scatterplot(data=fish_num_climate_2, x='temperature', y='fish_num',
hue='weather', palette=colors, ax=ax)
# GLMMの予測値の太線の描画(ランダム切片を除く平均予測値)
ax.plot(x_line, y_pred_cloudy, color='tomato', lw=2)
ax.plot(x_line, y_pred_sunny, color='tab:blue', lw=2)
# GLMMの予測値の細線の描画（ランダム切片を含む個々のデータの予測値）
for i, row in r_mean_df.iterrows():
weather, r, _ = row
weather_num = 1 if weather=='sunny' else 0
y_line = np.exp(beta0_hat + beta1_hat * weather_num + beta2_hat * x_line + r)
ax.plot(x_line, y_line, color=colors[weather], lw=0.5, alpha=0.3)
# 修飾
ax.set_xlabel('気温', fontsize=12)
ax.set_ylabel('釣獲尾数', fontsize=12)
ax.set_ylim(top=16);【実行結果】
赤太線は曇りの予測分布平均値、青太線は晴れの予測分布平均値です。
これらにはランダム切片を含めていません。
薄色の線はランダム切片を含めた個々のデータの予測分布です。
平均はそこそこ立派な線を描いているものの、１つ１つのデータに目を向けると、平均では表現しきれない、多様なばらつきを持つ、個性的で元気な分布を目の当たりにします。
ランダム切片（ランダム効果）をモデルに追加するだけで、モデルの表現力が格段に上がるのだなと実感いたします。
🔵🔵🔵
今回の記事は以上です。
楽しかったですね！
シリーズの記事
次の記事
前の記事
PyMC版
目次
ブログの紹介
note で８つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．統計・データ分析とつながる
シリーズ「統計・データ分析とつながる」は、統計・データ分析との「つながり」を発掘して、コラム風に仕立てたブログシリーズです。
生成 AI の力を借りながら、統計・データ分析の入り口をイメージして、自由気ままに書きました。
たとえば…
・日常生活と統計のつながり
・統計検定２級からその先へのつながり
気楽にお読みいただけたら嬉しいです🍀
３．実験！たのしいベイズモデリング１＆２を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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！
