私の座右の銘は「限界突破」でした。それは、自分が想像できる範囲の自分で決めた限界を乗り越えることに喜びを覚えていたわけです。まるで、判別境界を少しでも超えたらその瞬間に異なるクラスに振り分けられるように、設定した”限界”を超えた瞬間に突破と考えていたようです。
前置きが長くなりましたが、今回は”限界”という名が付くものの、全くそのニューアンスのない限界効果（Maginal Effect）について書いきます。限界効果とは、経済学でよく使われており、財を1単位追加して消費することによる効用の増加分のとことです（限界効用（Wikipedia））。線形回帰の話でよく聞くかと思います。限界効果がわかると、各変数が1単位増加したときに被説名変数がどれだけ増減するかがわかり施策を行いやすいです。
線形回帰の限界効果に関しては、巷に溢れているので今回は扱いません。今回はロジスティック回帰の限界効果に関して数式の導出とpythonを使っての実践の両面で書いていきます。線形回帰では、偏回帰係数が各変数の限界効果ですが、ロジスティック回帰では同じようにいきません。それは、ロジスティック回帰が非線形モデルだからです。ロジスティック回帰で限界効果を求めるためには、何かしらの処理が新たに必要です。今回のコードはこちらにあります。
ロジスティック回帰の限界効果の数式
ロジスティック回帰の限界効果は、モデルの偏回帰係数ではないです。なので、新たに計算をしなければいけません。そこで、限界効果が何だったかを思い出すと、変数が1単位増加したときに被説明変数の増減量のことでした。簡単な単回帰モデルで考えると、変数の傾きが限界効果に当たります。単回帰モデルを \(y(x)=ax+b\) とすると、下記の式変形で得られます。
y(x+1)-y(x) &= a(x+1)+b - (ax+b) \\
&= a
\end{align*}
これをロジスティック回帰で行いたいのですが、ロジスティック回帰は非線形モデルなので、傾きは各変数の値で変わってしまい、線形回帰モデルのように簡単には求まりません。ではどうするかというと、変数の値で傾きは変わってしまうのならば、傾きを平均すればええやんてな感じになります。
関数の傾きは微分により求めることができます。そこで、ここからは数式によりロジスティック回帰の限界効果（平均限界効果）を求めていきます。データ \(i\) 番目の \(j\) 番目の変数を \(x_{ij}\) とし、被説名変数 \(Y_i\) が1となる確率を \(\mathrm Pr (Y_i=1)\) とすると、限界確率効果（Marginal Probability Effect）\(MPE_{ij}\) は次式となります。
ここでロジスティック回帰モデルにおいて、\( \mathrm Pr (Y_i = 1) \) は、 \( \mathrm Pr (Y_i = 1) = f( \mathrm X_i \mathrm{ \beta } ) = \frac {e^x} {1+e^x} \) となります。
変数の数を \(k\) とすると、\( \mathrm X_i = [x_{i1}, x_{i2}, ..., x_{ik}] \) であり、\( \mathrm X_i \) はデータ \( i \) の各変数の値のベクトルです。 \( \mathrm{\beta} \) はロジスティック回帰モデルの偏回帰係数のベクトルであり、 \( \mathrm{\beta = [\beta_1, \beta_2, ..., \beta_j ]} \) です。これらをもとに \( MPE_{ij} \) は次式となります。
このときに注目されるのは、分子の \( X_i \) はベクトル変数であり、分母の \( x_{ij} \) はスカラー変数です。よって、このままでは偏微分できず、式変形を必要とします。
上式の右側の \( \frac { \partial X_i \beta} { \partial x_{ij} } \) は \( \beta_j\) となります。
\( \frac { \partial f ( \mathrm X_i \mathrm{ \beta })} { \partial \mathrm X_i \mathrm{ \beta }} \) はロジスティック回帰の偏微分です。
f'(x) &= { e^x }' \cdot (1 + e^x )^{-1} + e^x \cdot ( (1 + e^x )^{-1} )' \\
&= e^x \cdot (1 + e^x )^{-1} + e^x \cdot ( -(1 + e^x )^{-2} )e^x \\
&= \frac { e^x } { 1 + e^x } - \frac { e^{2x} } {(1+e^x)^2} \\
&= \frac { e^x } { (1 + e^x)^2 } \\
&= \frac {e^x} {1+e^x} \cdot (1 - \frac { e^x } { 1+e^x }) \\
&= f(x) \cdot ( 1 - f(x) )
\end{align*}
上式を見てもらえば分かる通り、ロジスティック回帰の偏微分の結果は、ロジスティック回帰 \(f(x)\) で表すことができ、それはつまり、ロジスティック回帰モデルの推測値だけを使って計算できるというとても便利なわけです。これまでの導出結果を踏まえて、限界確率効果 \(MPE_{ij}\) は以下で求められます。
今回の本題である限界効果は、限界確率効果の平均をとることで求まります。これは、変数の値で傾きが変わるのだから、平均をとってしまえばよいということです。この方法により求まる限界効果は、平均限界効果（Average Marginal Probability Effect）と呼ばれます。
ロジスティック回帰の限界効果（平均限界効果）を求める数式がわかりました。結果を見ると、大層おおがかりなものではなく、ロジスティック回帰の予測値と係数だけを使う非常にシンプルなものとなりました。せっかくなので、実際にデータを使って求めてみます。
実際にロジスティック回帰の限界効果を求める
これまでの話で、ロジスティック回帰における限界効果（平均限界効果）を求めるための数式がわかりました。なので、その数式通り愚直に計算をすればよいわけです。それでは、実践します。
実践のための準備
必要なモジュールをインポートし、結果を保存するためのディレクトリを作成しています。今回は、scikit-learnのロジスティック回帰を使います。
from pathlib import Path import pandas as pd import numpy as np import seaborn as sns import matplotlib.pyplot as plt from sklearn.datasets import load_breast_cancer from sklearn.model_selection import train_test_split from sklearn.linear_model import LogisticRegression import scipy.stats as ss result_dir_path = Path('result') if not result_dir_path.exists(): result_dir_path.mkdir(parents=True)
使用データの準備
今回使用するデータは乳がんの診断データです。分類目的で使用されるデータであり、各ケースにおいて良性腫瘍か悪性腫瘍かの診断結果がついています。レコード数は569、カラム数は30です。
cancer = load_breast_cancer() cancer_data = pd.DataFrame(cancer.data, columns=cancer.feature_names) cancer_data['y'] = cancer.target X_data = cancer_data.drop(columns='y') y_data = cancer_data['y'].tolist() X_train, X_test, y_train, y_test = train_test_split( X_data, y_data, test_size=0.3, random_state=42 )
ロジスティック回帰モデルの構築
乳がん診断データを分類するロジスティック回帰モデルを構築します。ここらへんは流れ作業ですね。
model = LogisticRegression(max_iter=10000) model.fit(X_train, y_train) print('正答率:{:.0f}%'.format( model.score(X_test, y_test)* 100 ))
正答率は下記の図の通り、テストデータにおいて98%です。まぁまぁですね。
モデルの係数の確認
上記でモデルを構築することができたので、すかさず限界効果を求めようとなりますが、各係数の状況を先に把握しておきます。
モデルの係数の信頼性を確認するためにp値を見たいのですが、scikit-learnのロジスティック回帰ではp値を求める関数はないです。そこで、下記の方のコードを参考にさせてもらいました。
P values for sklearn logistic regression · GitHub
def calc_logistic_p_value(model, X, y, alpha=0.05): denom = (2.0 * (1.0 + np.cosh(model.decision_function(X)))) denom = np.tile(denom,(X.shape[1],1)).T F_ij = np.dot((X/denom).T,X) ## Fisher Information Matrix Cramer_Rao = np.linalg.inv(F_ij) ## Inverse Information Matrix sigma_estimates = np.sqrt(np.diagonal(Cramer_Rao)) coefs = model.coef_[0] n = len(X) z_scores = coefs/sigma_estimates # z-score for eaach model coefficient p_values = [ss.norm.sf(abs(x))*2 for x in z_scores] ### two tailed test for p-values z_ = ss.norm.ppf(1-alpha/2) uppers = coefs + z_ * sigma_estimates lowers = coefs - z_ * sigma_estimates return pd.DataFrame({ 'colname' : X.columns.tolist(), 'coef' : coefs, 'sigma' : sigma_estimates, 'z_score' : z_scores, 'p_values' : p_values, 'lower_{:.3f}'.format(alpha/2) : lowers, 'upper_{:.3f}'.format(1-alpha/2) : uppers, 'significant' : ['*' if x <= alpha else '' for x in p_values] }) result = calc_logistic_p_value( model=model, X=X_train, y=y_train, alpha=0.05 ) result.to_csv(result_dir_path.joinpath('logit_summary.csv'), index=False)
下記の図は、学習済みのロジスティック回帰モデルの係数とp値です。p値が有意水準5%を下回った変数は「worst texture」だけです。「worst texture」の値が乳がんの診断に影響していることがわかります。「worst texture」の係数はマイナスになっており、値が小さくなると悪性腫瘍と診断される確率が上がるようです。どの変数が影響するかはわかりましたが、モデルの係数は限界効果ではないため、悪性腫瘍への診断への効果がいまいちつかめません。それでは、今回のメインテーマであるロジスティック回帰の平均限界効果を実際に求めてみます。
ロジスティック回帰の限界効果（平均限界効果）
ここまででロジスティック回帰のモデルを構築し、どの変数が被説明変数に影響しているかをモデルの係数により把握できましたが、どれだけ被説明変数に影響するかはわかりませんでした。なので、上記で導出した数式を使い、平均限界効果を求めます。実際のコードは大したことをしておらず、シンプルに書けます。
def calc_marginal_probability_effect(model, X, alpha=0.05): predict_values = model.predict_proba(X)[:, 1] predict_values = predict_values * (1 - predict_values) mpe = np.tile(predict_values.reshape((-1, 1)), X.shape[1]) * model.coef_[0] n = len(X) ampes = mpe.mean(axis=0) stds = mpe.std(axis=0) t_values = ampes / (stds/np.sqrt(n-1)) p_values = [1 - ss.t.cdf(x, df=n-1) if x > 0 else ss.t.cdf(x, df=n-1) for x in t_values] t_ = ss.t.ppf(1-alpha/2, df=n-1) uppers = ampes + t_ * (stds/np.sqrt(n-1)) lowers = ampes - t_ * (stds/np.sqrt(n-1)) return pd.DataFrame({ 'colname' : X_data.columns.tolist(), 'coef' : model.coef_[0], 'ampe' : ampes, 'mpe_std' : stds, 't_value' : t_values, 'p_value' : p_values, 'lower_{:.3f}'.format(alpha/2) : lowers, 'upper_{:.3f}'.format(1-alpha/2) : uppers, 'significant' : ['*' if x <= alpha else '' for x in p_values] }) result = calc_marginal_probability_effect( model=model, X=X_train, alpha=0.05 ) result = result.sort_values('ampe') result.to_csv(result_dir_path.joinpath('logit_ampe.csv'), index=False)
下記の図では、平均限界効果の計算結果を示しています。先程、ロジスティック回帰モデルの係数で優位なのは「worst texture」であることがわかりました。「worst texture」の平均限界効果を見てみると、-1.9%です。つまり、「worst texture」がわずかに上がると悪性腫瘍と診断される確率が1.9%減少することがわかりました。
グラフにより確認
実際に「worst texture」の値が被説明変数に影響がありそうかを、グラフにより確認してみます。確かに「worst texture」は、良性と悪性で差があるように見えます。ただ他の変数も同じような傾向がある気もします。こちらは今度詳しく見てみたいです。
終わりに
最近は因果推定に興味があり、傾向スコアでよく使うロジスティック回帰に関して、勉強がてら今回のテーマを書いてみようと思いました。因果推定を勉強しているので計量経済学の本をいくつか読んでいます。その中でも下記の本は、実際の分析結果の具体例があり理解しやすかったです。
