「ベイズ統計モデリングによるデータ分析入門」をPythonとPyMCで写経 ～ Vol.12 ダミー変数と分散分析モデル

書籍の著者　馬場真哉 先生

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

書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
今回は説明変数がダミー変数の「分散分析モデル」を３つのツール・手法で取り組みます。

• 回帰分析による分散分析モデル

• Bambiのベイズ分散分析モデル

• PyMCのベイズ分散分析モデル

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

はじめに

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

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

準備

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

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

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

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

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

# インポート

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

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

# 統計モデリング
import statsmodels.formula.api as smf

# デザイン行列
from patsy import dmatrices

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

第６章 ダミー変数と分散分析モデル

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

テキスト 6.4 節に相当します。
テキストの仮想のビールの売り上げデータを引用いたします。
csv ファイルを pandas のデータフレーム形式で変数 sales_weather に読み込みます。

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

# ファイルの読み込み
sales_weather = pd.read_csv('./data/3-6-1-beer-sales-3.csv')

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

【実行結果】
標本サイズ 150、変数 sales は 売り上げ（単位：万円）、weather は天気の種類（晴れ、曇り、雨）です。
天気と売り上げの関係を分散分析モデルで分析します。

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

# 天気ごとのデータ件数
sales_weather['weather'].value_counts().to_frame()

【実行結果】
曇り、晴れ、雨それぞれ 50 件のデータがあります。

データの要約統計量を確認します。
まずは全体です。

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

【実行結果】
売り上げは平均 70、最小値 25、最大値 116 です。範囲が広い感じ。

次は天気別売り上げの要約統計量です。

# 天気別の要約統計量
sales_weather.groupby(['weather'])['sales'].describe().round(2)

【実行結果】
平均 mean を見ると、曇りと雨の売り上げは約 63 万円、晴れのときは＋20 万円の約 83 万円です。

データを可視化しましょう。
テキスト 図 3.6.1 に相当するバイオリンプロットです。
seaborn ライブラリを利用します。

# p.203 図3.6.1 ビールの売上と天気のバイオリンプロット

# バイオリンプロット
sns.violinplot(data=sales_weather, x='weather', y='sales', hue='weather',
               fill=False)
# 散布図（ストリッププロット）
sns.swarmplot(data=sales_weather, x='weather', y='sales', hue='weather')
# 修飾
plt.title('観測データ：天気別のバイオリンプロット');

【実行結果】
晴れの売り上げが大きい感じを掴めると思います。

天気と売り上げの相関を見てみましょう。
質的変数 vs 量的変数なので「相関比」を使ってみます。
相関比算出関数を定義します。

# 相関比算出関数の定義
# ※引数: x=量的変数（numpy配列）, gp=質的変数（numpy配列）

def correlation_ratio(x, gp):
    # 全変動: 全データの偏差平方和を計算
    ST = ((x - x.mean())**2).sum()
    # 群内変動: 群gごとに偏差平方和を計算して総合計する
    SW = sum(((x[gp==g] - x[gp==g].mean())**2).sum() for g in np.unique(gp))
    # 相関比 = 1 - 群内変動 / 全変動 = 群間変動 / 全変動
    eta2 = 1 - SW / ST
    # 戻り値　　相関比:eta2, 総変動:ST, 群内変動:SW, 群間変動:SB
    return {'eta2': eta2, 'ST': ST, 'SW': SW, 'SB': ST-SW}

データの相関比を計算しましょう。

# 相関比の算出 ※決定係数 R² と一致する
result = correlation_ratio(sales_weather['sales'], sales_weather['weather'])
print(f'相関比η² = {result["eta2"]:.3f}')

【実行結果】
相関比は 0.247。
天気と売り上げとの間には弱い正の相関があるようです。

🚀🚀🚀

回帰分析による分散分析モデル

① モデルの概要
まず分散分析モデルを回帰分析で確認しておき、あとでベイズ統計モデルの結果と比べてみましょう。
Python の統計ライブラリ statsmodels を利用して、次のモデルを実装します！

$${\varepsilon}$$ は誤差です。
曇りを基準要素とし、ダミー変数から除外します。

この分散分析モデルを statsmodels の最小二乗法 olsにあたえる formula 構文に変換します。
$${\sim}$$ を挟んで、左辺は目的変数、右辺は説明変数です。

質的変数（Categorical variable）の天気 weather を C() で括る書き方もあります。

② 回帰分析の実行
では回帰分析を実行します！

# 回帰分析 by statsmodels
formula = 'sales ~ weather'
res_sm = smf.ols(formula=formula, data=sales_weather).fit()
res_sm.summary()

【実行結果】
こちらは（お馴染みの？）回帰分析のサマリーです。
決定係数（R-squared）の 0.247 は相関比と一致しています。
一元配置分散分析モデルの場合、データの相関比と回帰モデル（切片付き）の決定係数は一致するようですね！

③ 推定値の確認
係数の推定値にフォーカスしてみます。

# 要約表から係数の推定値の部分を取り出し
res_sm.summary().tables[1]

【実行結果】

Intercept：切片の推定値（coef）の 63.0 は基準要素「曇り（cloudy）」の平均売り上げです。
雨の平均売り上げは、「weather[T.rainy]」の推定値 - 0.3 と曇りの推定値を合算した 62.7 です。
晴れの平均売り上げは、「weather[T.sunny]」の推定値 20.0 と曇りの推定値を合算した 83.0 です。

誤差 $${\varepsilon}$$ の標準偏差を確認しましょう。
回帰分析の結果 res_sm の scale 属性で誤差分散を取り出して正の平方根をとります。

# 誤差の標準偏差
print('誤差の標準偏差:', np.sqrt(res_sm.scale))

【実行結果】
誤差の標準偏差は 16.78 です。
ベイズ統計モデルの観測値 sales の標準偏差と照らし合わせましょう。

④ 天気別平均売上の可視化
最後に天気別平均売上の95%信頼区間を可視化しましょう。
回帰分析の結果 res_sm に対して get_prediction メソッドを適用するなどして、予測値を得ます。

# 天気別平均売上の95%信頼区間をエラーバーで描画

# 予測用データの作成（各グループ1つずつ）
test_data = pd.DataFrame({'weather': ['cloudy', 'rainy', 'sunny']})

# 予測値を取得
prediction = res_sm.get_prediction(test_data)
pred_summary = prediction.summary_frame(alpha=0.05) # alpha=0.05で95%信頼区間

# 95%信頼区間付きエラーバーの描画
fig, ax = plt.subplots()
ax.errorbar(test_data['weather'], pred_summary['mean'], 
            yerr=[pred_summary['mean'] - pred_summary['mean_ci_lower'], 
                  pred_summary['mean_ci_upper'] - pred_summary['mean']],
            fmt='o', capsize=5)
plt.ylabel('sales')
plt.title('天気別平均売上の95%信頼区間')
ax.margins(x=0.3)
plt.show()

【実行結果】
「曇り・雨」と「晴れ」の信頼区間に重なりはなく、大きな差があります。
晴れると売り上げが上がる感覚をイメージできます。
ではベイズ分散分析モデルではどうなるでしょう？

ベイズ流の分散分析モデルへ進みます。

🚀🚀🚀

ベイズモデリング by Bambi

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

① モデルの概要
次のモデルを実装します。
先ほどの statsmodels の formula と一緒です。

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

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

# モデリング

# 無情報事前分布（想定）の設定
uninformed_prior = bmb.Prior('Normal', mu=0, sigma=1e5)
sigma_prior = bmb.Prior('HalfNormal', sigma=1e5)
# 事前分布を辞書にとりまとめ
priors = {'Intercept': uninformed_prior,
          'weather': uninformed_prior,
          'sigma': sigma_prior}

# モデルの定義
model_bmb = bmb.Model(
    formula='sales ~ weather',       # フォーミュラ式
    data=sales_weather,              # データ
    priors=priors,                   # 事前分布（辞書）
)

【実行結果】なし

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

# モデルの表示
model_bmb

【実行結果】
Intercept, temperarure, sigma の事前分布には先ほど設定した内容が表示されています。

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

# モデルの可視化
model_bmb.build()
model_bmb.graph()

【実行結果】
weather 係数の次元は２です（weather_dim (2) の 2 ）。
ダミー変数化するときに要素の１つを除外しています。

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

%%time
# MCMCの実行
idata_bmb = model_bmb.fit(
    draws=1000, tune=1000, chains=4, random_seed=123, nuts_sampler='nutpie'
)

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

④ 収束確認
収束の確認をします。
MCMC サンプルの要約表を表示します。

# p.203 要約統計量の表示
var_names = ['Intercept', 'weather', '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 はほぼ重なっており、かつ１峰です。
右側のチャートがゲジゲジしています。

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

⑤ 天気別の平均売上の可視化
テキストの図 3.6.2 に相当します。
Bambi モデル model_bmb に対して predict メソッドを適用することで平均売上（mu の事後分布に相当）のMCMCサンプルを取得します。

# 曇り、雨、晴れの平均売上の予測値の算出

# 予測用データの作成（データフレーム）
weathers = ['cloudy', 'rainy', 'sunny']
new_data = pd.DataFrame({'weather': weathers})

# 予測の実行
model_bmb.predict(idata=idata_bmb, data=new_data, kind='response_params')

# 平均売上の事後分布サンプルデータのデータフレーム化
sales_weather_result = pd.DataFrame(
    az.extract(idata_bmb.posterior).mu.data.T, columns=new_data['weather'])
sales_weather_result

【実行結果】
天気ごとに 4000 個の平均売上のMCMCサンプルを取得しました。

では可視化を実行します。
seaborn の pointplot を利用します。

# p.204 図3.5.3 回帰直線95%信用区間付き

# 天気別平均売上のエラーバープロットの描画
sns.pointplot(data=sales_weather_result, errorbar=('pi', 95), capsize=0.2,
              linestyle='none')
# 修飾
plt.title('事後分布サンプル：天気別平均売上の95%信用区間')
plt.xlabel('天気', fontsize=12)
plt.ylabel('売上金額', fontsize=12);

【実行結果】
晴れの平均売上は他の天気と比べて 20 万円くらい大きいですね！

⑥ 平均売上の事後分布プロットの描画
天気別平均売上の分布を可視化しましょう。

# 平均売上の事後分布プロット

# 天気名を付与
mu_plot_data = (idata_bmb.posterior.mu
                .assign_coords({'__obs__': new_data['weather']}))

# 事後予測のプロット
az.plot_posterior(mu_plot_data, hdi_prob=0.95);

【実行結果】

🚀🚀🚀

ベイズモデリング by PyMC

テキスト 6.6、6.7 節に相当します。
デザイン行列を作成して PyMC でモデリングします。

① モデルの概要
次のモデルを実装します。

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

# p.205 デザイン行列の作成

# formula構文の設定
formula_anova = 'sales ~ weather'
# デザイン行列（説明変数）Xの作成
Y_dm, X_dm = dmatrices(formula_anova, sales_weather, return_type='dataframe')
# デザイン行列の先頭5行の表示
X_dm.head()

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

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

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

③ モデル定義
上式のモデルを PyMC で記述します。

# モデリング

# coordsの設定
coords = {'id': sales_weather.index.values,          # 観測データの識別子
          'coefs': ['intercept', 'rainy', 'sunny']}  # 係数ベクトルbの識別子

# モデルの定義
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'))

    ## 事前分布: 無情報事前分布的な分布
    # 係数ベクトル
    b = pm.Normal('b', mu=0, sigma=1e5, dims='coefs')
    # 標準偏差
    sigma = pm.HalfNormal('sigma', sigma=1e5)

    ## 尤度関数: 正規分布を仮定
    mu = pm.Deterministic('mu', X @ b, dims='id')
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='id')

【実行結果】なし

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

# モデルの表示
model_pm

【実行結果】

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

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

【実行結果】

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

%%time
# MCMCの実行 ※高速化したいときは nuts_sampler='nutpie'
with model_pm:
    idata_pm = pm.sample(draws=1000, tune=1000, chains=4, random_seed=1,
                         nuts_sampler='nutpie')

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

⑤ 収束確認
収束の確認をします。
MCMC サンプルの要約表を表示します。

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

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

（参考：Bambiモデルの要約統計量）

パラメータの事後分布はBambiと似た値になっています。
ただ…有効サンプル数（ess_bulk、ess_tail）がBambiよりも小さくなっている点は気になりますね…

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

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

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

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

⑥ 天気別の平均売上の可視化
テキストの図 3.6.2 に相当します。
sample_posterior_predictive 関数を用いて、平均売上（mu の事後分布に相当）のMCMCサンプルを取得します。

# 曇り、雨、晴れの平均売上の予測値の算出

# 予測用データの作成
# 説明変数Xの予測用データの設定
X_new = [[1, 0, 0],  # cloudy
         [1, 1, 0],  # rainy
         [1, 0, 1]]  # sunny
# idの追加値の設定
id_position = range(len(X_dm), len(X_dm) + len(X_new))  # range(150, 153)

# 平均売上muの予測の実行
with model_pm:
    # Xに新しいデータの設定
    pm.set_data({'X': X_new}, coords={'id': id_position})
    # muのout-of-sample 予測の実行
    predictions = pm.sample_posterior_predictive(
        idata_pm, var_names=['mu'], random_seed=123, predictions=True)

# 平均売上の事後分布サンプルデータのデータフレーム化
sales_weather_result = pd.DataFrame(
    az.extract(predictions.predictions).mu.data.T, columns=weathers
)
sales_weather_result

【実行結果】
天気ごとに 4000 個の平均売上のMCMCサンプルを取得しました。

では可視化を実行します。

# p.204 図3.5.3 回帰直線95%信用区間付き

# 天気別平均売上のエラーバープロットの描画
sns.pointplot(data=sales_weather_result, errorbar=('pi', 95), capsize=0.2,
              linestyle='none')
# 修飾
plt.title('事後分布サンプル：天気別平均売上の95%信用区間')
plt.xlabel('天気', fontsize=12)
plt.ylabel('売上金額', fontsize=12);

【実行結果】
晴れの平均売上は他の天気と比べて 20 万円くらい大きいですね！

⑦ 平均売上の事後分布プロットの描画
天気別平均売上の分布を可視化しましょう。

# 平均売上の事後分布プロット

# 天気名を付与
mu_plot_data = (predictions.predictions.mu
                .assign_coords({'id': new_data['weather']}))

# 事後予測のプロット
az.plot_posterior(mu_plot_data, hdi_prob=0.95);

【実行結果】

🚀🚀🚀

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

シリーズの記事

次の記事

前の記事

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

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