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

書籍の著者　馬場真哉 先生

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

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

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

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

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

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

はじめに

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

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

準備

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

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

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

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

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

# インポート

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

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

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

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

第10章 交互作用 by Bambi

交互作用とは？

テキスト 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 Bambi

テキスト 10.5、10.6、10.7 節に相当します。
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, 'publicity': uninformed_prior,
          'bargen': uninformed_prior, 'publicity:bargen':uninformed_prior,
          'sigma': sigma_prior}

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

【実行結果】なし

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

# モデルの表示
model_bmb1

【実行結果】
交互作用項「publicity:bargen」も他のパラメータと同じように無情報的な事前分布を設定しています。

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

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

【実行結果】
右上に交互作用項「publicity:bargen」が表示されています。

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

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

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

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

# 要約統計量の表示
var_names = ['Intercept', 'publicity', 'bargen', 'publicity:bargen', 'sigma']
fit1_summary = az.summary(idata_bmb1, var_names=var_names, hdi_prob=0.95)
fit1_summary

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

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

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

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

# トレースプロットの描画
az.plot_trace(idata_bmb1, compact=False, var_names=var_names,
              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)

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

b. 係数推定値から平均売上を計算
テキスト 表 3.10.2 に相当します。
係数推定値に対して、カテゴリ変数「宣伝あり／なし、安売りあり／なし」の値に応じて「足し算」して、平均売上を計算します。

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

【実行結果】

c. 平均売上 $${\mu}$$ の予測値の要約統計量の計算
まず、model_bmb に対して predict メソッドを適用して $${\mu}$$ の予測値の MCMC サンプルを取得します。
次に予測値のMCMCサンプルを arviz の summary 関数に与えて、要約統計量を得ます。

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

# MCMCサンプルに基づく平均売上 mu の予測値の算出
mu_pred1_idata = model_bmb1.predict(
    idata_bmb1, kind='response_params', data=newdata_1, inplace=False
)

# 要約統計量の算出 az.summary()を利用
stats_df1 = az.summary(
    mu_pred1_idata, var_names=['mu'], kind='stats', hdi_prob=0.95
).reset_index(drop=True)

# newdataと要約統計量を結合
stats_df1 = pd.concat([newdata_1, stats_df1], axis=1)

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

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

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

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

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

# muの予測値のnumpy配列化 shape=(4, 4000)
mu_pred1_np = az.extract(mu_pred1_idata).mu.to_numpy()
# 説明変数とmuの予測値配列の結合 shape=(16000, 3)
mu_pred1_tmp = np.column_stack(
    [np.tile(newdata_1, (4000, 1)), mu_pred1_np.T.flatten()]
)
# pandasデータフレーム化
mu_pred1_df = pd.DataFrame(
    mu_pred1_tmp, columns=['publicity', 'bargen', 'sales']
)
# 結果の表示
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 Bambi

テキスト 10.8、10.9、10.10 節に相当します。
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,
          'publicity': uninformed_prior,
          'temperature': uninformed_prior,
          'publicity:temperature':uninformed_prior,
          'sigma': sigma_prior}

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

【実行結果】なし

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

# モデルの表示
model_bmb2

【実行結果】
交互作用項「publicity:temperature」も他のパラメータと同じように無情報的な事前分布を設定しています。

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

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

【実行結果】
右上に交互作用項「publicity:temperature」が表示されています。

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

%%time
# MCMCの実行 ※nutpieのvar_namesワーニングは無視する...
idata_bmb2 = model_bmb2.fit(
    draws=1000, tune=1000, chains=4, random_seed=123, nuts_sampler='nutpie'
)

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

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

# 要約統計量の表示
var_names = [
    'Intercept', 'publicity', 'temperature', 'publicity:temperature', 'sigma'
]
fit2_summary = az.summary(idata_bmb2, var_names=var_names, hdi_prob=0.95)
fit2_summary

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

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

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

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

# トレースプロットの描画
az.plot_trace(idata_bmb2, compact=False, var_names=var_names,
              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)

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

b. 係数推定値から平均売上を計算
テキスト 表 3.10.3 に相当します。
係数推定値に対して、カテゴリ変数「宣伝あり／なし、安売りあり／なし」の値に応じて「足し算」して、平均売上を計算します。

# p.234 表 3.10.3 カテゴリ×数量の交互作用があるときの予測値の変化のパターン
def get_mean(pub, temp):
    multiply = np.array([1, pub, temp, pub*temp])
    return sum(fit2_summary.iloc[:4, 0] * multiply)
pd.concat([
    newdata_2,
    pd.DataFrame(
        [get_mean(*param) for param in [(0, 0), (0, 10), (1, 0), (1, 10)]],
    columns=['mean'])
], axis=1)

【実行結果】

c. 平均売上 $${\mu}$$ の予測値の要約統計量の計算
まず、model_bmb に対して predict メソッドを適用して $${\mu}$$ の予測値の MCMC サンプルを取得します。
次に予測値のMCMCサンプルを arviz の summary 関数に与えて、要約統計量を得ます。

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

# MCMCサンプルに基づく平均売上 mu の予測値の算出
mu_pred2_idata = model_bmb2.predict(
    idata_bmb2, kind='response_params', data=newdata_2, inplace=False
)

# 要約統計量の算出 az.summary()を利用
stats_df2 = az.summary(
    mu_pred2_idata, var_names=['mu'], kind='stats', hdi_prob=0.95
).reset_index(drop=True)

# newdataと要約統計量を結合
stats_df2 = pd.concat([newdata_2, stats_df2], axis=1)

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

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

⑥ 宣伝の有無別・気温別の平均売上の可視化
テキストの図 3.10.2 に相当します。

最初に平均売上の予測値を算出します。
Bambi の interpret.predictions 関数を使って、予測値の平均値と 95% HDI 区間を取得します。
予測用の説明変数データを作らなくても取得可能です。

# 平均売上の予測値の取得 ※Bambiの予測計算関数を利用

pred_df2b = bmb.interpret.predictions(
    model_bmb2,                               # Bambiモデル
    idata_bmb2,                               # idata
    conditional=['temperature', 'publicity'], # 条件付けする共変量（変化させる変数）
    target='mean',                            # muの事後分布
    use_hdi=True,                             # True: HDI、False: 分位数（信用区間）
    prob=0.95,                                # 区間の確率
)
pred_df2b

【実行結果】
気温と宣伝の有無の組み合わせ $${50×2=100}$$ に対する平均売上の平均値と 95% HDI 区間のデータを取得しました。

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

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

# 色の設定
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, color in colors.items():
    # 描画データの取得
    temp, mean_sales, lower, upper = (
        pred_df2b[pred_df2b['publicity']==pub].iloc[:, [0, 2, 3, 4]].to_numpy().T
    )
    # 平均売上の平均値の描画
    plt.plot(temp, mean_sales, color=color)
    # 平均売上の95%HDI区間の塗りつぶし描画
    plt.fill_between(temp, lower, upper, color=color, 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', hue_order=clerk_order, 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 Bambi

テキスト 10.11、10.12、10.13 節に相当します。
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, 'product': uninformed_prior,
          'clerk': uninformed_prior, 'product:clerk':uninformed_prior,
          'sigma': sigma_prior}

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

【実行結果】なし

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

# モデルの表示
model_bmb3

【実行結果】
交互作用項「product:clerk」も他のパラメータと同じように無情報的な事前分布を設定しています。

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

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

【実行結果】
右上に交互作用項「product:clerk」が表示されています。

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

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

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

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

# 要約統計量の表示
var_names = ['Intercept', 'product', 'clerk', 'product:clerk', 'sigma']
fit3_summary = az.summary(idata_bmb3, var_names=var_names, hdi_prob=0.95)
fit3_summary

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

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

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

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

# トレースプロットの描画
az.plot_trace(idata_bmb3, compact=False, var_names=var_names,
              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)

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

b. 係数推定値から平均売上を計算
テキスト 表 3.10.4 に相当します。
係数推定値に対して、カテゴリ変数「宣伝あり／なし、安売りあり／なし」の値に応じて「足し算」して、平均売上を計算します。

# p.238 表 3.10.4 数量×数量の交互作用があるときの予測値の変化のパターン
def get_mean(product, clerk):
    multiply = np.array([1, product, clerk, product*clerk])
    return sum(fit3_summary.iloc[:4, 0] * multiply)
pd.concat([
    newdata_3,
    pd.DataFrame(
        [get_mean(*param) for param in [(0, 0), (10, 0), (0, 10), (10, 10)]],
    columns=['mean'])
], axis=1)

【実行結果】

c. 平均売上 $${\mu}$$ の予測値の要約統計量の計算
まず、model_bmb に対して predict メソッドを適用して $${\mu}$$ の予測値の MCMC サンプルを取得します。
次に予測値のMCMCサンプルを arviz の summary 関数に与えて、要約統計量を得ます。

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

# MCMCサンプルに基づく平均売上 mu の予測値の算出
mu_pred3_idata = model_bmb3.predict(
    idata_bmb3, kind='response_params', data=newdata_3, inplace=False
)

# 要約統計量の算出 az.summary()を利用
stats_df3 = az.summary(
    mu_pred3_idata, var_names=['mu'], kind='stats', hdi_prob=0.95
).reset_index(drop=True)

# newdataと要約統計量を結合
stats_df3 = pd.concat([newdata_3, stats_df3], axis=1)

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

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

⑥ 製品の種類数別・店員数別の平均売上の可視化
テキストの図 3.10.4 に相当します。

最初に平均売上の予測値を算出します。
Bambi の interpret.predictions 関数を使って、予測値の平均値と 95% HDI 区間を取得します。
予測用の説明変数データを作らなくても取得可能です。

# 平均売上の予測値の取得 ※Bambiの予測計算関数を利用

pred_df3b = bmb.interpret.predictions(
    model_bmb3,                           # Bambiモデル
    idata_bmb3,                           # idata
    conditional=['product', 'clerk'],     # 条件付けする共変量（変化させる変数）
    target='mean',                        # muの事後分布
    use_hdi=True,                         # True: HDI、False: 分位数（信用区間）
    prob=0.95,                            # 区間の確率
)
pred_df3b

【実行結果】
製品の種類数と店員数の組み合わせ $${41×9=100}$$ に対する平均売上の平均値と 95% HDI 区間のデータを取得しました。

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

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

# 色の設定
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, clerk in enumerate(clerk_order):
    # 描画データの取得
    product, mean_sales, lower, upper = (
        pred_df3b[pred_df3b['clerk']==clerk].iloc[:, [0, 2, 3, 4]].to_numpy().T
    )
    # 平均売上の平均値の描画
    plt.plot(product, mean_sales, color=cmap(i))
    # 平均売上の95%HDI区間の塗りつぶし描画
    plt.fill_between(product, lower, upper, 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.274 が効いていますね。
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, ax in zip(clerk_order, axes.flat):

    # 観測値の散布図の描画
    clerk_df = interaction_3[interaction_3['clerk']==clerk]
    sns.scatterplot(data=clerk_df, x='product', y='sales', s=30, ax=ax)
    # 平均売上の平均値の描画
    clerk_df2 = pred_df3b[pred_df3b['clerk'] == clerk]
    ax.plot(clerk_df2['product'], clerk_df2['estimate'], lw=1, color='tomato')
    # 平均売上の95%HDI区間の塗りつぶし描画
    ax.fill_between(
        clerk_df2['product'], clerk_df2['lower_3.0%'], clerk_df2['upper_97.0%'],
        color='tomato', alpha=0.2
    )
    # タイトルの表示、x,y軸ラベルの非表示
    ax.set(title=f'店員数: {clerk}', xlabel='', ylabel='')

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

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

🦌🦌🦌

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

シリーズの記事

次の記事

前の記事

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

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