こんにちは、えびかずきです。
今回はPythonでロジスティック回帰を使う方法について説明していきます。
こんな人におすすめ：
・2値のクラス分け問題をモデル化したい
・ロジスティック回帰とは何かを勉強したい
結論として、統計的な目的の場合「statsmodels」、機械学習的な目的の場合「scikit-learn」を使って実装するのが良いです。
では、順を追って解説していきましょう。
開発環境
Python 3.7.3
StatsModels 0.12.2
scikit-learn 0.24.2
NumPy 1.20.2
Pandas 1.2.3
IDE：Jupyter Notebook
回帰分析に使う生データ
ロジスティック分類は、2値分類の用途で使います。
ということで今回は、ある資格試験の合否に関する2値問題を題材として扱います。
ここでは、試験の合否[result]とその勉強時間[hours]で1セットのデータを100人分用意しました。（私が勝手に作ったものです）
これを使って資格試験の合否を[0:不合格,1:合格]の2値分類に見立て、ロジスティック回帰のモデルにフィッティングしていきましょう。
データはcsvで貼り付けてますので、自由に使ってください。
実装手順
statsmodelsを使う場合
まず、「statsmotels」というライブラリを使って実装する方法を説明します。
statsmodelsは、以下のようにpipでインストールできます。NumPy、Pandasと共に事前にインストールしておいてください。
$ pip install statsmodelsでは、Pythonによる実装に入っていきます。
まずライブラリのインストールです。
# ライブラリのインポート
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
import statsmodels.api as sm
import seaborn as sns
sns.set()
# jupyter Notebookでグラフ表示する設定
%matplotlib inline続いて、csvの生データをPandasで読み込みます。
データをグラフ化したものが下図です。
# データの読み込み
exam_data = pd.read_csv("exam_data.csv")
# データのグラフ化
# 縦軸は平均値、ciは信頼区間（ここでは表示しない）
sns.barplot(x = "hours",y = "result",
data = exam_data, ci=None, palette='Blues' )ここではデータの表示にbarplotを使いました。
10きざみの勉強時間に対して、結果が0~1の範囲で示されています。
生データにはresultの値として、0（不合格）か1（合格）のどちらかしかありませんが、ここで示されているのは各勉強時間における平均値です。
つまり縦軸は合格率を表していることになります。
グラフをみると、長く勉強するほど合格率が高まることがわかりますね。
信頼区間表示について：
seabornのbarplotは引数の「ci」に値を指定しない場合、デフォルトでブートストラップ法に基づく信頼区間が表示されます。ci=”sd”とすれば、標準偏差が表示されます。詳細は以下参照。
seaborn.barplot/公式ドキュメント
ブートストラップ法/Wikipedia
ブートストラップ法で信頼区間を求める時の注意点
それでは、ロジスティック回帰のモデルを作っていきましょう。
ロジスティック回帰は「statsmodels.formula.api」の一般線形モデル（glm）を使って実装します。
# モデル化
logistic_model = smf.glm(formula = "result ~ hours",
data = exam_data,
family=sm.families.Binomial(link=sm.genmod.families.links.logit())
)
logistic_result = logistic_model.fit() 一般線形モデル（glm）の引数の説明：
formula：”目的変数~説明変数”の形式で指定。
data：データセットをPandasのDataFrameで指定。
family：データの分布関数とリンク関数を指定。
※分布関数に二項分布（Binomial）を指定すればロジスティック回帰になる。
※Binomialを指定すると自動でリンク関数をlogitとしてモデルが作られるが、ここではあえて明示的にリンク関数を指定した。
さてそれでは回帰の結果をみていきましょう。
結果は下のようにsummaryメソッドで確認できます。
# 結果の出力
logistic_result.summary()これを実行すると、以下のようなデータが出力されます。
いろいろ出てますが、赤枠に回帰の係数などの情報が記載されています。
フィッティングしたモデルを使って、予測をやってみましょう。
予測はpredictメソッドを使って以下のようにします。
ここでは10,20,30,・・・,90という10きざみのデータをテストデータとして予測を実施していきます。
# モデルを使った予測のやり方
# 0~90の10間隔のテストデータを生成
test_data = pd.DataFrame({
"hours": np.arange(0, 100, 10)
})
# 成功確率の予測値
pred =logistic_result.predict(test_data)
pred結果は下のように合格確率として表示されます。
scikit-learnを使う場合
続いて、scikit-learnを使ってロジスティック回帰を実装する方法を説明していきます。
scikit-learnはpipで以下のようにしてインストールできます。
$ pip install scikit-learnではPythonコードを書いていきましょう。
まずライブラリのインポートとデータの読み込みです。
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
# データの読み込み
test_result = pd.read_csv("exam_data.csv")
hours = test_result['hours']
result = test_result['result']
X_train = hours.values.reshape(-1, 1)
Y_train = result.valuesロジスティックモデルはsklearn.linear_modelのLogisticRegressionクラスを使って実装します。
# ロジスティック回帰モデルの作成（インスタンス化）
lr = LogisticRegression()
# ロジスティック回帰モデルを学習
lr.fit(X_train, Y_train) 回帰結果は以下のようにして取得できます。
# 結果
print("coefficient = ", lr.coef_)
print("intercept = ", lr.intercept_)回帰係数として以下のアウトプットが表示されました。
statsmodelsで実行した結果と同じですね。
scikit-learnでフィッティングしたモデルを使って予測をする場合は、以下のようにします。
ここでも10,20,30,・・・,90という10きざみのデータをテストデータとして予測を実施していきます。
# 予測
test_data = np.arange(0, 100, 10)
X_test = test_data.reshape(-1,1)
probs = lr.predict_proba(X_test)
print(probs)結果として以下のアウトプットを得ました。
ここでは不合格の確率（左）と合格の確率（右）がセットで表示されています。
表示のされ方は異なりますが、結果自体はstatsmodelsの時と同じですね。
ロジスティック回帰とは何か？
最後にロジスティック回帰とは何かについて簡単に説明しておきます。
ロジスティック回帰とは一般線形回帰という手法の一つで、コインの表裏みたいな二項分布に従う事象の確率を、ロジット関数を使って線形回帰モデルに当てはめる方法のことを言います。
これだけだと頭に入ってこないと思うので、具体的に説明しましょう。
今回紹介した資格試験の合否をロジスティック回帰でモデル化する場合は、資格試験で合格する確率pを、下のようにロジット関数で線形問題に当てはめます。
ここで左辺の関数がロジット関数です。
\(\displaystyle log(\frac{p}{1-p})=\beta_0+\beta_1x\)
\(\beta\)は回帰係数,\(x\)は勉強時間を表す
これを変形すると、
\(\displaystyle p=\frac{1}{1+e^{-(\beta_0+\beta_1x)}}\)
となります。
この時の右辺がロジスティック関数と呼ばれます。
今回の記事で紹介したロジスティック回帰を図示すると、以下のようになります。試験の結果は0（不合格）と1（合格）の2値しかないわけですが、その確率が勉強時間とが、ロジット関数を使うことで線形で結びつけられた結果、このような曲線になります。
ロジット関数の中身のp/1-pの部分は、オッズ比とも呼ばれています。
したがってロジスティック回帰は、対数オッズ比の線形回帰であると説明されることもあります。
もっと詳しく知りたい方は、記事下の参考をご参照ください。
statsmodelsとscikit-learnはどちらを使えば良い？
今回ロジスティック回帰の実装について二つ方法を説明しました。
statsmodelsとscikit-learn。
さてどちらを使えば良いでしょうか？
極論やっていることは同じなので、どちらを使っても良いのですが、データ解析の目的に応じて選択するのが良いです。
要するに、目的が統計的ならstatsmodels、か機械学習的ならscikit-learnということです。
データの構造を考察するという目的で解析する統計処理的な用途ならはstatmodelsを使う方が自然ですし、未来のデータを予測するというような機械学習的なようとであればscikit-learnを使う方が自然です。
自分の用途に応じて、どちらを使うか選択してみてください。
ちなみに、これはロジスティック回帰に限った話ではありません。
統計学と機械学習はかなり密接に関係していますから内容としてオーバーラップする部分は多いです。
つまり何らかのツールを使ってデータ解析をする場合は、内容は似ていても統計的な文脈のものと、機械学習的な文脈のものとで使いやすさが違ってくるケースがありますので、作業を始める前に一度確認することをお勧めします。
まとめ
今回は、Pythonでロジスティック回帰を使う方法について説明しました。statsmodelsとscikit-learnの２種類の実装方法を解説しましたが、どちらを使うかは解析の目的に応じておのおの考えてみてください。
いろいろ書きましたが、そもそもPython自体が機械学習的な文脈で使われることが多いので、迷ったらscikit-learnを使う方が自然でしょう。
統計ガチ勢はPythonでなくRを好む場合が多いですからね。
参考
記事作成にあたって、以下の情報を参考にしました。
URL
統計処理的な目的の場合
機械学習的な目的の場合
Scikit-learn でロジスティック回帰（確率予測編）/Qiita
参考書籍
統計処理的な目的の場合
第6部3章にロジスティック回帰の実装例があります。
実装例だけでなく、理論面が割と易しめに解説されているのでおすすめ。
機械学習的な目的の場合
2.3.3.5 クラス分類のための線形モデルにてロジスティック回帰の実装例があります。言わずと知れた機械学習実務の入門書です。
