「ベイズ統計モデリングによるデータ分析入門」をPythonとStanで写経 ～ Vol.16 交互作用

書籍の著者　馬場真哉 先生

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

書籍の第３部は「一般化線形モデル」のベイズ統計モデリングです。
今回は３種類の交互作用項を含む線形回帰モデルを Stan で取り組みます。

• カテゴリ変数 × カテゴリ変数の交互作用

• カテゴリ変数 × 量的変数の交互作用

• 量的変数 × 量的変数の交互作用

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

はじめに

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

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

準備

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

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

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

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

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

# インポート

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

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

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

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

# ユーティリティ
import os

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

第10章 交互作用 by Stan

交互作用とは？

テキスト p.228 によると、交互作用は「説明変数同士が交互に影響を与え合う」ことです。
線形予測子において、説明変数同士の積で表現します。

テキストが正規線形モデル（第７章）と呼ぶ「線形回帰モデル」（重回帰モデル）で交互作用を３つのケースで体感いたします。

線形回帰モデルとは？

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

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

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

• 確率分布：正規分布

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

🔵🔵🔵

ケース１：カテゴリ変数 × カテゴリ変数

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

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

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

# ファイルの読み込み
interaction_1 = pd.read_csv('./data/3-10-1-interaction-1.csv')

# 結果の表示
# sales:売上, publicity: 宣伝の有無, bargen: 安売りの有無
print('interaction_1.shape: ', interaction_1.shape)
interaction_1.head(3)

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

• sales：売上金額 （単位：万円）

• publicity：宣伝の有無（あり：to_implement、なし：not）

• bargen：安売りの有無（あり：to_implement、なし：not）

publicity と bargen は二値のカテゴリ変数です。
sales は 量的変数です。

「定数項 + 宣伝の有無 + 安売りの有無 + 宣伝の有無 × 安売りの有無」を説明変数にした線形回帰モデルで分析します。

カテゴリ変数の要素ごとのデータ件数をカウントします。
pandas の value_counts メソッドを利用します。

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

【実行結果】
宣伝なし／あり、安売りなし／ありで、それぞれ 50 件のデータがあります。

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

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

【実行結果】
売上は平均 127.2、最小値 55.7、最大値 191.7 です。

次は売上平均値のクロス集計表です。

# クロス集計表：売上の平均値
cross_table = interaction_1.pivot_table(
    index='bargen', columns='publicity', values='sales', aggfunc='mean')
cross_table

【実行結果】
宣伝あり＆安売りあり（右下セル）が最も売上平均値が大きいです。
安売りなし→安売りありの売上平均値の増加は大きい感じがします。

クロス集計表を可視化しましょう。
そうです、ヒートマップです。
seaborn の heatmap を利用します。

# 可視化
sns.heatmap(cross_table, annot=True, fmt='.2f', annot_kws={'fontsize': 16},
            cmap='Greens', vmin=0)
plt.xlabel('宣伝の有無', fontsize=12)
plt.ylabel('安売りの有無', fontsize=12);

【実行結果】
色の濃さで売上平均値の大きさが直感的に分かりますね！

どんどんデータを可視化しましょう。
次はポイントプロットで、売上の平均値と95%信頼区間を比べてみます。
seaborn の pointplot を利用します。

# ポイントプロット：平均と95%信頼区間

# 色の設定
colors = {'not': 'tomato', 'to_implement': 'tab:blue'}
# ポイントプロットの描画
sns.pointplot(data=interaction_1, x='publicity', y='sales',
              hue='bargen', palette=colors, capsize=0.1)
# 修飾
plt.xlabel('宣伝の有無', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.legend(title='安売りの有無');

【実行結果】
横軸：宣伝の有無、縦軸：売上金額、第３の軸：安売りの有無です。
中央の点は売上平均値、ヒゲは売上平均値の95%信頼区間です。
宣伝と安売りを併用すると売上増分（傾き）は大きくなっています。
「宣伝 × 安売り」の交互作用を予感させますね！

続いて、箱ひげ図で分布を確認します。
seaborn の boxplot を利用します。

# 箱ひげ図＋スウォームプロット

# 色の設定
colors = {'not': 'tomato', 'to_implement': 'tab:blue'}
# 箱ひげ図の描画
sns.boxplot(data=interaction_1, x='publicity', y='sales',
            hue='bargen', palette=colors, fill=False,
            showmeans=True, meanprops={'ms': 10})
# スウォームプロットの描画
sns.swarmplot(data=interaction_1, x='publicity', y='sales', 
              hue='bargen', palette=colors, dodge=True, legend=False)
# 修飾
plt.xlabel('宣伝の有無', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.legend(title='安売りの有無');

【実行結果】
宣伝の有無・安売りの有無が売上金額の大きさに影響している感じがします。
また、「宣伝なし＋安売りあり」「 宣伝あり＋安売りなし」は分布がほぼ被っていることも分かります。

🔵🔵🔵

線形回帰モデル by statsmodels

① モデルの概要
statsmodels の 線形回帰 OLS で次のモデルを実装します。

$${\varepsilon}$$ は誤差です。
各ダミー変数は、あり：$${1}$$、なし：$${0}$$ の値を取ります。

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

この formula 文は省略形です。次の formula 文と同等です。

② 線形回帰の実行
では線形回帰を実行します！

# 線形回帰モデル by statsmodels

# formulaの設定
formula = 'sales ~ publicity * bargen'
# 回帰分析の実行
res_sm1 = smf.ols(formula=formula, data=interaction_1).fit()
res_sm1.summary()

【実行結果】
こちらは線形回帰のサマリーです。
自由度調整済み決定係数（Adj. R-squared）は 0.592 です。
まあまあの当てはまり具合でしょうか。

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

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

【実行結果】
最下行が交互作用項です。5%水準で有意です（$${p = 0.005}$$）。
宣伝あり・安売りありの場合、平均売上が追加的に 20.8 万円大きくなる、と解釈できます。
宣伝なし・安売りなしの平均売上は切片の 103.4 万円と推定されます。

④ 宣伝の有無別・安売りの有無別の平均売上の可視化
最後に宣伝の有無別・安売りの有無別の平均売上の95%信頼区間を可視化しましょう。

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

# 平均売上の予測値の算出

# 予測用データの作成
test_data1 = pd.DataFrame({
    'publicity': np.repeat(['not', 'to_implement'], 2),   # 宣伝の有無
    'bargen': np.tile(['not', 'to_implement'], 2)         # 安売りの有無
})
# 予測値の取得
prediction1 = res_sm1.get_prediction(test_data1)
pred_df1 = prediction1.summary_frame(alpha=0.05) # alpha=0.05で95%信頼区間
# 予測値データフレームの作成
pred_df1 = pd.concat([test_data1, pred_df1], axis=1)
pred_df1

【実行結果】

それでは描画します。

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

# 設定
colors = {'not': 'tomato', 'to_implement': 'tab:blue'}  # 安売り有無の色
x = np.array([-0.15, 0.15, 0.85, 1.15])                 # errorbarのx軸位置

# 95%信頼区間付きエラーバーの描画
plt.figure(figsize=(8, 5))
for i, row in pred_df1.iterrows():
    plt.errorbar(x[i], row['mean'], 
                 yerr=[[row['mean'] - row['mean_ci_lower']], 
                       [row['mean_ci_upper'] - row['mean']]],
                 fmt='o', ms=8, capsize=12, elinewidth=2, capthick=2,
                 color=colors[row['bargen']])
# スウォームプロットの描画
sns.swarmplot(data=interaction_1, x='publicity', y='sales',
              hue='bargen', palette=colors);
# 凡例処理
plt.legend(title='安売りの有無', bbox_to_anchor=(1, 1))
# 修飾
plt.title('平均売上の95%信頼区間', fontsize=14)
plt.xlabel('宣伝の有無', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.margins(x=0.3);

【実行結果】
宣伝あり・安売りありのエラーバーの位置が断然高い＝平均売上が大きいことが一目瞭然です。
交互作用の効果を実感できますね！

ベイズ流の線形回帰モデルへ進みます。

🔵🔵🔵

ベイズモデリング by Stan

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

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

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

② Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
テキストの Stan コードを引用いたします。
📑ファイル名：3-10-1-glm-intract-1-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.230 デザイン行列の作成

# formula構文の設定
formula_interaction1 = 'sales ~ publicity * bargen'
# 目的変数Y, デザイン行列（説明変数）Xの作成
Y, X = dmatrices(formula_interaction1, interaction_1, return_type='dataframe')
# int型の設定
X = X.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['sales'], X=X)

【実行結果】なし

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

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

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

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

【実行結果】

MCMC の準備が整いました！

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

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

【実行結果】

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

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

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

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

# p.231 結果
fit1_summary = fit1.summary(percentiles=[2.5, 50, 97.5])
fit1_summary.round(2)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（N_Eff）に問題はなさそうです。
５行目が交互作用項です。宣伝あり・安売りありの場合、平均売上が追加的に 20.4 万円大きくなる、と解釈できます。
１行目の切片から、宣伝なし・安売りなしの平均売上は切片の 103.2 万円と推定されます。

statsmodels の線形回帰の推定値の表を並べましょう。
（参考：statsmodels の推定値）

手法間で係数推定値は近い値になっています。

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

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

【実行結果】

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

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

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

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

⑦ 係数の解釈
テキスト 10.6 節の計算に準拠して進めます。

a. 説明変数を作る
売上（平均）$${\mu}$$ の予測値を算出するために、予測用の説明変数およびデザイン行列を作成します。

# p.232 説明変数を作る

# publicityとbargenのデータフレームを作成
newdata_1 = pd.DataFrame({
    'publicity': ['not', 'to_implement', 'not', 'to_implement'],
    'bargen': ['not', 'not', 'to_implement', 'to_implement'],
})
print('【説明変数】')
display(newdata_1)

# デザイン行列化
newdata_1_dm = dmatrix('publicity * bargen', newdata_1).base
print('【デザイン行列】')
print(newdata_1_dm)

【実行結果】
宣伝の有無と安売りの有無の組み合わせです。

b. 係数推定値から平均売上を計算
テキスト 表 3.10.2 に相当します。
予測用のデザイン行列と係数推定値（平均値）の積で平均売上を計算します。

# p.231 表 3.10.2 カテゴリ×カテゴリの交互作用があるときの予測値の変化のパターン
pd.concat(
    [newdata_1,
     pd.Series(newdata_1_dm @ fit1_summary.iloc[[1, 2, 3, 4], 0], name='mean')],
    axis=1
).round(3)

【実行結果】

c. 平均売上 $${\mu}$$ の予測値の要約統計量の計算
まず、予測用のデザイン行列とパラメータの MCMC サンプルの行列積で $${\mu}$$ の予測値の MCMC サンプルを取得します。
次に予測値のMCMCサンプルを arviz の summary 関数に与えて、要約統計量を得ます。

# p.232 売上の予測値の算出

# モデルの係数を取得
b_samples1 = az.extract(idata1.posterior).b.values

# MCMCサンプルに基づく売上平均の予測値の算出 shape=(4, 4000)
linear_fit1 = (newdata_1_dm @ b_samples1)

# 統計量の算出 az.summary()を利用
for i in range(len(linear_fit1)):
    tmp_df = az.summary(linear_fit1[i, :], kind='stats', hdi_prob=0.95)
    tmp_df.index = [i]
    if i==0:
        stats_df1 = tmp_df
    else:
        stats_df1 = pd.concat([stats_df1, tmp_df], axis=0)
# newdataと結合
stats_df1 = pd.concat([newdata_1, stats_df1], axis=1)

# 結果の表示
stats_df1.round(3)

【実行結果】
mean が平均売上 $${\mu}$$ の予測値です。
b. の結果とほぼ同じになりました。

⑧ 宣伝の有無別・安売りの有無別の平均売上の可視化
テキストの図 3.10.1 に相当します。

最初に描画用データを作成します。
先ほど生成した平均売上 $${\mu}$$ のMCMCサンプルのデータ形式を縦持ちに変更します。

# 事後分布の描画用データの作成

# 線形予測子 mu の予測値に宣伝有無・安売り有無を付与して１つのデータフレームにまとめる
mu_pred1_df = (
    # 縦持ちに変換
    pd.melt(
        pd.concat([newdata_1, pd.DataFrame(linear_fit1)], axis=1),
        id_vars=['publicity', 'bargen'], var_name='sample', value_name='sales')
    # インデックスを再発番
    .reset_index(drop=True)
)

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

【実行結果】
パラメータのMCMCサンプル数 4000 × 予測用の説明変数ケース 4 = 16000 個のサンプルです。

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

# p.232 図3.10.1 カテゴリ×カテゴリの交差作用がもたらした売上の変化

# 色の設定
colors = {'not': 'tomato', 'to_implement': 'tab:blue'}

# 描画領域の設定
plt.figure(figsize=(8, 5))

# 観測値のスウォームプロットの描画
sns.swarmplot(data=interaction_1, x='publicity', y='sales', hue='bargen',
              palette=colors)

# 事後分布：平均売上の95%信用区間のエラーバープロットの描画
sns.pointplot(data=mu_pred1_df, x='publicity', y='sales',
              hue='bargen', palette=colors, legend=False, ls='none', 
              errorbar=('pi', 95), capsize=0.1, dodge=0.3, err_kws={'lw': 2})
# 修飾
plt.title('平均売上の95%HDI区間', fontsize=14)
plt.xlabel('宣伝の有無', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.legend(bbox_to_anchor=(1, 1), title='安売りの有無');

【実行結果】
宣伝あり・安売りありのエラーバーの位置が断然高い＝平均売上が大きいことが一目瞭然です。
交互作用の効果を実感できますね！
statsmodels の 95% 信頼区間とも似ている感じです。

（参考：statsmodels の95%信頼区間）

🔵🔵🔵

ケース２：カテゴリ変数 × 量的変数

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

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

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

# ファイルの読み込み
interaction_2 = pd.read_csv('./data/3-10-2-interaction-2.csv')

# 結果の表示
# sales:売上, publicity: 宣伝の有無, temperature: 気温℃
print('interaction_2.shape: ', interaction_2.shape)
interaction_2.head(3)

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

• sales：売上金額 （単位：万円）

• publicity：宣伝の有無（あり：to_implement、なし：not）

• temperature：気温

publicity は二値のカテゴリ変数、temperature は量的変数です。
sales は 量的変数です。

「定数項 + 宣伝の有無 + 気温 + 宣伝の有無 × 気温」を説明変数にした線形回帰モデルで分析します。

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

# 宣伝有無ごとのデータ件数
for col in interaction_2.columns[interaction_2.dtypes=='object']:
    display(interaction_2[col].value_counts().to_frame())

【実行結果】
宣伝なし／ありで、それぞれ 50 件のデータがあります。

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

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

【実行結果】
売上は平均 123.5、最小値 25.9、最大値 263.3 です。
気温は平均 15.5、最小値 0.4、最大値 29.8 です。

次は宣伝の有無別の売上の要約統計量です。

# 宣伝有無別の売上の要約統計量
interaction_2.groupby(['publicity'])['sales'].describe().round(2)

【実行結果】
宣伝あり（to_implement）の売上の方が大きいです！

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

# 宣伝有無別の気温の要約統計量
interaction_2.groupby(['publicity'])['temperature'].describe().round(2)

【実行結果】
大きな違いは見られません。

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

# ジョイントプロットの描画
colors = {'not': 'tomato', 'to_implement': 'tab:blue'}
sns.jointplot(data=interaction_2, x='temperature', y='sales',
              hue='publicity', palette=colors)
# 修飾
plt.xlabel('気温', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.legend(title='宣伝の有無');

【実行結果】
散布図では気温の増加と売上の増加の関係が見られ、しかも、宣伝ありの方が気温上昇につれて売上増分（傾き）が大きくなる傾向が見られます。
気温のKDEプロット（上）では、宣伝の有無でほぼ同じ分布です。
売上のKDEプロット（右）では、宣伝なしのばらつきが小さく、宣伝ありのばらつきが大きいことが分かります。

続いて、回帰係数付きの散布図です。
seaborn の lmplot を利用します。

# 売上と気温・宣伝有無の回帰直線付き散布図

# 回帰直線付き散布図の描画
colors = {'not': 'tomato', 'to_implement': 'tab:blue'}
sns.lmplot(data=interaction_2, x='temperature', y='sales', height=4, aspect=1.8,
           hue='publicity', palette=colors, legend=False,
           scatter_kws={'s': 20})
# 修飾
plt.title('気温・宣伝有無と売上の散布図', fontsize=14)
plt.xlabel('気温', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.xticks(range(0, 31, 2))
plt.legend(title='宣伝の有無');

【実行結果】
宣伝あり・なしの回帰直線は平行になっていなくて、宣伝ありの傾きが明らかに大きくなっています。
「宣伝 × 気温」の交互作用を予感させますね！

🔵🔵🔵

線形回帰モデル by statsmodels

① モデルの概要
statsmodels の 線形回帰 OLS で次のモデルを実装します。

$${\varepsilon}$$ は誤差です。
各ダミー変数は、あり：$${1}$$、なし：$${0}$$ の値を取ります。

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

この formula 文は省略形です。次の formula 文と同等です。

② 線形回帰の実行
では線形回帰を実行します！

# 線形回帰モデル by statsmodels

# formulaの設定
formula = 'sales ~ publicity * temperature'
# 回帰分析の実行
res_sm2 = smf.ols(formula=formula, data=interaction_2).fit()
res_sm2.summary()

【実行結果】
こちらは線形回帰のサマリーです。
自由度調整済み決定係数（Adj. R-squared）は 0.903 です。
かなり当てはまりがよいですね！

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

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

【実行結果】
最下行が交互作用項です。5%水準で有意です（$${p = 0.000}}$）。
宣伝ありの場合、気温が１度上昇すると平均売上が追加的に 4.2 万円大きくなる、と解釈できます。
宣伝なし・気温 0 の平均売上は切片の 43.1 万円と推定されます。

④ 宣伝の有無別・気温別の平均売上の可視化
最後に宣伝の有無別・気温別の平均売上の95%信頼区間を可視化しましょう。

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

# 平均売上の予測値の算出

# 予測用データの作成
n_samples = 100
test_data = pd.DataFrame({
    'publicity': np.repeat(['not', 'to_implement'], n_samples),   # 宣伝の有無
    'temperature': np.tile(np.linspace(0, 30, n_samples), 2)      # 気温
})
# 予測値の取得
prediction2 = res_sm2.get_prediction(test_data)
pred_df2 = prediction2.summary_frame(alpha=0.05) # alpha=0.05で95%信頼区間
# 予測値データフレームの作成
pred_df2 = pd.concat([test_data, pred_df2], axis=1)
pred_df2

【実行結果】
宣伝の有無 2 × 気温 100 刻み に対して、平均売上 mean、95% 信頼区間 mean_ci_lower, mean_ci_upper を取得しました。

それでは描画します。

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

# 色の設定
colors = {'not': 'tomato', 'to_implement': 'tab:blue'}
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=interaction_2, x='temperature', y='sales',
                hue='publicity', palette=colors)
# 宣伝の有無ごとに平均売上の予測値と95%信頼区間の描画を繰り返し処理
for pub, color in colors.items():
    # 描画データの取得
    temp, mean_sales, lower, upper = (
        pred_df2[pred_df2['publicity']==pub].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%信頼区間', fontsize=14)
plt.xlabel('気温', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.xticks(range(0, 31, 2))
plt.legend(title='宣伝の有無');

【実行結果】
宣伝あり・なしの回帰直線は平行になっていなくて、宣伝ありの傾きが明らかに大きくなっています。
「宣伝 × 気温」の交互作用を予感させますね！

ベイズ流の線形回帰モデルへ進みます。

🔵🔵🔵

ベイズモデリング by Stan

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

① モデルの概要
次のモデルを実装します。
ケース１と同じモデルです（説明変数の数も同じです）。

② Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
テキストの Stan コードを引用いたします。
📑ファイル名：3-10-2-glm-intract-2-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 ライブラリを用いてデザイン行列等を作成します。

# デザイン行列の作成

# formula構文の設定
formula_interaction2 = 'sales ~ publicity * temperature'

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

# int型の設定
X = X.astype({'Intercept': int, 'publicity[T.to_implement]': 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['sales'], X=X)

【実行結果】なし

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

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

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

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

【実行結果】

MCMC の準備が整いました！

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

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

【実行結果】

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

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

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

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

# p.234 MCMCの結果の確認
fit2_summary = fit2.summary(percentiles=[2.5, 50, 97.5])
fit2_summary.round(2)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（N_Eff）に問題はなさそうです。
４行目が交互作用項です。宣伝ありの場合、気温が１度上昇すると平均売上が追加的に 4.2 万円大きくなる、と解釈できます。
宣伝なし・気温 0 の平均売上は切片の 43.0 万円と推定されます。

statsmodels の線形回帰の推定値の表を並べましょう。
（参考：statsmodels の推定値）

手法間で係数推定値は近い値になっています。

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

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

【実行結果】

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

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

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

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

⑦ 係数の解釈
テキスト 10.9 節の計算に準拠して進めます。

a. 説明変数を作る
売上（平均）$${\mu}$$ の予測値を算出するために、予測用の説明変数およびデザイン行列を作成します。

# p.235 説明変数を作る

# publicityとtemperatureのデータフレームを作成
newdata_2 = pd.DataFrame({
    'publicity': ['not', 'not', 'to_implement', 'to_implement'],
    'temperature': [0, 10, 0, 10],
})
print('【説明変数】')
display(newdata_2)

# デザイン行列化
newdata_2_dm = dmatrix('publicity * temperature', newdata_2).base
print('【デザイン行列】')
print(newdata_2_dm)

【実行結果】
宣伝の有無と気温 $${\{0, 10\}}$$ の組み合わせです。

b. 係数推定値から平均売上を計算
テキスト 表 3.10.3 に相当します。
予測用のデザイン行列と係数推定値（平均値）の積で平均売上を計算します。

# p.234 表 3.10.3 カテゴリ×数量の交互作用があるときの予測値の変化のパターン
pd.concat(
    [newdata_2,
     pd.Series(newdata_2_dm @ fit2_summary.iloc[[1, 2, 3, 4], 0], name='mean')],
    axis=1
).round(3)

【実行結果】

c. 平均売上 $${\mu}$$ の予測値の要約統計量の計算
まず、予測用のデザイン行列とパラメータの MCMC サンプルの行列積で $${\mu}$$ の予測値の MCMC サンプルを取得します。
次に予測値のMCMCサンプルを arviz の summary 関数に与えて、要約統計量を得ます。

# p.235 売上の予測値の算出

# モデルの係数を取得
b_samples2 = az.extract(idata2.posterior).b.values

# 線形予測子の予測値
linear_fit2 = (newdata_2_dm @ b_samples2)

# 統計量の算出 az.summary()を利用
for i in range(len(linear_fit2)):
    tmp_df = az.summary(linear_fit2[i, :], kind='stats', hdi_prob=0.95)
    tmp_df.index = [i]
    if i==0:
        stats_df2 = tmp_df
    else:
        stats_df2 = pd.concat([stats_df2, tmp_df], axis=0)
# newdataと結合
stats_df2 = pd.concat([newdata_2, stats_df2], axis=1)

# 結果の表示
stats_df2.round(3)

【実行結果】
mean が平均売上 $${\mu}$$ の予測値です。
b. の結果と同じになりました。

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

# p.235 図3.10.2 カテゴリ×数量の交互作用を加えたときの回帰直線

## MCMCサンプルに基づく平均売上の予測値の作成
# MCMCサンプルからβ0, β1, β2, β4を取り出し
intercept, pub_yes, temp, pub_x_temp = az.extract(idata2.posterior).b.values
# x軸の値
x_val = np.linspace(X['temperature'].min(), X['temperature'].max(), 1001)
# MCMCサンプルから平均売上を算出
not_pub_sales = (intercept + pub_yes*0 + np.outer(x_val, temp)).T
imp_pub_sales = (intercept + pub_yes*1 + np.outer(x_val, (temp + pub_x_temp))).T
# 平均売上データと描画色を宣伝有無別の辞書に格納
all_sales = {'not': not_pub_sales, 'to_implement': imp_pub_sales}
colors = {'not': 'tomato', 'to_implement': 'tab:blue'}

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

【実行結果】
宣伝あり・なしの回帰直線は平行になっていなくて、宣伝ありの傾きが明らかに大きくなっています。
交互作用の効果を実感できますね！
statsmodels の 95% 信頼区間とも似ている感じです。

（参考：statsmodels の95%信頼区間）

🔵🔵🔵

ケース３：量的変数 × 量的変数

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

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

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

# ファイルの読み込み
interaction_3 = pd.read_csv('./data/3-10-3-interaction-3.csv')

# 結果の表示
# sales:売上, product: 製品の種類数, clerk: 店員の数
print('interaction_3.shape: ', interaction_3.shape)
interaction_3.head(3)

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

• sales：売上金額 （単位：万円）

• product：製品の種類数（整数）

• clerk：店員数（整数）

「定数項 + 製品の種類数 + 店員数 + 製品の種類数 × 店員数」を説明変数にした線形回帰モデルで分析します。

データの要約統計量を確認します。

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

【実行結果】
売上は平均 204.0、最小値 35.5、最大値 487.2 です。
製品の種類数は平均 29.7、最小値 10、最大値 50 です。
店員数は平均 4.9、最小値 1、最大値 9 です。

後続処理で活用する目的で、店員数の要素のリストを作ります。

# チャートの店員数の並び順を取得
clerk_order = sorted(interaction_3['clerk'].unique())
clerk_order

【実行結果】
店員数は 1 ～ 9 までですね。

データを可視化しましょう。
最初に変化球のヒートマップで、製品の種類数と店員数の組み合わせの平均売上を可視化してみます。
seaborn の heatmap を利用します。

# ヒートマップの描画

# 行：製品の種類数、列：店員数のピボットテーブルを作成
pivot_df = interaction_3.pivot_table(
    index='product', columns='clerk', values='sales'
)
# ヒートマップの描画
sns.heatmap(data=pivot_df, cmap='Reds')
plt.xlabel('店員数', fontsize=12)
plt.ylabel('製品の種類数', fontsize=12)
plt.gca().invert_yaxis();               # 製品の種類数

【実行結果】
赤が濃くなるにつれて平均売上が大きくなります。
店員数が増えるにつれて平均売上が大きくなります。
店員数４以下では製品の種類数にかかわらず、店員数と平均売上はほぼ同じに見えます。
店員数５以上では製品の種類数が大きくなるにつれて、平均売上が大きくなります。

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

# p.237 図3.10.3 売上と製品数・従業員数の散布図

# 描画領域の設定
plt.figure(figsize=(8, 4))
# 散布図
sns.scatterplot(data=interaction_3, x='product', y='sales', hue='clerk',
                palette='tab10')
# 修飾
plt.title('製品の種類数・店員数と売上の散布図')
plt.xlabel('製品の種類数', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.legend(bbox_to_anchor=(1, 1), title='店員数');

【実行結果】
店員数が大きいほど売上が大きい傾向です。
製品の種類数が大きくなると、店員数が大きいほどに売上増分が大きい傾向です。

続いて、回帰係数付きの散布図です。
seaborn の lmplot を利用します。

# 売上と製品の種類数・店員数の散布図（回帰直線付き）

# 散布図
sns.lmplot(data=interaction_3, x='product', y='sales',
           hue='clerk', hue_order=clerk_order, palette='tab10',
           height=4, aspect=1.8, legend=False,
           scatter_kws={'s': 20})
# 修飾
plt.title('製品の種類数・店員数と売上の散布図')
plt.xlabel('製品の種類数', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.xticks(range(10, 51, 10))
plt.legend(bbox_to_anchor=(1, 1), title='店員数');

【実行結果】
店員数別に着色しています。
店員数が多いほど、製品種類数と売上の回帰直線の傾きが大きい印象です。
「製品の種類数 × 店員数」の交互作用を予感させますね！

🔵🔵🔵

線形回帰モデル by statsmodels

① モデルの概要
statsmodels の 線形回帰 OLS で次のモデルを実装します。

$${\varepsilon}$$ は誤差です。

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

この formula 文は省略形です。次の formula 文と同等です。

② 線形回帰の実行
では線形回帰を実行します！

# 線形回帰モデル by statsmodels

# formulaの設定
formula = 'sales ~ product * clerk'
# 回帰分析の実行
res_sm3 = smf.ols(formula=formula, data=interaction_3).fit()
res_sm3.summary()

【実行結果】
こちらは線形回帰のサマリーです。
自由度調整済み決定係数（Adj. R-squared）は 0.967 です。
めちゃくちゃ当てはまりがよいですね！

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

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

【実行結果】
最下行が交互作用項です。5%水準で有意です（$${p = 0.000}}$）。
交互作用項にフォーカスすると、製品種類数 × 店員数 × 1.1 万円 が追加的な売上増加になる、と解釈できます。

④ 製品の種類数別・店員数別の平均売上の可視化
最後に製品の種類数別・店員数別の平均売上の95%信頼区間を可視化しましょう。

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

# 平均売上の予測値の算出

# 予測用データの作成
n_samples = 100
test_data = pd.DataFrame({
    'clerk': np.repeat(clerk_order, n_samples),                          # 店員数
    'product': np.tile(np.linspace(0, 50, n_samples), clerk_order[-1]),  # 種類数
})
# 予測値の取得
prediction3 = res_sm3.get_prediction(test_data)
pred_df3 = prediction3.summary_frame(alpha=0.05) # alpha=0.05で95%信頼区間
# 予測値データフレームの作成
pred_df3 = pd.concat([test_data, pred_df3], axis=1)
pred_df3

【実行結果】
店員数 9 × 製品種類数 100 刻み に対して、平均売上 mean、95% 信頼区間 mean_ci_lower, mean_ci_upper を取得しました。

それでは描画します。

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

# 色の設定
cmap = plt.get_cmap('tab10')

# 描画領域の設定
plt.figure(figsize=(8, 4))

# 観測値の散布図の描画
sns.scatterplot(data=interaction_3, x='product', y='sales',
                hue='clerk', hue_order=clerk_order, palette='tab10')

# 店員数ごとに平均売上の予測値と95%信頼区間の描画を繰り返し処理
for i, clerk in enumerate(clerk_order):
    # 描画データの取得
    product, mean_sales, lower, upper = (
        pred_df3[pred_df3['clerk']==clerk].iloc[:, [1, 2, 4, 5]].to_numpy().T
    )
    # 平均売上の予測値の描画
    plt.plot(product, mean_sales, color=cmap(i))
    # 平均売上の95%信頼区間の塗りつぶし描画
    plt.fill_between(product, lower, upper, color=cmap(i), alpha=0.2)

# 修飾
plt.title('平均売上の95%信頼区間')
plt.xlabel('製品の種類数', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.xticks(range(10, 51, 10))
plt.legend(bbox_to_anchor=(1, 1), title='店員数');

【実行結果】
店員数別の回帰直線は平行になっていなくて、店員数が増えるにつれて傾きが大きくなっています。
「製品の種類数 × 店員数」の交互作用の効果を確認できます。
また、店員数が２以下の場合、回帰直線の傾きが負になっています。
製品の種類数の係数推定値 -2.2742 が効いていますね。

ベイズ流の線形回帰モデルへ進みます。

🔵🔵🔵

ベイズモデリング by Stan

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

① モデルの概要
次のモデルを実装します。
ケース１・ケース２と同じモデルです（説明変数の数も同じです）。

② Stan のモデル設定
Stan ファイル（Stan コード）を作成します。
テキストの Stan コードを引用いたします。
📑ファイル名：3-10-3-glm-intract-2-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 ライブラリを用いてデザイン行列等を作成します。

# デザイン行列の作成

# formula構文の設定
formula_interaction3 = 'sales ~ product * clerk'

# 目的変数Y, デザイン行列（説明変数）Xの作成
Y, X = dmatrices(formula_interaction3, interaction_3, 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-10-3-glm-intract-3-design-matrix.stan'
current_dir = os.path.abspath(os.getcwd())
stan_path = os.path.join(current_dir, 'stan', stan_file)

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

【実行結果】

MCMC の準備が整いました！

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

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

)

【実行結果】

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

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

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

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

# p.238 MCMCの結果の確認
fit3_summary = fit3.summary(percentiles=[2.5, 50, 97.5])
fit3_summary.round(2)

【実行結果】
$${\widehat{R}}$$（r_hat）、有効サンプル数（N_Eff）に問題はなさそうです。
交互作用項にフォーカスすると、製品種類数 × 店員数 × 1.1 万円 が追加的な売上増加になる、と解釈できます。

statsmodels の線形回帰の推定値の表を並べましょう。
（参考：statsmodels の推定値）

手法間で係数推定値は近い値になっています。

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

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

【実行結果】

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

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

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

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

⑦ 係数の解釈
テキスト 10.12 節の計算に準拠して進めます。

a. 説明変数を作る
売上（平均）$${\mu}$$ の予測値を算出するために、予測用の説明変数およびデザイン行列を作成します。

# p.239 説明変数を作る

# productとclerkのデータフレームを作成
newdata_3 = pd.DataFrame({
    'product': [0, 10, 0, 10],
    'clerk': [0, 0, 10, 10],
})
print('【説明変数】')
display(newdata_3)

# デザイン行列化
newdata_3_dm = dmatrix('product * clerk', newdata_3).base
print('【デザイン行列】')
print(newdata_3_dm)

【実行結果】
製品の種類数 $${\{0, 10\}}$$ と店員数 $${\{0, 10\}}$$ の組み合わせです。

b. 係数推定値から平均売上を計算
テキスト 表 3.10.4 に相当します。
予測用のデザイン行列と係数推定値（平均値）の積で平均売上を計算します。

# p.238 表 3.10.4 数量×数量の交互作用があるときの予測値の変化のパターン
pd.concat(
    [newdata_3,
     pd.Series(newdata_3_dm @ fit3_summary.iloc[[1, 2, 3, 4], 0], name='mean')],
    axis=1
).round(3)

【実行結果】

c. 平均売上 $${\mu}$$ の予測値の要約統計量の計算
まず、予測用のデザイン行列とパラメータの MCMC サンプルの行列積で $${\mu}$$ の予測値の MCMC サンプルを取得します。
次に予測値のMCMCサンプルを arviz の summary 関数に与えて、要約統計量を得ます。

# p.239 売上の予測値の算出

# モデルの係数を取得
b_samples3 = az.extract(idata3.posterior).b.values

# 線形予測子の予測値
linear_fit3 = (newdata_3_dm @ b_samples3)

# 統計量の算出 az.summary()を利用
for i in range(len(linear_fit3)):
    tmp_df = az.summary(linear_fit3[i, :], kind='stats', hdi_prob=0.95)
    tmp_df.index = [i]
    if i==0:
        stats_df3 = tmp_df
    else:
        stats_df3 = pd.concat([stats_df3, tmp_df], axis=0)
# newdataと結合
stats_df3 = pd.concat([newdata_3, stats_df3], axis=1)

# 結果の表示
stats_df3.round(3)

【実行結果】
mean が平均売上 $${\mu}$$ の予測値です。
b. の結果とほぼ同じになりました。

⑧ 製品の種類数別・店員数別の平均売上の可視化
テキストの図 3.10.4 に相当します。
パラメータのMCMCサンプルを使って、PyMCの外で平均売上の予測値を計算し、arviz の plot_hdi 関数で 95%HDI 区間を描画します。

# p.240 図3.10.4 数量×数量の交互作用を加えたときの回帰曲線

## MCMCサンプルに基づく平均売上の予測値の作成
# MCMCサンプルからβ0, β1, β2, β4を取り出し
intercept, product, clerk, pro_x_clerk = az.extract(idata3.posterior).b.values
# x軸の値
x_val = np.linspace(X['product'].min(), X['product'].max(), 1001)
# MCMCサンプルから店員数ごとの平均売上を算出
all_sales = []
for i in clerk_order:
    sales = (intercept + np.outer(x_val, product + pro_x_clerk*i) +  clerk*i).T
    all_sales.append(sales)

## 描画
# 色の設定
cmap = plt.get_cmap('tab10')
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 観測値の散布図の描画
sns.scatterplot(data=interaction_3, x='product', y='sales',
                hue='clerk', hue_order=clerk_order, palette='tab10')
# 店員数ごとに平均売上の平均値と95%HDI信用区間の描画を繰り返し処理
for i, sales in enumerate(all_sales):
    # 平均売上の平均値の描画
    plt.plot(x_val, sales.mean(axis=0), color=cmap(i))
    # 平均売上の95%HDI信用区間の塗りつぶし描画
    az.plot_hdi(x_val, sales.reshape(4, 1000, -1), hdi_prob=0.95,
                fill_kwargs={'color': cmap(i), 'alpha': 0.2})
# 修飾
plt.title('平均売上の95%HDI区間')
plt.xlabel('製品の種類数', fontsize=12)
plt.ylabel('売上金額', fontsize=12)
plt.xticks(range(10, 51, 10))
plt.legend(bbox_to_anchor=(1, 1), title='店員数');

【実行結果】
店員数別の回帰直線は平行になっていなくて、店員数が増えるにつれて傾きが大きくなっています。
「製品の種類数 × 店員数」の交互作用の効果を確認できます。
また、店員数が２以下の場合、回帰直線の傾きが負になっています。
製品の種類数の係数推定値 -2.28 が効いていますね。
statsmodels の 95% 信頼区間とも似ている感じです。

（参考：statsmodels の95%信頼区間）

⑨ 店員数ごとに分割した製品の種類数別の平均売上の可視化
テキストの図 3.10.5 に相当します。

# p.241 図3.10.5 数量×数量の交互作用を加えたときの回帰曲線（グラフ分割）

## 描画
# 描画領域の設定
fig, axes = plt.subplots(3, 3, figsize=(10, 8), sharex=True, sharey=True)

# 店員数ごとにチャート描画を繰り返し処理
for clerk, sales, ax in zip(clerk_order, all_sales, axes.flat):

    # 観測値の散布図の描画
    tmp = interaction_3[interaction_3['clerk']==clerk]
    sns.scatterplot(data=tmp, x='product', y='sales', s=30, ax=ax)
    # 平均売上の平均値の描画
    ax.plot(x_val, sales.mean(axis=0), lw=1, color='tomato')
    # 平均売上の95%HDI区間の塗りつぶし描画
    az.plot_hdi(x_val, sales.reshape(4, 1000, -1), hdi_prob=0.95,
                color='tomato', fill_kwargs={'alpha': 0.2}, ax=ax)
    # タイトルの表示、x,y軸ラベルの非表示
    ax.set(title=f'clerk: {clerk} 人', xlabel='', ylabel='')

# チャート全体の修飾
fig.suptitle('店員数ごとの製品の種類数別平均売上')
fig.supxlabel('製品の種類数', fontsize=12)
fig.supylabel('売上金額', fontsize=12)
plt.tight_layout()
plt.show()

【実行結果】
店員数が大きくなるにつれて、回帰直線の傾きが大きくなることがよく分かるチャートです。
傾きの変化は店員数と製品の種類数の交互作用の影響です。

🔵🔵🔵

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

シリーズの記事

次の記事

前の記事

PyMC版

Bambi版

目次

ブログの紹介

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

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