医療統計などでは「生存 / 死亡」のような結果に各因子がどの程度寄与しているのかをロジスティック回帰分析を用いて検討することがよくあります。ここではPythonを用いてロジスティック回帰分析を行う方法を見てみましょう。
開発環境
- pandas 1.3.3
- Python 3.7.11
- seaborn 0.11.2
- statsmodels 0.12.2
Contents
ロジスティック回帰分析
ロジスティック回帰分析の数学的背景
ある事象の確率pのロジットlogit(p)は
$$logit(p)=log(\frac{p}{1-p})$$
であり、説明変数\(x_1, x_2, x_3, … ,x_n\)を用いて
$$logit(p)=a_1*x_1+a_2*x_2+a_3*x_3+…+a_n*x_n+b$$
というように、説明変数と回帰係数\(a_1, a_2, …, a_n\)の１次結合で\(logit(p)\)を表すのがロジスティック回帰モデルです。これを\(p\)について解くと、標準ロジスティック関数で表せるので、
$$p=\frac{1}{1+e^{-logit(p)}}=\frac{1}{1+e^{-(a_1*x_1+a_2*x_2+a_3*x_3+…+a_n*x_n+b)}}$$
となり、説明変数\(x_1, x_2, x_3, … ,x_n\)によって目的変数\(p\)(その事象が起こる確率)が表されていることが分かります。
このように、ロジスティック回帰モデルではある事象の起こる確率を説明変数の1次結合を用いて表します。これにより「生存 / 死亡」「陽性 / 陰性」「悪性 / 良性」のような2値の質的変数(目的変数)を説明変数(「年齢・性別・検査データなど」)で表すことができ、このモデルを用いてそれぞれの説明変数の結果への寄与を推定する分析法をロジスティック回帰分析と言います。
詳細は以下をご覧ください。
Pythonを用いたロジスティック回帰分析
Pythonにおけるロジスティック回帰モデルの実装はscikit-learnが有名です。
ただし、scikit-learnはあくまでも機械学習のライブラリなので、説明変数(\(X_1, X_2, X_3, … \))を用いて目的変数(\(y\))をモデル化することは可能ですが、説明変数(\(X_1, X_2, X_3, … \))が統計的に有意差を持っているのかどうかやそのp値を求めることはできません。
実際にロジスティック回帰分析を用いて統計処理を行う場合は、目的変数(「生存 / 死亡」など)を説明変数(「年齢・性別・検査データなど」)の数式で表すことが目的ではなく、p値を求めてどの説明変数の項目がその結果に関わっている(有意差を持っている)かの検定を行ったり、そのオッズ比を求めることが目的となります。そのような用途のためにPythonではstatsmodelsという統計ライブラリが用意されているので、ここではstatsmodelsを用いてロジスティック回帰分析を行ってみましょう。
サンプルデータの取得と概観
ここでは有名なタイタニック号沈没事故の乗客の生存者/死亡者のデータセットをサンプルデータとして用います。タイタニックのデータセットを取得する方法はいくつかありますが、seabornのサンプルデータとしてのタイタニックのデータセットを読み込んでみましょう。
import seaborn as sns
titanic_data = sns.load_dataset('titanic')
print(titanic_data)
- survival : 生存 or 死亡 (0=死亡、1=生存)
- pclass : チケットクラス (1=上層クラス、2=中級クラス、3=下層クラス)
- sex : 性別 (male=男性、female=女性)
- age : 年齢
- sibsp : タイタニックに同乗している兄弟/配偶者の数
- parch : タイタニックに同乗している親/子供の数
- fare : 料金
- embarked : 出港地 (C=Cherbourg、Q=Queenstown、S=Southampton)
- who : 男性 or 女性
- adult_male : 成人男性であるかどうか
- deck : 乗船していたデッキ
- embark_town : 出港地
- alive : 生存したかどうか
- alone : 一人であったかどうか
ここでは説明変数として生死に関係のありそうな「pclass / sex / age」を選択します。なお、目的変数は「生存 / 死亡」を表すsurvivedです。それでは取得したデータセットから説明変数と目的変数のみを抽出しましょう。
titanic_data = titanic_data[['pclass', 'sex', 'age', 'survived']]
このデータセットの中身を表示してみます。
print(titanic_data)
pclass sex age survived 0 3 male 22.0 0 1 1 female 38.0 1 2 3 female 26.0 1 3 1 female 35.0 1 4 3 male 35.0 0 .. ... ... ... ... 886 2 male 27.0 0 887 1 female 19.0 1 888 3 female NaN 0 889 1 male 26.0 1 890 3 male 32.0 0 [891 rows x 4 columns]
sexの項目は「male / female」と文字列のカテゴリカル変数になっているので、解析前に数字に変換しておく必要があることが分かります。またageにはNaNが含まれているのが見えるので、欠損値もありそうですね。そこで、それぞれの項目のデータ数を表示して欠損値の有無を確認してみましょう。DataFrameのisnaメソッドで欠損値かどうかを判定し、その合計を求めましょう。
print(titanic_data.isna().sum())
pclass 0 sex 0 age 177 survived 0 dtype: int64
これによりageのみ欠損値があることが分かりました。
データの前処理
先ほど取得したデータ(titanic_data)には欠損値があり、またsexがカテゴリカル変数となっていることが分かりました。そこでこの2点について解析に適した形に前処理を行いましょう。
欠損値を含むデータを除外する
欠損値の扱い方はいろいろな流儀があり、特に機械学習の分野で予測モデルを作成する目的であれば中央値で埋めたりという方法がとられる場合が多いようです。それは欠損値を含んだデータの予測もせざるを得ないので、欠損値を含んだデータで予測モデルを作る必要があるからです。しかし、ここでは予測モデルを作成することではなく、それぞれの因子が結果にどの程度寄与しているのかを明らかにすることが目的なので、シンプルに欠損値を含むデータを除外してしまいましょう(ただし、欠損値を持つデータに偏りがないと仮定します)。
DataFrameから欠損値を含むデータを除外するのはdropnaメソッドで簡単に行うことができます。
titanic_data = titanic_data.dropna()
カテゴリカル変数を数字に変換する
続いて、「male / female」のようなどのカテゴリーに属するかを表す値も数字に変換します。ロジスティック回帰分析では説明変数の値を数式に代入して、その結果として目的変数を求めるので、数式に代入できる「0 / 1」のような値にしておく必要があるのです。
ここではDataFrameのカテゴリカル変数を数字に変換する方法として、get_dummies関数を用います。
titanic_data = pd.get_dummies(titanic_data, drop_first=True)
なお、ここでは多重共線性の問題を回避するためにdrop_first引数をTrueに設定しています。
これでもう一度データを概観してみましょう。
print(titanic_data)
pclass age survived sex_male 0 3 22.0 0 1 1 1 38.0 1 0 2 3 26.0 1 0 3 1 35.0 1 0 4 3 35.0 0 1 .. ... ... ... ... 885 3 39.0 0 0 886 2 27.0 0 1 887 1 19.0 1 0 889 1 26.0 1 1 890 3 32.0 0 1 [714 rows x 4 columns]
sexはsex_maleの「0 / 1」に変換されて「男性か否か」で表現されているのが分かります。また、行数は「714 rows」と表示されているので、ageが欠損値となっていたデータが除外されたことが分かります。
ロジスティック回帰分析を行う
データの前処理が完了したので、説明変数(X)と目的変数(y)のデータセットに分離しましょう。
X = titanic_data[['pclass', 'sex_male', 'age']]
y = titanic_data['survived']
statsmodelsでロジスティック回帰分析を行うにはロジットモデルを用いますが、そもそもロジットモデルは一般化線形モデル(GLM)のリンク関数にロジット関数を指定したものなので、一般化線形モデルを用いてロジスティック回帰分析を行うことも可能です。
ここではその2つの方法について見ていきましょう。
ロジットモデルを用いる方法
statsmodelsではロジットモデルはLogitクラスで用意されています。Logitクラスに目的変数・説明変数を指定してインスタンス化することで、目的変数と説明変数を登録したロジットモデルのインスタンスを生成します。そのインスタンスのfitメソッドを用いることで最尤法でそのロジットモデルの係数や切片を求めて、LogitResultsインスタンスのラッパーとして結果を返します。
それでは、先ほどのデータセットを用いて、ロジスティック回帰分析を行ってみましょう。結果として得られるLogitResultsインスタンスのsummaryメソッドを実行することで、その結果を要約して表示することができます。
import statsmodels.api as sm
logit_model = sm.Logit(y, X)
logit_model_result = logit_model.fit()
print(logit_model_result.summary())
Optimization terminated successfully.
Current function value: 0.548600
Iterations 5
Logit Regression Results
==============================================================================
Dep. Variable: survived No. Observations: 714
Model: Logit Df Residuals: 711
Method: MLE Df Model: 2
Date: Thu, 21 Oct 2021 Pseudo R-squ.: 0.1878
Time: 13:55:34 Log-Likelihood: -391.70
converged: True LL-Null: -482.26
Covariance Type: nonrobust LLR p-value: 4.691e-40
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
pclass -0.0402 0.063 -0.635 0.525 -0.164 0.084
sex_male -2.0429 0.185 -11.049 0.000 -2.405 -1.681
age 0.0242 0.005 5.335 0.000 0.015 0.033
==============================================================================
一般化線形モデル(GLM)にロジット関数を指定する方法
statsmodelsでは一般化線形モデルはGLMクラスで用意されています。GLMクラスに目的変数・説明変数を指定してインスタンス化することで、目的変数と説明変数を登録したGLMモデルのインスタンスを生成しますが、その際に使用する分布をfamily引数でBinomialとします。これでリンク関数はBinomialのデフォルトのLogitになりました。あとは先ほどと同様に、そのインスタンスのfitメソッドを用いることで最尤法でそのGLMの係数や切片を求めて、GLMResultsインスタンスのラッパーとして結果を返します。
それでは、先ほどのデータセットを用いて、ロジスティック回帰分析を行ってみましょう。結果として得られるGLMResultsインスタンスのsummaryメソッドを実行することで、その結果を要約して表示することができます。
import statsmodels.api as sm
glm_model = sm.GLM(y, X, family=sm.families.Binomial())
glm_model_result = glm_model.fit()
print(glm_model_result.summary())
Generalized Linear Model Regression Results
==============================================================================
Dep. Variable: survived No. Observations: 714
Model: GLM Df Residuals: 711
Model Family: Binomial Df Model: 2
Link Function: logit Scale: 1.0000
Method: IRLS Log-Likelihood: -391.70
Date: Thu, 21 Oct 2021 Deviance: 783.40
Time: 14:03:12 Pearson chi2: 756.
No. Iterations: 4
Covariance Type: nonrobust
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
pclass -0.0402 0.063 -0.635 0.525 -0.164 0.084
sex_male -2.0429 0.185 -11.049 0.000 -2.405 -1.681
age 0.0242 0.005 5.335 0.000 0.015 0.033
==============================================================================
ロジスティック回帰分析の結果の解釈
ロジットモデルを用いる場合と一般化線形モデルを用いる場合のいずれの方法でも上記の通りロジスティック回帰分析の結果が得られました。結果として得られるインスタンスはLogitResultsかGLMResultsかの違いがあるので、summaryメソッドで出力される項目に若干の違いがありますが、以下の下段の部分については全く同じであることが分かります。
==============================================================================
coef std err z P>|z| [0.025 0.975]
------------------------------------------------------------------------------
pclass -0.0402 0.063 -0.635 0.525 -0.164 0.084
sex_male -2.0429 0.185 -11.049 0.000 -2.405 -1.681
age 0.0242 0.005 5.335 0.000 0.015 0.033
==============================================================================
この部分はそれぞれ以下の情報を表します。
| ラベル | 説明 |
|---|---|
| coef | 回帰係数 |
| std err | 回帰係数の標準誤差 |
| z | zスコア |
| P>|z| | p値 |
| [0.025 0.975] | 回帰係数の95%信頼区間 |
つまり、上記の結果から性別(男性)と年齢がp値 < 0.05となり、有意であることが分かりました。この時の回帰係数が対数オッズ比に相当し、それぞれ-2.0429と0.0242となります。対数オッズ比はオッズ比の自然対数をとったものなので、オッズ比は\(e^{対数オッズ比}\)となり、それぞれ0.1296522と1.024495です。
ここで求められたオッズ比の概念については以下をご覧ください。
なお、statsmodelsではp値が小さいと0.000と表示されてしまいますが、正確な値を知りたい場合は以下のようにしてpvalues属性から取得することが可能です。
ロジットモデルを用いる場合
print(logit_model_result.pvalues)
pclass 5.252232e-01 sex_male 2.212353e-28 age 9.578283e-08 dtype: float64
一般化線形モデルを用いる場合
print(glm_model_result.pvalues)
pclass 5.252232e-01 sex_male 2.212344e-28 age 9.578276e-08 dtype: float64
ちなみにこの結果から分かるように、ロジットモデルを用いる場合と一般化線形モデルを用いる場合とで、その計算結果は完全に一致するわけではないようです。(もちろん、5桁目までは一致しているので実用上は全く問題ないと思われますが)
コードまとめ
最後に上記のコードをまとめてみましょう。ここではロジットモデルを用いる場合で、ロジスティック回帰分析のサマリーを表示させるコードをお示しします。
import pandas as pd
import seaborn as sns
import statsmodels.api as sm
# データの取得・確認
titanic_data = sns.load_dataset('titanic')
titanic_data = titanic_data[['pclass', 'sex', 'age', 'survived']]
print(titanic_data)
# データの前処理
titanic_data = titanic_data.dropna()
titanic_data = pd.get_dummies(titanic_data, drop_first=True)
print(titanic_data)
# 説明変数と目的変数に分ける
X = titanic_data[['pclass', 'sex_male', 'age']]
y = titanic_data['survived']
# ロジスティック回帰分析
logit_model = sm.Logit(y, X)
logit_model_result = logit_model.fit()
print(logit_model_result.summary())