【Python】ロジスティック回帰で限界効果を求めたい

はじめに

2値のバイナリ変数を被説明変数とした回帰分析を行う場合、ロジスティック回帰を行うことが多いと思います。

ロジスティック回帰で推定された係数$\beta$の解釈を行う際には、ネイピア数を推定された係数乗する($e^{\beta}$)ことでオッズ比を求めることができます。
例えば、とある試験の合格率に効く変数を調べるために、試験前の勉強時間(時間)を説明変数としたロジスティック回帰を実行し、係数が0.32と推定された場合を考えます。
この場合 $e^{0.32} =1.37$ となり、「勉強時間が1時間増えると合格するオッズが1.37倍になる」というような解釈になります。これは少し解釈がしにくいのかなと個人的には思っています。
線形回帰のように「勉強時間が1時間増えると合格率がY%増える」というような解釈ができる方が、直感的で分かりやすいですし、ビジネスでも活用しやすいのではないでしょうか。

限界効果

上述のような説明変数が1単位変化した際の被説明変数の変化量は、限界効果(marginal effect) と呼ばれます。
$y = \beta_0 + \beta_1x + \epsilon$ というような一般的な線形回帰モデルであれば、推定された係数の値($\beta_1$)が限界効果となります。
ではロジスティック回帰のような非線形モデルの場合はどうなるでしょうか？

ロジスティック回帰の限界効果

ロジスティック回帰における係数の限界効果を求めてみます。以下のようなモデルを仮定します。

$$P(y_i = 1 | x_i) = \frac{1}{1 + e^{-\beta_0 + \beta_1x_i + \epsilon_i}}$$

先ほどの試験の合格率の話で言えば、$i$が生徒を表すインデックス、$y_i$が生徒$i$が合格する確率、$x_i$が生徒$i$の勉強時間となります。

この式を $x_i$ で偏微分して限界効果を算出します。

$$\begin{aligned} \frac{\partial}{\partial x_i} P(y_i = 1 | x_i) &= \frac{\partial}{\partial x_i} (\frac{1}{1 + e^{-\beta_0 + \beta_1x_i + \epsilon_i}}) \\\\ &= \frac{-(1 + e^{-\beta_0 + \beta_1x_i + \epsilon_i})'}{(1 + e^{-\beta_0 + \beta_1x_i + \epsilon_i})^2} \\\\ &= \frac{e^{-\beta_0 + \beta_1x_i + \epsilon_i}}{(1 + e^{-\beta_0 + \beta_1x_i + \epsilon_i})^2} \cdot \beta_1 \\\\ &= \frac{1}{1 + e^{-\beta_0 + \beta_1x_i + \epsilon_i}} \cdot (1 - \frac{1}{1 + e^{-\beta_0 + \beta_1x_i + \epsilon_i}}) \cdot \beta_1 \\\\ &= P(y_i = 1 | x_i) \cdot (1 - P(y_i = 1 | x_i)) \cdot \beta_1 \end{aligned}$$

このように、ロジスティック回帰の限界効果は説明変数に依存しており$P(y_i = 1 | x_i)$の値で算出できること、限界効果はレコードごと($i$ごと)に異なることが分かります。
線形な回帰の場合、限界効果は説明変数に依存することはありませんし、全てのレコードで共通の値となります。

限界効果は説明変数に依存するため、予めどのような特徴を持つレコードにおける効果を算出したいのかを決める必要があります。
試験の例で言えば、単に母集団における平均的な限界効果を知りたければ、全ての生徒の限界効果を平均するのが良いでしょう。
日頃から勉強熱心な生徒における限界効果を知りたければ、勉強時間が大きい生徒に限定して限界効果を平均する、もしくは、上の式に$x_i = 100$等の具体的な値を代入して得られる限界効果を利用するのも良いでしょう。
分析の目的に沿って算出するのが望ましいかと思います。

Pythonによる実装

Pythonでロジスティック回帰を行う場合、statsmodels と scikit-learn が候補に上がるかと思います。
statsmodels では限界効果だけでなくその信頼区間等も算出してくれるメソッドが実装されていて至れり尽くせりなので、限界効果を算出したい際はstatsmodelsを利用するのが推奨されます。

回帰の結果である statsmodels.discrete.discrete_model.LogitResults というクラスに get_margeff というメソッドが生えているので、これを利用します。
https://www.statsmodels.org/devel/generated/statsmodels.discrete.discrete_model.LogitResults.get_margeff.html

今回はタイタニックのデータに適用してみます。生存率を被説明変数、性別と年齢を説明変数としてロジスティック回帰を実行します。

import pandas as pd
import statsmodels.api as sm

# データ読み込み
df = pd.read_csv('https://raw.githubusercontent.com/pandas-dev/pandas/master/doc/data/titanic.csv')

# 前処理
df_ = df[df['Age'].notnull()]
y = df_['Survived'].astype('float')
X = pd.concat(
    [
        pd.get_dummies(df_['Sex'], drop_first=True)
        , df_['Age']
    ]
    , axis=1
).astype('float')

# ロジスティック回帰の実行
model = sm.Logit(y, X)
res = model.fit()

# 限界効果の算出
print(res.get_margeff(at='overall').summary())

# =====================================
# Dep. Variable:               Survived
# Method:                          dydx
# At:                           overall
# ==============================================================================
#                 dy/dx    std err          z      P>|z|      [0.025      0.975]
# ------------------------------------------------------------------------------
# male          -0.3808      0.017    -21.771      0.000      -0.415      -0.347
# Age            0.0041      0.001      6.515      0.000       0.003       0.005
# ==============================================================================

出力結果を見ると、年齢が1増えるごとに生存率が0.41%増加することが分かります。

引数の at で限界効果の算出方法を変更することができます。今回指定したat='overall' では、全レコードに対して算出された限界効果を平均したものになります。

実行環境

Last updated: Fri Mar 01 2024

Python implementation: CPython
Python version       : 3.10.12
IPython version      : 7.34.0

statsmodels: 0.14.1
pandas     : 1.5.3

Watermark: 2.4.3