第８章「階層モデル」

書籍の著者　松浦健太郎 先生

この記事は、テキスト第８章「階層モデル」の8.4節「ロジスティック回帰の階層モデル」の PyMC5写経 を取り扱います。
難易度高めのグラフ描画が出現します。

はじめに

StanとRでベイズ統計モデリングの紹介

この記事は書籍「StanとRでベイズ統計モデリング」（共立出版、「テキスト」と呼びます）のベイズモデルを用いて、PyMC Ver.5で「実験的」に写経する翻訳的ドキュメンタリーです。

テキストは、2016年10月に発売され、ベイズモデリングのモデル式とプログラミングに関する丁寧な解説とモデリングの改善ポイントを網羅するチュートリアル「実践解説書」です。もちろん素晴らしいです！
「アヒル本」の愛称で多くのベイジアンに愛されてきた書籍です！

テキストに従ってStanとRで実践する予定でしたが、RのStan環境を整えることができませんでした（泣）
そこでこのシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「StanとRでベイズ統計モデリング」初版第13刷、著者 松浦健太郎、共立出版

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

PyMC環境の準備

Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。

8.4 ロジスティック回帰の階層モデル

インポート

### インポート

# 数値・確率計算
import pandas as pd
import numpy as np

# PyMC
import pymc as pm
import arviz as az

# ROC-AUC
from sklearn.metrics import roc_curve, roc_auc_score

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

データの読み込み

サンプルコードのデータを読み込みます。
学生個人ごとの属性データです。

### データの読み込み ◆データファイル8.4 data-attendance-4-1.txt
# PersonID:学生ID, A:バイト好き区分(1:好き), Score:学問の興味の強さ（0~200）

data1 = pd.read_csv('./data/data-attendance-4-1.txt')
print('data1.shape: ', data1.shape)
display(data1.head())

【実行結果】

もう１ファイル、データを読み込みます。
説明変数の候補である学生個人のID、科目ID、天気、そして、目的変数となる授業の出欠実績のデータです。

### データの読み込み ◆データファイル8.5 data-attendance-4-2.txt
# PersonID:学生ID, CourseID:科目ID, Weather:天気（A:晴れ、B:曇り、C:雨）
# Y:出欠区分：授業に出席したかどうか（0:欠席、1:出席）

data2 = pd.read_csv('./data/data-attendance-4-2.txt')
print('data2.shape: ', data2.shape)
display(data2.head())

【実行結果】

8.4.4 Stanで実装

PyMCでモデル式8-8を実装します。
データの前処理から。

### データの前処理：天気の重み係数を設定
Weather_dict2 = {'A': 0, 'B': 0.2, 'C': 1}
data2['Weather_w'] = data2['Weather'].map(Weather_dict2)
data2.head()

【実行結果】

続いてモデルの定義です。

### モデルの定義 ◆モデル式8-8 model8-8.stan

with pm.Model() as model:
    
    ### データ関連定義
    ## coordの定義
    model.add_coord('data', values=data2.index, mutable=True)
    model.add_coord('person', values=data1['PersonID'].values, mutable=True)
    model.add_coord('course', values=sorted(data2['CourseID'].unique()),
                    mutable=True)
    model.add_coord('coef', values=[1, 2, 3, 4], mutable=True)

    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data2['Y'].values, dims='data')
    # 説明変数 A
    A = pm.ConstantData('A', value=data1['A'].values, dims='person')
    # 説明変数 Score
    Score = pm.ConstantData('Score', value=data1['Score'].values / 200,
                            dims='person')
    # 説明変数 Weather
    Weather = pm.ConstantData('Weather', value=data2['Weather_w'].values,
                             dims='data')
    # インデックス PersonIdx data2のPersonIDを0始まりに変換
    PersonIdx = pm.ConstantData('PersonIdx',
                                value=data2['PersonID'].values - 1,
                                dims='data')
    # インデックス CourseIdx data2のCourseIDを0始まりに変換
    CourseIdx = pm.ConstantData('CourseIdx',
                                value=data2['CourseID'].values - 1,
                                dims='data')
    
    ### 事前分布
    # 回帰係数
    b = pm.Uniform('b', lower=-100, upper=100, dims='coef')
    # 学生に依存する項
    sigmaP = pm.Uniform('sigmaP', lower=0, upper=100)
    bP = pm.Normal('bP', mu=0, sigma=sigmaP, dims='person')
    xP = pm.Deterministic('xP', b[1]*A + b[2]*Score + bP, dims='person')
    # 科目に依存する項
    sigmaC = pm.Uniform('sigmaC', lower=0, upper=100)
    bC = pm.Normal('bC', mu=0, sigma=sigmaC, dims='course')
    xC = pm.Deterministic('xC', bC, dims='course')
    # 授業に依存する項
    xJ = pm.Deterministic('xJ', b[3]*Weather, dims='data')

    ### 線形予測子 x  ※切片はこちらに集約
    x = pm.Deterministic('x', b[0] + xP[PersonIdx] + xC[CourseIdx] + xJ,
                         dims='data')
    ### 確率パラメータq ※線形予測子を逆ロジット変換
    q = pm.Deterministic('q', pm.invlogit(x), dims='data')

    ### 尤度関数
    obs = pm.Bernoulli('obs', p=q, observed=Y, dims='data')

モデルの定義内容を見ます。

### モデルの表示
model

【実行結果】

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

【実行結果】
かなり巨大なモデルになりました。
学生個人の層（person）と科目の層（course）が分かれていて、それぞれが事前分布の階層構造をもっています。

MCMCを実行します。

### 事後分布からのサンプリング 2分30秒
with model:
    idata  = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                       nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

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

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

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['b', 'bP', 'bC', 'sigmaP', 'sigmaC']
pm.summary(idata, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

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

### トレースプロットの表示
pm.plot_trace(idata, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

8.4.5 推定結果の解釈

事後分布の要約統計量を算出します。
算出関数を定義します。

### mean,sd,2.5%,25%,50%,75%,97.5%パーセンタイル点をデータフレーム化する関数の定義
def make_stats_df(y):
    probs = [2.5, 25, 50, 75, 97.5]
    columns = ['mean', 'sd'] + [str(s) + '%' for s in probs]
    quantiles = pd.DataFrame(np.percentile(y, probs, axis=0).T, index=y.columns)
    tmp_df = pd.concat([y.mean(axis=0), y.std(axis=0), quantiles], axis=1)
    tmp_df.columns=columns
    return tmp_df

要約統計量を算出します。

### 要約統計量の算出・表示 ※図8.9左のパラメータ
vars = ['sigmaP', 'sigmaC']
param_samples = idata.posterior[vars].to_dataframe().reset_index(drop=True)
b_samples_df = pd.DataFrame(
            idata.posterior.b.stack(sample=('chain', 'draw')).data.T,
            columns=[f'b[{i+1}]' for i in range(4)])
param_samples = pd.concat([b_samples_df, param_samples], axis=1)
display(make_stats_df(param_samples).round(2))

【実行結果】

変数の内容は次のとおりです。

### 要約統計量の算出・表示 ※図8.9右のパラメータ
bC_samples_df = pd.DataFrame(
            idata.posterior.bC.stack(sample=('chain', 'draw')).data.T,
            columns=[f'bC[{i+1}]' for i in range(10)])
display(make_stats_df(bC_samples_df).round(2))

【実行結果】
bCは科目別の変数であり、線形予測子の科目差（xC）と同値です。

パラメータ$${b,\ \sigma_P,\ \sigma_C}$$の事後分布をヴァイオリンプロットで描画します。
テキスト図8.9左に相当します。

### パラメータb,sigmaP,sigmaCの事後分布のヴァイオリンプロットの描画 ◆図8.9左

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
b_samples = idata.posterior.b.stack(sample=('chain', 'draw')).data
s_P_samples = idata.posterior.sigmaP.stack(sample=('chain', 'draw')).data
s_C_samples = idata.posterior.sigmaC.stack(sample=('chain', 'draw')).data
# データフレーム化
params_df = pd.DataFrame({
    'b1': b_samples[0],
    'b2': b_samples[1],
    'b3': b_samples[2],
    'b4': b_samples[3],
    's_P': s_P_samples,
    's_C': s_C_samples
})
display(params_df)

## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 5))
# ヴァイオリンプロットの描画
sns.violinplot(data=params_df, orient='h', fill=False, inner=None, linewidth=5,
               alpha=0.5, ax=ax)
# エラーバー付きポイントプロットの描画（ヴァイオリンプロットの中）
# estimatorで中央値、errorbarで95パーセンタイルインターバル(pi)を指定
sns.pointplot(data=params_df, estimator='median', errorbar=('pi', 95),
              orient='h', linestyle='none', color='tab:red', ms=10, ax=ax)
# x=0の垂直点線の描画
ax.axvline(0, color='black', ls='--')
# 修飾
ax.set(xlabel='value', ylabel='Parameter')
plt.grid(lw=0.5);

【実行結果】

続いてパラメータ$${b_C}$$の事後分布をKDEプロットで描画します。
テキスト図8.9右に相当します。

### パラメータb_Cの事後分布のKDEプロットの描画 ◆図8.9右

## 描画用データの作成
# 推論データからパラメータのサンプリングデータを取得
b_C_samples = idata.posterior.bC.stack(sample=('chain', 'draw')).data

## 描画処理
# カラーマップの取得（線の色分けに利用）
cmap = plt.get_cmap('tab10')
# 描画領域の指定
plt.figure(figsize=(15, 5))
# 科目ごとのb_C[n]のKDEプロットの描画
for i, b_C in enumerate(b_C_samples):
    # 色の設定
    color = cmap(i / 10)
    # x軸の値を設定
    x_ticks = np.linspace(b_C.min(), b_C.max(), 1001)
    # KDEプロットの描画
    kde_plot = sns.kdeplot(b_C, color=color)
    # KDEプロットで描画した線から、x軸,y軸の値を取得
    kde_data = kde_plot.get_lines()[i].get_data()
    # KDEプロットのyの最大値のインデックスを取得
    kde_max_idx = kde_data[1].argmax()
    # KDEプロットのyの最大値≒MAP推定値の垂直線の描画
    plt.vlines(x=kde_data[0][kde_max_idx], ymin=0, ymax=kde_data[1][kde_max_idx],
               color=color, ls='--')
    # KDEプロットのyの最大値≒MAP推定値のy=0付近のx点の描画
    plt.scatter((kde_data[0][kde_max_idx]), (-0.05), marker='s', s=20,
                color=color)

# 修飾
plt.xlabel('value')
plt.title(r'科目差のパラメータ $\beta_C$ の事後分布：KDEプロット')
# plt.ylim(-3, 45)
plt.grid(lw=0.5)

【実行結果】

蛇足でROC曲線を描きます。
描画用データの作成から。

### ROC曲線の描画

## 描画用データの作成
# 確率qの事後分布サンプルデータの取り出し
q_samples = idata.posterior.q.stack(s=('chain', 'draw')).data
# 描画項目 Y, q_median をデータフレーム化
plot_data = data2[['Y']]
# qの事後分布の中央値をデータフレームに追加
plot_data['q'] = np.median(q_samples, axis=1)
# qの80%区間をデータフレームに追加
plot_data[['q10%', 'q90%']] = np.quantile(q_samples, q=[0.1, 0.9], axis=1).T
display(plot_data)

【実行結果】

描画処理です。

## ROC曲線描画用のデータ作成
# ROC曲線のx軸：fpr、y軸：tprの算出
fpr, tpr, threshold = roc_curve(plot_data['Y'], plot_data['q'])
fpr10, tpr10, threshold10 = roc_curve(plot_data['Y'], plot_data['q10%'])
fpr90, tpr90, threshold90 = roc_curve(plot_data['Y'], plot_data['q90%'])
# ROC-AUCスコアの算出
q_score = roc_auc_score(plot_data['Y'], plot_data['q'])
q10_score = roc_auc_score(plot_data['Y'], plot_data['q10%'])
q90_score = roc_auc_score(plot_data['Y'], plot_data['q90%'])

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 5))
ax = plt.subplot()
# ROC曲線の青線の描画
ax.plot(fpr, tpr)
# 80%ベイズ信頼区間の描画（ただし、10%と90%のy軸が一致していない為、概算です）
ax.fill_betweenx(tpr, fpr10[:-1], fpr90, alpha=0.3)
# 対角線の赤点線の描画
ax.plot([0, 1], [0, 1], color='tab:red', ls='--')
# 修飾
ax.set(xlabel='False Positive: FPR', ylabel='True Positive: TPR',
       title=f'ROC曲線\nAUC: '
             f'{q_score:.4f} [{q90_score:.4f}, {q10_score:.4f}]')
plt.grid(lw=0.5);

【実行結果】
AUCの値はテキストと微妙に違っています（汗）

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

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