pythonを使ったデータ探索と回帰【kaggle, EDA,randomforest】

今回は[コレ](https://www.kaggle.com/hely333/eda-regression)を紹介しながら写経します。

#データについて
試用するデータは保険料のデータ(と思われる)。

2019/07/30追記 : データは健康保険から請求される費用とのこと。

データには住所や年齢、喫煙、子供の数などと、保険料が入っている。

このデータをEDAで色々な角度から見ていってモデルの変数選択をして、
最終的に回帰で保険料を算出しようという試み。

#早速はじめます

```python
import numpy as np 
import pandas as pd 
import os
import matplotlib.pyplot as pl
import seaborn as sns
import warnings
warnings.filterwarnings('ignore')
data = pd.read_csv('C:/・・・/insurance.csv')
```

いつもの

```py
data.head()
data.isnull().sum()
```

```
 	age 	sex 	bmi 	children 	smoker 	region 	charges
0 	19 	female 	27.900 	0 	yes 	southwest 	16884.92400
1 	18 	male 	33.770 	1 	no 	southeast 	1725.55230
2 	28 	male 	33.000 	3 	no 	southeast 	4449.46200
3 	33 	male 	22.705 	0 	no 	northwest 	21984.47061
4 	32 	male 	28.880 	0 	no 	northwest 	3866.85520
```

```
age         0
sex         0
bmi         0
children    0
smoker      0
region      0
charges     0
dtype: int64
```

データにNaNが含まれてない日はついている日です。
カテゴリカルデータをエンコードして、変数の相関を計算します。

```py
from sklearn.preprocessing import LabelEncoder
#sex
le = LabelEncoder()
le.fit(data.sex.drop_duplicates()) 
data.sex = le.transform(data.sex)
# smoker or not
le.fit(data.smoker.drop_duplicates()) 
data.smoker = le.transform(data.smoker)
#region
le.fit(data.region.drop_duplicates()) 
data.region = le.transform(data.region)
```

カテゴリカルデータにフィルターをかけて重複しているものをファクターに変換しています。
regionのような変数にはonehotencorderなどを使用すべきなのでしょうが今回は怠慢にもlabelencorderで行くようです。

```py
data.corr()['charges'].sort_values()
```

```
region     -0.006208
sex         0.057292
children    0.067998
bmi         0.198341
age         0.299008
smoker      0.787251
charges     1.000000
```


```py
f, ax = pl.subplots(figsize=(10, 8))
corr = data.corr()
sns.heatmap(corr, mask=np.zeros_like(corr, dtype=np.bool), cmap=sns.diverging_palette(240,10,as_cmap=True),
            square=True, ax=ax)
```

![01.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/4947d549-8a06-8797-e948-8fd1ce06814c.png)

相関を計算した結果、喫煙者であることと、保険料に正の相関がありそうです。
筆者はbmiと相関があると想定していたようです。

費用の分布を確認します。

```py
from bokeh.io import output_notebook, show
from bokeh.plotting import figure
output_notebook()
import scipy.special
from bokeh.layouts import gridplot
from bokeh.plotting import figure, show, output_file
p = figure(title="Distribution of charges",tools="save",background_fill_color="#E8DDCB")
hist, edges = np.histogram(data.charges)
p.quad(top=hist, bottom=0, left=edges[:-1], right=edges[1:],fill_color="#036564", line_color="#033649")
p.xaxis.axis_label = 'x'
p.yaxis.axis_label = 'Pr(x)'
```

元のソースコードでは

```
show(gridplot(p,ncols = 2, plot_width=400, plot_height=400, toolbar_location=None))
```

なのですが、動かなかったので

```
show(p)
```

![bokeh_plot.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/65d20f78-2e8b-b657-cf1a-a49c630d960f.png)

費用の分布をみるとべらぼうに高い値は少ないようです。
大体10000くらい。

```py
f= pl.figure(figsize=(12,5))

ax=f.add_subplot(121)
sns.distplot(data[(data.smoker == 1)]["charges"],color='c',ax=ax)
ax.set_title('Distribution of charges for smokers')

ax=f.add_subplot(122)
sns.distplot(data[(data.smoker == 0)]['charges'],color='b',ax=ax)
ax.set_title('Distribution of charges for non-smokers')
```

喫煙者と非喫煙者での費用分布の違いを見ると、
喫煙者はふたこぶの分布で、さらに費用は高い傾向にある。
喫煙していなければ10000くらいで大体のデータが収まっています。


![02.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/794f9b6f-3392-7ca6-c143-c359708333ae.png)


```py
sns.catplot(x="smoker", kind="count",hue = 'sex', palette="pink", data=data)
```


![03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/709426a7-ce6e-f9f6-c7d0-54bc4101924d.png)



女性は1,男性は0に相当します
喫煙者を確認してみると、男性のほうが喫煙者でありそうです。
この関係から、男性全体の治療費は女性よりも多くなりそうですね。

```py
sns.catplot(x="sex", y="charges", hue="smoker",
            kind="violin", data=data, palette = 'magma')
```

![04.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/77a4eff4-aba3-bfec-8358-74e0fcdcbf53.png)


バイオリンplotで治療費と人数の密度を見てみましょう
こっちのほうが分かりやすい視覚化ですよね


```py
pl.figure(figsize=(12,5))
pl.title("Box plot for charges of women")
sns.boxplot(y="smoker", x="charges", data =  data[(data.sex == 1)] , orient="h", palette = 'magma')
```

![05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/9771992e-1449-0ac9-3e07-4cd9dbf22117.png)


喫煙者のboxplot 女性から

```py
pl.figure(figsize=(12,5))
pl.title("Box plot for charges of men")
sns.boxplot(y="smoker", x="charges", data =  data[(data.sex == 0)] , orient="h", palette = 'rainbow')
```

![06.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/f0823e18-be4e-536e-87a8-1736690e6d69.png)


男性も。



年齢との関係を見てみましょう。

```py
pl.figure(figsize=(12,5))
pl.title("Distribution of age")
ax = sns.distplot(data["age"], color = 'g')
```

![07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/39bb3278-f73f-beab-badb-e343a6e032fc.png)


患者の最もヤングな年齢が18
シニアは64です
18歳で吸ってる人はいるのでしょうか？？

```py
sns.catplot(x="smoker", kind="count",hue = 'sex', palette="rainbow", data=data[(data.age == 18)])
pl.title("The number of smokers and non-smokers (18 years old)")
```

![08.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/dc23d1c7-c4e1-2dab-d6ca-533bc9e8024b.png)


あらま
18歳でも吸ってる人がいるんですね。
やはり18歳でも治療費は高くなるのでしょうか。

```py
pl.figure(figsize=(12,5))
pl.title("Box plot for charges 18 years old smokers")
sns.boxplot(y="smoker", x="charges", data = data[(data.age == 18)] , orient="h", palette = 'pink')
```

![09.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/1ad00956-3a30-3292-8192-83b76c60b387.png)


18歳の金額分布。
18歳でも喫煙者は治療費が高くなっているようですね。

```py
g = sns.jointplot(x="age", y="charges", data = data[(data.smoker == 0)],kind="kde", color="m")
g.plot_joint(pl.scatter, c="w", s=30, linewidth=1, marker="+")
g.ax_joint.collections[0].set_alpha(0)
g.set_axis_labels("$X$", "$Y$")
ax.set_title('Distribution of charges and age for non-smokers')
```

![10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/63efd7e7-f706-90c8-e8e6-6d3f5de3d666.png)

喫煙している人の金額分布
(若い人ほど単純にやすい)


```py
g = sns.jointplot(x="age", y="charges", data = data[(data.smoker == 1)],kind="kde", color="c")
g.plot_joint(pl.scatter, c="w", s=30, linewidth=1, marker="+")
g.ax_joint.collections[0].set_alpha(0)
g.set_axis_labels("$X$", "$Y$")
ax.set_title('Distribution of charges and age for smokers')
```


![11.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/34325864-06da-3c58-9b7f-665730fad488.png)


喫煙者の金額分布
年齢と金額をみると金額は二こぶになっている。


```py
#non - smokers
p = figure(plot_width=500, plot_height=450)
p.circle(x=data[(data.smoker == 0)].age,y=data[(data.smoker == 0)].charges, size=7, line_color="navy", fill_color="pink", fill_alpha=0.9)

show(p)
```

![bokeh_plot(1).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/2a3fe716-e16c-cbf1-d1f3-6443ce8e7ef1.png)



```py
#smokers
p = figure(plot_width=500, plot_height=450)
p.circle(x=data[(data.smoker == 1)].age,y=data[(data.smoker == 1)].charges, size=7, line_color="navy", fill_color="red", fill_alpha=0.9)
show(p)
```

![bokeh_plot(2).png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/cfac30e8-2ec1-05c0-3fed-f8a4fc1c54e5.png)



```py
sns.lmplot(x="age", y="charges", hue="smoker", data=data, palette = 'inferno_r', size = 7)
ax.set_title('Smokers and non-smokers')
```

![12.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/f984425b-579c-3a8a-b51e-4ab7865d2b94.png)


重ね合わせてplotみるとこんな感じ

喫煙者の2グループについて、それぞれのグループに分けて回帰をしている。

非喫煙者は直線が年齢と比例して増加していっている。
自然の摂理ですよね。健康第一!

喫煙者ではばらつきが大変大きく、単純に年齢を重ねたら高くなるのか？と聞かれると説明できない部分が多いデータになっている。

Bmiと比べてみましょう
治療費とbmiは関係があるのでしょうか？

```py
pl.figure(figsize=(12,5))
pl.title("Distribution of bmi")
ax = sns.distplot(data["bmi"], color = 'm')
```

![13.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/adbf1966-1f95-ffeb-67b4-d0274c2d242a.png)


Bmiの分布は平均30のきれいな分布です。
30ってどんなもんなのかグーグルに聞いてみましょう
30は肥満の始まりらしいです。Bmiを30で切り分け、30以上以下で費用分布を見ましょう

```py
pl.figure(figsize=(12,5))
pl.title("Distribution of charges for patients with BMI greater than 30")
ax = sns.distplot(data[(data.bmi >= 30)]['charges'], color = 'm')
```

![14.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/62628db3-4e9c-e11e-4c29-858e44bcecbe.png)


30以上の分布


```py
pl.figure(figsize=(12,5))
pl.title("Distribution of charges for patients with BMI less than 30")
ax = sns.distplot(data[(data.bmi < 30)]['charges'], color = 'b')
```

![15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/198e7085-c81a-d179-65b8-7f4df7a9fbfb.png)


30以下の分布

Bmiが30を超えるほど高額を支払っています


```py
g = sns.jointplot(x="bmi", y="charges", data = data,kind="kde", color="r")
g.plot_joint(pl.scatter, c="w", s=30, linewidth=1, marker="+")
g.ax_joint.collections[0].set_alpha(0)
g.set_axis_labels("$X$", "$Y$")
ax.set_title('Distribution of bmi and charges')
```

![16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/62277062-ea02-464d-ee4f-81cf22a8de25.png)


Bmiに関して治療費と合わせてみてみましょ

さらに喫煙者とbmiではどうなるでしょうか

```py
pl.figure(figsize=(10,6))
ax = sns.scatterplot(x='bmi',y='charges',data=data,palette='magma',hue='smoker')
ax.set_title('Scatter plot of charges and bmi')

sns.lmplot(x="bmi", y="charges", hue="smoker", data=data, palette = 'magma', size = 8)
```

![17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/7c38cea1-aee0-adb2-94e3-bef45d85ab96.png)


![18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/0e97df3c-a85b-2b17-cc2d-7e7665eb6f2a.png)


喫煙者であればbmiと治療費に関する直線の傾きが上がっていますね。


子供を持っていたら何か特徴はみられるのでしょうか？


```py
sns.catplot(x="children", kind="count", palette="ch:.25", data=data, size = 6)
```

![19.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/a4a7f16e-9164-1138-1365-2b27aae634b3.png)


ほとんど子供を持っていないようです
子供を抱えていたら喫煙したりするのでしょうか？

```py
sns.catplot(x="smoker", kind="count", palette="rainbow",hue = "sex",
            data=data[(data.children > 0)], size = 6)
ax.set_title('Smokers and non-smokers who have childrens')
```

![20.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/c81c9fbe-4016-cd72-e1b2-d335d4cfee0a.png)


子持ちに喫煙者がいますね。
でも喫煙してない親のほうがかなり多いことがわかります

#それぞれの変数から回帰

```py
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import PolynomialFeatures
from sklearn.metrics import r2_score,mean_squared_error
from sklearn.ensemble import RandomForestRegressor
```

```py
x = data.drop(['charges'], axis = 1)
y = data.charges

x_train,x_test,y_train,y_test = train_test_split(x,y, random_state = 0)
lr = LinearRegression().fit(x_train,y_train)

y_train_pred = lr.predict(x_train)
y_test_pred = lr.predict(x_test)

print(lr.score(x_test,y_test))
```

```
0.7962732059725786
```

線形回帰してみました
予測値と実際の値を比較して精度を測定すると0.79でした
データが必ずしも綺麗とは限らないので参考程度にしましょう

```py
X = data.drop(['charges','region'], axis = 1)
Y = data.charges

quad = PolynomialFeatures (degree = 2)
x_quad = quad.fit_transform(X)

X_train,X_test,Y_train,Y_test = train_test_split(x_quad,Y, random_state = 0)

plr = LinearRegression().fit(X_train,Y_train)

Y_train_pred = plr.predict(X_train)
Y_test_pred = plr.predict(X_test)

print(plr.score(X_test,Y_test))
```

```
0.8849197344147237
```

変数を減らしてシンプルにしています。
Regionを消すと、予測の精度は0.88まで上がりました。

ランダムフォレストを使ってみます

```py
forest = RandomForestRegressor(n_estimators = 100,
                              criterion = 'mse',
                              random_state = 1,
                              n_jobs = -1)
forest.fit(x_train,y_train)
forest_train_pred = forest.predict(x_train)
forest_test_pred = forest.predict(x_test)

print('MSE train data: %.3f, MSE test data: %.3f' % (
mean_squared_error(y_train,forest_train_pred),
mean_squared_error(y_test,forest_test_pred)))
print('R2 train data: %.3f, R2 test data: %.3f' % (
r2_score(y_train,forest_train_pred),
r2_score(y_test,forest_test_pred)))
```

```
MSE train data: 3729086.094, MSE test data: 19933823.142
R2 train data: 0.974, R2 test data: 0.873
```


ランダムフォレストといえばclassifierをよく見る気がしますが、
こんかいはregressorを使用していきます。

二条平均平方根MSEとR二乗値を確認します


```py
pl.figure(figsize=(10,6))

pl.scatter(forest_train_pred,forest_train_pred - y_train,
          c = 'black', marker = 'o', s = 35, alpha = 0.5,
          label = 'Train data')
pl.scatter(forest_test_pred,forest_test_pred - y_test,
          c = 'c', marker = 'o', s = 35, alpha = 0.7,
          label = 'Test data')
pl.xlabel('Predicted values')
pl.ylabel('Tailings')
pl.legend(loc = 'upper left')
pl.hlines(y = 0, xmin = 0, xmax = 60000, lw = 2, color = 'red')
pl.show()
```

![22.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/275572/6ff2762c-d696-8252-ef01-b2c5c39ee9a8.png)



X軸に予測値
Y軸に予測値と実際の値との差をplotしている
Trainとtestどちらもplotしている


#以上

seabornやbokehはplotがとても分かりやすくなって素敵な資料が作成できそうです。
そこに時間を割くかは置いといても、見やすいplotは見にくいplotよりは理解を深めてくれるのでマスターしていきたいです。

個人的にはlabelencorderの使用例が見られたのがうれしかった。

こんな感じでランダムフォレストで回帰がうごくのかって勉強になりました。