回帰分析フルパイプライン_欠損・外れ値・相関/VIF・エンコーディング・Lasso/RF選択・RF/XGB比較まで一気通貫（Python/ scikit-learn）

本記事では、回帰分析のフルパイプラインを構築します。
具体的には、欠損値の削除、補完外れ値の確認（IQR/3σ）、相関係数・VIFによる多重共線性の整理、カテゴリ変数のエンコーディング（One-Hot, Rare統合）、LassoCV / ランダムフォレストによる特徴量選択、ランダムフォレストとXGBoostのモデル比較 & ハイパーパラメータチューニング、クロスバリデーション＋ホールドアウトによる精度評価までを一気通貫で実装します。
可視化やログ出力を組み込み、再現性と解釈性の高い分析ワークフローを再現します。

使用するデータセットは、UCバークレー大学の「UCI Machine Leaning Repository」で公開されている、自動車価格のデータセットです。
データセットの構成は下記となっていて目的変数はPriceで回帰分析を想定したデータ分析を実施します。

ID	属性名	説明
0	symboling	シンボル（-3〜+3）。車の危険度を示す指標で、+3が最も危険度が高い。保険会社の視点で設定。
1	normalized-losses	正規化損失（65〜256）。保険会社が支払う平均損失額を正規化したもの。値が大きいほどリスク高。
2	make	メーカー（アルファロメロ、アウディ、BMW、…、フォルクスワーゲン、ボルボ）
3	fuel-type	燃料タイプ（diesel, gas）
4	aspiration	自然吸気方式（std, turbo）
5	num-of-doors	ドア数（two, four）
6	body-style	ボディスタイル（hardtop, wagon, sedan, hatchback, convertible）
7	drive-wheels	駆動方式（4wd, fwd, rwd）
8	engine-location	エンジン位置（front, rear）
9	wheel-base	ホイールベース（86.6〜120.9）
10	length	車長（141.1〜208.1）
11	width	車幅（60.3〜72.3）
12	height	車高（47.8〜59.8）
13	curb-weight	車両重量（1488〜4066）
14	engine-type	エンジンタイプ（dohc, ohc, ohcv, rotor など）
15	num-of-cylinders	気筒数（two, three, four, five, six, eight, twelve）
16	engine-size	エンジンサイズ（61〜326）
17	fuel-system	燃料システム（1bbl, 2bbl, 4bbl, mpfi, spfi, など）
18	bore	ボア径（2.54〜3.94）
19	stroke	ストローク長（2.07〜4.17）
20	compression-ratio	圧縮比（7〜23）
21	horsepower	馬力（48〜288）
22	peak-rpm	ピーク回転数（4150〜6600）
23	city-mpg	市街地燃費（13〜49 mpg）
24	highway-mpg	高速道路燃費（16〜54 mpg）
25	price	価格（5118〜45400）

本記事で実施する回帰分析は、以下のステップに沿って進めます。

データセットの読み込み
欠損値処理
外れ値処理
特徴量選択（エンコーディング前）
カテゴリ変数エンコーディング
特徴量選択（エンコード後）＆次元爆発対策
モデル候補とハイパーパラメータ探索
評価指標と検証
モデル解釈と可視化
モデル保存・運用

※今回のポイント
今回のデータセットはカテゴリ変数が多いため、ステップ4で外れ値の影響が大きい、または不要な特徴量を整理した後、ステップ6でカテゴリ変数をエンコードし、さらにステップ7で次元爆発対策として特徴量選択を再実施しています。
これにより、学習効率を高めつつ予測性能を維持するバランスの取れたモデル構築が可能となります。

1. データ読み込み

import warnings
warnings.filterwarnings('ignore')
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import OneHotEncoder
from sklearn.impute import SimpleImputer
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LassoCV
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
from statsmodels.stats.outliers_influence import variance_inflation_factor
import joblib
# 日本語フォントの設定
plt.rcParams['font.family'] = 'Meiryo'  # Windowsの場合
# Macの場合は 'Hiragino Sans' など
plt.rcParams['axes.unicode_minus'] = False # 負の符号の文字化けを防ぐ

# 「?」を欠損として読み込む
df = pd.read_csv('../dataset/Automobile/data.csv', sep=',', na_values='?')
print(df.info())

target = 'price'
# 目的変数の欠損行を削除
df.dropna(subset=[target], inplace=True)
df[target] = df[target].astype(float)

print(f"データの形状: {df.shape}")
print(f"列名: {list(df.columns)}")
print(f"目的変数 '{target}' の欠損件数: {df[target].isnull().sum()} → 欠損行を削除しました。")

実行ログ

RangeIndex: 205 entries, 0 to 204
Data columns (total 26 columns):

Column Non-Null Count Dtype

0 symboling 205 non-null int64
1 normalized_losses 164 non-null float64
2 make 205 non-null object
3 fuel_type 205 non-null object
4 aspiration 205 non-null object
5 num_of_doors 203 non-null object
6 body_style 205 non-null object
7 drive_wheels 205 non-null object
8 engine_location 205 non-null object
9 wheel_base 205 non-null float64
10 length 205 non-null float64
11 width 205 non-null float64
12 height 205 non-null float64
13 curb_weight 205 non-null int64
14 engine_type 205 non-null object
15 num_of_cylinders 205 non-null object
16 engine_size 205 non-null int64
17 fuel_system 205 non-null object
18 bore 201 non-null float64
19 stroke 201 non-null float64
20 compression_ratio 205 non-null float64
21 horsepower 203 non-null float64
22 peak_rpm 203 non-null float64
23 city_mpg 205 non-null int64
24 highway_mpg 205 non-null int64
25 price 201 non-null float64
dtypes: float64(11), int64(5), object(10)
memory usage: 41.8+ KB
None
データの形状: (201, 26)
列名: ['symboling', 'normalized_losses', 'make', 'fuel_type', 'aspiration', 'num_of_doors', 'body_style', 'drive_wheels', 'engine_location', 'wheel_base', 'length', 'width', 'height', 'curb_weight', 'engine_type', 'num_of_cylinders', 'engine_size', 'fuel_system', 'bore', 'stroke', 'compression_ratio', 'horsepower', 'peak_rpm', 'city_mpg', 'highway_mpg', 'price']
目的変数 'price' の欠損件数: 0 → 欠損行を削除しました。

2. 目的変数（ターゲット）の分布確認と変換

目的変数priceの分布をヒストグラムで確認し、歪度と尖度を計算します。分布の歪みが大きいため、対数変換を適用して正規性に近づけます。

正規性に近づける2つの主なメリット

モデルの性能向上
多くの回帰モデル（例：線形回帰）は、誤差項が正規分布に従うという前提に立っています。目的変数を正規分布に近づけることで、この前提が満たされやすくなり、モデルの予測性能が向上します。また、回帰係数の推定値が安定し、信頼性が高まります。
モデル解釈の改善
変換後のデータはより対称的な分布になるため、モデルの予測結果や特徴量重要度の解釈がしやすくなります。外れ値の影響が緩和され、データ内の関係性をより正確に捉えることができます。また、予測区間や信頼区間の計算も、正規性の前提に立つことで、より信頼できるものになります。

from scipy.stats import skew, kurtosis

# 変換前のヒストグラムと歪度・尖度
print(f"変換前の歪度: {skew(df[target])}")
print(f"変換前の尖度: {kurtosis(df[target])}")

plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
sns.histplot(df[target], kde=True, bins=30)
plt.title('変換前のpriceの分布')

# 対数変換の実行
df[target] = np.log1p(df[target])

# 変換後のヒストグラムと歪度・尖度
print(f"変換後の歪度: {skew(df[target])}")
print(f"変換後の尖度: {kurtosis(df[target])}")
plt.subplot(1, 2, 2)
sns.histplot(df[target], kde=True, bins=30)
plt.title('変換後のpriceの分布')
plt.tight_layout()
plt.show()

実行ログ

変換前の歪度: 1.7961422058437058
変換前の尖度: 3.1220053568033848
変換後の歪度: 0.6736525970922048
変換後の尖度: -0.20235110336701823

html_img_001.png)

3. 欠損値処理

normalized_lossesなど、一部の変数に存在する欠損値を補完します。数値変数には中央値、カテゴリ変数には最頻値を適用します。

import missingno as msno

# 欠損値の可視化
msno.bar(df)
plt.title('欠損値の棒グラフ (補完前)')
plt.show()

# 数値変数とカテゴリ変数の分類
numerical_cols = df.select_dtypes(include=np.number).columns.tolist()
categorical_cols = df.select_dtypes(include='object').columns.tolist()
numerical_cols.remove(target)

# 数値変数の欠損値を中央値で補完
for col in numerical_cols:
    if df[col].isnull().any():
        df[col].fillna(df[col].median(), inplace=True)

# カテゴリ変数の欠損値を最頻値で補完
for col in categorical_cols:
    if df[col].isnull().any():
        df[col].fillna(df[col].mode()[0], inplace=True)

# 補完後の確認
print("--- 欠損値処理後の確認 ---")
print(f"データの形状: {df.shape}")
print("各列の欠損値数:")
print(df.isnull().sum())
print("\n主要な数値列の補完前後ヒストグラムを比較")
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
sns.histplot(df['horsepower'], kde=True, bins=20)
plt.title('horsepower 補完後の分布')
plt.subplot(1, 2, 2)
sns.histplot(df['normalized_losses'], kde=True, bins=20)
plt.title('normalized_losses 補完後の分布')
plt.tight_layout()
plt.show()

実行ログ

--- 欠損値処理後の確認 --- データの形状: (201, 26) 各列の欠損値数: symboling 0 normalized_losses 0 make 0 fuel_type 0 aspiration 0 num_of_doors 0 body_style 0 drive_wheels 0 engine_location 0 wheel_base 0 length 0 width 0 height 0 curb_weight 0 engine_type 0 num_of_cylinders 0 engine_size 0 fuel_system 0 bore 0 stroke 0 compression_ratio 0 horsepower 0 peak_rpm 0 city_mpg 0 highway_mpg 0 price 0 dtype: int64
主要な数値列の補完前後ヒストグラムを比較

4. 外れ値処理

ツリー系モデルは外れ値に頑健なため、今回は外れ値の除去は行わず、そのまま学習データとして使用します。

ツリー系モデルとは、データを木（ツリー）のような構造で分割し、予測や分類を行う機械学習モデルの総称です。主に決定木、ランダムフォレスト、勾配ブースティング（XGBoostなど）がこれに該当します。

外れ値に頑健な理由
多くの統計モデル（特に線形回帰モデル）は、外れ値の影響を強く受けます。外れ値が一つ存在するだけで、回帰直線が大きく歪み、予測精度が著しく低下することがあります。
一方、ツリー系モデルはデータの分割に基づいて学習するため、外れ値の影響を受けにくいという特性があります。

分割による学習：決定木はデータを特定の閾値で分割していきます。外れ値は、その値が非常に大きくても小さくても、単一のリーフ（木の末端）ノードに隔離されることが多く、モデル全体の学習に与える影響が限定されます。

多数決による安定性：ランダムフォレストや勾配ブースティングは、複数の決定木を組み合わせるアンサンブル学習です。個々の木が外れ値の影響を受けても、全体の予測は多数の木の予測の平均や加重平均で決まるため、結果として頑健な予測が可能になります。

この頑健性から、前処理の段階で外れ値を除去する手間を省くことができ、分析プロセスを効率化できるというメリットがあります。

print("外れ値は除去せず、そのまま学習データとして使用します。")

# 外れ値の可視化 (参考)
plt.figure(figsize=(15, 5))
plt.subplot(1, 3, 1)
sns.boxplot(y=df['horsepower'])
plt.title('horsepower boxplot')
plt.subplot(1, 3, 2)
sns.boxplot(y=df['curb_weight'])
plt.title('curb_weight boxplot')
plt.subplot(1, 3, 3)
sns.scatterplot(x=df['engine_size'], y=df[target])
plt.title('engine_size vs price')
plt.tight_layout()
plt.show()

print("処理前後の行数差: 0行")
print("除去件数: 0件")

実行ログ

外れ値は除去せず、そのまま学習データとして使用します。
処理前後の行数差: 0行除去件数: 0件

5. 特徴量選択（エンコーディング前）

このコードは、回帰分析で問題となる多重共線性（特徴量同士が非常に強い相関を持つ状態）を自動的に検出し、整理するための包括的なロジックを実装しています。

主な目的は、相関ヒートマップと**VIF（分散拡大係数）**という2つの指標を組み合わせて、冗長な特徴量を特定し、機械学習モデルの安定性を高めることです。

多重共線性を排除するため、相関係数が0.9以上、またはVIFが20以上の特徴量を整理します。

コードの主要なステップ

このコードは、以下の5つの主要なステップで構成されています。

1. 設定値と変数準備

分析で使用する相関の閾値（CORR_THRESH = 0.9）やVIFの閾値（VIF_THRESH = 20.0）を定義します。また、数値型の特徴量のみを抽出し、VIF計算で問題となる定数列（値がすべて同じ列）を事前に除外します。

2. 相関分析と可視化

すべての数値特徴量間の相関ヒートマップを作成し、視覚的に相関の強さを確認します。その後、設定した閾値（0.9）を超える高相関な特徴量ペアを抽出し、リストとして出力します。これにより、多重共線性の原因となる組み合わせを特定できます。

3. VIFの計算と可視化

VIF（分散拡大係数）を各特徴量について計算し、その結果を降順にソートして表示します。VIFは、他の特徴量によってその特徴量がどの程度説明できるかを示す指標です。VIFが高い（一般的に10以上）ほど、多重共線性の影響を強く受けていることを意味します。このステップでは、VIFの結果を棒グラフでも可視化し、閾値を超える特徴量を一目で把握できるようにします。

4. 削除候補の決定ロジック

このセクションは、このコードの中核となる部分です。

特定列の優先削除: 事前に多重共線性の影響が大きいと分かっている列を優先的に削除候補に加えます。
高相関ペアからの削除: 相関の高いペア（例：lengthとcurb_weight）から、どちらか一方を削除します。削除する特徴量は、VIFがより高い方、または他の列との平均的な相関がより高い方を選びます。これは、「より多くの情報が重複している」特徴量を削除することで、情報の損失を最小限に抑えるためです。
VIF閾値の反復削除: VIFは他の列が削除されると値が変わるため、VIFが20.0を超える特徴量がなくなるまで、最もVIFが高い特徴量を繰り返し削除するロジックを実装しています。これにより、多重共線性の問題を体系的に解消します。

5. 最終的な削除の実行

ステップ4で決定したすべての削除候補を、元のデータフレームから実際に削除します。これにより、多重共線性の問題を解決した、より安定した特徴量セットが次の分析ステップに引き渡されます。

import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from statsmodels.stats.outliers_influence import variance_inflation_factor

# --- 1. 設定値と変数準備 ---
CORR_THRESH = 0.9  # 高相関と判断する閾値
VIF_THRESH = 20.0  # VIFが高すぎると判断する閾値（VIF > 10が一般的だが、ここではより厳しく設定）

# 目的変数を定義
target = 'price'

# 数値カラムのみを抽出し、目的変数は除外
num_cols_all = df.select_dtypes(include=np.number).columns.tolist()
num_feats = [c for c in num_cols_all if c != target]

# 定数列（分散=0）やほぼ定数列を除外
# これらの列はVIFが無限大になるため、計算から除外
def _is_near_constant(s, tol=1e-12):
    return (s.max() - s.min()) < tol

near_const_cols = [c for c in num_feats if _is_near_constant(df[c])]
if near_const_cols:
    print("[LOG] Near-constant numeric columns (drop from VIF calc):", near_const_cols)
# 該当列をVIF計算用のデータフレームから削除
X_num = df[num_feats].drop(columns=near_const_cols, errors="ignore")

# --- 2. 相関分析と可視化 ---
print("--- 相関分析 ---")
# 相関ヒートマップ（数値特徴量同士）
plt.figure(figsize=(10, 8))
corr = X_num.corr()
sns.heatmap(corr, cmap="coolwarm", center=0, annot=False, fmt=".2f", square=True)
plt.title("相関ヒートマップ（数値特徴量間；目的変数除外）")
plt.tight_layout()
plt.show()

# 高相関ペアの抽出（|r| >= 0.90）
high_corr_pairs = []
cols = corr.columns
for i in range(len(cols)):
    for j in range(i+1, len(cols)):
        r = corr.iloc[i, j]
        if np.abs(r) >= CORR_THRESH:
            high_corr_pairs.append((cols[i], cols[j], float(r)))
high_corr_pairs_sorted = sorted(high_corr_pairs, key=lambda x: -abs(x[2]))
print(f"[LOG] 高相関ペア（|r|>={CORR_THRESH}）件数: {len(high_corr_pairs_sorted)}")
if high_corr_pairs_sorted:
    for a, b, r in high_corr_pairs_sorted[:20]:
        print(f"  - {a} vs {b}: r={r:.3f}")
else:
    print("  - 該当なし")


# --- 3. VIFの計算と可視化 ---
print("\n--- VIF分析 ---")
# VIFは欠損があると計算できないため、事前に補完が必要
X_for_vif = X_num.copy()
# statsmodelsのVIFはnp.arrayを渡すため、列順を固定
vif_cols = X_for_vif.columns.tolist()
X_vals = X_for_vif.values
vif_list = []
for i in range(X_vals.shape[1]):
    try:
        vif_val = variance_inflation_factor(X_vals, i)
    except Exception as e:
        vif_val = np.nan # 計算エラーの場合はNaNを代入
    vif_list.append(vif_val)
# 結果をDataFrameにまとめ、VIF降順でソート
vif_df = pd.DataFrame({"feature": vif_cols, "VIF": vif_list}).sort_values("VIF", ascending=False).reset_index(drop=True)

# VIFの棒グラフを可視化
plt.figure(figsize=(8, max(3, 0.3*len(vif_df))))
sns.barplot(data=vif_df.sort_values("VIF", ascending=True), x="VIF", y="feature")
plt.axvline(VIF_THRESH, color="r", linestyle="--", linewidth=1) # 閾値に赤の点線を引く
plt.title("VIF（数値特徴量）")
plt.tight_layout()
plt.show()

print("[LOG] VIF 上位（抜粋）")
print(vif_df.head(20))


# --- 4. 削除候補の決定ロジック ---
to_drop = set()

# 1) 特定列の優先削除
# 'vehicle_length'など、事前に削除を決めた列があれば追加
if "vehicle_length" in X_num.columns:
    to_drop.add("vehicle_length")
    print("\n[LOG] 'vehicle_length' を高相関/VIF対策として削除候補に追加")

# 2) 高相関ペアからの削除方針
# 相関が高いペアから、VIFがより高い方、または他列との平均相関が高い方を削除候補に追加
corr_abs = corr.abs()
mean_abs_corr = corr_abs.mean(axis=1)  # 各列の平均絶対相関を計算
vif_map = dict(zip(vif_df["feature"], vif_df["VIF"])) # VIFを辞書化して高速参照

for a, b, r in high_corr_pairs_sorted:
    if a in to_drop or b in to_drop: # 既に削除候補ならスキップ
        continue
    a_vif = vif_map.get(a, np.nan)
    b_vif = vif_map.get(b, np.nan)
    # まずVIFで比較、同程度なら平均|相関|で比較
    if (not np.isnan(a_vif) and not np.isnan(b_vif)) and (a_vif != b_vif):
        drop_col = a if a_vif >= b_vif else b
    else:
        drop_col = a if mean_abs_corr.get(a, 0) >= mean_abs_corr.get(b, 0) else b
    to_drop.add(drop_col)

# 3) VIF > 閾値 の列を反復的に削除
# VIFは他の列が削除されると値が変わるため、VIFが閾値を超える限り反復的に削除
def compute_vif_frame(df_num):
    # VIFを再計算する関数
    cols_ = df_num.columns.tolist()
    Xv = df_num.values
    out = []
    for i in range(Xv.shape[1]):
        try:
            v = variance_inflation_factor(Xv, i)
        except Exception:
            v = np.nan
        out.append(v)
    return pd.DataFrame({"feature": cols_, "VIF": out}).sort_values("VIF", ascending=False).reset_index(drop=True)

X_iter = X_num.drop(columns=list(to_drop), errors="ignore").copy()
changed = True
while changed and len(X_iter.columns) >= 2:
    changed = False
    vif_now = compute_vif_frame(X_iter)
    worst = vif_now.iloc[0]
    if pd.notnull(worst["VIF"]) and worst["VIF"] > VIF_THRESH:
        to_drop.add(worst["feature"])
        X_iter = X_iter.drop(columns=[worst["feature"]])
        changed = True

# --- 5. 最終的な削除の実行 ---
# 決定した列を元のデータフレームから削除
print(f"\n[LOG] 最終的に削除する数値列: {sorted(to_drop) if to_drop else 'なし'}")
df_preselect = df.drop(columns=list(to_drop), errors="ignore").copy()

print("\n[LOG] 前段選択（相関/VIF）完了")
print(f"[LOG] 形状: {df.shape} -> {df_preselect.shape}")

# 処理結果を次のステップに引き渡すために df を更新
df = df_preselect

実行ログ

--- 相関分析 ---
[LOG] 高相関ペア（|r|>=0.9）件数: 1 - city_mpg vs highway_mpg: r=0.972
--- VIF分析 ---
[LOG] VIF 上位（抜粋） feature VIF 0 width 2899.176689 1 wheel_base 2344.652133 2 length 1944.530402 3 height 1009.086363 4 highway_mpg 533.288219 5 city_mpg 447.973627 6 curb_weight 404.078590 7 bore 294.530819 8 peak_rpm 229.807582 9 stroke 125.880805 10 horsepower 75.951180 11 engine_size 73.980656 12 normalized_losses 23.810968 13 compression_ratio 15.890320 14 symboling 3.206431
[LOG] 最終的に削除する数値列: ['bore', 'curb_weight', 'engine_size', 'height', 'highway_mpg', 'length', 'peak_rpm', 'stroke', 'wheel_base', 'width']

[LOG] 前段選択（相関/VIF）完了
[LOG] 形状: (201, 26) -> (201, 16)

6. カテゴリ変数エンコーディング

概要
カテゴリ変数をOne-Hotエンコーディングで数値化します。出現率1%未満のカテゴリは"Other"に統合し、高カーディナリティを抑えます。

# 'num_of_doors'の欠損値を最頻値で補完
if df['num_of_doors'].isnull().any():
    df['num_of_doors'].fillna(df['num_of_doors'].mode()[0], inplace=True)

# 低頻度カテゴリの統合
for col in df.select_dtypes(include='object').columns:
    counts = df[col].value_counts(normalize=True)
    rare_categories = counts[counts < 0.01].index
    if not rare_categories.empty:
        df[col] = df[col].apply(lambda x: 'Other' if x in rare_categories else x)
        print(f"'{col}' の低頻度カテゴリを 'Other' に統合しました。")

# One-Hotエンコーディング
df_encoded = pd.get_dummies(df, columns=df.select_dtypes(include='object').columns, drop_first=True)

print(f"エンコーディング前の列数: {len(df.columns)}")
print(f"エンコーディング後の列数: {len(df_encoded.columns)}")
print(f"エンコード後のデータ形状: {df_encoded.shape}")

実行ログ

'make' の低頻度カテゴリを 'Other' に統合しました。
'num_of_cylinders' の低頻度カテゴリを 'Other' に統合しました。
'fuel_system' の低頻度カテゴリを 'Other' に統合しました。
エンコーディング前の列数: 16
エンコーディング後の列数: 51
エンコード後のデータ形状: (201, 51)

7. 特徴量選択（エンコード後）＆次元爆発対策

このコードは、機械学習モデルの予測性能を向上させ、モデルの解釈性を高めるために、不要な特徴量を自動的に絞り込むプロセスを実装しています。
LassoCVで係数0の特徴量を削除し、RandomForestのfeature_importances_で重要度上位の特徴量を選択し、最終的な特徴量を決定します。

概要

このコードは、2つの異なる機械学習モデル（LassoCVとRandom Forest）の特徴量選択能力を利用し、最も予測に貢献する特徴量のセットを特定します。LassoCVは不要な特徴量の係数をゼロにすることで特徴量を選択し、Random Forestは各特徴量の予測に対する貢献度（重要度）を計算します。最終的には、両方のモデルが重要だと判断した特徴量を統合して、最終的な特徴量セットを決定します。

詳細解説

1. LassoCVによる特徴量選択

LassoCVは、Lasso（L1正則化）という手法をクロスバリデーション（CV=5）と組み合わせて、最適な正則化の強さ（alpha）を自動的に見つけ出します。Lassoの特徴は、重要でない特徴量の係数を強制的にゼロにすることです。

lasso_model.coef_ != 0：この部分で、Lassoによって係数がゼロにならなかった、つまりモデルが重要だと判断した特徴量だけを抽出しています。
このアプローチは、特に線形モデルにおいて、冗長な特徴量を効率的に排除し、過学習を抑制するのに役立ちます。

2. Random Forestによる特徴量重要度

RandomForestRegressorは、多数の決定木を組み合わせたアンサンブルモデルです。このモデルは、予測を行う際にどの特徴量が最も貢献したかを数値化する**特徴量重要度（Feature Importances）**を算出できます。

rf_model.feature_importances_：各特徴量の重要度を数値として取得します。
importances > 0.005：このコードでは、重要度が0.005という任意の閾値を超える特徴量を選択しています。閾値を用いることで、予測への貢献がごくわずかな特徴量を効率的に除外できます。

Random Forestは非線形な関係も捉えることができるため、LassoCVとは異なる視点から特徴量の重要性を評価します。

3. 最終的な特徴量の決定

最終ステップでは、LassoCVとRandom Forest、両方の手法で選ばれた特徴量を組み合わせます。

set(selected_by_lasso).union(set(selected_by_rf))：union（和集合）を取ることで、どちらか一方でも重要だと判断された特徴量をすべて含めるようにしています。
このハイブリッドなアプローチは、特定のモデルに依存することなく、より包括的で頑健な特徴量セットを構築するための戦略です。

このようにして選ばれたfinal_featuresは、多重共線性の影響を受けにくく、モデルの予測性能を最大化する可能性が高い、最も重要な特徴量集合となります。

X = df_encoded.drop(columns=target)
y = df_encoded[target]

# LassoCVによる特徴量選択
lasso_model = LassoCV(cv=5, random_state=42).fit(X, y)
selected_by_lasso = X.columns[lasso_model.coef_ != 0]

print(f"LassoCVで係数0になった列数: {len(X.columns) - len(selected_by_lasso)}列")

# RandomForestによる特徴量重要度
rf_model = RandomForestRegressor(n_estimators=100, random_state=42).fit(X, y)
importances = pd.Series(rf_model.feature_importances_, index=X.columns).sort_values(ascending=False)
selected_by_rf = importances[importances > 0.005].index # 閾値0.005以上の特徴量を選択

# 最終的な特徴量の決定（両方の手法で選択された特徴量の和集合）
final_features = list(set(selected_by_lasso).union(set(selected_by_rf)))
X_final = df_encoded[final_features]

print(f"最終的に残った特徴量数: {len(final_features)}")
print(f"最終的な特徴量リスト: {final_features}")

実行ログ

LassoCVで係数0になった列数: 44列
最終的に残った特徴量数: 14
最終的な特徴量リスト: ['body_style_hatchback', 'drive_wheels_rwd', 'num_of_cylinders_four', 'fuel_system_2bbl', 'compression_ratio', 'city_mpg', 'aspiration_turbo', 'make_mercedes-benz', 'symboling', 'drive_wheels_fwd', 'fuel_system_mpfi', 'normalized_losses', 'horsepower', 'make_bmw']

8. モデル候補とハイパーパラメータ探索

概要

このコードは、機械学習モデルの予測性能を最大限に引き出すために、最も適切なハイパーパラメータを自動的に見つけ出すプロセスを実装しています。
具体的には、データセットを訓練用とテスト用に分割し、GridSearchCVという手法を用いて、2つの異なるモデル（RandomForestRegressorとXGBoost）の最適なハイパーパラメータを探索します。これにより、手動でパラメータを調整する手間を省き、モデルが過学習や未学習に陥るのを防ぎます。

詳細解説

1. データ分割

まず、train_test_split関数を使って、データセットを**訓練データ（X_train, y_train）とテストデータ（X_test, y_test）**に分割します。

test_size=0.2：データ全体の20%をテストデータとして確保します。
random_state=42：このパラメータを設定することで、何度実行しても同じ方法で分割されるため、結果の再現性を確保します。

モデルは訓練データで学習し、テストデータでその性能を評価することで、未知のデータに対する予測能力を測ることができます。

2. ハイパーパラメータ探索（`GridSearchCV`）

GridSearchCVは、指定されたハイパーパラメータの組み合わせをすべて試すことで、最も性能の良い組み合わせを自動的に特定する手法です。

rf_param_gridとxgb_param_grid：これらは、探索したいハイパーパラメータと、それぞれ試したい値のリストを辞書形式で定義したものです。
- 'n_estimators'：モデルが構築する決定木の数。
- 'max_depth'：個々の決定木の最大の深さ。
- 'learning_rate'：XGBoostにおける学習の進み具合。
GridSearchCV(...)
- cv=5：5-分割交差検証（Cross-Validation）を意味します。データを5つのグループに分け、1つのグループをテストデータ、残りを訓練データとして5回繰り返し学習と評価を行います。これにより、特定のデータ分割に依存しない、より信頼性の高い評価が可能になります。
- scoring='r2'：モデルの評価指標としてR²スコアを使用します。R²スコアは、モデルが目的変数の変動をどの程度説明できているかを示す指標で、1に近いほど良いモデルとされます。
- n_jobs=-1：利用可能なすべてのCPUコアを使って計算を並列化します。これにより、探索時間を大幅に短縮できます。

rf_grid.fit()とxgb_grid.fit()を実行すると、GridSearchCVはそれぞれのモデルについて、設定されたすべてのパラメータの組み合わせを交差検証で評価し、最も高いR²スコアを出したパラメータの組み合わせをbest_params_として返します。

from sklearn.model_selection import GridSearchCV

# データを訓練用とテスト用に分割
X_train, X_test, y_train, y_test = train_test_split(X_final, y, test_size=0.2, random_state=42)

# RandomForestRegressorのハイパーパラメータ探索
rf_param_grid = {'n_estimators': [100, 300], 'max_depth': [5, 10]}
rf_grid = GridSearchCV(RandomForestRegressor(random_state=42), rf_param_grid, cv=5, scoring='r2', n_jobs=-1)
rf_grid.fit(X_train, y_train)

print("RandomForestRegressor 最良パラメータ:", rf_grid.best_params_)
print("RandomForestRegressor 最良CVスコア:", rf_grid.best_score_)

# XGBoostのハイパーパラメータ探索
xgb_param_grid = {'n_estimators': [100, 300], 'learning_rate': [0.05, 0.1]}
xgb_grid = GridSearchCV(XGBRegressor(random_state=42), xgb_param_grid, cv=5, scoring='r2', n_jobs=-1)
xgb_grid.fit(X_train, y_train)

print("XGBoost 最良パラメータ:", xgb_grid.best_params_)
print("XGBoost 最良CVスコア:", xgb_grid.best_score_)

実行ログ

RandomForestRegressor 最良パラメータ: {'max_depth': 10, 'n_estimators': 300}
RandomForestRegressor 最良CVスコア: 0.8208416925448392
XGBoost 最良パラメータ: {'learning_rate': 0.1, 'n_estimators': 100}
XGBoost 最良CVスコア: 0.8008819025419406

9. 評価指標と検証

概要このプロセスは、モデルがまだ学習していない未知のデータに対して、どの程度の精度で予測できるかを客観的に測定するために不可欠です。複数の主要な評価指標（R², RMSE, MAE）を計算し、予測値と実際の値の関係を散布図で可視化することで、モデルのパフォーマンスを包括的に検証します。

詳細解説

1. 最良モデルによる予測

best_model = xgb_grid.best_estimator_：前のステップ（ハイパーパラメータ探索）でGridSearchCVが特定した、最も性能の高いモデル（この場合はXGBoost）を選択します。
y_pred = best_model.predict(X_test)：選択した最良モデルを使って、学習には使われなかったテストデータ（X_test）の目的変数（price）を予測します。このy_predが、モデルがどれだけ正確かを測るための予測値となります。

2. 評価指標の計算

以下の3つの主要な指標を計算し、モデルの性能を数値で示します。

R²（決定係数）：モデルが目的変数の変動をどの程度説明できているかを示す指標です。0から1までの値を取り、1に近いほど良いモデルとされます。
RMSE（平均二乗誤差の平方根）：予測値と実測値の差（誤差）の大きさを評価します。単位は目的変数と同じ（この場合はpriceの対数変換値）になり、値が小さいほど予測精度が高いことを意味します。
MAE（平均絶対誤差）：予測値と実測値の差の絶対値の平均です。RMSEと同様に、値が小さいほど良いことを示し、外れ値の影響をRMSEより受けにくいという特徴があります。

3. 予測値 vs 実測値の可視化

散布図の作成：plt.scatter(y_test, y_pred, ...)で、実測値（y_test）を横軸に、それに対するモデルの予測値（y_pred）を縦軸にプロットします。
理想線の描画：plt.plot(...)で、y = xの理想線（赤の破線）を引きます。もしモデルの予測が完璧であれば、すべての点がこの理想線上に並びます。

この散布図を見ることで、モデルが全体として予測を過大評価しているのか、過小評価しているのか、あるいは特定の範囲で予測誤差が大きいのかといった傾向を視覚的に把握できます。点群が理想線に密接に沿っているほど、モデルの予測精度が高いことを示します。

from sklearn.metrics import mean_squared_error, mean_absolute_error, r2_score

# 最良モデルの選択
best_model = xgb_grid.best_estimator_
y_pred = best_model.predict(X_test)

# 評価指標の計算
r2 = r2_score(y_test, y_pred)
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
mae = mean_absolute_error(y_test, y_pred)

print("--- ホールドアウト検証結果 ---")
print(f"R²: {r2:.4f}")
print(f"RMSE: {rmse:.4f}")
print(f"MAE: {mae:.4f}")

# 予測値 vs 実測値の散布図
plt.figure(figsize=(8, 6))
plt.scatter(y_test, y_pred, alpha=0.7)
plt.plot([min(y_test), max(y_test)], [min(y_test), max(y_test)], 'r--')
plt.title('予測値 vs 実測値')
plt.xlabel('実測値')
plt.ylabel('予測値')
plt.show()

実行ログ

--- ホールドアウト検証結果 ---
R²: 0.8644
RMSE: 0.2294
MAE: 0.1681

10. モデル解釈と可視化

概要

このコードは、以下の3つの主要な分析を通じて、モデルの最終的な検証を行います。

特徴量重要度: どの特徴量がモデルの予測に最も貢献したかを可視化します。
残差分析: モデルの予測誤差（残差）を分析し、モデルの予測が偏っていないか、特定の傾向がないかを確認します。
残差の分布: 誤差が正規分布に従っているかを視覚的に確認します。これは、モデルの前提条件が満たされているかを知る上で重要です。

詳細解説

1. 特徴量重要度 (Feature Importances)

feature_importances_：XGBoostのようなツリー系モデルは、各特徴量が予測の決定にどれだけ影響を与えたかを数値（重要度）として算出できます。
棒グラフ: この数値を使って棒グラフを作成することで、priceの予測に最も寄与した上位20個の特徴量を一目で把握できます。例えば、engine_sizeやhorsepowerが重要度の上位に来る場合、モデルが妥当な判断をしていると解釈できます。

2. 残差プロットと残差ヒストグラム

残差：モデルの予測値と実際の値との差（y_test - y_pred）です。この誤差を分析することで、モデルの弱点を見つけます。
残差プロット：横軸に予測値、縦軸に残差をプロットします。
- 理想的な状態は、残差が予測値の大小にかかわらず**均等にばらついている（等分散性）**ことです。もし、残差が特定のパターン（例：予測値が大きくなると残差も大きくなる）を示している場合、モデルが特定の範囲でうまく予測できていないことを示唆します。
残差ヒストグラム：残差の分布をヒストグラムで可視化します。
- 理想的な状態は、分布が**平均0を中心とした釣鐘型（正規分布）**になることです。誤差が正規分布に従うことは、多くの統計モデルや機械学習モデルの前提であり、モデルの予測が体系的なバイアスを持っていないことの証明になります。

3. Q-Qプロット (Quantile-Quantile Plot)

stats.probplot(...)：残差の分布が正規分布にどれだけ近いかをより厳密に確認するためのプロットです。
プロットの解釈:
- プロット上の点が赤い対角線（理想的な正規分布）に沿って並んでいる場合、残差はほぼ正規分布に従っていると判断できます。
- 点が対角線から離れている場合、分布に歪みや外れ値が含まれている可能性があり、モデルの予測誤差に何らかの偏りがあることを示唆します。

これらの可視化と分析を通じて、モデルの予測結果を単に数値で評価するだけでなく、その背後にある理由や潜在的な問題を深く理解することができます。

import scipy.stats as stats

# 特徴量重要度Top20
feature_importances = pd.Series(best_model.feature_importances_, index=X_final.columns).sort_values(ascending=False)
plt.figure(figsize=(10, 8))
sns.barplot(x=feature_importances.head(20).values, y=feature_importances.head(20).index)
plt.title('特徴量重要度 Top 20')
plt.show()

# 残差プロット
residuals = y_test - y_pred
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
sns.scatterplot(x=y_pred, y=residuals)
plt.axhline(y=0, color='r', linestyle='--')
plt.title('残差プロット')
plt.xlabel('予測値')
plt.ylabel('残差')

# 残差ヒストグラム
plt.subplot(1, 2, 2)
sns.histplot(residuals, kde=True)
plt.title('残差ヒストグラム')
plt.show()

# Q-Qプロット
plt.figure(figsize=(8, 6))
stats.probplot(residuals, dist='norm', plot=plt)
plt.title('Q-Qプロット')
plt.show()

実行ログ

11. モデル保存・運用

前処理から学習器までをPipelineにまとめ、pickleファイルとして保存します。これにより、モデルを再利用する際に、前処理の手順を繰り返す必要がなくなります。

import joblib
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer

# パイプラインのための元のデータセットを再読み込み
df_pipe = pd.read_csv('../dataset/Automobile/data.csv', na_values='?')
df_pipe.dropna(subset=[target], inplace=True)
df_pipe[target] = df_pipe[target].astype(float)
df_pipe[target] = np.log1p(df_pipe[target])

# 数値列とカテゴリ列を再度分類（元の列構成でパイプラインを構築するため）
numerical_features = df_pipe.select_dtypes(include=np.number).columns.tolist()
numerical_features.remove(target)
categorical_features = df_pipe.select_dtypes(include='object').columns.tolist()

# 前処理パイプラインの構築
preprocessor = ColumnTransformer(
    transformers=[
        ('num', SimpleImputer(strategy='median'), numerical_features),
        ('cat', OneHotEncoder(handle_unknown='ignore'), categorical_features)
    ])

# 最終的なパイプラインに LassoCV と XGBoost を組み込む
final_pipeline = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('regressor', XGBRegressor(**xgb_grid.best_params_))
])

# 訓練データ全体でパイプラインを学習
X_full = df_pipe.drop(columns=target)
y_full = df_pipe[target]
final_pipeline.fit(X_full, y_full)

# パイプラインを保存
joblib.dump(final_pipeline, './final_regression_model.pkl')

print("--- モデル保存完了 ---")
print("保存パス: 'models/final_regression_model.pkl'")
print("Pipeline 構成:")
print(final_pipeline)

実行ログ

--- モデル保存完了 ---
保存パス: 'models/final_regression_model.pkl'
Pipeline 構成:
Pipeline(steps=[('preprocessor',
ColumnTransformer(transformers=[('num',
SimpleImputer(strategy='median'),
['symboling',
'normalized_losses',
'wheel_base', 'length',
'width', 'height',
'curb_weight', 'engine_size',
'bore', 'stroke',
'compression_ratio',
'horsepower', 'peak_rpm',
'city_mpg', 'highway_mpg']),
('cat',
OneHotEncoder(handle_unknown='ignore'),
['make', 'fuel...
feature_types=None, feature_weights=None,
gamma=None, grow_policy=None,
importance_type=None,
interaction_constraints=None, learning_rate=0.1,
max_bin=None, max_cat_threshold=None,
max_cat_to_onehot=None, max_delta_step=None,
max_depth=None, max_leaves=None,
min_child_weight=None, missing=nan,
monotone_constraints=None, multi_strategy=None,
n_estimators=100, n_jobs=None,
num_parallel_tree=None, ...))])

まとめ

前処理から保存までの“実務テンプレ”として再利用できます。
非線形性は木系モデルで、線形の読み解きはLasso係数で補完。
実運用では実験管理（MLflow等）とあわせて使うと便利です。