見出し画像

「ベイズ統計モデリングによるデータ分析入門」をPythonとPyMCで写経 ～ Vol.18 ランダム切片モデル

ネイピア DS

書籍の著者　馬場真哉 先生

この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第４部第２章「ランダム切片モデル」の Python 写経活動記録です。

第４部は一般化線形混合モデル（通称：GLMM）の階層ベイズモデルです。
今回は「ランダム切片モデル」を３つのツール・手法で取り組みます。

  • Bambiのベイズランダム切片モデル

  • PyMCのベイズランダム切片モデル

  • アディショナルタイム：GPBoostのランダム切片モデル

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


はじめに

このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。

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

準備

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

2.2 分析の準備
2.5 brms によるランダム切片モデルの推定
2.6 回帰曲線の図示

ベイズ統計モデリングの理論面が気になった場合には、ぜひテキストの第１部および本章の理論面の記載をご覧いただき、本記事との繋がりをご確認下さいませ。

■ コード記述法
Jupyter Notebook 形式でコードを記述します。

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。

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

# インポート

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

# ベイズ統計モデリング
import pymc as pm                      # pymc
import bambi as bmb                    # bambi
import arviz as az                     # 分析・可視化

# デザイン行列
from patsy import dmatrices, dmatrix

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme()                        # ggplot風のスタイル
plt.rcParams['font.family'] = 'Meiryo'

第２章 ランダム切片モデル

前回記事から引き続き、ポアソンGLMMの階層ベイズモデルである「ランダム切片モデル」を実装します。
GLMMおよび階層ベイズモデルの概要はぜひ前回記事をお読み下さい！

🚀🚀🚀

データの読み込みと外観の確認

テキスト 2.2 節に相当します。
テキストの仮想の魚の釣獲尾数データを引用いたします。
csv ファイルを pandas のデータフレーム形式で変数 fish_num_climate_3 に読み込みます。

# p.255 分析対象データの読み込み

# ファイルの読み込み
fish_num_climate_3 = pd.read_csv('./data/4-2-1-fish-num-3.csv')

# 結果の表示
print('fish_num_climate_3.shape: ', fish_num_climate_3.shape)
fish_num_climate_3.head(3)

【実行結果】
標本サイズ 100のデータです。

  • fish_num：釣った魚の数（釣獲尾数）

  • weather：天気の種類（晴れ：sunny、曇り：cloudy）

  • temperature：気温

  • human：釣り人の名前（A ～J までの 10 人）

のちのち使う目的で共通データを作成します。
１つ目は釣り人の名前リストです。

# 釣り人の名前リスト
humans = sorted(fish_num_climate_3['human'].unique())
humans

【実行結果】
A ～J までの 10 人の釣り人のデータです。

２つ目は天気ごとの色設定です。

# 天気ごとの色の設定
colors = {'cloudy': 'tomato', 'sunny': 'tab:blue'}
colors

【実行結果】
曇りがトマト色（オレンジ寄りの赤）、晴れが青色です。

🚀

天気ごとおよび釣り人ごとのデータ件数をカウントします。
pandas の value_counts メソッドを利用します。

# カテゴリ変数の要素ごとのデータ件数
for col in fish_num_climate_3.columns[fish_num_climate_3.dtypes == 'object']:
    display(fish_num_climate_3[col].value_counts().to_frame())

【実行結果】
曇り・晴れで 50 件ずつ、釣り人ごとに 10 件ずつのデータです。

データの要約統計量を確認します。
まずは全体（量的変数のみ）です。

# データの要約統計量
fish_num_climate_3.describe().T.round(2)

【実行結果】
釣獲尾数は平均 2.5、最小値 0、最大値 15 です。
気温は平均 15.6、最小値 0.4、最大値 29.7 です。

釣獲尾数にポアソン分布を仮定する際には、平均 2.5 と分散 9 （標準偏差の二乗）の乖離、つまり過分散の考慮が必要かもしれません。

次は釣り人別・天気別の釣獲尾数の要約統計量です。

# 釣り人 × 天気 ごとの要約統計量
fish_num_climate_3.groupby(['human', 'weather'])['fish_num'].describe().round(2)

【実行結果】
さまざまにばらついている感じです。
数字を追うのは大変💦ですので、後で可視化しましょう。

データを可視化しましょう。
最初に散布図で気温と釣獲尾数の関係を可視化します。
seaborn の scatterplot を利用します。

# 天気別散布図の描画
sns.scatterplot(
    data=fish_num_climate_3, x='temperature', y='fish_num',
    hue='weather', palette=colors
)
# 修飾
plt.xlabel('気温 [℃]', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12)
plt.legend(title='天気');

【実行結果】
散布図では気温が上昇するにつれて、釣獲尾数の増加傾向とばらつき増加傾向の関係が見られます。

続いて、箱ひげ図で釣り人別・天気別の釣獲尾数のばらつきを確認します。
seaborn の boxplot を利用します。

# 箱ひげ図の描画

# 描画領域の設定
plt.figure(figsize=(10, 4))
# 箱ひげ図の描画
sns.boxplot(data=fish_num_climate_3, x='human', y='fish_num',
            hue='weather', palette=colors, fill=False, gap=0.2)
# スウォームプロットの描画
sns.swarmplot(data=fish_num_climate_3, x='human', y='fish_num',
              hue='weather', palette=colors, dodge=True)
# 修飾
plt.xlabel('釣り人', fontsize=12)
plt.ylabel('釣獲尾数', fontsize=12)
plt.legend(title='天気');

【実行結果】
釣り人によって釣れた魚の数のばらつきが大きく異なっている印象です。
また、曇りの方が釣果がいい人、晴れの方が釣果がいい人に分かれている感じもします。

最後に釣り人・気温と釣獲尾数の散布図を描画します。
seaborn の relplot で釣り人別の小さなチャートを配置します。

# 釣り人別散布図の描画
g = sns.relplot(
    data=fish_num_climate_3, x='temperature', y='fish_num',  # データ, x軸, y軸
    col='human', col_wrap=3,                                 # 列の変数, 列の数
    hue='weather', palette=colors,                           # hueの変数、色
    kind='scatter', s=80, height=3, aspect=1.2
)
# 修飾
g.set_titles(col_template='釣り人 {col_name}')
g.set_axis_labels('気温 [℃]', '釣獲尾数')
g._legend.set_title('天気');

【実行結果】

A さん、C さんがよく釣れている感じ、D ～ F さんはあまり釣れていない感じ、天気と釣果の関係は…今ひとつ分かりません…。
ベイズモデルに期待しましょう！

なお、この relplot を用いた散布図の描画は、以後の分析でも頻出するので、お楽しみに🍀

🚀🚀🚀

ベイズモデリング by Bambi

テキスト 2.5、2.6 節に相当します。
brms の代わりに Bambi を利用します。

① モデルの数式
今回は次の数式で示されるランダム切片モデルに取り組みます。
$${r_k}$$ が釣り人 $${k}$$ ごとのランダム効果（ランダム切片）です。

$$
\begin{align*}
r_k &\sim \text{Normal}(0, \sigma_r^2) \\
\log(\lambda_i) &= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_k \\
y_i &\sim \text{Poisson}(\lambda_i)
\end{align*}
$$

テキスト p.255 式（4.3）を一部改変して引用

💡 ランダム切片の気持ちを探る
ランダム切片 $${r_k}$$ は平均 $${0}$$、標準偏差 $${\sigma_r}$$（分散 $${\sigma_r^2}$$）の正規分布に従うとしています。
つまり、釣り人のランダム切片の効果は１つの正規分布からサンプリングされて決まるという関係になっています。
そして、標準偏差 $${\sigma_r}$$ が小さい場合は各釣り人のランダム効果は平均 $${0}$$ に近く、みなの釣果は似たりよったりになります。
標準偏差 $${\sigma_r}$$ が大きい場合は各釣り人のランダム効果はばらつきが大きく、個性的な釣果が生まれる可能性があります。

一見、ランダム切片は「ダミー変数の係数」に似ているように見えます。
しかし、ダミー変数はカテゴリーごとに独立して推定されるため、共通の確率分布の設計（または縛り）を持っていません。
Gemini はこう言ってます。

ダミー変数のようには、少ないデータで釣り人の実力を決めつけない。この『正規分布という設計図』の範囲内に収まるはずだと考えることで、極端な予測ミスを防ぐ。これこそが、ランダム効果モデルが持つ『賢さ』の正体なのです。

② モデルの概要
確率分布をポアソン分布、リンク関数を対数リンク関数（ポアソン分布のデフォルト）、線形予測子を次の formula で実装します。

$$
\texttt{fish\_num} \sim \texttt{weather} + \texttt{temperature} + \texttt{(1|human)}
$$

$${\texttt{(1|human)}}$$ がランダム切片の表現です。
テキストの $${r_k}$$ のことです。
ランダム効果は $${|}$$ で変数を挟む次の書き方で表現します。

$$
(\ \texttt{ランダム効果が与えられる対象}\ |\ \texttt{ランダム効果変数}\ )
$$

今回、ランダム効果が与えられる対象は「切片」（つまり $${1}$$）、ランダム効果の変数は釣り人を識別する「$${\texttt{human}}$$」です。

テキストの無情報事前分布に合わせるべく、Bambi に以下の事前分布情報を与えます。

$$
\begin{align*}
intercept &\sim \text{Normal}\ (0, (1e5)^2) \\
weather &\sim \text{Normal}\ (0, (1e5)^2) \\
temperature &\sim \text{Normal}\ (0, (1e5)^2) \\
(1 | human) &\sim \text{Normal}\ (0, \sigma^2) \\
\sigma &\sim \text{HalfNormal}(10^2) \ ^\dagger
\end{align*}
$$

$${^\dagger}$$：標準偏差パラメータには弱情報に相当する $${10}$$ を与えています。

③ モデル定義
上式のモデルを Bambi で記述します。
bmb.prior() で事前分布を定義します。複数ある場合は辞書でまとめます。
bmb.Model() において、prior 引数で定義した事前分布を与えます。

# モデリング

# 無情報事前分布（想定）の設定
uninformed_prior = bmb.Prior('Normal', mu=0, sigma=1e5)
# ランダム効果用の"Hyper prior"の設定
Hyper_prior = bmb.Prior('Normal', mu=0, sigma=bmb.Prior('HalfNormal', sigma=10))
# 事前分布を辞書にとりまとめ
priors = {'Intercept': uninformed_prior, 'weather': uninformed_prior,
          'temperature': uninformed_prior, '1|human': Hyper_prior}

# モデルの定義
model_bmb = bmb.Model(
    formula='fish_num ~ weather + temperature + (1|human)',  # フォーミュラ式
    data=fish_num_climate_3,                              # データ
    family='poisson',                                     # 確率分布：ポアソン分布
    priors=priors,                                        # 事前分布（辞書）
)

【実行結果】なし

モデルの内容を表示します。

# モデルの表示
model_bmb

【実行結果】
確率分布 Family に ポアソン分布 poisson、リンク関数 Link に対数リンク関数 mu = log が設定されました。

ランダム切片は「Group-level effects」（グループ効果）に設定されました。
その他の線形予測子のパラメータは「Common-level effects」（固定効果）に在ります。

モデルをグラフィカルモデルで描画します。

# モデルの可視化 ※1|human_offsetはNormal(0, 1)。非中心化パラメータ化しているらしい
# 1|human = 1|human_offset × 1|human_sigma
model_bmb.build()
model_bmb.graph()

【実行結果】
左上の３つの変数がランダム切片に関わるものです。
「$${\texttt{1|human}}$$」がランダム切片です。

④ MCMC の実行
MCMCを実行しましょう。
NUTS サンプラーに nutpie を利用します。
また、ダイバージェンス対策として、tune（バーンイン）を 10000 に増量しています。

%%time
# p.257 MCMCの実行 ※divargence回避のためtuneを増やす

idata_bmb = model_bmb.fit(
    draws=1000, tune=10000, chains=4, random_seed=123, nuts_sampler='nutpie',
)

【実行結果】
Divergences（ダイバージェンス）は０件です。

⑤ 収束確認
arviz の rhat 関数を用いて、$${\widehat{R} > 1.01}$$ のパラメータがないことを数値で確かめます。

# r_hat>1.01の確認

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

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

【実行結果】
4B の右側の数値が $${\widehat{R} > 1.01}$$ のパラメータの個数です。
すべて $${0}$$ ですので、$${\widehat{R} \leq 1.01}$$ と言えます。

MCMC サンプルの要約表を表示します。

# p.257 要約統計量の表示
var_names = ['Intercept', 'weather', 'temperature', '1|human_sigma']
az.summary(idata_bmb, var_names=var_names, hdi_prob=0.95)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。

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

# トレースプロットの描画
az.plot_trace(idata_bmb, compact=False, var_names=var_names,
              backend_kwargs={'tight_layout': True});

【実行結果】
左側のチャートの４本の Chain はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。

以上のチェックに基づいて、収束していると考えましょう。

⑥ 推定されたモデルの解釈
線形予測子の数式表現を再掲します。

$$
\log(\lambda_i) = \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_k \\
$$

両辺の指数 $${\exp}$$ をとると次のようになります。

$$
\begin{align*}
\lambda_i &= \exp(\beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + r_k) \\
&= \exp(\beta_0) \times \exp(\beta_1 x_{i1}) \times \exp(\beta_2 x_{i2}) \times \exp(r_k)
\end{align*}
$$

係数（および変数の値）の変動によって、平均パラメータ $${\lambda_i}$$ は指数 $${\exp(\cdot)}$$ 倍変動します。
これを踏まえて、GLMM モデルの事後分布要約統計量を確認しましょう。

前回記事の Gemini 謹製分析テンプレを利用して読み解きを進めます。

1. 気温の影響は「確実」にある
$${\texttt{temperature}}$$ の 95%HDIが 0.080 〜 0.115 となっており、0 を含まずにプラス側に振り切っています。
気温が上がると釣果が伸びるという傾向は、この過分散なデータの中でも統計的にしっかりと検出できています。

2. ランダム効果の標準偏差 $${\texttt{1|human\_sigma}}$$ の存在感
ここが今回の肝です。ランダム切片の標準偏差の平均が 0.643 です。
対数リンク関数を使っているので、この数値は「個体ごとのゆらぎによって、釣果が平均的に（1 標準偏差的に）$${\exp(0.643) \approx 1.9}$$ 倍程度、あるいは $${\exp(-0.643) \approx 0.53}$$ 程度まで変動する」という、ノイズの大きさを物語っています。

3. 天気の効果 $${\texttt{weather[sunny]}}$$
$${\texttt{weather[sunny]}}$$ が -0.521 です。指数をとると $${\exp(-0.521) \approx 0.59}$$。
つまり、「晴れ」は「曇り」に比べて、気温などの条件が同じでも釣果が半分近くまで落ち込む傾向がある、と読み解けます。

⑦ 釣り人の影響力（ランダム切片）の確認
テキスト p.257 の ranef 関数と同等のことを az.summary 関数で実行します。95%HDI 区間を表示します。

# p.257 釣り人の影響の大きさ（ランダム切片の要約統計量）
az.summary(idata_bmb, var_names=['1|human'], kind='stats', hdi_prob=0.95)

【実行結果】
A ～ H の釣り人のランダム切片の推定値です。
指数変換前ですので、$${\log(\lambda_i)}$$ への加減算の効果です。
テキストのとおり、Aさんは釣りが得意、D・Fさんは釣りが苦手なのかもしれません。

ランダム効果を可視化しましょう。
フォレストプロットで表現します。
arviz の plot_forest 関数を利用します。
指数変換を行って、スケールを「倍率」にしてみましょう。

# フォレストプロット（expで「倍率」スケールに変換後）

# 値を指数変換 (exp) する
# これにより、単位が「対数」から「倍率」に変わります
ds_ratio = idata_bmb.posterior[['1|human']].copy()
ds_ratio['1|human'] = np.exp(ds_ratio['1|human'])

# フォレストプロットの描画
az.plot_forest(ds_ratio, var_names=['1|human'], combined=True, hdi_prob=0.95)
plt.axvline(1, color='tab:red');

【実行結果】
値が１倍より大きい釣り人 A、C さんは正の倍率効果がありそうです。
一方で１倍より小さい釣り人 D、F さんは負の倍率効果がありそうです。
その他の釣り人は 95%HDI 区間が１倍を含んでいるので、正・負の判断をしかねます。

⑧ 釣り人別・天気別の釣獲尾数の可視化
テキストの図 4.2.1 に相当します。
釣り人別にチャートを分割して、天気別平均獲釣尾数の平均値と 95%HDI 区間を描画します。
まず、平均獲釣尾数の平均値と 95%HDI 区間を算出します。
Bambi の interpret.predictions 関数を利用します。

# 平均釣獲尾数の平均値と95%HDI区間の取得 ※Bambiの予測計算関数を利用

mean_pred = bmb.interpret.predictions(
    model_bmb,                                        # Bambiモデル
    idata_bmb,                                        # idata
    conditional=['temperature', 'weather', 'human'],  # 条件付けする共変量
    target='mean',                                    # muの事後分布
    use_hdi=True,                                     # True: HDI, False:分位数
    prob=0.95,                                        # 区間の確率
)
mean_pred

【実行結果】
モデルの学習データの値を使って、平均値 estimate、2.5% HDI、97.5 HDI を計算しました（HDI 列名の％は気にしないで下さい）。

描画します。
seaborn の relplot の axes を取得して、その上に平均値・95%HDI 区間を重ね描きするのがこの記事のクライマックスです！（当社比）

# p.259 図4.2.1 釣り人別の回帰曲線

# 釣り人別散布図の描画
g = sns.relplot(
    data=fish_num_climate_3, x='temperature', y='fish_num',  # データ, x軸, y軸
    col='human', col_wrap=3,                                 # 列の変数, 列の数
    hue='weather', palette=colors,                           # hueの変数、色
    kind='scatter', s=70, height=3, aspect=1.2
)

# 釣り人別の平均釣獲尾数の平均値と95% HDI区間の描画
# sns.relplotのaxesを取り出して、釣り人ごとに描画を繰り返し処理
for human, ax in g.axes_dict.items():
    # 天気ごとに描画を繰り返し処理
    for weather in colors.keys():
        # 対象の釣り人と天気をデータフレームから取得
        query = (mean_pred['human']==human) & (mean_pred['weather']==weather)
        subset = mean_pred[query]
        # 平均釣獲尾数の平均値の描画
        ax.plot(subset['temperature'], subset['estimate'], color=colors[weather])
        # 平均釣獲尾数の95% HDI区間の塗りつぶし描画
        ax.fill_between(subset['temperature'],
                        subset['lower_3.0%'], subset['upper_97.0%'],
                        color=colors[weather], alpha=0.2)
        # 修飾
        ax.set(title=f'釣り人 {human}', xlabel='', ylabel='')

# チャート全体のx,y軸ラベルの表示
g.figure.supxlabel('気温 [℃]', fontsize=14)
g.figure.supylabel('釣獲尾数', fontsize=14)
g._legend.set_title('天気');

【実行結果】

気温が高くなるにつれて釣獲尾数が増える傾向と、曇りの方が晴れよりも釣獲尾数が多い傾向が可視化されています。
先ほど確認した「気温と晴れダミーのパラメータ推定値」を釣り人共通（ランダム効果を除く線形予測子）で適用して、予測値を算出しているので、当然と言えば当然です。
あとは、釣り人の効果としてランダム切片を加味することで、各釣り人の回帰曲線が大きくなったり小さくなったりしています。

🚀🚀🚀

ベイズモデリング by PyMC

テキスト 2.5、2.6 節を PyMC で実装します。
デザイン行列を作成して PyMC でモデリングします。

① モデルの概要
次のモデルを実装します。
ランダム切片は $${\bm r}$$（次元 10：釣り人の人数）です。

$$
\begin{align*}
\bm Y &\sim \text{Poisson}(\bm \lambda) \\
\log(\bm \lambda) &= \bm{X \beta + \bm r}[human] \\
\bm \beta &\sim \text{Normal}(0, (1e5)^2, \text{dims}=3) \\
\bm r &\sim \text{Normal}(0, \sigma_r^2, \text{dims}=10) \\
\sigma_r &\sim \text{HalfNormal}((1e5)^2) \ ^{\dagger}\\
\end{align*}
$$

$${^\dagger}$$：標準偏差パラメータには無情報に相当する $${1e5}$$ を与えています。Bambi モデルと変えています。

② データセットの作成
patsy ライブラリを用いてデザイン行列等を作成します。

# デザイン行列の作成

# formula構文の設定
formula_pois_glmm = 'fish_num ~ weather + temperature'

# 目的変数Y, デザイン行列（説明変数）Xの作成
Y_dm, X_dm = formula_lm = dmatrices(formula_pois_glmm, fish_num_climate_3,
                                    return_type='dataframe')

# デザイン行列の先頭5行の表示
X_dm.head()

【実行結果】
デザイン行列は定数項、天気のダミー変数（晴れ）、気温で構成されます。

# 目的変数の先頭5行の表示
Y_dm.head()

【実行結果】
目的変数も作ってくれました。

③ モデル定義
上式のモデルを PyMC で記述します。
データフレームの human 列の値を数値化するために、pandas の factorize 関数を利用して、数値化後の変数 human_ids を得て、Data クラスで変数 $${\texttt{human}}$$ を作ります。
次元 10 のランダム切片 $${r}$$ を配列に見立てて、線形予測子に $${r[\texttt{human}]}$$ を足しています。

# モデリング

# 釣り人のコード変換 ※A,B,C... から 0,1,2...へ変換
human_vars, human_ids = pd.factorize(fish_num_climate_3['human'], sort=True)

# coordsの設定
coords = {'id': fish_num_climate_3.index.values,           # 観測データの識別子
          'coefs': ['intercept', 'sunny', 'temperature'],  # 係数ベクトルbの識別子
          'human_id': human_ids}                           # 釣り人のID

# モデルの定義
with pm.Model(coords=coords) as model_pm:
    
    ## dataの設定
    # 目的変数: 釣獲尾数データ
    Y = pm.Data('Y', value=Y_dm.values.flatten(), dims='id')
    # 説明変数: デザイン行列
    X = pm.Data('X', value=X_dm.values, dims=('id', 'coefs'))
    # 説明変数: 釣り人（ランダム効果の要素）
    human = pm.Data('human', value=human_vars, dims='id')

    ## 事前分布: 1e5は無情報事前分布的な分布
    # 係数ベクトル
    b = pm.Normal('b', mu=0, sigma=1e5, dims='coefs')
    # ランダム効果
    sigma_r = pm.HalfNormal('sigma_r', sigma=1e5)
    r = pm.Normal('r', mu=0, sigma=sigma_r, dims='human_id')

    ## 尤度関数: ポアソン分布を仮定
    lam = pm.Deterministic('lam', pm.math.exp(X @ b + r[human]), dims='id')
    obs = pm.Poisson('obs', mu=lam, observed=Y, dims='id')

【実行結果】なし

モデルを数式ライクに表示します。

# モデルの表示
model_pm

【実行結果】

モデルをグラフィカルモデルで描画します。

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

【実行結果】
右側の $${\texttt{sigma\_r}}$$ と $${\texttt{r}}$$ が階層を形成しています。

④ MCMC の実行
MCMCを実行しましょう。
NUTS サンプラーに nutpie を利用します。

%%time
# p.257 MCMCの実行
with model_pm:
    idata_pm = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1,
                         nuts_sampler='nutpie')

【実行結果】
Divergences（ダイバージェンス）は０件です。

⑤ 収束確認
arviz の rhat 関数を用いて、$${\widehat{R} > 1.01}$$ のパラメータがないことを数値で確かめます。

# r_hat>1.01の確認

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

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

【実行結果】
4B の右側の数値が $${\widehat{R} > 1.01}$$ のパラメータの個数です。
すべて $${0}$$ ですので、$${\widehat{R} \leq 1.01}$$ と言えます。

MCMC サンプルの要約表を表示します。

# p.257 要約統計量の表示
var_names = ['b', 'sigma_r']
az.summary(idata_pm, var_names=var_names, hdi_prob=0.95)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（ess_bulk、ess_tail）に問題はなさそうです。

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

# トレースプロットの描画
az.plot_trace(idata_pm, compact=False, var_names=var_names,
              backend_kwargs={'tight_layout': True});

【実行結果】
左側のチャートの４本の Chain はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。

以上のチェックに基づいて、収束していると考えましょう。

⑥ 釣り人の影響力（ランダム切片）の確認
テキスト p.257 の ranef 関数と同等のことを az.summary 関数で実行します。95%HDI 区間を表示します。

# p.257 釣り人の影響の大きさ（ランダム切片の要約統計量）
az.summary(idata_pm, var_names=['r'], kind='stats', hdi_prob=0.95)

【実行結果】
A ～ H の釣り人のランダム切片の推定値です。
指数変換前ですので、$${\log(\lambda_i)}$$ への加減算の効果です。
テキストのとおり、Aさんは釣りが得意、D・Fさんは釣りが苦手なのかもしれません。

なお、Bambi の推定値とちょっとズレている点が気になりますね…

ランダム効果を可視化しましょう。
フォレストプロットで表現します。
arviz の plot_forest 関数を利用します。
指数変換を行って、スケールを「倍率」にしてみましょう。

# フォレストプロット（expで「倍率」スケールに変換後）

# 事後分布のMCMCサンプルを指数変換 (exp) する
# これにより、単位が「対数」から「倍率」に変わります
ds_ratio = idata_pm.posterior[['r']].copy()
ds_ratio['r'] = np.exp(ds_ratio['r'])

# フォレストプロットの描画
az.plot_forest(ds_ratio, var_names=['r'], combined=True, hdi_prob=0.95)
plt.axvline(1, color='tab:red');

【実行結果】
値が１倍より大きい釣り人 A、C さんは正の倍率効果がありそうです。
一方で１倍より小さい釣り人 D、F さんは負の倍率効果がありそうです。
その他の釣り人は 95%HDI 区間が１倍を含んでいるので、正・負の判断をしかねます。

⑦ 釣り人別・天気別の釣獲尾数の可視化
テキストの図 4.2.1 に相当します。
釣り人別にチャートを分割して、天気別平均獲釣尾数の平均値と 95%HDI 区間を描画します。

作成するデータは Bambi と異なっています。
平均獲釣尾数のMCMCサンプルを生成します。
せっかくですので次の２つの計算方法を試してみます！
- その１：PyMC の外で計算する方法
- その２：PyMC の予測機能で計算する方法

【その１】PyMC の外で計算する方法
PyMC で生成したMCMCサンプルを使い、GLMMの数式に沿って $${\bm \lambda}$$ を計算します。

# 平均釣獲尾数のMCMCサンプルの生成 ※PyMCの外で算出

# MCMCサンプルから取り出し
# 係数β0, β1, β2 各shape=(4000,)
intercept, sunny, temper = az.extract(idata_pm.posterior).b.to_numpy()
# ランダム切片r shape=(10, 4000)
rs = az.extract(idata_pm.posterior).r.to_numpy()

# x軸の値（気温）の設定
x_val = np.linspace(fish_num_climate_3['temperature'].min(), 
                    fish_num_climate_3['temperature'].max(), 100)

# 平均釣獲尾数のMCMCサンプルの生成
# ※numpy多次元配列：shape=(10, 2, 100, 4000)⇒(釣り人, 晴れ, 気温, MCMCサンプル)
pred_samples = np.exp(
    intercept                                    # 切片 shape=(4000)
    + np.outer([0, 1], sunny)[None, :, None, :]  # 晴れ shape=(2,4000)
    + np.outer(x_val, temper)[None, None, :, :]  # 気温 shape=(100, 4000)
    + rs[:, None, None, :]                       # ランダム切片 shape=(10, 4000)
)
print('pred_samples.shape:', pred_samples.shape)

【実行結果】

４軸の numpy 配列にしています。

  • 第１軸：釣り人 10

  • 第２軸：晴れダミー 2（曇りと晴れ）

  • 第３軸：気温 100 個

  • 第４軸：MCMCサンプル 4000 個

描画します。
後利用できるように描画ヘルパー関数化しておきます。

# p.259 図4.2.1 釣り人別の回帰曲線

# 釣り人別散布図の描画
g = sns.relplot(
    data=fish_num_climate_3, x='temperature', y='fish_num',  # データ, x軸, y軸
    col='human', col_wrap=3,                                 # 列の変数, 列の数
    hue='weather', palette=colors,                           # hueの変数、色
    kind='scatter', s=70, height=3, aspect=1.2
)

# 釣り人別の平均釣獲尾数の平均値と95% HDI区間の描画
# sns.relplotのaxesを取り出して、釣り人ごとに描画を繰り返し処理
for human, ax in g.axes_dict.items():
    # 天気ごとに描画を繰り返し処理
    for weather in colors.keys():
        # 対象の釣り人と天気をデータフレームから取得
        query = (mean_pred['human']==human) & (mean_pred['weather']==weather)
        subset = mean_pred[query]
        # 平均釣獲尾数の平均値の描画
        ax.plot(subset['temperature'], subset['estimate'], color=colors[weather])
        # 平均釣獲尾数の95% HDI区間の塗りつぶし描画
        ax.fill_between(subset['temperature'],
                        subset['lower_3.0%'], subset['upper_97.0%'],
                        color=colors[weather], alpha=0.2)
        # 修飾
        ax.set(title=f'釣り人 {human}', xlabel='', ylabel='')

# チャート全体のx,y軸ラベルの表示
g.figure.supxlabel('気温 [℃]', fontsize=14)
g.figure.supylabel('釣獲尾数', fontsize=14)
g._legend.set_title('天気');

【実行結果】
Bambi の描画結果とほぼ同じのチャートです（当たり前 ^^;;;）。

【その２】PyMC の予測機能で計算する方法
予測用の説明変数データ（デザイン行列等）を作成し、PyMC の set_data メソッドでデータをセットして、sample_posterior_predictive 関数で平均釣獲尾数の予測値（MCMCサンプル）を生成します。

# 平均釣獲尾数のMCMCサンプルの生成 ※PyMCのsample_posterior_predictive利用

# 追加インポート
from itertools import product

# 予測条件の作成 (10人 × 2天気 × 100気温 = 2000行)
x_val = np.linspace(*np.sort(fish_num_climate_3['temperature'])[[0, -1]], 100)

# デザイン行列と釣り人インデックスを作成
# - 釣り人 × 天気 × 気温 の全組み合わせ（直積）の作成
prod = list(product(range(len(humans)), range(2), x_val))
# - デザイン行列の作成：[定数項, 天気, 気温]
X_new = [[1.0, p[1], p[2]] for p in prod]
# - 釣り人インデックスの作成
human_new = [p[0] for p in prod]

# PyMCで期待値 λ のMCMCサンプルを生成
with model_pm:
    # 予測用データのセット
    pm.set_data(
        {
            'X': X_new,                 # 説明変数のデザイン行列
            'human': human_new,         # 釣り人のインデックス
            'Y': range(2000)            # Yも長さを拡張する
        },
        coords={'id': np.arange(2000)}  # coords の id も長さを拡張する
    )
    # 事後予測サンプリングの実行
    post = pm.sample_posterior_predictive(
        idata_pm,
        var_names=['lam'],              # λ の事後分布をサンプリング
        predictions=True                # predictions グループに格納
    )

# λを抽出してNumPy配列として抽出し、[10人, 2天気, 100気温, 4000サンプル] に整形
pred_samples_pm = (
    post.predictions['lam'].to_numpy()
    .reshape(-1, len(humans), 2, 100)
    .transpose(1, 2, 3, 0)
)
print('pred_samples_pm.shape:', pred_samples_pm.shape)

【実行結果】

実は、MCMCサンプル取得後、データを numpy 配列に変換しています。
先ほどの描画ヘルパー関数を利用したいからです。
では描画します。

# p.259 図4.2.1 釣り人別の回帰曲線

# 描画処理
plot_helper_function(pred_samples_pm)

【実行結果】
Bambi や１つ前の描画結果とほぼ同じのチャートです（当たり前 ^^;;;）。

🚀🚀🚀

アディショナルタイム：GPBoost 実装例

非ベイズの GLMM をやってみましょう。
前回記事で利用した GPBoost で今回のランダム切片モデルを実装します。

① 追加インポート

# ライブラリの追加インポート
import gpboost as gpb
gpb.__version__

【実行結果】
今回利用する GPBoost のバージョンは 1.6.1 です（最新版では無いです）。

② GLMM のモデリング
ランダム切片の変数 human を 引数 group_data に設定します。
ランダム切片と説明変数に pandas の DataFrame、目的変数に pandas の Series を渡すと、分析結果に変数名が表示されるのでおすすめです。
このモデリングでは先ほど作成したデザイン行列を再利用しています。

# ポアソンGLMMモデルによる最尤推定：ラプラス近似を用いた（制限付き）最尤法

# 1. データセットの作成
# ランダム効果（グループ）を作成 ※pandas DataFrame
group_data = fish_num_climate_3[['human']]

# 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())

【実行結果】

③ 予測の実行
GPBoost のモデルで予測を行います。
予測用の説明変数データ（デザイン行列等）を作成し、モデル gp_model に対して predict メソッドを適用して予測し、平均値と95%信頼区間を算出してデータフレーム化します。

# 1. 予測用データの作成
# 予測条件の作成 (10人 × 2天気 × 100気温 = 2000行)
x_val = np.linspace(*np.sort(fish_num_climate_3['temperature'])[[0, -1]], 100)
# 定数項 × 釣り人 × 天気 × 気温 の全組み合わせ（直積）の作成
prod = list(product(range(1,2), humans, range(2), x_val))
# データフレーム化
pred_gp_df = pd.DataFrame(
    prod, columns=['Intercept', 'human', 'weather_idx', 'temperature']
)

# 2. GPModel.predict による潜在変数の予測
pred_gp = gp_model.predict(
    X_pred=pred_gp_df[['Intercept', 'weather_idx', 'temperature']],
    group_data_pred=pred_gp_df['human'],
    predict_response=False,               # 線形予測子スケール（log(λ)）を取得
    predict_var=True                      # 分散を取得
)
# 'mu' と 'var' を取得
mu_latent = pred_gp['mu']
var_latent = pred_gp['var']

# 3. 期待値と95%信頼区間の算出 (指数変換)
pred_gp_df['mean_fish'] = np.exp(mu_latent)
pred_gp_df['lower'] = np.exp(mu_latent - 1.96 * np.sqrt(var_latent))
pred_gp_df['upper'] = np.exp(mu_latent + 1.96 * np.sqrt(var_latent))
pred_gp_df['weather'] = pred_gp_df['weather_idx'].map({0: 'cloudy', 1: 'sunny'})

pred_gp_df

【実行結果】
Intercept ～ temperature までが予測に用いた説明変数、mean_fish が予測平均値、lower と upper が予測平均の 95% 信頼区間です。

④ 回帰曲線の可視化
テキスト 図 4.2.1 のような釣り人別・天気別の回帰曲線を描画します。

# 釣り人別の回帰曲線の描画（図4.2.1 ライク）

# 釣り人別散布図の描画
g = sns.relplot(
    data=fish_num_climate_3, x='temperature', y='fish_num',  # データ, x軸, y軸
    col='human', col_wrap=3,                                 # 列の変数, 列の数
    hue='weather', palette=colors,                           # hueの変数、色
    kind='scatter', s=70, height=3, aspect=1.2
)

# 釣り人別の平均釣獲尾数の平均値と95%信頼区間の描画
# sns.relplotのaxesを取り出して、釣り人ごとに描画を繰り返し処理
for human, ax in g.axes_dict.items():
    # 天気ごとに描画を繰り返し処理
    for weather in colors.keys():
        # 対象の釣り人と天気をデータフレームから取得
        query = (pred_gp_df['human']==human) & (pred_gp_df['weather']==weather)
        subset = pred_gp_df[query]
        # 平均釣獲尾数の平均値の描画
        ax.plot(subset['temperature'], subset['mean_fish'], color=colors[weather])
        # 平均釣獲尾数の95%信頼区間の塗りつぶし描画
        ax.fill_between(subset['temperature'], subset['lower'], subset['upper'],
                        color=colors[weather], alpha=0.2)
        # 修飾
        ax.set(title=f'釣り人 {human}', xlabel='', ylabel='')

# チャート全体のx,y軸ラベルの表示
g.figure.supxlabel('気温 [℃]', fontsize=14)
g.figure.supylabel('釣獲尾数', fontsize=14)
g._legend.set_title('天気');

【実行結果】
ベイズモデルとよく似た描画結果になりました。

🚀🚀🚀

今回の記事は以上です。
楽しかったですね！

シリーズの記事

次の記事

前の記事

Stan版

目次

ブログの紹介

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

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

いいなと思ったら応援しよう！

ネイピア DS 応援ありがとうございます。これからもがんばって記事を作成します！

この記事が参加している募集

#夏の1コマ

2.9万件

まで

#

推し活

11万件

#

最近の学び

31万件