2. Pythonで綴る多変量解析 5-1. ロジスティック回帰分析(scikit-learn)

• 重回帰分析は、説明変数が複数ある回帰モデルで、目的変数は数値です。これに対して、ロジスティック回帰分析は、目的変数がカテゴリデータになります。
• ある商品を買うか買わないか、ある会社が倒産するかしないか、といったカテゴリ（属性）を扱います。しかも２分類だけでなく３分類以上にもつかえます。
• たとえば、ユーザーグループを分けるのに効果的な変数を見つけたい、個人特性やライフスタイル特性から自社・他社どちらの製品を買いそうか予測したい、といった目的にも利用されます。
• 分類を扱うアルゴリズムなのに、回帰という名がついているのも、なんとなく頷けます。
• さて、重回帰分析では、予測値と実測値との誤差がなるべく小さくなるような回帰係数を見つけてきます。ロジスティック回帰分析では、実測値に対して予測値が外れる確率（誤答率）の最小化を目指します。これを交差エントロピー誤差といいます。

####⑴ ライブラリのインポート

# データ加工・処理・分析ライブラリ
import numpy as np
import numpy.random as random
import scipy as sp
from pandas import Series, DataFrame
import pandas as pd

# 可視化ライブラリ
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
%matplotlib inline

# 機械学習ライブラリ
import sklearn

# 小数第3位まで表示
%precision 3

####⑵ データの取得と読み込み

# データを取得
url = 'https://raw.githubusercontent.com/yumi-ito/datasets/master/datasets_adult.csv'

# 取得したデータをDataFrameオブジェクトとして読み込み
df = pd.read_csv(url, header=None)

# 列ラベルを設定
df.columns = ['age', 'workclass', 'fnlwgt', 'education', 'education-num', 
              'marital-status', 'occupation', 'relationship', 'race', 'sex', 
              'capital-gain', 'capital-loss', 'hours-per-week', 'native-country', 'flg-50K']

• 職業、学歴、婚姻歴、人種などの属性情報をもとに、年収が50万ドルを超えるかどうかを予測するデータセットです。
• アメリカの1994 年国勢調査データベースから抽出されたもので、詳細はこちら。https://archive.ics.uci.edu/ml/datasets/adult
• 概要を以下に整理しました。説明変数は、数値型・属性型を合わせて全14個あります。目的変数は一つ、年収が「50万ドル超か否か」であり、標本数は32,561人となっています。

	変数名	意訳	項目（解説）	データ型
0	age	年齢	連続値	int
1	workclass	職業階級	Self-emp-not-inc(個人事業主), Self-emp-inc(自営業法人), State-gov(州政府), など	varchar(可変長文字列)
2	fnlwgt	最終的な重み	連続値（回答したユニットが表す対象母集団のユニット数）	int
3	education	学歴	Bachelors(学士), HS-grad(高校卒), Prof-school(専門職系大学院), など	varchar
4	education-num	就学年数	連続値	int
5	marital-status	婚姻歴	Married-civ-spouse(既婚の配偶者), Divorced(離別), Never-married(未婚), など	varchar
6	occupation	職業	Sales(営業職), Exec-managerial(経営幹部), Adm-clerical(事務職), など	varchar
7	relationship	同居家族	Wife(妻), Own-child(実子), Husband(夫), など	varchar
8	race	人種	White(白人), Asian-Pac-Islander(アジア系), Black(黒人), など	varchar
9	sex	性別	Female(女性), Male(男性).	varchar
10	capital-gain	固定資産売却益	連続値	int
11	capital-loss	固定資産売却損	連続値	int
12	hours-per-week	週当たり稼働時間	連続値	int
13	native-country	国籍	United-States(アメリカ), Japan(日本), Cuba(キューバ), など	varchar
target	flg-50K	年収が50Kドル以上か否か	50K, <=50K.	varchar

# データの形式と欠損数を出力
print(df.shape)
print('欠損の数：{}'.format(df.isnull().sum().sum()))

# データの先頭5行を出力
df.head()

####⑶ データの前処理（目的変数を変換）

# 目的変数の各項目の個数を確認
df.groupby('flg-50K').size()

• 年収50万ドル以上は7841人で、全体の約24％となっています。
• ちなみに、50万ドルは日本円で53,062,750.00 円（2020年7月24日 22:30 UTC）、つまり年収5000万円超の高額所得者が四分の一近くいます。
• 目的変数のデータの値は「<=50K」と「>50K」で、このままでは扱いにくいので、ダミー変数（該当=1、非該当=0）に変換します。

# カラムを追加し、1または0のフラグを設定
df['fin-flg'] = df['flg-50K'].map(lambda x: 1 if x ==' >50K' else 0)
print(df)

# 各項目の個数を再確認
df.groupby('fin-flg').size()

####⑷ モデル構築

# ロジスティック回帰のモデル構築のためのインポート
from sklearn.linear_model import LogisticRegression

# データ分割（訓練データとテストデータ）のためのインポート
from sklearn.model_selection import train_test_split

• 説明変数として、数値型のage、fnlwgt、education-num、capital-gain、capital-lossの5変数を採用します。
• 目的変数は、フラグを立てたfin-flgです。

# 説明変数と目的変数の設定
x = df[['age', 'fnlwgt', 'education-num', 'capital-gain', 'capital-loss']]
y = df['fin-flg']

• 説明変数(X)と目的変数(Y)を、それぞれ訓練データ(train)とテストデータ(test)に分けます。
• 訓練データをつかった学習でモデルを生成し、テストデータをつかってそのモデルを評価します。

# 訓練データとテストデータに分割
X_train, X_test, Y_train, Y_test = train_test_split(x, y, test_size=0.5, random_state=0)

• 引数random_state=0は、乱数の生成を固定し（この場合は0に設定）、何度実行しても同じように分類させます。ある行が、実行するたびに訓練用に分類されたりテスト用に分類されたりするのを制御して結果を一定し、モデルに再現性を持たせます。
• ここではtest_sizeを0.5として訓練用・テスト用を半々に分けています。

# LogisticRegressionクラスの初期化
model = LogisticRegression()

# 学習の実行
model.fit(X_train, Y_train)

####⑸ モデル評価

# 訓練データの正解率
train_score = format(model.score(X_train, Y_train))
print('正解率(train):', train_score)

# テストデータの正解率
test_score = format(model.score(X_test, Y_test))
print('正解率(test):', test_score)

• 訓練データ、テストデータともに約80％の正解率となっており、過学習は起きていないと判断できます。

# 学習済みモデルの各変数の係数を取得
regression_coefficient = model.coef_

# 行ラベル･列ラベルを付与してDataFrameに変換
pd.DataFrame(regression_coefficient.T,
             index = [x.columns.values],
             columns = ['回帰係数'])

• 説明変数の各々が目的変数に及ぼす影響度を表わす係数はモデル名.coef_で取得します。
• 固定資産の売却損益(capital-gain、capital-loss)は年収に影響する、というわかりやすい結果です。
• 次いで、オッズ比を算出してみます。

# 学習済みモデルの各変数のオッズ比を取得
odds_ratio = np.exp(model.coef_)

# 行ラベル･列ラベルを付与してDataFrameに変換
pd.DataFrame(np.exp(model.coef_).T,
             index = [x.columns.values],
             columns = ['オッズ比'])

• オッズ比は、それぞれの係数が1増加したとき、予測確率にどの程度影響があるかを示します。影響がないときは1.0となります。

---

追って、交差エントロピー誤差とオッズ比について、もう少し詳しく見ておきたいと思います。