はじめに
- pythonによるベイズ統計モデリング入門ではベイズ統計学の基本的な考え方とpythonによるベイズ線形回帰の例を紹介しました
- 本記事では、解釈性も高く実務上もよく利用される一般化線形モデルについて、引き続きベイズ統計学的な観点から紹介をしてみようと思います
- また、pythonの実装例では、ロジスティック回帰分析をベイズ統計学のアプローチでパラメータ推定を行います
- アウトプットイメージからscikit-leanやstatsmodelsとの違い、ベイズアプローチの特徴を掴んで頂ければと思います
- また、前回記事と同様に分かりやすさ・イメージを中心に紹介しておりますので、厳密な各種理論の定義や解説は書籍をご覧ください
一般化線形モデルとは
線形モデルとは
- 「線形モデル」とは、目的変数yと説明変数Xの関係性が"線形"で表現できるモデルを指します
- 線形モデルの例
- 単回帰分析: $y = ax + ε$
- 重回帰分析 $y = ax^2 + bx + ε$
- "線形"とはざっくり言うと、説明変数が「足し算」の形で表現されたモデルです
- 数学的な定義では、加法性(additivity)や斉次性(homogeneity)といった難しい言葉もあるのですが、ここでは割愛します
- また、逆の表現では"非線形"な関係が含まれていない（"非線形"な関係とはlog、exp、三角関数など）とも言えます
- ここで重要な点は、線形モデルでは誤差εは正規分布と仮定している点です
"一般化"へのモチベーション
- 世の中のデータは誤差に正規分布を仮定できないデータが多く存在します
- 例えば、目的変数yがカテゴリーのデータ( 例: クイズの正解 or 不正解)や正の値しかとり得ないデータ（例: 動物園のペンギンの数）などが該当します
- こういったデータに線形モデルを当てはめることは、平均を中心に-∞から+∞までの推定値が起き得ることを暗に仮定してしまいます
- そこで、誤差を正規分布に仮定できないデータにも線形モデルを適用したい、というモチベーションから"一般化"への研究が進められました
一般化線形モデルの考え方
- 線形モデルの一般化には3つの要素があります
- 一般化線形モデルはこの3要素の組み合わせなので、ここだけ抑えておけば、上述の様々なデータにも線形モデルを適用できる様になります
-
1. 確率分布
- 誤差εが仮定する確率分布
- 例：売上データなら正規分布、コインの裏表なら二項分布
- 確率分布には様々なパターンがある
- 参考：代表的な確率分布の特徴まとめ
-
2. 線形予測子
- 説明変数の線形結合部分
- y = Xβ であれば Xβ の部分
- デザイン行列とも呼ばれる
-
3. リンク関数
- 目的変数yと線形予測子の関係性を表す関数
- 連結関数とも呼ばれる
ロジスティック回帰の例
- 顧客の購買有無の確率を予測する様なモデルとして、ロジスティック回帰の例で考えてみます
- ロジスティック回帰では、以下の様な前提を置きます
- 確率分布 = 二項分布（目的変数が2つの値しかとらない為）
- 線形予測子 = 観測されたデータx1とx2の線形結合
- リンク関数 = ロジット関数（確率値は0~1なので、その範囲の出力にしたい）
- この様に確率分布、線形予測子、リンク関数の組み合わせで様々なデータを線形モデルで表現することができます
- このフレームワークで考えると、回帰分析（例：y=ax+b）は確率分布は正規分布、線形予測子をax+b、リンク関数を恒等関数（何も変換しない関数）と捉えることができます
確率分布とリンク関数の組み合わせ
- 非常に沢山の確率分布があるので、3つの組み合わせだと自由度が高すぎて、どの様にモデリングに手をつけたらいいのか悩むかもしれません
- しかし、実際上は確率分布ごとによく使われるリンク関数は決まっています
| データタイプ | 確率分布 | リンク関数 |
|---|---|---|
| 離散値 | ベルヌーイ分布 | Logit |
| 二項分布 | Logit | |
| 連続値 | ポアソン分布 | log |
| ガンマ分布 | log | |
| 正規分布 | - |
※データ解析のための統計モデリングより作成
一般化線形モデルの注意点
- 一般化線形モデルは”一般化”により対応できるデータ範囲を拡張させたとはいえ、線形モデルをベースにしていることに変わりません
- そのため、線形モデルにおける注意事項（多重共線性など）は同様に当てはまります
- 機械的にモデリングをする前に、データの基本統計量や変数間の相関の確認は必須です
- 以前、新規データを受領したら最初にすべき10ステップとい記事を作成しましたので、この辺りの方法についても興味がある方はご参照ください
Pythonによる一般化線形モデル
- モデル構築には、確率モデル、線形予測子、リンク関数を設定した後、手元のデータからパラメータを推定する必要があります
- パラメータの推定方法としては、以下の2パターンがあります
- パラメータの値を直接推定する(statsmodels, sklean)
- パラメータの確率分布を推定する(pystan, pymc3)
- 1の方法はQiitaを含めて多くの記事があるので、ここでは記事の紹介だけに留め、本記事では2の方法について説明します
1. statsmodelsやskleanによるモデル構築
- Pythonの場合、統計モデリングのアプローチではstatsmodels、機械学習的なアプローチではscikit-leanをを利用することが多いです
- statsmodelsでは、pandas.dataframeに対して、確率分布、線形予測子、リンク関数を指定すれば簡単にモデルの構築が可能です
- scikit-leanでは配列として準備した特徴量Xと目的変数yを用意し、fit（パラメータ学習）させるだけでモデルの学習が完了します
2. pystanによるモデル構築（ベイジアンアプローチ）
- ベイズ統計学のアプローチでは、パラメータの確率分布を推定することで、予測結果の確率分布を取得することができます（パラメータの点推定でなく、パラメーターの分布推定）
- これは、予測の確信度をセットで入手することができるので、実務上の意思決定に役立つことに繋がります
- 本アプローチはベイズ統計学の考え方を前提にしていますが、ここではpythonによる実装部分のみ紹介しよう思います
- また、本記事ではpyMCを利用してロジスティック回帰を解いています
デモ: pyMCによるベイズロジスティック回帰
- ここではirisのデータセット（2クラス分類へデータを修正）を利用して、ベイズロジスティック回帰を試します
- pyMCの使い方は前回記事の方が詳しいので、詳細が気になる方はご参照ください
- はじめに必要なライブラリをインポートします
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns
import pymc3 as pm
サンプルデータ作成
- seabornのデータセットからirisデータをロードします
- 今回は2クラス分類問題へロジスティック回帰分析を行いたいので、正解ラベル（species）を修正します
- 散布図の可視化をしていますが、グラフをみるだけでも綺麗に分離できそうです
# サンプルデータの作成（irisのデータを2クラスに変更）
iris_df = sns.load_dataset('iris') # データセットの読み込み
# 2クラス分類の問題へ変更
label = (iris_df['species']=='versicolor') | (iris_df['species']=='virginica')
iris_df['species'] = label * 1
# データの確認: species(0 or 1)が目的変数y、それ以外を特徴量Xとする
iris_df.head()
# データの確認: 分布と散布図の可視化
sns.pairplot(iris_df, hue='species')
モデル構築
- はじめにpd.DataFramde型だとデータが取り扱えないので、配列形式(np.array)に変換します
- その際に、モデル設定に必要なパラメータを取得しておきます
# y, Xを配列で取得
y = iris_df['species'].values
X = iris_df.drop('species', axis=1).values
# 回帰モデルの定数項を表現するために1の列を追加（線形基底関数モデル）
X = np.hstack((np.ones((X.shape[0], 1)), X))
# モデル設定に必要なパラメータを取得
n, k = X.shape # 説明変数Xの行列サイズ
b0 = np.zeros(k) # 初期値0のβ
A0 = 0.01 * np.eye(k) # 分散共分散行列の逆行列（精度行列）
- つぎにpyMCを用いてモデルの定義を行います
- ここでは回帰係数$β_i$の事前分布に正規分布を仮定しており、回帰係数は複数あるので多変量正規分布を用いています
- 数式で記述すると以下の様になります
$logit(p) = β_1sepal_length + β_2sepal_width + β_3petal_length + β_4petal_width + Constant $
# モデルの設定
logit_model = pm.Model()
with logit_model as model:
# 回帰係数βの事前分布として多変量正規分布を設定
b = pm.MvNormal('b', mu=b0, tau=A0, shape=k)
# 線形予測子を内積（説明変数Xと係数β）で表現
y_hat = pm.math.dot(X, b)
# 尤度の算出
## ロジスティック回帰なので、ベルヌーイ分布を仮定
## リンク関数としてロジット（logit_p）を指定
likelihood = pm.Bernoulli('y', logit_p=y_hat, observed=y)
パラメータの推定（事後分布のサンプリング）
- モデルの準備だできたので、 MCMCを利用したサンプリングを行います
- MCMCのよるサンプリングの引数の説明は前回記事をご参照ください
# 事後分布からのサンプリング
n_draws = 5000
n_chains = 4
n_tune = 1000
with logit_model:
trace = pm.sample(draws=n_draws, chains=n_chains, tune=n_tune, random_seed=123)
# 評価に利用するサンプルの取得
trace_n = trace[1000:]# "バーン・イン"として初めの1000件を捨てる
推移結果を利用した予測
- ここまででロジスティック回帰分析の回帰係数のパラメータの分布の推定が完了しました
- この推定結果を用いて、上手く分類できるモデルが構築できているか、正解データと比べてみます
# 予測分布の取得
ppc = pm.sample_posterior_predictive(trace_n, samples=1000, model=logit_model)
# 予測分布の平均値
y_pred = ppc['y'].mean(axis=0)
# 0.5を閾値にバイナリラベルに変換
y_pred_binary = np.where(y_pred>0.5, 1, 0)
# 正解データyとの比較
print("Accuracy: ", np.sum(y == y_pred_binary)/len(y))
Accuracy: 1.0
- 綺麗に分離するモデルが構築されている様です
回帰係数の評価
- ここからがベイズアプローチの特徴が活かされるパートになります
- 逆に言えば、予測結果だけ欲しい場合は、skleanの方がコードも簡易で高速です
- 推定した回帰係数の分布を可視化してみます
# 推定したロジスティック回帰係数の分布の可視化
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
sns.histplot(trace_n['b'][:,0], ax=axes[0,0], kde=True, alpha=.7).set_title('Constant')
sns.histplot(trace_n['b'][:,1], ax=axes[0,1], kde=True, alpha=.7).set_title('b_1 (sepal_length)')
sns.histplot(trace_n['b'][:,2], ax=axes[0,2], kde=True, alpha=.7).set_title('b_2 (sepal_width)')
sns.histplot(trace_n['b'][:,3], ax=axes[1,0], kde=True, alpha=.7).set_title('b_3 (petal_length)')
sns.histplot(trace_n['b'][:,4], ax=axes[1,1], kde=True, alpha=.7).set_title('b_4 (petal_width)')
- このように回帰係数ごとの分布を取得することができます
- ロジスティック回帰を用いる際は回帰係数から各変数ごとの影響度を考察するかと思います
- この様に回帰係数の分布も可視化することで、回帰係数ごとの確信度（分布のばらつき）や歪みを視覚的に捉えることができる様になります
さいごに
- 本記事では一般化線形モデルの概要とpythonによるベイズ統計モデリングを用いたロジスティック回帰分析の紹介を行いました
- 今回はロジスティック回帰分析のみの紹介に留めましたが、同様の方法でポアソン回帰や時系列モデルへの拡張として状態空間モデルも取り扱うことができます
- 昨今はニューラルネットワークや機械学習を用いた予測精度重視の話題が多いのですが、データそのものの深い理解や意思決定に役立つ統計モデリングという点でベイズ統計学は非常に有用かと思います
- 本記事はベイズモデリングの導入の導入...程度の位置付けなので、興味を持った方は専門書でより深く学習をお勧めします