Pythonで前処理から始めるロジスティック回帰
ロジスティック回帰の初め方
ロジスティック回帰をpythonでとにかく初めてみたい人向けの記事です。
前処理から評価まで一連の流れが解説されている記事ってあまり無いと思ったので書きました。
実務やコンペでロジスティック回帰を使いたい人（予測を重視）向けに、一通りの使い方を書いています。とりあえず使ってみたい人向けで、研究者とか仕組みを理解したい人向けでは無いかもしれません。
理論については下記書籍が大変分かりやすいです。
ロジスティック回帰とは？
分類予測タスクに使われる回帰モデルです。
ある病気に罹る/罹らない、ある顧客が商品を買うか/買わないか等の２値分類に利用されます。
また、オッズ比によって各説明変数が目的変数へ与える影響を確認できます。
教師データが少なかったり、予測の根拠が欲しいケースで使用されることが多いです。
ここでは、タイタニックデータをサンプルにして解説します。
なお、複数の分類を行う応用モデルもありますが、ここでは単純な2値分類タスクを対象にします。
基本的な前処理について
変数の特性によって前処理の方法はもちろん変わりますが、ここではロジスティック回帰利用時に用いる基本的な前処理について解説します。
欠損値への対処
ロジスティック回帰では変数に欠損があると処理できません。そのため、欠損値の除去や欠損値埋めを行う必要があります。
pandasの場合、isnull().sum()やinfo()で確認できます。下の例では欠損値が1個以上の変数のみを表示していますが、df.isnull().sum()だけでも問題ないです。
import seaborn as sns
df = sns.load_dataset('titanic')
df.isnull().sum()[df.isnull().sum() > 0]
実行結果は下記の通り、ageやembarked、deck、embark_townで欠損があることがわかります。
age 177
embarked 2
deck 688
embark_town 2
dtype: int64
欠損値除去
まずは欠損値の除去を説明します。欠損値があるレコードを削除する手法です。
上述のembarkedなど欠損レコードが少ない場合などに有効です。
下記例ではdropnaでembarked列について欠損値除去しています。
import seaborn as sns
df = sns.load_dataset('titanic')
df = df.dropna(subset = ['embarked'])
df.isnull().sum()[df.isnull().sum() > 0]
実行結果は下記の通りです。embarked列の欠損が消えていることが分かります。
その他の手法として、レコードの中で一つでも欠損があれば削除するリストワイズや複数の組み合わせで削除するレコードを決定するペアワイズなどもありますが、予測モデル作成ではあまり使われている印象はなく説明は他に譲ります。
age 177
deck 688
dtype: int64
欠損値埋め
欠損値除去を行うと簡単に欠損値処理を行うことができるのですが、一方でデータ数が減ってしまいます。そのため、欠損値を埋めてデータ数を減らさない工夫が求められます。
下記例ではfillnaとmeanを利用して平均値で欠損値埋めを行っています。
import seaborn as sns
df = sns.load_dataset('titanic')
df['age'] = df['age'].fillna(df['age'].mean())
df.isnull().sum()[df.isnull().sum() > 0]
実行結果は下記の通りです。age列の欠損が消えていることが分かります。
embarked 2
deck 688
embark_town 2
dtype: int64
その他にも、特定のカテゴリ変数と相関がある場合は、そのカテゴリごとの平均の代入も有効です。
今回の場合ですとclassによってageの平均が異なるため、classごとの平均を代入します。
df.groupby('class').mean()['age']
class
First 37.048118
Second 29.866958
Third 26.403259
Name: age, dtype: float64
下記の様にするとclassごとの平均の値を代入できます。
import seaborn as sns
df = sns.load_dataset('titanic')
df.loc[df.age.isnull(),'age'] = df.groupby('class')['age'].transform('mean')
予測モデルによる欠損値補完
欠損値のある変数と他の説明変数に相関がある場合などは、予測モデル作成し欠損値を予測するといった手法もあります。
下記は決定木で欠損値のある変数を目的変数とするモデルを作成し、欠損値補完を行っています。
import seaborn as sns
from sklearn.tree import DecisionTreeClassifier
#データロード,前処理
df = sns.load_dataset('titanic')
df = df.dropna(subset = ['embarked'])
df.loc[df.age.isnull(),'age'] = df.groupby('class')['age'].transform('mean')
X_train = df.loc[~df.deck.isnull(),['pclass','fare','age']]
y_train = df.loc[~df.deck.isnull(),'deck']
X_test = df.loc[df.deck.isnull(),['pclass','fare','age']]
#決定木による予測
clf = DecisionTreeClassifier(max_depth = 3)
clf = clf.fit(X_train,y_train)
df.loc[df.deck.isnull(),'deck'] = clf.predict(X_test)[0]
カテゴリ変数の処理
カテゴリ変数はそのままでは扱えないため、ダミー変数化を行い数値化します。
yes,noなど2値のカテゴリ変数であれば、[0,1]に変換することが一般的です。
Sランク,Aランク,Bランクなどの3値以上の水準があれば、それぞれに対してその水準かどうかを[0,1]で表現することが一般的です。ただし、後述の多重共線性のため複数の水準のうち一つを省略します。
| 変換前 | rank |
|---|---|
| RankSの場合 | S |
| RankAの場合 | A |
| RankBの場合 | B |
| 変換後 | rank_A | rank_B |
|---|---|---|
| RankSの場合 | 0 | 0 |
| RankAの場合 | 1 | 0 |
| RankBの場合 | 0 | 1 |
pandasのget_dummies関数を利用すると楽です。drop_firstをTrueにすることで複数の水準の内、最初の水準による変数が省略されます。
import seaborn as sns
df = sns.load_dataset('titanic')
display(df.loc[:5,'sex':'embarked'])
df = pd.get_dummies(df,drop_first = True,)
display(df.loc[:5,'sex_male':'embarked_S'])
実行結果は下記の通りです。ダミー変数化されていることが分かります。
sex age sibsp parch fare embarked
0 male 22.0 1 0 7.2500 S
1 female 38.0 1 0 71.2833 C
2 female 26.0 0 0 7.9250 S
3 female 35.0 1 0 53.1000 S
4 male 35.0 0 0 8.0500 S
5 male NaN 0 0 8.4583 Q
sex_male embarked_Q embarked_S
0 1 0 1
1 0 0 0
2 0 0 1
3 0 0 1
4 1 0 1
5 1 1 0
標準化について
線形回帰などでは標準化などを行う必要がありますが、ロジスティック回帰においては必須ではありません。モデリングの処理速度などが若干早くなると言う話がありますが、通常はあまり気にする必要はないかと思います。
変数削減について
多重回帰においても同様ですが、関連性が高い説明変数のペアが存在するとモデリングが不安定になることが知られています。
また、変数の数が多いと過学習と言う問題も発生する可能性があります。過学習とはモデルが学習データに適合しすぎてしまう事象で、未知のデータに対して高い精度の予測ができないという問題です。
（機械学習系の書籍で良く述べられている事象なので詳細は割愛します）
ここでは回避方法をいくつか解説しますが、単純に予測を目的とするのであれば、まずは正則化を行っていただくのが良いのでは無いかと思います。
相関係数から削除する変数を決定
変数間の相関係数を確認し、高いペアを見つけてどちらかを削除すると言う手法がデータ分析では良く使われている様です。
相関係数は下記の様にcorrメソッドにて求めることができます。1や-1に近いペアを確認します。
ただし、前提として相関係数はどの様な変数に対しても有効ではなく、各変数の分布などを確認する必要があります。
また、相関係数がいくつ以上であれば削除するかの定義が難しいです。
import seaborn as sns
df = sns.load_dataset('titanic')
df = pd.get_dummies(df,drop_first = True)
df.corr(method='pearson')
ステップワイズ法について
変数の選択を機械的に行う手法としてステップワイズ法があります。ステップワイズ法は説明変数の組み合わせを色々と試して最も良さそうな説明変数を導出する手法です。
pythonの有名どころのライブラリにはステップワイズの処理が無いため実装する必要があります。下記などが参考になります。
http://tomoshige-n.hatenablog.com/entry/2019/12/24/060514
正則化
pythonでサクッとロジスティック回帰を試すのであれば、正則化が一番簡単だと思います。
前処理として変数を削減するのではなく、（誤解を恐れずに言うのであれば）モデリング時に変数の選択を行う手法です。
理論などは下記が大変分かりやすいです。
https://toeming.hatenablog.com/entry/2020/04/03/000925
実装方法はロジスティック回帰実行時のパラメータ指定で実現できるのため、次章を確認ください。
次元削減
次元削減とはいくつかの変数をまとめて新たな別の変数として定義し直す手法です。因子分析と主成分分析が良く使われる手法ですが、予測モデルの作成であれば主成分分析の方が実行が簡単なのでまずはこちらをお勧めします。
なお、次元削減の前処理として標準化の処理が必要だったり、どの変数を次元削減の対象とするかの選択が必要なのでデータをしっかりと確認する必要がある。
主成分分析の手法としてPCAやt-SNEなどありますが、以下の例ではUMAPという比較的新しい手法を使用しています。
import umap
#標準化
target = df[['sibsp','parch','adult_male']]
target = (target - target.mean()) / target.std()
#UMAPでの時限削減
um = umap.UMAP()
df[['umap_1','umap_2']] = um.fit_transform(target[['sibsp','parch','adult_male']])
ロジスティック回帰の実行
下記の通りsklearnのロジスティック回帰のクラスを利用して実行する。
penalty='l1'にすることでL1正則化となり、いくつかの変数が削減される。
import pandas as pd
import seaborn as sns
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import confusion_matrix
#データロード
df = sns.load_dataset('titanic')
del df['alive']
#欠損値除去と欠損値埋め
df = df.dropna(subset = ['embarked'])
df.loc[df.age.isnull(),'age'] = df.groupby('class')['age'].transform('mean')
#決定木による欠損値埋め
clf = DecisionTreeClassifier(max_depth = 3)
clf = clf.fit(df.loc[~df.deck.isnull(),['pclass','fare','age']],
df.loc[~df.deck.isnull(),'deck'])
df.loc[df.deck.isnull(),'deck'] = clf.predict(df.loc[df.deck.isnull(),['pclass','fare','age']])[0]
#ダミー変数化
df = pd.get_dummies(df,drop_first = True)
#データセット作成
X = df.drop('survived',axis=1)
y = df['survived']
X_train,X_test, y_train,y_test, = train_test_split(X,y,test_size = 0.2)
#モデル作成
lr = LogisticRegression(penalty='l1', solver='liblinear')
lr.fit(X_train, y_train)
予測について
予測と正解率の算出を同時に行いたい場合は、scoreメソッドで可能。
0,1の予測結果を出したい場合はpredict、確率を算出したい場合はpredict_probaを利用する。
下記は正解率と、predictを利用した混同行列、predict_probaを利用したAUCの算出例。
from sklearn.metrics import roc_auc_score
#正解率
print('Accuracy={}'.format(lr.score(X_test,y_test)))
#混同行列
y_pred = lr.predict(X_test)
print('ConfusionMatrix\n{}'.format(confusion_matrix(y_test, y_pred)))
#AUC
y_pred_prob = lr.predict_proba(X_test)[:,1]
print('AUC={}'.format(roc_auc_score(y_test, y_pred_prob)))
実行結果
Accuracy=0.8146067415730337
ConfusionMatrix
[[97 15]
[18 48]]
AUC=0.8566693722943723
調整オッズについて
ロジスティック回帰の各係数や調整オッズ比は下記で求められます。
調整オッズ比は各変数が１上昇した際のオッズ比で、説明変数が目的変数へどの程度影響を与えるかの指標になります。1であれば影響は無しとなります。
import numpy as np
#ロジスティック回帰の各係数
display(pd.DataFrame(lr.coef_,columns=X.columns))
#調整オッズ比
display(pd.DataFrame(np.exp(lr.coef_),columns=X.columns))
各係数は下記の通り。
正則化によっていくつかの変数が0となり、変数削減されていることが分かります。
||pclass|age|sibsp|parch|fare|adult_male|alone|sex_male|embarked_Q|embarked_S|...|who_man|who_woman|deck_B|deck_C|deck_D|deck_E|deck_F|deck_G|embark_town_Queenstown|embark_town_Southampton|
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
|0|-0.201031|-0.014119|-0.582273|-0.31315|0.004862|-2.198662|-0.32373|0.0|0.0|-0.181672|...|-0.794377|0.06919|0.0|-0.22423|0.474059|1.060814|-0.53949|0.0|0.0|-0.003317|
調整オッズ比は下記の通り。
who_manやdeck_Eなどが目的変数への影響が大きいことが分かります。
||pclass|age|sibsp|parch|fare|adult_male|alone|sex_male|embarked_Q|embarked_S|...|who_man|who_woman|deck_B|deck_C|deck_D|deck_E|deck_F|deck_G|embark_town_Queenstown|embark_town_Southampton|
| ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- | ---- |
|0|0.817887|0.98598|0.558627|0.73114|1.004874|0.110951|0.723445|1.0|1.0|0.833875|...|0.451862|1.07164|1.0|0.799131|1.606501|2.888721|0.583046|1.0|1.0|0.996688|
Discussion