pythonによる回帰分析〜リッジ回帰モデル

こんにちは、デジタルボーイです。今回は回帰モデルの中の1モデル、リッジ回帰モデルについて、Pythonとscikit-learnを使い、解説したいと思います!

リッジ回帰は、線形回帰の一種であり、特に多重共線性(特徴量同士が相関しすぎている状態)が問題になる場合に有効な手法です。本記事では、Scikit-learnを用いてカリフォルニアの住宅価格データを分析し、モデルの結果をグラフで可視化しながら解釈していきます。

ちなみにこれまで、単回帰分析、重回帰分析についても解説しています。リッジ回帰モデルはこれらのモデルの応用モデルと言えるので、必要があれば以下の解説を先におさらいしてみてください。

記事を書いた人

デジタルボーイです。
データサイエンス歴20年以上のおっさんです。中小企業診断士として、データサイエンス、WEBマーケティング、SEOに関するデータ分析、コンサルティングの仕事をしています。自己紹介の詳細はコチラ

目次

リッジ回帰とは?

リッジ回帰(Ridge regression)は、通常の線形回帰に L2正則化(ペナルティ)を加えたものです。通常の線形回帰では、説明変数が多すぎると過学習(学習データに過剰に適合してしまい、汎用性が低下する)を起こすことがあります。リッジ回帰では、回帰係数の大きさを抑えることで、モデルの安定性を向上させることができます。

リッジ回帰の直感的な理解

リッジ回帰の重要なポイントは 回帰係数の大きさを制限する ことです。例えば、通常の線形回帰では以下のような回帰直線が求められます。

$$
y = w_1 x_1 + w_2 x_2 + \dots + w_n x_n + b
$$

リッジ回帰では、これに 罰則項(L2正則化)を加え、極端に大きな回帰係数を抑制します。具体的には、以下のような形で損失関数を定義します。

\[
\text{損失関数} = \sum (y_i – \hat{y_i})^2 + \alpha \sum w_j^2
\]

ここで

  • 最初の項は通常の回帰の誤差(二乗誤差)
  • α(アルファ) は正則化の強さを決めるハイパーパラメータ
  • 罰則項(L2ノルム)によって、回帰係数 \( w_j​ \) が大きくなりすぎるのを防ぐ

αが大きいほど正則化が強くなり、係数は小さくなります。逆にαが0に近づくと通常の線形回帰に近づきます。

ちなみに、リッジ回帰と似たような回帰である、ラッソ回帰(Lasso回帰)は、L1正則化を適用した線形回帰モデルであり、以下の損失関数を最小化します。

\[
\text{損失関数} = \sum (y_i – \hat{y}_i)^2 + \alpha \sum |w_j|
\]

ラッソ回帰は、L1正則化項を加えることで、不要な特徴量の係数をゼロにし、モデルのスパース性を向上させます。 過学習を防ぎつつ、重要な特徴のみを選択するため、高次元データの特徴選択に有効です。 リッジ回帰との違いは、ラッソ回帰は一部の係数をゼロにするのに対し、リッジ回帰は係数を小さく抑える点です。

現実での利用例

回帰分析において複数の特徴量が強い相関を持つ状態のことを多重共線性と言います。これが発生すると、回帰係数が不安定になり、予測の解釈や精度に悪影響を及ぼします。リッジ回帰は、以下のような多重共線性のあるデータ場面で活用されています。

  • 経済指標の分析(GDP、失業率、物価指数などが相関している) 特徴量が多いデータセット
  • テキストデータのベクトル化(TF-IDFなど)や画像のピクセルデータ ノイズの影響を受けやすいデータ
  • 株価予測、医療データの予測(ノイズを抑えつつ全体の傾向を掴む)

必要なライブラリのインストール

以下のコマンドで必要なライブラリをインストールできます。

pip install scikit-learn matplotlib pandas numpy

pipを使ったインストール方法はこちらになります。

まずは、必要なライブラリをインポートしましょう。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error
from sklearn.datasets import load_diabetes

データの概要(糖尿病データセット)

サンプルデータとして、sklearn.datasets に含まれている 糖尿病データセット(diabetes dataset) を使用します。
このデータセットは、10個の特徴量(年齢、BMI、血圧など)を使って糖尿病の進行度を予測するためのデータです。
load_diabetes() を用いてデータを読み込み、概要を確認しましょう。

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

# データの特徴量とターゲット変数
X = data.data
y = data.target

# データフレーム化
df = pd.DataFrame(X, columns=data.feature_names)
df["target"] = y

# 先頭5行を表示
print(df.head())

コード解説

  • データの読み込み:load_diabetes() を使って、糖尿病データセットを取得。
  • 特徴量とターゲット変数の取得data.data には特徴量(説明変数)が、data.target には予測対象(目的変数)が格納されている。
  • データフレーム化pd.DataFrame() を使って、特徴量をPandasのデータフレームに変換。
    列名には data.feature_names を設定し、最後に target 列を追加。
  • データの先頭5行を表示df.head() でデータの最初の5行を確認し、データの概要を把握。

アウトプットはこんな感じです。

このデータセットには 10個の特徴量 があり、それぞれが糖尿病の進行度に関連しています。また、target 列が予測対象(目的変数)になります。

それぞれのデータの意味は以下となります。

特徴量の説明

特徴量説明
age年齢(標準化済み)
sex性別(標準化済み、具体的な値の意味は明示されていない)
bmiBMI(体格指数) = 体重(kg) / 身長(m)^2
bp平均血圧(標準化済み)
s1血清総コレステロール(TC: Total Cholesterol)
s2低密度リポタンパクコレステロール(LDL: Low-Density Lipoprotein Cholesterol)
s3高密度リポタンパクコレステロール(HDL: High-Density Lipoprotein Cholesterol)
s4TCH(総コレステロール / HDL コレステロール比)
s5血清トリグリセリドレベル(Serum Triglycerides)
s6血糖値(Blood Sugar Level)

分析の目的とゴール

分析の目的とゴールは以下とします。

  • この分析では、糖尿病データセットを用いて 患者の糖尿病進行度を予測 することを目的とする。
  • モデルの評価には、MSE(平均二乗誤差)、RMSE(ルート平均二乗誤差)、決定係数 \( R^2 \) を使用し、
    予測の精度を確認する。
  • 通常の線形回帰モデルと リッジ回帰(Ridge Regression) を比較し、
    正則化がモデルの性能に与える影響を評価することをゴールとする。

モデル実装

以下のコードでモデルを構築します。

# 必要なライブラリをインポート
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import Ridge
from sklearn.metrics import mean_squared_error, r2_score
import numpy as np

# データの準備
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# 特徴量の標準化(リッジ回帰ではスケーリングが重要)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# リッジ回帰モデルの学習(α=1.0)
ridge = Ridge(alpha=1.0)
ridge.fit(X_train_scaled, y_train)

# 予測
y_pred = ridge.predict(X_test_scaled)

コードの解説

  • train_test_split() を使用し、データを 80%(学習用)・20%(テスト用) に分割。
  • StandardScaler()特徴量を標準化(平均0、標準偏差1に変換)。
    • リッジ回帰は特徴量のスケールに敏感なため、標準化が重要。
  • Ridge(alpha=1.0) を使用してリッジ回帰モデルを作成し、学習データで fit() させる。
  • predict() を用いて、テストデータの予測値を算出。

モデル評価

では、構築したモデルを評価してみましょう。

# MSE(平均二乗誤差)
mse = mean_squared_error(y_test, y_pred)

# RMSE(ルート平均二乗誤差)
rmse = np.sqrt(mse)

# R²(決定係数)
r2 = r2_score(y_test, y_pred)

# 結果の表示
print(f"MSE: {mse:.2f}")
print(f"RMSE: {rmse:.2f}")
print(f"R²: {r2:.2f}")

コードの解説

  • mean_squared_error()MSE(誤差の2乗の平均) を算出。
  • np.sqrt(mse)RMSE(MSEの平方根) を計算。
    • 予測誤差の平均的な大きさを直感的に理解しやすくするため。
  • r2_score()決定係数 \( R^2 \) を算出。
    • 1.0に近いほど良いモデル(0.0ならランダムな予測と同じ)。
  • print() で結果を表示し、モデルの予測精度を確認。

アウトプットはこんな感じでした。

モデルの評価の解説

  • MSE(Mean Squared Error, 平均二乗誤差)
    • 誤差の2乗の平均。値が小さいほど良いモデル。
    • ただし、単位が元のデータの単位の 二乗 なので、直感的に解釈しにくい。
  • RMSE(Root Mean Squared Error, ルート平均二乗誤差)
    • MSEの平方根を取った値。
    • 単位が元のデータと同じになるため、誤差の大きさを直感的に把握しやすい。
  • 決定係数 \( R^2 \) (Coefficient of Determination)
    • 予測値がどれだけ実際の値を説明できるかを表す。
    • 1に近いほど良い(1.0: 完全予測、0.0: ランダム予測、負: モデルが悪すぎる)。

今回のモデルの決定係数は0.45でした。現実場面ではこの数値だと、「使い物にならないくらい精度の低い」モデルですかね(汗)。今回はモデルの紹介が目的なので、こんなもんにしときましょう(焦)

アウトプットのグラフ化

予測結果を 実際の値(正解データ)と比較 し、どの程度正確に予測できているかを可視化します。さらに、 相関係数(ピアソンの相関) を計算し、予測値と実測値の関連性を数値で示しましょう。

import matplotlib.pyplot as plt
import numpy as np

# フォント設定(Macの場合)
plt.rcParams["font.family"] = "Hiragino sans"

# 実測値 vs 予測値の相関係数を計算
correlation = np.corrcoef(y_test, y_pred)[0, 1]

# 実測値 vs 予測値の散布図
plt.figure(figsize=(8, 6))
plt.scatter(y_test, y_pred, alpha=0.7, edgecolors="k")
plt.plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], "r--")  
plt.xlabel("実測値(y_test)")
plt.ylabel("予測値(y_pred)")
plt.title(f"実測値 vs 予測値(相関係数: {correlation:.2f})")
plt.grid(True)
plt.show()

コードの解説

  • np.corrcoef(y_test, y_pred)[0, 1]実測値と予測値の相関係数 を計算し、どの程度の相関があるか確認。
  • scatter(y_test, y_pred)実測値と予測値の対応関係を散布図でプロット
  • plot([y_test.min(), y_test.max()], [y_test.min(), y_test.max()], "r--")基準線 を描画。
  • plt.rcParams["font.family"] = "Hiragino sans"Mac用の日本語フォントを指定
    • Windowsの場合は plt.rcParams["font.family"] = "Meiryo" などに変更

結果はこんな感じでした。

相関係数が0.68ということで、概ね、右肩上がりの関係性で、「実測値が上がれば予測値も直線的に上がっている」がとなっていますね!

モデルの回帰係数の確認

リッジ回帰では、通常の線形回帰と異なり、係数(重み)が 正則化によって抑制 されます。ここでは、モデルが学習した 回帰係数(weight, β) を可視化し、各特徴量の影響度を確認します。

# 各特徴量の回帰係数を取得
ridge_coefficients = ridge.coef_

# 特徴量名と回帰係数をデータフレーム化
coef_df = pd.DataFrame({"特徴量": data.feature_names, "回帰係数": ridge_coefficients})

# 回帰係数を降順に並べる
coef_df = coef_df.sort_values(by="回帰係数", ascending=False)

# グラフ化
plt.figure(figsize=(10, 6))
plt.barh(coef_df["特徴量"], coef_df["回帰係数"], color="skyblue")
plt.xlabel("回帰係数の値")
plt.ylabel("特徴量")
plt.title("リッジ回帰の回帰係数")
plt.grid(True)
plt.show()

コードの解説

  • ridge.coef_学習済みモデルの回帰係数 を取得。
  • pandas.DataFrame()特徴量と回帰係数を対応付けた表を作成
  • sort_values()影響度が大きい特徴量順に並び替え
  • barh()特徴量ごとの回帰係数を横棒グラフで可視化
    • 係数が大きいほど、モデルに与える影響が強い特徴量
    • 正の係数はターゲット変数を増加させる方向に働き、負の係数は減少させる方向に働く

アウトプットはこんな感じでした。

回帰係数のグラフを見ると、リッジ回帰がどの特徴量を重要視しているかがわかります。

  • 係数が大きい特徴量 → 予測に大きく寄与している
  • 係数が小さい(0に近い)特徴量 → 正則化の影響で影響が抑制されている

S1(血清総コレステロール)が正の方向に、S5(血清トリグリセリドレベル)が負の方向に、それぞれモデルの説明力に大きく貢献していることがわかります。つまり、血清総コレステロールが上がると糖尿病進行度が上がり、血清トリグリセリドレベルが上がると糖尿病進行度が下がる、ということがモデルから言えます。

ちなみに、通常、総コレステロール(S1)が高いことは動脈硬化や心血管疾患のリスクを高める ため、糖尿病の進行とも関連している可能性がます。また、トリグリセリド(S5)が高いことは糖尿病のリスクを上げると考えられており、この場合、「逆じゃねーか!」と思われます。ここら辺は僕も、医学の専門家でないので、全くわかりませんが、回帰係数は、「他の変数が一定だった場合に、当該変数が変化した場合の予測値の変化量」なので、通常、トリグリセリドが上がれば(両者に相関があるため)コレステロール値も上がるが、モデルではコレステロール値などの変数が一定だった場合のトリグリセリドの上昇に対する糖尿病リスクであり、その場合は糖尿病リスクは低いと言える、ということなのかもしれません。

リッジ回帰の特性について

以下ではリッジ回帰の特性について、実データから簡単に実験しながら見ていきましょう

通常の線形回帰との比較

リッジ回帰と通常の 線形回帰(最小二乗法) の違いを確認するために、両者を比較します。
通常の線形回帰では、特徴量のスケールが異なると 過学習しやすく、回帰係数が極端に大きくなる ことがあります。
リッジ回帰では、正則化によって 回帰係数が制御され、より安定した予測が可能 になります。

from sklearn.linear_model import LinearRegression

# 線形回帰モデルの学習
lin_reg = LinearRegression()
lin_reg.fit(X_train_scaled, y_train)

# 線形回帰の予測
y_pred_lin = lin_reg.predict(X_test_scaled)

# MSE, RMSE, R2の比較
mse_lin = mean_squared_error(y_test, y_pred_lin)
rmse_lin = np.sqrt(mse_lin)
r2_lin = r2_score(y_test, y_pred_lin)

mse_ridge = mean_squared_error(y_test, y_pred)
rmse_ridge = np.sqrt(mse_ridge)
r2_ridge = r2_score(y_test, y_pred)

# 結果の表示
print(f"線形回帰 - MSE: {mse_lin:.2f}, RMSE: {rmse_lin:.2f}, R²: {r2_lin:.2f}")
print(f"リッジ回帰 - MSE: {mse_ridge:.2f}, RMSE: {rmse_ridge:.2f}, R²: {r2_ridge:.2f}")

コードの解説

  • LinearRegression() を使い、通常の線形回帰モデルを学習。
  • predict() を用いて、テストデータの予測値を取得。
  • MSE, RMSE, R²を比較 し、リッジ回帰のほうが過学習しにくいことを確認。

結果はこうでした。

ぬぬぬ。。。想定ではリッジ回帰の方が良くなるところをお見せしたかったんですが、この程度の精度だったら、どっちもどっちですね。。。もっとたくさんの説明変数がある場合は、結構リッジ回帰の方が精度が高くなることがあります。

α(正則化の強さ)を変化させた場合の影響の可視化

リッジ回帰の α(正則化パラメータ) を変えると、回帰係数の大きさや予測精度が変化します。αが 小さい ほど通常の線形回帰に近づき、αが 大きい ほど正則化が強くなり、回帰係数が抑制されます。

# フォント設定(Macの場合)
plt.rcParams["font.family"] = "Hiragino sans"

# αの値を変えてリッジ回帰を実行
alpha_values = [0.001, 0.01, 0.1, 1.0, 10, 100]
mse_values = []

for alpha in alpha_values:
    ridge_model = Ridge(alpha=alpha)
    ridge_model.fit(X_train_scaled, y_train)
    y_pred_ridge = ridge_model.predict(X_test_scaled)
    mse_values.append(mean_squared_error(y_test, y_pred_ridge))

# グラフ化
plt.figure(figsize=(8, 6))
plt.plot(alpha_values, mse_values, marker="o", linestyle="-")
plt.xscale("log")
plt.xlabel("正則化パラメータ α(対数スケール)")
plt.ylabel("MSE(平均二乗誤差)")
plt.title("α の変化によるMSEの変動")
plt.grid(True)
plt.show()

コードの解説

  • alpha_values に異なる αの値(0.001 〜 100) を設定。
  • 各αに対してリッジ回帰を学習し、MSE(平均二乗誤差) を計算。
  • x軸を対数スケール にして αの増加とMSEの関係 を可視化。
    • αが 小さすぎると過学習 し、MSEが高くなる。
    • αが 大きすぎると、過剰に正則化されてしまい、精度が低下 する。

結果はこんな感じです。

グラフでは100 (=1.0)のところでガクッと下がっていますね。実際の分析場面でもこのようにざっくりとグラフを見ながら、現実場面では多くの場合、1未満で調整するといいでしょう。

特徴量ごとの係数の変化の確認

リッジ回帰では αが増えると、回帰係数の大きさが抑制される ことが重要なポイントです。以下のコードで、 αを変化させたときの特徴量ごとの回帰係数の変動 を可視化します。

# フォント設定(Macの場合)
plt.rcParams["font.family"] = "Hiragino sans"

# αの値リスト
alpha_values = [0.001, 0.01, 0.1, 1.0, 10, 100]
coefs = []

# 各αで学習し、回帰係数を取得
for alpha in alpha_values:
    ridge_model = Ridge(alpha=alpha)
    ridge_model.fit(X_train_scaled, y_train)
    coefs.append(ridge_model.coef_)

# グラフ化
plt.figure(figsize=(10, 6))
for i in range(len(data.feature_names)):
    plt.plot(alpha_values, [coef[i] for coef in coefs], label=data.feature_names[i])

plt.xscale("log")
plt.xlabel("正則化パラメータ α(対数スケール)")
plt.ylabel("回帰係数の値")
plt.title("α の変化による回帰係数の変動")
plt.legend(loc="best", fontsize=10)
plt.grid(True)
plt.show()

コードの解説

  • 異なるαの値 でリッジ回帰を学習し、各特徴量の回帰係数 を取得。
  • x軸をαの値(対数スケール)y軸を回帰係数の値 としてプロット。
  • αが増加すると回帰係数の絶対値が小さくなる(特徴量の影響が均等化される)。

結果はこんな感じです。

αを変化させた時、回帰係数の大きな特徴量ほど、ガクッと0に近くなっていますね。この正則化の影響で、大きな回帰係数ほど強くペナルティを受け、αの増加に伴い急激に0に近づく ことが見て取れます。この抑制の仕組みによって、過学習を防ぐ効果が期待できるんですね!

まとめ

以上が、リッジ回帰の解説でした。リッジ回帰を使うことで、特徴量の影響を抑制する仕組みを、コードを用いて実感できたと思います!

目次