Zenn
🧁

【Python】scikit-learnのdiabetes datasetを使って回帰分析をしてみる

2023/07/05に公開
Python
scikit-learn
機械学習
tech


scikit-learnサイトにはToy dataと言われる様々な公開データがあります。今回はその中の「糖尿病患者の診療データセット」を使って単回帰分析や重回帰分析をしてみます。

目次[clickで開く/閉じる]
  • データの概要
  • データの中身の確認
  • 単回帰分析の流れ
  • 単回帰分析のソースコード
  • 詳細
  • 重回帰分析の流れ
  • 重回帰分析のソースコード
  • 詳細
  • まとめ

scikit-learnとは、誰もが無料で使えるPythonのオープンソース機械学習ライブラリです。今回はその中でも「diabetes dataset」(糖尿病患者の診療データセット)を用いて回帰分析に取り組みます。データセットに含まれる様々な変数から住宅価格を予測することが最終目標です。
　では、実際のデータの中身についてみていきましょう。

  • データ数：442
  • カラム数：11
変数名 詳細
age 年齢
sex 性別
bmi BMI
bp 平均血圧
s1 tc : 総血清コレステロール
s2 ldl : 密度リボタンパク質
s3 hdl : 高密度リボタンパク質
s4 tch : 総コレステロール/HDL
s5 ltg : 結成トリグリセリドベルの対数
s6 glu : 血糖値
(target) （目的変数）

説明変数が10列で目的変数のtargetが1列あります。

各データの詳細については以下のscikit-learnのページもご参考ください。
https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_diabetes.html

データの中身の確認

まずは必要ライブラリをまとめてインポートしておきます。

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_diabetes

インストールしたデータを見てみます。この時についでにデータの大きさやデータの型、欠損値の有無についても確認しておきましょう。
変数名をわかりやすくするために変数名s1~s6を具体的な名前に変更しておきます。

# データの読み込み
dia = load_diabetes()

# 目的変数
exp_data = pd.DataFrame(dia.data, columns=dia.feature_names)
# 説明変数
tar_data = pd.DataFrame(dia.target, columns=['target'])
# データを結合
df = pd.concat([exp_data, tar_data], axis=1)

df.rename({"s1":"tc", "s2":"ldl", "s3":"hdl", "s4":"tch", "s5":"ltg", "s6":"glu"}, axis=1, inplace=True)
df.dropna(inplace = True)
print(df.shape)
display(df.head())
print(df.dtypes)
print(df.isnull().sum())

では次にデータ同士の相関関係を確認するためんいヒートマップを作成します。

# ヒートマップで相関係数を確認を表示
corr = df.corr()
fig = plt.figure(figsize=(6, 5))
sns.heatmap(corr, cmap="coolwarm", vmin=-1, vmax=1, annot=True)

作成したヒートマップがこちらです。

目的変数のtargetに注目してヒートマップを見てみるとbmiとの相関が一番大きいようなので単回帰分析はbmi, targetで行います。

単回帰分析の流れ

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

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

今回は説明変数をbmi、目的変数をtargetとして進めていきます。以下が単回帰分析のソースコードになります。

# 目的変数と説明変数に分割
X = df[['bmi']]
y = df[['target']]

# 学習
from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X, y)

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

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

# 回帰直線と散布図を表示
plt.figure(figsize=(5, 5))
plt.scatter(df['bmi'], df['target'])
plt.plot(X, model.predict(X), color='red')
plt.xlabel('bmi')
plt.ylabel('target')
plt.title('bmi vs target (Linear Regression)')

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

exp_var = 'bmi' # 説明変数
tar_var = 'target' # 目的変数

データの内容の確認

2変数に絞って散布図と統計量を見てみます。

# bmiとtargetの散布図を表示
fig = plt.figure(figsize=(5, 5))
plt.scatter(df["bmi"], df["target"])
plt.xlabel('bmi')
plt.ylabel('target')
plt.title('bmi vs target')
plt.savefig('bmi_vs_target.png')

# 統計量を表示
df[["bmi", "target"]].describe()

散布図を見てみると正の相関があるようです。念の為外れ値の除去を行います。

外れ値の除去

# 外れ値を削除
q_95 = df['bmi'].quantile(0.95)
print('95%点の分位数', q_95)
df = df[df['bmi'] < q_95]

データの内容の確認

# 散布図を表示
plt.figure(figsize=(5, 5))
plt.scatter(df['bmi'], df['target'])
plt.xlabel('bmi')
plt.ylabel('target')
plt.title('bmi vs target2')
plt.savefig('bmi cs target2.png')
plt.show()

# 記述統計量を確認
df[['bmi', 'target']].describe()

外れ値を除去したところ、bmiの標準偏差がほとんどなくなったため今後の分析もこのデータを使用して行います。それでは単回帰分析に取り掛かります。

説明変数と目的変数に分割する

# 目的変数と説明変数に分割
X = df[['bmi']]
y = df[['target']]

単回帰モデルを学習する

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

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(X, y)

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

精度の確認

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

結果はこちらになります。

決定係数が約0.2586であることから予測精度は悪いことがわかります。やはり1変数だけでは予測が難しそうです。

回帰直線とデータを図示する

# 回帰直線と散布図を表示
plt.figure(figsize=(5, 5))
plt.scatter(df['bmi'], df['target'])
plt.plot(X, model.predict(X), color='red')
plt.xlabel('bmi')
plt.ylabel('target')
plt.title('bmi vs target (Linear Regression)')
plt.savefig('bmi vs target (Linear Regression).png')

散布図と回帰直線を見比べてみると、回帰直線が若干下の方にずれているように思います。
では、続いて重回帰分析に取り組んでみます。単回帰分析の時と同様に分析の流れについて説明し、その後ソースコードを指名sてからコードの詳細に触れていきます。
今回は目的変数がtarget、説明変数がその他の変数になっています。

重回帰分析の流れ

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

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

必要に応じてコメントアウトしてください。

# 説明変数
exp_vars = ['bmi', 'bp', 'tc', 'hdl', 'tch', 'ltg', 'glu']
# 目的変数
tar_var = 'target'

# 外れ値を除去
for exp_var in exp_vars:
    q_95 = df[exp_var].quantile(0.95)
df = df[df[exp_var] < q_95]

# 訓練データとテストデータに分割
X = df[exp_vars]
y = df[[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)
print('X_train', X_train.shape)
print('X_test', X_test.shape)
print('y_train', y_train.shape)
print('y_test', y_test.shape)

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

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

# 予測値を計算
y_pred = model.predict(X_train_scaled)

# 偏回帰係数を計算
for xi, wi in zip(exp_vars, model.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, float(wi)))
    
# 決定係数の確認
r2 = model.score(X_train_scaled, y_train)
print('決定係数', r2)

# MSE
from sklearn.metrics import 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)
# 標準化
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)
X_train_scaled = pd.DataFrame(X_train_scaled, columns=exp_vars)
# 学習
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'])
for xi, wi in zip(exp_vars, ridge.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, wi))
print('L2ノルム: {:.3f}'.format(np.linalg.norm(ridge.coef_)))
# 訓練データに対するMSE
ridge_mse_train = mean_squared_error(y_train, ridge_y_pred)
print("訓練データに対するMSE: {:.3f}".format(ridge_mse_train))
# テストデータに対する MSE
X_test_scaled = scaler.transform(X_test)
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))
# 決定係数
ridge_r2 = ridge.score(X_train_scaled, y_train)
print('決定係数: {:.3f}'.format(ridge_r2))


# Lasso回帰
from sklearn.linear_model import Lasso
lasso = Lasso(alpha=1.0)
lasso.fit(X_train_scaled, y_train)
lasso_y_pred = lasso.predict(X_train_scaled)
# 偏回帰係数の確認
lasso_w = pd.Series(index=exp_vars, data=lasso.coef_)
for xi, wi in zip(exp_vars, lasso.coef_): # coef_は回帰直線の傾き
    # print(f'{xi, wi}')
    print('{0:7s}: {1:6.3f}'.format(xi, wi))
print('L2ノルム', np.linalg.norm(lasso_w))
# MSE
lasso_mse_train = mean_squared_error(y_train, lasso_y_pred)
print("訓練データに対するMSE: {:.3f}".format(lasso_mse_train))
lasso_X_test_scaled = scaler.transform(X_test)
lasso_y_pred_test = lasso.predict(lasso_X_test_scaled)
lasso_mse_test = mean_squared_error(y_test, lasso_y_pred_test)
print("テストデータに対するMSE: {:.3f}".format(lasso_mse_test))
# 決定係数
lasso_r2 = lasso.score(X_train_scaled, y_train)
print('決定係数: {:.3f}'.format(lasso_r2))
result_data = {'訓練データMSE':[mse_train, ridge_mse_train, lasso_mse_train], 
        'テストデータMSE':[mse_test, ridge_mse_test, lasso_mse_test],
        '決定係数':[r2, ridge_r2, lasso_r2]}
df_result = pd.DataFrame(data = result_data, index=['通常の重回帰', 'Ridge回帰', 'Lasso回帰'])
display(df_result)

# 特徴量選択
# ラッパー法
import statsmodels.api as sm
# 切片追加
X_ = sm.add_constant(X)
# 重回帰学習
lr = sm.OLS(y,X_).fit()
# 各係数のp値
lr.pvalues
# 関数定義（変数減少法）
## 引数colsは特徴量、引数cは基準
def backward_elimination(cols,c):
    
    while (len(cols)>0):
        
        # p値を格納するハコ
        p = []
        
        # モデル構築（重回帰）
        X_ = sm.add_constant(X_train[cols])
        lr = sm.OLS(y_train,X_).fit()
        
        # p値の抽出
        p = pd.Series(lr.pvalues.values[1:],index = cols)      
        
        # ｐ値が最大の特徴量がc以上の場合に除外
        if(max(p)>c):  
            cols.remove(p.idxmax())
        else:
            break
    return cols
# 特徴量選択（変数選択）の実施
X_selected = backward_elimination(cols=list(X.columns),c=0.05)
print(X_selected)
# 学習
model.fit(X_train[X_selected], y_train)
# MSE
rapper_mse_train = mean_squared_error(y_train, model.predict(X_train[X_selected]))
print('MSE(train)', rapper_mse_train)
rapper_mse_test = mean_squared_error(y_test, model.predict(X_test[X_selected]))
print('MSE(test)', rapper_mse_test)
# 決定係数の確認
# rapper_r2_train = model.score(X_train[X_selected], y_train)
# print('決定係数(train)', rapper_r2_train)
rapper_r2_test = model.score(X_test[X_selected], y_test)
print('決定係数(test)', rapper_r2_test)
# 偏回帰係数を計算
for xi, wi in zip(X_selected, model.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, float(wi)))
    
# フィルター法
# 相関係数（学習データ）
cor = pd.concat([X_train, y_train], axis=1).corr()
# 目的変数との相関係数の絶対値
target_cor = abs(cor['target'])
print(target_cor)
# 基準
c = 0.2
# 選択の実施
X_selected = target_cor[target_cor>c]
X_selected = X_selected.drop('target').index
# 選択した特徴量（説明変数）
print(X_selected)
# 選択した説明変数同士の相関係数
display(X_train[X_selected].corr())
# 学習（学習データ利用）
model.fit(X_train[X_selected], y_train)
# MSE
filter_mse_train = mean_squared_error(y_train, model.predict(X_train[X_selected]))
print('MSE（学習データ）:', filter_mse_train)
filter_mse_test = mean_squared_error(y_test, model.predict(X_test[X_selected]))
print('MSE（テストデータ）:', filter_mse_test)
# 精度検証（決定係数R2）
# filter_r2_train = model.score(X_train[X_selected], y_train)
# print('決定係数R2（学習データ）:', filter_r2_train)
filter_r2_test = model.score(X_test[X_selected], y_test)
print('決定係数R2（テストデータ）:', filter_r2_test)
# 偏回帰係数を計算
for xi, wi in zip(X_selected, model.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, float(wi)))
    
# RFE法
from sklearn.feature_selection import RFE
rfe = RFE(LinearRegression())
# 特徴量選定（変数選択）の実施
X_rfe = rfe.fit(X_train,y_train)  
X_selected = X.columns[X_rfe.support_]
# 選択した特徴量
print(X_selected)
# 学習（学習データ利用）
model.fit(X_train[X_selected], y_train)
# MSE
rfe_mse_train = mean_squared_error(y_train, model.predict(X_train[X_selected]))
print('MSE（学習データ）:', rfe_mse_train)
rfe_mse_test = mean_squared_error(y_test, model.predict(X_test[X_selected]))
print('MSE（テストデータ）:', rfe_mse_test)
# 精度検証（決定係数R2）
rfe_r2_train = model.score(X_train[X_selected], y_train)
print('決定係数R2（学習データ）:', rfe_r2_train)
rfe_r2_test = model.score(X_test[X_selected], y_test)
print('決定係数R2（テストデータ）:', rfe_r2_test)
# 偏回帰係数を計算
for xi, wi in zip(X_selected, model.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, float(wi)))
result_data2 = {'訓練データMSE':[mse_train, ridge_mse_train, lasso_mse_train, rapper_mse_train, filter_mse_train, rfe_mse_train], 
        'テストデータMSE':[mse_test, ridge_mse_test, lasso_mse_test, rapper_mse_test, filter_mse_test, rfe_mse_test],
        '決定係数':[r2, ridge_r2, lasso_r2, rapper_r2_test, filter_r2_test, rfe_r2_test]}
df_result2 = pd.DataFrame(data = result_data2, index=['通常の重回帰', 'Ridge回帰', 'Lasso回帰', 'ラッパー法', 'フィルタ法', 'RFE法'])
display(df_result2)

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

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

目的変数をwxp_vars, 説明変数をtar_varとし、データを分割します。

# 説明変数
exp_vars = ['bmi', 'bp', 'tc', 'hdl', 'tch', 'ltg', 'glu']
# 目的変数
tar_var = 'target'

先ほどと同様に外れ値を除去し、記述統計量を見てみます。

外れ値の除去

# 外れ値を除去
for exp_var in exp_vars:
    q_95 = df[exp_var].quantile(0.95)
df = df[df[exp_var] < q_95]
display(df.describe())

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

学習用に訓練データとテストデータを8:2に分割します。

# 訓練データとテストデータに分割
X = df[exp_vars]
y = df[[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)
print('X_train', X_train.shape)
print('X_test', X_test.shape)
print('y_train', y_train.shape)
print('y_test', y_test.shape)

標準化

標準化をすることで異なる単位の変数を同一尺度で表せるようになります。

# 標準化
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
scaler.fit(X_train)
X_train_scaled = scaler.transform(X_train)

X_train_scaled = pd.DataFrame(X_train_scaled, columns=exp_vars)

モデルの学習

モデルは単回帰分析の時と同様にLinearRegressionを使用します。

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

# 予測値を計算
y_pred = model.predict(X_train_scaled)

偏回帰係数を計算

# 偏回帰係数を計算
for xi, wi in zip(exp_vars, model.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, float(wi)))

学習と予測ができたところでどのくらいの精度稼働かをみるため、モデルの評価を行います。評価尺度には決定係数とMSEを使用します。

モデルの評価

# 決定係数の確認
r2 = model.score(X_train_scaled, y_train)
print('決定係数', r2)

from sklearn.metrics import 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.4958095621701042
  • 訓練データに対するMSE：2775.980
  • テストデータに対するMSE：3749.166

決定係数は約0.5となりあまり良くない結果になりました。また、MSEについても訓練データとテストデータで値に差があるためやや過学習ぎみと言えるでしょう。そのため正則化を行って精度の比較を行います。

正則化

正則化とは過学習を抑えるための仕組みのひとつで、これを行うことで学習データにモデルが過剰に適合してしまうことで未知のデータに対する予測能力が低下するのを防ぐ役割を持っています。回帰分析における目的は実データに対する予測を行うのであり既知のデータを正確に予測できるだけでは意味がありません。
今回は正則化としてL1ノルムとL2ノルムを用いてLasso回帰とRidge回帰を行います。Ridge回帰はパラメータの値が0となり疎なベクトルになる箇所が多くなるという特徴を持ち、Lasso回帰はパラメータ（今回は偏回帰係数）が小さくなるように学習されるという特徴を持ちます。それでは実際に正則化をやってみます。

まずはRidge回帰から、

# Ridge回帰
from sklearn.linear_model import Ridge
ridge = Ridge(alpha=1.0)

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

# 学習
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'])
for xi, wi in zip(exp_vars, ridge.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, wi))
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
X_test_scaled = scaler.transform(X_test)
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))

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

次にLasso回帰を、

# Lasso回帰
from sklearn.linear_model import Lasso
lasso = Lasso(alpha=1.0)
lasso.fit(X_train_scaled, y_train)
lasso_y_pred = lasso.predict(X_train_scaled)

# 偏回帰係数の確認
lasso_w = pd.Series(index=exp_vars, data=lasso.coef_)
for xi, wi in zip(exp_vars, lasso.coef_): # coef_は回帰直線の傾き
    # print(f'{xi, wi}')
    print('{0:7s}: {1:6.3f}'.format(xi, wi))
print('L2ノルム', np.linalg.norm(lasso_w))

lasso_mse_train = mean_squared_error(y_train, lasso_y_pred)
print("訓練データに対するMSE: {:.3f}".format(lasso_mse_train))

lasso_X_test_scaled = scaler.transform(X_test)
lasso_y_pred_test = lasso.predict(lasso_X_test_scaled)
lasso_mse_test = mean_squared_error(y_test, lasso_y_pred_test)
print("テストデータに対するMSE: {:.3f}".format(lasso_mse_test))

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

2種類の正則化を行ったところで、それぞれの精度の確認を行います。
今回は比較しやすいよう

  • 訓練データのMSE
  • テストデータのMSE
  • 決定係数
    の3項目でDataFrameを作成して比較します。
result_data = {'訓練データMSE':[mse_train, ridge_mse_train, lasso_mse_train], 
        'テストデータMSE':[mse_test, ridge_mse_test, lasso_mse_test],
        '決定係数':[r2, ridge_r2, lasso_r2]}
df_result = pd.DataFrame(data = result_data, index=['通常の重回帰', 'Ridge回帰', 'Lasso回帰'])
df_result

下記に実際に出力されたものを表示します。

MSE, 決定係数ともにほとんど変わりませんでした。正則化で精度が向上しなかったので、特徴量選択を工夫することで精度の向上を図りたいと思います。

特徴量選択

今回は特徴量選択として以下の3つに取り組みました。
※どのような原理なのかなどの詳細については触れていません。

  • ラッパー法：p-value >0.05の説明変数を削除する
  • フィルタ法：|相関係数|>cの説明変数を削除する
  • RFE法

それぞれの方法のコードを示し、最後に精度をまとめて確認します。

ラッパー法
# p値>0.05の説明変数を削除する（特徴量選択）
import statsmodels.api as sm

# 切片追加
X_ = sm.add_constant(X)
# 重回帰学習
lr = sm.OLS(y,X_).fit()
# 各係数のp値
lr.pvalues

# 関数定義（変数減少法）
## 引数colsは特徴量、引数cは基準
def backward_elimination(cols,c):
    
    while (len(cols)>0):
        
        # p値を格納するハコ
        p = []
        
        # モデル構築（重回帰）
        X_ = sm.add_constant(X_train[cols])
        lr = sm.OLS(y_train,X_).fit()
        
        # p値の抽出
        p = pd.Series(lr.pvalues.values[1:],index = cols)      
        
        # ｐ値が最大の特徴量がc以上の場合に除外
        if(max(p)>c):  
            cols.remove(p.idxmax())
        else:
            break
    return cols
# 特徴量選択（変数選択）の実施
X_selected = backward_elimination(cols=list(X.columns),c=0.05)
print(X_selected)

# 学習
model.fit(X_train[X_selected], y_train)

# MSE
rapper_mse_train = mean_squared_error(y_train, model.predict(X_train[X_selected]))
print('MSE(train)', rapper_mse_train)
rapper_mse_test = mean_squared_error(y_test, model.predict(X_test[X_selected]))
print('MSE(test)', rapper_mse_test)

# 決定係数の確認
# rapper_r2_train = model.score(X_train[X_selected], y_train)
# print('決定係数(train)', rapper_r2_train)
rapper_r2_test = model.score(X_test[X_selected], y_test)
print('決定係数(test)', rapper_r2_test)

# 偏回帰係数を計算
for xi, wi in zip(X_selected, model.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, float(wi)))
フィルタ法
# 相関係数（学習データ）
cor = pd.concat([X_train, y_train], axis=1).corr()
# 目的変数との相関係数の絶対値
target_cor = abs(cor['target'])
print(target_cor)

# 基準
c = 0.2
# 選択の実施
X_selected = target_cor[target_cor>c]
X_selected = X_selected.drop('target').index
# 選択した特徴量（説明変数）
print(X_selected)

# 選択した説明変数同士の相関係数
display(X_train[X_selected].corr())

# 学習（学習データ利用）
model.fit(X_train[X_selected], y_train)

# MSE
filter_mse_train = mean_squared_error(y_train, model.predict(X_train[X_selected]))
print('MSE（学習データ）:', filter_mse_train)
filter_mse_test = mean_squared_error(y_test, model.predict(X_test[X_selected]))
print('MSE（テストデータ）:', filter_mse_test)

# 精度検証（決定係数R2）
# filter_r2_train = model.score(X_train[X_selected], y_train)
# print('決定係数R2（学習データ）:', filter_r2_train)
filter_r2_test = model.score(X_test[X_selected], y_test)
print('決定係数R2（テストデータ）:', filter_r2_test)

# 偏回帰係数を計算
for xi, wi in zip(X_selected, model.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, float(wi)))
RFE法
# インスタンス
from sklearn.feature_selection import RFE
rfe = RFE(LinearRegression())
# 特徴量選定（変数選択）の実施
X_rfe = rfe.fit(X_train,y_train)  
X_selected = X.columns[X_rfe.support_]
# 選択した特徴量
print(X_selected)

# 学習（学習データ利用）
model.fit(X_train[X_selected], y_train)

# MSE
rfe_mse_train = mean_squared_error(y_train, model.predict(X_train[X_selected]))
print('MSE（学習データ）:', rfe_mse_train)
rfe_mse_test = mean_squared_error(y_test, model.predict(X_test[X_selected]))
print('MSE（テストデータ）:', rfe_mse_test)

# 精度検証（決定係数R2）
rfe_r2_train = model.score(X_train[X_selected], y_train)
print('決定係数R2（学習データ）:', rfe_r2_train)
rfe_r2_test = model.score(X_test[X_selected], y_test)
print('決定係数R2（テストデータ）:', rfe_r2_test)

# 偏回帰係数を計算
for xi, wi in zip(X_selected, model.coef_[0]):
    print('{0:7s}: {1:6.3f}'.format(xi, float(wi)))

正則化の時と同様にDataFrameでまとめて精度の確認を行います。
※正則化の精度も一緒に比較します。

result_data2 = {'訓練データMSE':[mse_train, ridge_mse_train, lasso_mse_train, rapper_mse_train, filter_mse_train, rfe_mse_train], 
        'テストデータMSE':[mse_test, ridge_mse_test, lasso_mse_test, rapper_mse_test, filter_mse_test, rfe_mse_test],
        '決定係数':[r2, ridge_r2, lasso_r2, rapper_r2_test, filter_r2_test, rfe_r2_test]}
df_result2 = pd.DataFrame(data = result_data2, index=['通常の重回帰', 'Ridge回帰', 'Lasso回帰', 'ラッパー法', 'フィルタ法', 'RFE法'])
df_result2

結果がこちらになります。

みていただけるとわかるようにMSE, 決定係数ともに大きな変化はありませんでした。

今回はここで分析を終わりますが、より良い精度を出せるようにパラメータや変数に対して工夫が必要なようです。

最後までご覧いただきありがとうございました。
内容に誤りがあった場合は教えていただけるとありがたいです。
（2023/07/05）

Discussion