【Python】seabornのpenguins datasetを使って回帰分析をしてみる

seabornには様々な公開データがあります。今回はその中の「prnguins dataset」を使って単回帰分析や重回帰分析をしてみます。

目次[clickで開く/閉じる]

• データの概要
• データの中身の確認
• 単回帰分析の流れ
• 単回帰分析のソースコード
• 詳細
• 重回帰分析の流れ
• 重回帰分析のソースコード
• 詳細
• まとめ

データの概要

seabornはPythonデータ視覚化ライブラリでいくつかデータセットが用意されています。今回はその中でも「Paalmer Penguins」というペンギンの測定データを含んでいる「penguins dataset」(ペンギンデータセット)を用いて回帰分析に取り組みます。データセットに含まれる様々な変数に対して目的変数と説明変数を設定して分析を進めます。目的変数と説明変数を自由に設定できますがまずは目的変数をペンギンの種類として進めていこうと思います。
では、実際のデータの中身についてみていきましょう。

• データ数：344
• カラム数：7

変数名	詳細
island	ペンギンが生息する島の名前('Torgersen', 'Biscoe', 'Dream')
bill_length_mm	ペンギンのくちばしの長さ(mm)
bill_depth_mm	ペンギンのくちばしの奥行き(mm)
flipper_length_mm	ペンギンのヒレの長さ(mm)
body_mass_g	ペンギンの体重(g)
sex	ペンギンの性別('Male', 'Female')
species	ペンギンの種類('Adelie', 'Chinstrap', 'Gentoo')
year	データが収集された年

今回はyearは除外して扱います

各データの詳細については以下のページもご参考ください。
https://github.com/mwaskom/seaborn-data/blob/master/penguins.csv
https://github.com/allisonhorst/palmerpenguins

データの中身の確認

まずは基本的な必要ライブラリをまとめてインポートしておきましょう。

import seaborn as sns
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

データを出力してみます。この時についでにデータの大きさやデータの型、欠損値の有無について確認しておきましょう。

# データセットの読み込み
data = sns.load_dataset("penguins")
print(data.shape)
display(data.head())
print(data.dtypes)
print(data.isnull().sum())

dataframeで表示した結果欠損値が含まれていることがわかりました。それでは欠損値を含む行を削除していきましょう。

# 欠損値削除
data = data.dropna()
print(data.shape)
display(data.head())
print(data.dtypes)
print(data.isnull().sum())

きちんと削除できたか確かめます。

欠損値が全て0になり、データの大きさも元々あったsexの欠損値11行分少なくなり(333, 7)に更新されています。これで欠損値のないきれいなデータを作成できました。
回帰分析をするにはもう1つデータを変更しなければならないところがあります。island, sex, speciesなどのように定性的なデータを回帰分析で使用する場合はsexのMale = 0, Female = 1のように数値データとして扱う必要があります。これをダミー変数と言います。このようにカテゴリカルなデータを使用するときはダミー変数に変換することでデータの特徴を適切に表すことができます。それではisland, sex, speciesをダミー変数に変換していきましょう。

# speciesをダミー変数にする
data['species'] = data['species'].replace('Adelie', 0)
data['species'] = data['species'].replace('Chinstrap', 1)
data['species'] = data['species'].replace('Gentoo', 2)

# islandをダミー変数にする
data['island'] = data['island'].replace('Torgersen', 0)
data['island'] = data['island'].replace('Dream', 1)
data['island'] = data['island'].replace('Biscoe', 2)

# sexをダミー変数にする
data['sex'] = data['sex'].replace('Male', 0)
data['sex'] = data['sex'].replace('Female', 1)

print(data.shape)
display(data.head())

それではデータを確認してみます。

island, sex, species列の値が0, 1, 2で表されています。これでようやく分析の準備が整いました。
ここで一度データをペアプロットとヒートマップで可視化してみます。

# データの可視化
# 10s位かかります
sns.pairplot(data, hue='species', palette = 'muted')
plt.legend()
plt.title('Pairplot of penguins')
plt.savefig('pairplot.png')

# ヒートマップを表示
# ヒートマップを表示
plt.figure(figsize=(12, 9))
sns.heatmap(data.corr(), annot=True, cmap='coolwarm')
plt.title('Correlation Matrix')
plt.savefig('heatmap.png')

ついでに相関係数も確認しておきます。今回はわかりやすいように相関係数の大きさで色を付けてDataFrameで表示します。

# 相関係数を計算
corr = data.corr()
corr.style.background_gradient(cmap='coolwarm', axis=None)

単回帰分析

相関係数を求めて一番相関が大きかったデータを使用して単回帰分析をしてみます。今回は相関係数 = 0.872979になっているbody_mass_g, flipper_length_mm で行います。説明変数をflipper_length_mm, 目的変数をbody_mass_gにしてみます。

単回帰分析の流れ

1. データのダウンロード（データの中身の確認で実行済）
2. データを説明変数と目的変数に分割・外れ値の除去
3. 学習
4. 予測精度の確認 - 決定係数や回帰直線
5. 結果の図示 - 散布図などを使用する

単回帰分析のソースコード

# 説明変数と目的変数を設定
exp_var = 'flipper_length_mm'
tar_var = 'body_mass_g'

# 散布図を表示
plt.figure(figsize=(8, 6))
plt.scatter(data[exp_var], data[tar_var])
plt.xlabel(exp_var)
plt.ylabel(tar_var)
plt.savefig('scatter.png')

# 記述統計量を確認 
data[['flipper_length_mm', 'body_mass_g']].describe()

# 外れ値を除去
q_95 = data['flipper_length_mm'].quantile(0.95)
print('95%点の分位数', q_95)

# 絞り込む
data = data[data['flipper_length_mm'] < q_95]

# 散布図を表示
plt.figure(figsize=(12, 9))
plt.scatter(data[exp_var], data[tar_var])
plt.xlabel(exp_var)
plt.ylabel(tar_var)
plt.savefig('Penguins2.png')
plt.show()

# 記述統計量を確認
data[['flipper_length_mm', 'body_mass_g']].describe()

# 説明変数と目的変数にデータを分割
X = data[[exp_var]]
y = data[[tar_var]]

# 学習
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(data[[exp_var]], data[tar_var])

# 精度の確認
print('回帰直線の切片', model.intercept_)
print('回帰係数', model.coef_[0])
print('決定係数', model.score(data[[exp_var]], data[tar_var]))
print('回帰直線', 'y = ', model.coef_[0], 'x + ', model.intercept_)

# 回帰直線と散布図を表示
plt.figure(figsize=(12, 9))
plt.scatter(data[exp_var], data[tar_var])
plt.plot(data[exp_var], model.predict(data[[exp_var]]), color='red')
plt.xlabel(exp_var)
plt.ylabel(tar_var)
plt.title('penguins_regression2')
plt.savefig('penguins_regression.png')

説明変数と目的変数を指定する

# 説明変数と目的変数を設定
exp_var = 'flipper_length_mm'
tar_var = 'body_mass_g'

先ほど全体の概形はつかみましたが、今回使用する2変数に絞ってデータを見てみましょう。散布図と記述統計量を見てみます。

データの中身を確認

# 散布図を表示
plt.figure(figsize=(8, 6))
plt.scatter(data[exp_var], data[tar_var])
plt.xlabel(exp_var)
plt.ylabel(tar_var)
plt.savefig('scatter.png')

# 記述統計量を確認 
data[['flipper_length_mm', 'body_mass_g']].describe()

念の為外れ値の影響がないかみてみます。

# 外れ値を除去

q_95 = data['flipper_length_mm'].quantile(0.95)
print('95%点の分位数', q_95)

# 絞り込む
data = data[data['flipper_length_mm'] < q_95]

# 散布図を表示
plt.figure(figsize=(12, 9))
plt.scatter(data[exp_var], data[tar_var])
plt.xlabel(exp_var)
plt.ylabel(tar_var)
plt.savefig('Penguins2.png')
plt.show()

# 記述統計量を確認
data[['flipper_length_mm', 'body_mass_g']].describe()

外れ値除去前後の散布図に大きな変化はありませんが、記述統計量の標準偏差(std)に注目すると、両変数とも小さくなっていることがわかります。どっちのデータを使ってもほとんど変わらなさそうですが、折角外れ値除去したので単回帰分析では外れ値除去後の値を使って分析します。

説明変数と目的変数を設定する

目的変数(exp_var), 説明変数(tar_var)を別の変数に代入しておきます。
今後、学習に使うためです。

# 説明変数と目的変数にデータを分割
X = data[[exp_var]]
y = data[[tar_var]]

単回帰モデルを学習

今回のモデルは線形回帰による予測を行うクラスのLinearRegressionとします。LinearRegressionの詳細はこちら。
https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LinearRegression.html

# 学習
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(data[[exp_var]], data[tar_var])

精度の確認

学習精度の確認のために回帰直線と決定係数の値を確認します。

# 学習
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(data[[exp_var]], data[tar_var])

# 精度の確認
print('回帰直線の切片', model.intercept_)
print('回帰係数', model.coef_[0])
print('決定係数', model.score(data[[exp_var]], data[tar_var]))
print('回帰直線', 'y = ', model.coef_[0], 'x + ', model.intercept_)

# 回帰直線と散布図を表示
plt.figure(figsize=(12, 9))
plt.scatter(data[exp_var], data[tar_var])
plt.plot(data[exp_var], model.predict(data[[exp_var]]), color='red')
plt.xlabel(exp_var)
plt.ylabel(tar_var)
plt.title('penguins_regression2')

私の場合は決定係数が0.712となり、まあまあの精度が出ました。相関関係があった割には悪いような気もしますが、、、

続いて、重回帰分析に取り組みます。単回帰分析の時と同様にまずは分析の流れについて説明し、その後ソースコードを示してからコードの詳細に触れていきます。今回は目的変数がspecies, 説明変数がその他の変数になります。

重回帰分析の流れ

1. 説明変数と目的変数の決定・分割・外れ値の除去
2. 説明変数と目的変数それぞれを訓練データとテストデータに分割
3. 標準化
4. 訓練データの学習
5. 学習から得られた予測値・偏回帰係数の確認
6. 説明変数と目的変数それぞれでモデルの評価
   　　決定係数や二乗平均誤差（MSE）などを用いる
7. 正則化の実施
8. 精度の比較

重回帰分析のソースコード

# 説明変数
exp_vars = ['island', 'bill_length_mm', 'bill_depth_mm','flipper_length_mm', 'body_mass_g', 'sex']
# 目的変数
tar_var = 'species'

# 外れ値を除去
for exp_var in exp_vars:
    q_95 = data[exp_var].quantile(0.95)
    print(exp_var, '95%点の分位数', q_95)
    data = data[data[exp_var] < q_95]
    
# 説明変数
X = data[exp_vars]
# 目的変数
y = data[tar_var]

# 訓練データとテスとデータに分割
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state= 0)

# 標準化
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)
# print('X_train_scaled', X_train_scaled.shape)
X_train_scaled = pd.DataFrame(X_train_scaled, columns = exp_vars)
# display(X_train_scaled.head())

# 学習
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X_train_scaled, y_train)

# 予測値を計算
y_pred = model.predict(X_train_scaled)
y_pred[:10]

# 偏回帰係数を計算
xi_list = []
wi_list = []
for xi, wi in zip(exp_vars, model.coef_):
    xi_list.append(xi)
    wi_list.append(wi)
    # print(xi,':', wi)
    df_w = pd.DataFrame(wi_list, index = xi_list, columns = ['重回帰'])
display(df_w)

# 決定係数
print('決定係数:{:.3f}'.format(model.score(X_train_scaled, y_train)))
from sklearn.metrics import mean_squared_error

# 訓練データに対する Mean Squared Error (MSE)
mse_train = mean_squared_error(y_train, y_pred)
print("訓練データに対するMSE: {:.3f}".format(mse_train))

# テストデータに対する MSE
X_test_scaled = scaler.transform(X_test) # テストデータを訓練データから得られた平均と標準偏差で標準化
y_test_pred = model.predict(X_test_scaled) # テストデータに対して予測する
mse_test = mean_squared_error(y_test, y_test_pred)
print("テストデータに対するMSE: {:.3f}".format(mse_test))

# Ridge回帰
from sklearn.linear_model import Ridge
ridge = Ridge(alpha=1.0)
ridge.fit(X_train_scaled, y_train)
ridge_y_pred = ridge.predict(X_train_scaled)

# 偏回帰係数の確認
ridge_w = pd.DataFrame(ridge.coef_.T, index=exp_vars, columns=['Ridge'])
xi_list = []
wi_list = []
for xi, wi in zip(exp_vars, ridge.coef_):
    xi_list.append(xi)
    wi_list.append(wi)
    # print(xi,':', wi)
df_w_ridge = pd.DataFrame(wi_list, index = xi_list, columns = ['Ridge回帰'])
display(df_w_ridge)
print('L2ノルム: {:.3f}'.format(np.linalg.norm(ridge.coef_)))

# 訓練データに対する Mean Squared Error (MSE)
ridge_mse_train = mean_squared_error(y_train, ridge_y_pred)
print("訓練データに対するMSE: {:.3f}".format(ridge_mse_train))

# テストデータに対する MSE
ridge_y_test_pred = ridge.predict(X_test_scaled) # テストデータに対して予測する
ridge_mse_test = mean_squared_error(y_test, ridge_y_test_pred)
print("テストデータに対するMSE: {:.3f}".format(ridge_mse_test))

# 決定係数
print('決定係数:{:.3f}'.format(ridge.score(X_train_scaled, y_train)))
# 偏回帰係数を比較
df_w_all = pd.concat([df_w, df_w_ridge], axis=1)
df_w_all

# MSEと決定係数を比較
data = {'訓練データMSE':[mse_train, ridge_mse_train], 
        'テストデータMSE':[mse_test, ridge_mse_test], 
        '決定係数':[model.score(X_test_scaled, y_test), ridge.score(X_test_scaled, y_test)]}
df_mse = pd.DataFrame(data=data, index=['重回帰', 'Ridge回帰'])
df_mse

では、ソースコードを分解してみていきましょう。流れは単回帰の時とほとんど同じです。

説明変数と目的変数を指定する

目的変数をexc_var、説明変数をtar_varとし、データを分割します。
display分やprint文は必要に応じてコメントアウトしてください。

# 説明変数
exp_vars = ['island', 'bill_length_mm', 'bill_depth_mm','flipper_length_mm', 'body_mass_g', 'sex']
# 目的変数
tar_var = 'species'

データの内容の確認

単回帰分析と同様に外れ値を除去します。

# 外れ値を除去
for exp_var in exp_vars:
    q_95 = data[exp_var].quantile(0.95)
    print(exp_var, '95%点の分位数', q_95)
    data = data[data[exp_var] < q_95]

目的変数と説明変数を設定

# 説明変数
X = data[exp_vars]
# 目的変数
y = data[tar_var]

訓練データとテストデータに分割

# 訓練データとテスとデータに分割
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size = 0.2, random_state= 0)

標準化

# 標準化
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)
# print('X_train_scaled', X_train_scaled.shape)
X_train_scaled = pd.DataFrame(X_train_scaled, columns = exp_vars)
# display(X_train_scaled.head())

モデルの学習

# 学習
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X_train_scaled, y_train)

予測

# 予測値を計算
y_pred = model.predict(X_train_scaled)
y_pred[:10]

偏回帰係数を計算

# 偏回帰係数を計算
xi_list = []
wi_list = []
for xi, wi in zip(exp_vars, model.coef_):
    xi_list.append(xi)
    wi_list.append(wi)
    # print(xi,':', wi)
    df_w = pd.DataFrame(wi_list, index = xi_list, columns = ['重回帰'])
display(df_w)

偏回帰係数を見るとbill_length_mmが関係していそうです。
予測ができたところで、実際のデータと比較してどのくらいの精度で予測ができているのかを決定係数と二乗平均誤差（MSE）で確かめます。

モデルの評価

# 決定係数
print('決定係数:{:.3f}'.format(model.score(X_train_scaled, y_train)))

from sklearn.metrics import mean_squared_error

# 訓練データに対する Mean Squared Error (MSE)
mse_train = mean_squared_error(y_train, y_pred)
print("訓練データに対するMSE: {:.3f}".format(mse_train))

# テストデータに対する MSE
X_test_scaled = scaler.transform(X_test) # テストデータを訓練データから得られた平均と標準偏差で標準化
y_test_pred = model.predict(X_test_scaled) # テストデータに対して予測する
mse_test = mean_squared_error(y_test, y_test_pred)
print("テストデータに対するMSE: {:.3f}".format(mse_test))

私の場合は、

• 決定係数：0.881
• 訓練データに対するMSE：0.028
• テストデータに対すMSE：0.022
  となりました。
  決定係数は1に近いほど精度が高いと言えるので、単回帰分析よりも精度が高いことがわかります。MSEについては訓練データの値とテストデータの値にそれほど差がないので過学習はしていないと考えられます。

では次に正則化もしてみましょう。

正則化

正則化とは過学習を抑える仕組みの人るで、これを行うことで訓練データにモデルが過剰に適合してしまうのを防ぎ、未知のデータに対する予測能力が低下するのを防いでくれます。今回はL2ノルムのRidge回帰を試してみます。

# Ridge回帰
from sklearn.linear_model import Ridge
ridge = Ridge(alpha=1.0)
ridge.fit(X_train_scaled, y_train)
ridge_y_pred = ridge.predict(X_train_scaled)

# 偏回帰係数の確認
ridge_w = pd.DataFrame(ridge.coef_.T, index=exp_vars, columns=['Ridge'])
xi_list = []
wi_list = []
for xi, wi in zip(exp_vars, ridge.coef_):
    xi_list.append(xi)
    wi_list.append(wi)
    # print(xi,':', wi)
df_w_ridge = pd.DataFrame(wi_list, index = xi_list, columns = ['Ridge回帰'])
display(df_w_ridge)
print('L2ノルム: {:.3f}'.format(np.linalg.norm(ridge.coef_)))

# 訓練データに対する Mean Squared Error (MSE)
ridge_mse_train = mean_squared_error(y_train, ridge_y_pred)
print("訓練データに対するMSE: {:.3f}".format(ridge_mse_train))

# テストデータに対する MSE
ridge_y_test_pred = ridge.predict(X_test_scaled) # テストデータに対して予測する
ridge_mse_test = mean_squared_error(y_test, ridge_y_test_pred)
print("テストデータに対するMSE: {:.3f}".format(ridge_mse_test))

# 決定係数
print('決定係数:{:.3f}'.format(ridge.score(X_train_scaled, y_train)))

正則化ができたところで、正則化なしの重回帰分析とRidge回帰でどのような変化があったのかみてみます。偏回帰係数と決定係数、MSEを比較します。

# 偏回帰係数を比較
df_w_all = pd.concat([df_w, df_w_ridge], axis=1)
df_w_all

# MSEと決定係数を比較
data = {'訓練データMSE':[mse_train, ridge_mse_train], 
        'テストデータMSE':[mse_test, ridge_mse_test], 
        '決定係数':[model.score(X_test_scaled, y_test), ridge.score(X_test_scaled, y_test)]}
df_mse = pd.DataFrame(data=data, index=['重回帰', 'Ridge回帰'])
df_mse

比較の結果両者ともにほとんど差がないことがわかりました。重回帰分析の時に過学習していなかったからですかね、、、
今回の分析はこれで以上になります。パラメータなどを変更することでより良い精度の分析ができると思いますのでぜひ試してみてください。

内容に不備があった場合はご連絡いただけると幸いです。