「ベイズ統計モデリングによるデータ分析入門」をPythonとStanで写経 ～ Vol.13 正規線形モデル

書籍の著者　馬場真哉 先生

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

書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
今回は「正規線形モデル」を２つのツール・手法で取り組みます。

• 回帰分析による正規線形モデル

• Stanのベイズ正規線形モデル

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

はじめに

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

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

準備

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

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

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

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

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

# インポート

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

# ベイズ統計モデリング
from cmdstanpy import CmdStanModel     # stan
import arviz as az                     # 分析・可視化

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

# デザイン行列
from patsy import dmatrices

# ユーティリティ
import os

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

第７章 正規線形モデル

正規線形モデルとは？

正規線形モデルは一般化線形モデルの文脈において、

• 線形予測子：複数の説明変数を用いる

• リンク関数：恒等関数（変換なし）

• 確率分布：正規分布

で構成される統計モデルです。
ざっくり、重回帰分析と同じだと思います！

🔵🔵🔵

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

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

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

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

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

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

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

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

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

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

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

【実行結果】
売り上げは平均 79.6、最小値 26.1、最大値 96.1 です。
気温は平均 19.9、最小値 10.1、最大値 29.8 です。

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

# 天気別の売り上げ要約統計量
sales_climate.groupby(['weather'])['sales'].describe().round(2)

【実行結果】
天気によって売り上げの平均値は異なっています。

次は天気別気温の要約統計量も念のため見ておきましょう。

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

【実行結果】
天気ごとの気温差はなさそうですね。

データを可視化しましょう。
最初にジョイントプロットで、要約統計量の感覚値を可視化してみます。
seaborn の jointplot を利用します。

# ジョイントプロットの描画
colors = {'cloudy': 'tomato', 'rainy': 'tab:green', 'sunny': 'tab:blue'}
sns.jointplot(
    data=sales_climate, x='temperature', y='sales',
    hue='weather', palette=colors
);

【実行結果】
散布図では気温の上昇と売り上げの増加の関係が見られます。
気温のKDEプロット（上）では、天気別気温の分布は似ています。
売り上げのKDEプロット（右）では、晴れの売り上げ分布が上側にシフトしています。晴れの売り上げは大きいようです。

続いて、テキスト 図 3.7.1 に相当する散布図です。
seaborn の scatterplot を利用します。

# p.209 図3.7.1 ビールの売上と天気・気温の散布図

# 描画領域の設定
plt.figure(figsize=(8, 4))
# 散布図の描画
colors = {'cloudy': 'tomato', 'rainy': 'tab:green', 'sunny': 'tab:blue'}
sns.scatterplot(data=sales_climate, x='temperature', y='sales',
                hue='weather', palette=colors)
# 修飾
plt.title('ビールの売上と気温・天気の関係', loc='left');

【実行結果】
ジョイントプロットで見てしまいましたが…気温の上昇と売り上げの増加の関係が見られます。

この散布図に天気ごとの回帰直線を重ねてみましょう。
seaborn の lmplot を利用します。

# ビールの売上と天気・気温の散布図（回帰直線付き）

# 回帰直線付き散布図
colors = {'cloudy': 'tomato', 'rainy': 'tab:green', 'sunny': 'tab:blue'}
sns.lmplot(data=sales_climate, x='temperature', y='sales',
           hue='weather', palette=colors, height=4, aspect=1.8)
# 修飾
plt.title('ビールの売上と気温・天気の関係（回帰直線付き）', loc='left');

【実行結果】
すべての天気で気温上昇と売り上げ増加の関係が見られます。
晴れの売り上げは他の天気よりも相対的に大きいです。
曇りと雨の気温と売り上げの関係は近い感じです。

🔵🔵🔵

回帰分析による正規線形モデル

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

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

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

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

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

【実行結果】
こちらは（お馴染みの？）回帰分析のサマリーです。
決定係数（R-squared）の 0.629 です。

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

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

【実行結果】

Intercept：切片の推定値（coef）の 20.2 は基準要素「曇り（cloudy）」の売り上げ増分です。
雨の売り上げ増分は「weather[T.rainy]」の推定値 - 3.5 と曇りの推定値を合算した 16.7 です。
晴れの売り上げ増分は「weather[T.sunny]」の推定値 29.5 と曇りの推定値を合算した 49.7 です。
さらに気温が１度上がるごとに売り上げは 2.5 増加します。

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

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

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

④ 天気別気温別平均売上の可視化
最後に天気別気温別平均売上の95%信頼区間を可視化しましょう。
先ほど見た lmplot と似た図になるはずです。

最初に平均売上の予測値を算出します。
回帰分析の結果 res_sm に対して get_prediction メソッドを適用するなどして、予測値を得ます。

## 平均売上予測値の算出
# 予測用データの作成
n_samples = 100
test_data = pd.DataFrame({
    'weather': np.repeat(['cloudy', 'rainy', 'sunny'], n_samples),  # 天気
    'temperature': np.tile(np.linspace(10, 30, n_samples), 3)       # 気温
})
# 予測値の取得
prediction = res_sm.get_prediction(test_data)
pred_df = prediction.summary_frame(alpha=0.05) # alpha=0.05で95%信頼区間
# 予測値データフレームの作成
pred_df = pd.concat([test_data, pred_df], axis=1)
pred_df

【実行結果】

それでは描画します。

## 平均売上の95%信頼区間の描画

# 色の設定
colors = {'cloudy': 'tomato', 'rainy': 'tab:green', 'sunny': 'tab:blue'}
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=sales_climate, x='temperature', y='sales',
                hue='weather', palette=colors)
# 天気ごとに平均売上の予測値と95%信頼区間の描画を繰り返し処理
for weather, color in colors.items():
    # 描画データの取得
    temp, mean_sales, lower, upper = (
        pred_df[pred_df['weather']==weather].iloc[:, [1, 2, 4, 5]].to_numpy().T
    )
    # 平均売上の予測値の描画
    plt.plot(temp, mean_sales, color=color)
    # 平均売上の95%信頼区間の塗りつぶし描画
    plt.fill_between(temp, lower, upper, color=color, alpha=0.2)
# 修飾
plt.title('平均売上の95%信頼区間')
plt.xlabel('気温', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.xticks(range(10, 31, 2));

【実行結果】
lmplot とほぼほぼ同様となりました。
ではベイズ正規線形モデルではどうなるでしょう？

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

🔵🔵🔵

ベイズモデリング

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

① モデルの概要
デザイン行列と係数ベクトルを用いて、次のモデルを実装します。

目的変数 $${\bm Y}$$（sales）は正規分布に従うと仮定しています。
正規分布の平均パラメータ $${\bm \mu}$$ は、デザイン行列 $${\bm X}$$ と係数ベクトル $${\bm \beta}$$ の積です。
係数ベクトル $${\bm \beta}$$、標準偏差パラメータ $${\sigma}$$ には事前分布を明示的に設定しません。

② Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
テキストの Stan コードを引用いたします。
📑ファイル名：3-7-1-lm-design-matrix.stan

data {
    int N;                  // 標本サイズ
    int K;                  // デザイン行列の列数（説明変数の数＋１）
    vector[N] Y;            // 目的変数
    matrix[N, K] X;         // 説明変数
}

parameters {
    vector[K] b;            // 切片を含む係数ベクトル
    real<lower=0> sigma;    // 標準偏差
}

model {
    vector[N] mu = X * b;
    Y ~ normal(mu, sigma);
}

【実行結果】なし

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

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

# formula構文の設定
formula_lm = 'sales ~ weather + temperature'
# 目的変数Y, デザイン行列（説明変数）Xの作成
Y, X = dmatrices(formula_lm, sales_climate, return_type='dataframe')
# デザイン行列の先頭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['sales'], X=X)

【実行結果】なし

④ モデルのコンパイル
モデルのコンパイルを実行します。

%%time
# モデルのコンパイル

# stanプログラムファイルのパス指定
stan_file = '3-7-1-lm-design-matrix.stan'
current_dir = os.path.abspath(os.getcwd())
stan_path = os.path.join(current_dir, 'stan', stan_file)

# モデルオブジェクトの作成（exeの作成）
model = CmdStanModel(stan_file=stan_path)  # stanのpathを設定

【実行結果】

MCMC の準備が整いました！

⑤ MCMC の実行
MCMCを実行しましょう。

%%time
# p.209 MCMCの実行
fit = model.sample(
    data=data_dict_design,   # 対象データ
    seed=1,                  # 乱数の種
    sig_figs=18,             # 出力CSV等に適用する数値精度
)

【実行結果】

⑥ 収束確認
収束の確認をします。
診断メソッド diagnose を利用します。

# 事後分布の診断
print(fit.diagnose())

【実行結果】（１行目のファイルパスは記載省略）
問題は検出されませんでした（no problems detected.）。

CmdStanPy の 出力結果 fit から MCMCサンプルの要約統計量を把握しましょう。

# p.210 結果
fit.summary(percentiles=[2.5, 50, 97.5]).round(2)

【実行結果】
$${\widehat{R}}$$（R_hat）、有効サンプル数（N_Eff）に問題はなさそうです。
b[1] は切片、b[2] は雨、b[3] は晴れ、b[4] は気温の係数に対応しています。

（参考：回帰分析の推定値と誤差の標準偏差）

平均値・係数推定値は近い値になっています。

トレースプロットを描画します。
fit を arviz の idata に変換します。

# arvizのidataに変換
idata = az.from_cmdstanpy(posterior=fit, log_likelihood='lp__')
idata

【実行結果】

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

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

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

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

⑦ 気温別天気別の平均売上の可視化
テキストの図 3.7.2 に相当します。
パラメータのMCMCサンプルを使って、PyMCの外で平均売上の予測値を計算し、arviz の plot_hdi 関数で 95%HDI 区間を描画します。

# p.210 図3.7.2 正規線形モデルから得られた回帰直線

## 平均売上データの作成
# MCMCサンプルからβ0, β1, β2, β3, σを取り出し
intercept, w_rainy, w_sunny, temper, sigmas = fit.draws_pd().iloc[:, -5:].values.T
# x軸の値
x_val = np.linspace(X['temperature'].min(), X['temperature'].max(), 1001)
# MCMCサンプルから曇り、雨、晴れの平均売上を算出
cloudy = (intercept + np.outer(x_val, temper)).T
rainy = (intercept + w_rainy + np.outer(x_val, temper)).T
sunny = (intercept + w_sunny + np.outer(x_val, temper)).T
# 平均売上データと描画色を天気別の辞書に格納
sales = {'cloudy': cloudy, 'rainy': rainy, 'sunny': sunny}
colors = {'cloudy': 'tomato', 'rainy': 'tab:green', 'sunny': 'tab:blue'}

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=sales_climate, x='temperature', y='sales', 
                hue='weather', palette=colors)
# 天気ごとに平均売上の平均値と95%HDI区間の描画を繰り返し処理
for weather in sales.keys():
    # 平均売上の平均値の描画
    plt.plot(x_val, sales[weather].mean(axis=0), color=colors[weather])
    # 平均売上の95%HDI区間の塗りつぶし描画
    az.plot_hdi(x_val, sales[weather].reshape(4, 1000, -1), hdi_prob=0.95,
                color=colors[weather], fill_kwargs={'alpha': 0.2})
# 修飾
plt.title('事後分布サンプル：平均売上の95%HDI区間')
plt.xlabel('気温', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.xticks(range(10, 31, 2));

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

（参考：回帰分析の 95% 信頼区間）

🔵🔵🔵

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

シリーズの記事

次の記事

前の記事

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

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