医療保険料のEDAと回帰問題

概要

会社で機械学習の知識がない課員にAIとかいうものの教育をする機会があるので、米国の医療保険料データを使ったEDA,クラスタリング,回帰問題の教材を作ることとした
データセットはKaggleのこれを使う
ざっと作ってみた感触は、ちょっと難しく感じられるかもしれなくって教材には使用しない
もっと簡単なものを別に用意することとしよう

• 実施期間: 2020年11月
• 環境：Google Colaboratory

データセット

insurance.csvの構成は下図の通り
医療保険料を示すchargesを今回の回帰による予測対象とする

各列の説明を転記しておく
Columns
age: age of primary beneficiary
sex: insurance contractor gender, female, male
bmi: Body mass index, providing an understanding of body, weights that are relatively high or low relative to height,
objective index of body weight (kg / m ^ 2) using the ratio of height to weight, ideally 18.5 to 24.9
children: Number of children covered by health insurance / Number of dependents
smoker: Smoking
region: the beneficiary's residential area in the US, northeast, southeast, southwest, northwest.
charges: Individual medical costs billed by health insurance

EDA(Explanatory Data Analysis)

馴染みのない言葉の方もおられるかもしれない
解析対象のデータセットの構造を確認し、どのように解析するか作戦を最初に立てるプロセスのこと（と思っている）
ツールの種類や慣れ、直感が必要
解析方法は人それぞれなので、以下は参考までにとどめていただきたい

まず、ファイルを開いて中身とデータに欠損がないか確認する

import numpy as np 
import pandas as pd 
import os
from sklearn.preprocessing import LabelEncoder
import matplotlib.pyplot as plt
import seaborn as sns

df = pd.read_csv('insurance.csv')
# Datasetの表示
print(df.tail())
# NaNの有無を確認
df.isnull().sum()

欠損がないことがわかる（あれば少し面倒なのでラッキー）
カテゴリデータが含まれているので、それらをLabelEncoderで数値化する

# カテゴリデータのエンコード
# sex
le = LabelEncoder()
le.fit(df.sex.drop_duplicates()) 
df.sex = le.transform(df.sex)
# smoker or not
le.fit(df.smoker.drop_duplicates()) 
df.smoker = le.transform(df.smoker)
# 　region
le.fit(df.region.drop_duplicates()) 
df.region = le.transform(df.region)
# Datasetの表示
print(df.tail()) 

整数に置き換わっている
全体の相関係数を眺めて、chargesに影響を与えている説明変数にあたりをつける

print(df.corr())

0.787251でsmokerがとても強い相関を示している
次が0.299008のage、肥満度を示すbmi(=体重/身長^2)の相関係数が思ったより小さいのが気になる
しかし、たった".corr()"だけで一発計算なので、Python使っててよかったと思うね

全体の分布を表示する

sns.pairplot(df.loc[:, ['age', 'sex', 'bmi', 'children', 'smoker', 'region', 'charges']], height=2);

特徴的な分布は下図なので、age,bmi,smokerに集中する
どうも3つのグループに分けられそうだが、何が原因で3つに分割されたのか解析する

age-charges分布がとても特徴的
相関が強かったsmokerについて掘ってみる

# smoker-chargesをハコヒゲで見てみる
sns.boxplot(x="smoker", y="charges", data=df)
# 数字で確認
df_s0 = df[df['smoker'] == 0]  # 非喫煙者
df_s1 = df[df['smoker'] == 1]  # 喫煙者
pd.set_option('display.max_rows', 20)
pd.set_option('display.max_columns', None)
print(' non smoker:\n' + str(df_s0.describe()))
print('\n smoker:\n' + str(df_s1.describe()))
print('\n 非喫煙者の医療保険料は、' + str(df_s1['charges'].mean() / df_s0['charges'].mean()) + '倍もかかっている。') 

箱ひげ図から非喫煙者はchargesヒストグラムの左の高い山に相当することがわかる
逆に喫煙者のオレンジはchargesヒストグラムの低い2つの山がまとまったものといえる

非喫煙者だけで分布を確認する

sns.pairplot(df_s0.loc[:, ['age', 'bmi', 'charges']], height=2); 

前述のage-charges分布が分離できた
あの特徴的な分布は喫煙か否かが影響してたらしい

非喫煙者でも年齢にかかわらず一定の割合で病院にかかるから、age-charges分布の上部散らばりがそれを表しているのか？
この散らばりとbmi-charges分布の散らばりがリンクしているようにみえるので、クラスタリングしてみる

from sklearn.cluster import KMeans
# from sklearn.cluster import MiniBatchKMeans

df_wk = df_s0.drop(['sex', 'children', 'smoker', 'region'], axis=1)
kmeans_model = KMeans(n_clusters=2, random_state=1, init='k-means++', n_init=10, max_iter=300, tol=0.0001, algorithm='auto', verbose=0).fit(df_wk)
# kmeans_model = MiniBatchKMeans(n_clusters=2, random_state=0, max_iter=300,  batch_size=100, verbose=0).fit(df_wk)

df_wk['cluster'] = kmeans_model.labels_
labels = kmeans_model.labels_
print(labels.sum())
color_codes = {0:'r', 1:'g', 2:'b'}
colors = [color_codes[x] for x in labels]
ax1 = df_wk.plot.scatter(x='age', y='charges', c=colors)
ax1 = df_wk.plot.scatter(x='bmi', y='charges', c=colors)

キレイに分離できなかったのは残念だが、前出の2つの分布図の散らばりはお互いに呼応しているといえる
ただ、BMIの大小にかかわらず全年齢でchargesが高い理由がわからない
childrenとregionはまだ調べてなく、どちらかというとchildrenのほうがbmiのばらつきに関係していそう
非喫煙者のchildren(=扶養家族数)ごとの分布を見れば、家族構成とbmiの関係がわかりそう

df_c0 = df_s0[df_s0['children'] == 0]  # children = 0
df_c1 = df_s0[df_s0['children'] == 1]  # children = 1
df_c2 = df_s0[df_s0['children'] == 2]  # children = 2
df_c3 = df_s0[df_s0['children'] == 3]  # children = 3
df_c4 = df_s0[df_s0['children'] == 4]  # children = 4
df_c5 = df_s0[df_s0['children'] == 5]  # children = 5
plt.figure(figsize=(12, 6))
plt.hist([df_c0['charges'], df_c1['charges'], df_c2['charges'], df_c3['charges'], df_c4['charges'], df_c5['charges']]

charges > 15000のヒストグラムに注目すると、chargesと子供の数に特徴は見えない
やっぱり、「非喫煙者でも年齢にかかわらず一定の割合で病院にかかるから」純粋にばらついているだけなのか？
charges > 15000においてbmiが正規分布でばらついているか検定してみる

df_b0 = df_s0[df_s0['charges'] > 15000]  # children = 0
plt.hist([df_b0['bmi']])

import scipy.stats as stats
stats.probplot(df_b0['bmi'], dist="norm", plot=plt)
stats.shapiro(df_b0['bmi'])

(0.9842953681945801, 0.3445127308368683)

ヒストグラムもQQプロットもシャピロ–ウィルク検定(p値=0.3445)も正規分布と言えないとは言えないことを示している
結局原因はなんかよくわからなかった
まぁ、正規分布ってことで、自然にばらついたことにしておこう…

今度は喫煙者だけで分布を確認する

sns.pairplot(df_s1.loc[:, ['age', 'bmi', 'charges']], height=2);

直感的に下図のようにageとbmiが呼応しているように見える

非喫煙者のときと同様にクラスタリングしてみる

df_wk = df_s1.drop(['sex', 'children', 'smoker', 'region'], axis=1)
kmeans_model = KMeans(n_clusters=2, random_state=1, init='k-means++', n_init=10, max_iter=300, tol=0.0001, algorithm='auto', verbose=0).fit(df_wk)
df_wk['cluster'] = kmeans_model.labels_
labels = kmeans_model.labels_
print(labels.sum())
color_codes = {0:'r', 1:'g', 2:'b'}
colors = [color_codes[x] for x in labels]
ax1 = df_wk.plot.scatter(x='age', y='charges', c=colors)
ax1 = df_wk.plot.scatter(x='bmi', y='charges', c=colors) 

思っていたとおりだが、今度はキレイに分離できた

EDAの結論

喫煙癖にかかわらず医療保険料が高額なのはBMIが高い人たちだが、年齢は関係ない（耳が痛い…）
喫煙者のデータからbmi=30を境に医療保険料が一段上昇することがわかった
この境がはっきりしすぎているので、例えば30を超えたら高額な検査項目が強制的に追加されるような行政上の縛りがあるのではないだろうか

sns.regplot(x="bmi", y="charges", data=df_wk)

つまり、例えば上図のような回帰直線を引いてはならないことが今回のEDAからわかる

単回帰、重回帰問題

非喫煙者について、chargesを回帰問題予測してみる
具体的には単回帰(ageのみ)と重回帰(age, bmi)でcharges予測精度の違いをRMSEで評価する
クラスタリング結果は使用しないので、ばらつきは含んだままとする

from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import mean_squared_error
from sklearn import metrics
from sklearn import linear_model

cv_alphas = [0.01, 0.1, 1, 10, 100]     #交差検証用のα
df_wk = df_s0.drop(['sex', 'children', 'smoker', 'region'], axis=1)

# Predict結果の評価とプロット
def regression_prid(clf_wk, df_wk, intSample=20, way='single'):
    arr_wk = df_wk.values[:intSample]
    arr_age, arr_agebmi, arr_charges = arr_wk[:,0], arr_wk[:,0:2], arr_wk[:,2]

    if way == 'single':
        arr_prid = clf_wk.predict(arr_age.reshape([intSample,1]))
        titel1='simple regression'
        titel2='RMSE(single) =\t'
    else:
        arr_prid = clf_wk.predict(arr_agebmi.reshape([intSample,2]))
        titel1='multiple regression'
        titel2='RMSE(multi) =\t'

    rmse = np.sqrt(mean_squared_error(arr_charges, arr_prid))
    print(titel2 + str(rmse))

    # プロット
    arr_chart = df_wk.values[:intSample].reshape([intSample,3])
    arr_chart = np.hstack([arr_chart, arr_prid.reshape([intSample,1])])

    fig = plt.figure(figsize=(8, 5))
    ax = fig.add_subplot(1,1,1)
    ax.scatter(arr_chart[:,0],arr_chart[:,2], c='red')
    ax.scatter(arr_chart[:,0],arr_chart[:,3], c='blue')
    ax.set_title(titel1)
    ax.set_xlabel('age')
    ax.set_ylabel('charges')
    fig.show() 

最小二乗法による単回帰

LinearRegression

clf = linear_model.LinearRegression()    

# 単回帰問題
clf.fit(df_wk[['age']], df_wk['charges'])
regression_prid(clf, df_wk, 50, 'single')

# 重回帰問題
clf.fit(df_wk[['age', 'bmi']], df_wk['charges'])
regression_prid(clf, df_wk, 50, 'multi')

上がageで予測した単回帰で下がageとbmiで予測した重回帰
RMSEは20点で予測させた時のもので、単回帰のほうがやや精度が高い
元データがほぼ線形なので単回帰でも重回帰でも予測は直線に乗っている

L2正則化による回帰

Ridge
RidgeCV
ちなみにCVはCloss Validation、交差検証のこと
交差検証の詳細はググってください

cv_flag = True   # True:交差検証

if cv_flag:
    clf = linear_model.RidgeCV(alphas=cv_alphas ,cv=3, normalize=False)

    # 単回帰問題
    clf.fit(df_wk[['age']], df_wk['charges'])
    print("alpha =\t", clf.alpha_)
    regression_prid(clf, df_wk, 50, 'single')

    # 重回帰問題
    clf.fit(df_wk[['age', 'bmi']], df_wk['charges'])
    print("alpha =\t", clf.alpha_)
    regression_prid(clf, df_wk, 50, 'multi')

else:
    clf = linear_model.Ridge(alpha=1.0)

    # 単回帰問題
    clf.fit(df_wk[['age']], df_wk['charges'])
    regression_prid(clf, df_wk, 50, 'single')

    # 重回帰問題
    clf.fit(df_wk[['age', 'bmi']], df_wk['charges'])
    regression_prid(clf, df_wk, 50, 'multi') 

交差検証を使って予測してみたが、やはり重回帰のほうがよろしくない

交差検証しなくてもRMSEはそれほど変わらない
やはりchargesがリニアに上昇しているから回帰手法による精度差が出にくなったと思う

L1正則化による回帰

Lasso
LassoCV
上のコード中のRidgeをLassoに置き換えるだけなので省略
RMSEはほとんど変わらなかった