筋肉で解決しないために。

日々出会うモノに対する考察をしたり、主に以下のテーマに関して書いています。 データサイエンス/人工知能/AI/機械学習/DeepLearning/Python//数学/統計学/統計処理

どんな妻が不倫しやすいかをロジスティック回帰してみた。

はじめに

このエントリでは、Statsmodelsの一部として提供されているサンプルデータを用いてロジスティック回帰を実装してみます。 サンプルデータは、1974年に行われた、女性に対して不倫の有無を聞いた調査です。 http://www.statsmodels.org/stable/datasets/generated/fair.html

おまじない

import numpy as np
import pandas as pd
from pandas import Series,DataFrame
import math

#プロット
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')
%matplotlib inline

#機械学習
from sklearn.linear_model import LogisticRegression
from sklearn.cross_validation import train_test_split
from sklearn import metrics
import statsmodels.api as sm

サンプルデータのインポート

df = sm.datasets.fair.load_pandas().data

データの概要を引用しておきます。

Number of observations: 6366 Number of variables: 9 Variable name definitions:

rate_marriage : How rate marriage, 1 = very poor, 2 = poor, 3 = fair, 4 = good, 5 = very good

age : Age

yrs_married : No. years married. Interval approximations. See original paper for detailed explanation.

children : No. children

religious : How relgious, 1 = not, 2 = mildly, 3 = fairly, 4 = strongly

educ : Level of education, 9 = grade school, 12 = high school, 14 = some college, 16 = college graduate, 17 = some graduate school, 20 = advanced degree

occupation : 1 = student, 2 = farming, agriculture; semi-skilled, or unskilled worker; 3 = white-colloar; 4 = teacher counselor social worker, nurse; artist, writers; technician, skilled worker, 5 = managerial, administrative, business, 6 = professional with advanced degree

occupation_husb : Husband's occupation. Same as occupation.

affairs : measure of time spent in extramarital affairs

df.head()
>>>

f:id:watarumon:20180709202346p:plain

データ概要と目的の整理(タスクの設定)

まず今回の目的は、ある諸データが与えられたとき、その人が不倫をしているかどうかを求めることのできるモデルの作成です。

そのため、データセットのaffairsは不倫をしている時間を表しているので、まずはこれを0と1で表現するように処理が必要です。これをタスク1と呼ぶことにします。

次に、数字が意味をもっていないものを処理してあげる必要があります。 そのカラムはoccupationとoccupation_husbです。 これは、数字を職業を表す記号として使用されているので、0と1で表現できるような何らかの処理が必要であります。 これをタスク2と呼ぶことにします。

タスク1の処理

まずは、タスク1から処理をしていきます。 0でない場合に1を返し、その他は0を返す関数をつくります。

def affair_check(x):
    if x != 0:
        return 1
    else:
        return 0

#Had_Affairという新しいカラムをつくります
df['Had_Affair'] = df['affairs'].apply(affair_check)

df.head()
>>>

f:id:watarumon:20180709202416p:plain

これで、affairsを0と1で表現することができました。

タスク2の処理

次にタスク2を処理していきます。

今回は、occupation(職業)のカテゴリが6種類あることより、列を6列作り、それぞれにマッチしていれば、1を返すことで、0と1で表現することとしたいと思います。 これは、pandas.get_dummies()関数によって簡単に処理することができます。 pandas.get_dummies()関数は、機械学習の前処理としてよく使われる手法です。

#ダミー変数に展開します。
occ_dummies = pd.get_dummies(df['occupation'])
hus_occ_dummies = pd.get_dummies(df['occupation_husb'])

#カラムに名前を付けます。
occ_dummies.columns = ['occ1','occ2','occ3','occ4','occ5','occ6']
hus_occ_dummies.columns = ['hocc1','hocc2','hocc3','hocc4','hocc5','hocc6']

#表示します
occ_dummies.head()

>>>

f:id:watarumon:20180709202520p:plain

状況整理(タスクの設定)

ここまでで、タスク1とタスク2を処理することができました。 次は、ダミー変数を用いたので、多重共線性を考慮しつつ、ここまで処理したdataframeを合体させたいと思います。 これを、タスク3と呼ぶことにします。

その次に、ロジスティックモデルにはめこめるように、説明変数をXというdataframeに、目的変数をYというnumpyarrayに処理しようと思います。 これをタスク4とと呼ぶことにします。

その後、目的のモデルを作成し、そのスコアを求めることをタスク5と呼ぶことにします。

最後に簡単にそのスコアの妥当性を検討します。 これをこれをタスク6と呼ぶことにします。

タスク3

ダミー変数を用いたので、多重共線性を考慮しつつ、ここまで処理したdataframeを合体させたいと思います。

多重共線性については、こちらのページが概要を理解するにはわかりやすかったです。 http://heartland.geocities.jp/ecodata222/ed/edj1-2-1-2-1.html

#occupationと、目的変数「Had_Affair」を削除します。
X = df.drop(['occupation','occupation_husb','Had_Affair','affairs'],axis=1)

#ダミー変数のDataFrameを繋げます。
dummies = pd.concat([occ_dummies,hus_occ_dummies],axis=1)

#Xとdummiesを繋げます。
X = pd.concat([X,dummies],axis=1)

#多重共線性を考慮して、occ1とhocc1を落とします。
X = X.drop('occ1',axis=1)
X = X.drop('hocc1',axis=1)
X.head()

>>>

f:id:watarumon:20180709202605p:plain

タスク4

タスク3で説明変数Xはやってしまいましたね。 目的変数Yをつくります。

Y = df.Had_Affair

#1次元のarrayにしておきます
Y = np.ravel(Y)
Y

タスク5

目的のモデルを作成し、そのスコアを求めます。

#LogisticRegressionクラスのインスタンスを作ります。
log_model = LogisticRegression()

#モデルを作ります。
log_model.fit(X,Y)

#モデルのスコア
log_model.score(X,Y)

>>>
0.7258875274897895

およそ73%一致のモデルを作成できました!

タスク6

簡単にこのスコアの妥当性を検討します。

まず、この問題は2値問題なので、例えば常に0と予想するモデルを作ったとして、その精度スコアと比べてみます。

Y.mean()

>>>
0.3224945020420987

68%の精度スコアということがわかります。この値よりは、先ほどのモデルの73%の方が高いことがわかります。

次に、sklearn.cross_validationで検証してみます。

#テストデータに分割
X_train, X_test, Y_train, Y_test = train_test_split(X, Y)

# 新しいモデルを作ります。
log_model2 = LogisticRegression()
log_model2.fit(X_train, Y_train)

# テスト用データを使って、予測。
class_predict = log_model2.predict(X_test)

# 精度スコアを出力。
print(metrics.accuracy_score(Y_test,class_predict))

>>>
0.7374371859296482

約73%の精度を得ました。すべてのデータを使ったときとおよそ同じ結果が得られました。

結論

どの変数が予測に寄与しているかを可視化してみます。

coeff_df = DataFrame([X.columns, log_model.coef_[0]]).T
plt.figure(figsize=(15,5))
sns.barplot(x=coeff_df[0], y=coeff_df[1])

f:id:watarumon:20180709201740p:plain

これを見ると、結婚生活が豊かなほど、宗教観が強いほど不倫しにくく、 逆に、結婚年月が長くなるほど、不倫率が上がる様です。

職業の面で言うと、夫の職業によって妻の不倫率はあまり変化なく、 3 = white-colloar、 5 = managerial,administrative, business, 6 = professional with advanced degree これらの、ホワイトカラー、経営者、実業家、高学歴専門職の様な人が不倫をしやすいことが示唆されます。

お金もステータスもある人が不倫しやすいって考えると、妥当かもしれませんね。

今回は、Statsmodelsの一部として提供されているサンプルデータを用いてロジスティック回帰を実装してみました。 お疲れ様でした!

それじゃー、また。