はじめに

機械学習において「特徴量エンジニアリング」は、モデルの性能を大きく左右する重要な要素です。今回は、LightGBMで高精度なモデルを構築するための特徴量エンジニアリング手法と、SHAP（SHapley Additive exPlanations）を使った最新の解釈可能性技術について詳しく解説します。

データセットの準備

今回は、ECサイトの売上予測を題材にします。時系列要素も含む実践的なデータセットです。

import pandas as pd
import numpy as np
import lightgbm as lgb
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import LabelEncoder
import matplotlib.pyplot as plt
import seaborn as sns
from datetime import datetime, timedelta
import warnings
warnings.filterwarnings('ignore')

# データ生成
np.random.seed(42)
n_days = 365 * 2  # 2年分のデータ
n_products = 100  # 商品数

# 日付の生成
start_date = datetime(2022, 1, 1)
dates = [start_date + timedelta(days=i) for i in range(n_days)]

# 商品マスタデータ
products = pd.DataFrame({
    'product_id': range(n_products),
    'category': np.random.choice(['Electronics', 'Clothing', 'Food', 'Books', 'Toys'], n_products),
    'price': np.random.uniform(10, 1000, n_products),
    'brand': np.random.choice(['A', 'B', 'C', 'D', 'E'], n_products),
    'is_seasonal': np.random.choice([0, 1], n_products, p=[0.7, 0.3])
})

# 売上データの生成
sales_data = []
for date in dates:
    for product_id in range(n_products):
        # 基本的な売上量
        base_sales = np.random.poisson(10)
        
        # 曜日効果
        weekday = date.weekday()
        if weekday in [5, 6]:  # 週末
            base_sales *= 1.5
        
        # 月次効果
        month = date.month
        if month in [11, 12]:  # 年末商戦
            base_sales *= 2
        
        # 商品特性による調整
        product = products.iloc[product_id]
        if product['is_seasonal'] and month in [6, 7, 8]:
            base_sales *= 1.3
        
        # ランダムノイズ
        sales = max(0, int(base_sales + np.random.normal(0, 3)))
        
        sales_data.append({
            'date': date,
            'product_id': product_id,
            'sales': sales,
            'weather': np.random.choice(['sunny', 'cloudy', 'rainy'], p=[0.5, 0.3, 0.2]),
            'temperature': np.random.normal(20 + 10 * np.sin((date.month - 1) * np.pi / 6), 5),
            'campaign': np.random.choice([0, 1], p=[0.8, 0.2]),
            'competitor_price': product['price'] * np.random.uniform(0.8, 1.2)
        })

# DataFrameに変換
df = pd.DataFrame(sales_data)
df = df.merge(products, on='product_id', how='left')

print(f"データセットのサイズ: {df.shape}")
print(f"期間: {df['date'].min()} から {df['date'].max()}")
print("\nデータの先頭:")
print(df.head())

基本的な特徴量エンジニアリング

1. 時系列特徴量の作成

def create_time_features(df):
    """時系列に関する特徴量を作成"""
    df = df.copy()
    
    # 基本的な時間特徴量
    df['year'] = df['date'].dt.year
    df['month'] = df['date'].dt.month
    df['day'] = df['date'].dt.day
    df['dayofweek'] = df['date'].dt.dayofweek
    df['weekofyear'] = df['date'].dt.isocalendar().week
    df['quarter'] = df['date'].dt.quarter
    
    # 週末・月末フラグ
    df['is_weekend'] = (df['dayofweek'] >= 5).astype(int)
    df['is_month_end'] = df['date'].dt.is_month_end.astype(int)
    df['is_month_start'] = df['date'].dt.is_month_start.astype(int)
    
    # 祝日・イベント（簡易版）
    df['is_year_end'] = ((df['month'] == 12) & (df['day'] >= 20)).astype(int)
    df['is_golden_week'] = ((df['month'] == 5) & (df['day'] <= 7)).astype(int)
    
    # 周期的特徴量（sin/cos変換）
    df['month_sin'] = np.sin(2 * np.pi * df['month'] / 12)
    df['month_cos'] = np.cos(2 * np.pi * df['month'] / 12)
    df['day_sin'] = np.sin(2 * np.pi * df['day'] / 31)
    df['day_cos'] = np.cos(2 * np.pi * df['day'] / 31)
    
    return df

df = create_time_features(df)
print("時系列特徴量を追加しました")

2. ラグ特徴量の作成

def create_lag_features(df, target_col='sales', lags=[1, 7, 14, 30]):
    """過去の売上データからラグ特徴量を作成"""
    df = df.sort_values(['product_id', 'date'])
    
    for lag in lags:
        # 単純なラグ
        df[f'{target_col}_lag_{lag}'] = df.groupby('product_id')[target_col].shift(lag)
        
        # 移動平均
        df[f'{target_col}_ma_{lag}'] = df.groupby('product_id')[target_col].transform(
            lambda x: x.rolling(window=lag, min_periods=1).mean()
        )
        
        # 移動標準偏差
        df[f'{target_col}_std_{lag}'] = df.groupby('product_id')[target_col].transform(
            lambda x: x.rolling(window=lag, min_periods=1).std()
        )
    
    # 前日比
    df[f'{target_col}_diff_1'] = df.groupby('product_id')[target_col].diff(1)
    
    # 前週比
    df[f'{target_col}_diff_7'] = df.groupby('product_id')[target_col].diff(7)
    
    return df

df = create_lag_features(df)
print("ラグ特徴量を追加しました")

3. 集約特徴量の作成

def create_aggregate_features(df):
    """カテゴリごとの集約特徴量を作成"""
    # カテゴリ別の統計量
    category_stats = df.groupby(['category', 'date'])['sales'].agg([
        'mean', 'std', 'min', 'max', 'sum'
    ]).reset_index()
    category_stats.columns = ['category', 'date'] + [f'category_sales_{col}' for col in category_stats.columns[2:]]
    
    # ブランド別の統計量
    brand_stats = df.groupby(['brand', 'date'])['sales'].agg([
        'mean', 'sum'
    ]).reset_index()
    brand_stats.columns = ['brand', 'date'] + [f'brand_sales_{col}' for col in brand_stats.columns[2:]]
    
    # マージ
    df = df.merge(category_stats, on=['category', 'date'], how='left')
    df = df.merge(brand_stats, on=['brand', 'date'], how='left')
    
    # 商品の相対的な売上位置
    df['sales_ratio_to_category'] = df['sales'] / (df['category_sales_mean'] + 1)
    df['sales_ratio_to_brand'] = df['sales'] / (df['brand_sales_mean'] + 1)
    
    return df

df = create_aggregate_features(df)
print("集約特徴量を追加しました")

4. 交互作用特徴量の作成

def create_interaction_features(df):
    """特徴量間の交互作用を作成"""
    # 価格と天気の交互作用
    df['price_weather_interaction'] = df['price'] * df['weather'].map({
        'sunny': 1.0, 'cloudy': 0.8, 'rainy': 0.6
    })
    
    # カテゴリと季節の交互作用
    df['category_season_interaction'] = (
        df['category'].map({'Electronics': 2, 'Clothing': 3, 'Food': 1, 'Books': 1, 'Toys': 2}) * 
        df['quarter']
    )
    
    # キャンペーンと週末の交互作用
    df['campaign_weekend_interaction'] = df['campaign'] * df['is_weekend']
    
    # 価格差と温度の交互作用
    df['price_diff'] = df['price'] - df['competitor_price']
    df['price_diff_temp_interaction'] = df['price_diff'] * df['temperature']
    
    return df

df = create_interaction_features(df)
print("交互作用特徴量を追加しました")

高度な特徴量エンジニアリング

1. ターゲットエンコーディング

from sklearn.model_selection import KFold

def target_encoding(df, categorical_cols, target_col='sales', n_splits=5):
    """クロスバリデーションを使ったターゲットエンコーディング"""
    df = df.copy()
    kf = KFold(n_splits=n_splits, shuffle=True, random_state=42)
    
    for col in categorical_cols:
        encoded_col = f'{col}_target_encoded'
        df[encoded_col] = 0
        
        for train_idx, val_idx in kf.split(df):
            # 訓練データでエンコーディングを学習
            encoding_map = df.iloc[train_idx].groupby(col)[target_col].mean()
            
            # 検証データに適用
            df.loc[val_idx, encoded_col] = df.iloc[val_idx][col].map(encoding_map)
        
        # 欠損値は全体平均で埋める
        df[encoded_col].fillna(df[target_col].mean(), inplace=True)
    
    return df

categorical_cols = ['category', 'brand', 'weather']
df = target_encoding(df, categorical_cols)
print("ターゲットエンコーディングを追加しました")

2. 統計的特徴量

def create_statistical_features(df):
    """統計的な特徴量を作成"""
    # 変動係数（標準偏差/平均）
    df['sales_cv_7'] = df['sales_std_7'] / (df['sales_ma_7'] + 1)
    df['sales_cv_30'] = df['sales_std_30'] / (df['sales_ma_30'] + 1)
    
    # トレンド（線形回帰の傾き）
    def calculate_trend(series, window=7):
        if len(series) < window:
            return np.nan
        x = np.arange(window)
        y = series.values[-window:]
        if np.isnan(y).any():
            return np.nan
        slope, _ = np.polyfit(x, y, 1)
        return slope
    
    df['sales_trend_7'] = df.groupby('product_id')['sales'].transform(
        lambda x: x.rolling(window=7, min_periods=7).apply(calculate_trend, raw=False)
    )
    
    # 季節性指標
    monthly_avg = df.groupby(['product_id', 'month'])['sales'].transform('mean')
    yearly_avg = df.groupby('product_id')['sales'].transform('mean')
    df['seasonality_index'] = monthly_avg / (yearly_avg + 1)
    
    return df

df = create_statistical_features(df)
print("統計的特徴量を追加しました")

データの準備とモデル構築

# カテゴリ変数のエンコーディング
label_encoders = {}
for col in ['category', 'brand', 'weather']:
    if col in df.columns:
        le = LabelEncoder()
        df[f'{col}_encoded'] = le.fit_transform(df[col])
        label_encoders[col] = le

# 欠損値を含む最初の期間を除外
df = df.dropna(subset=['sales_lag_30'])

# 特徴量とターゲットの準備
feature_cols = [col for col in df.columns if col not in [
    'date', 'product_id', 'sales', 'category', 'brand', 'weather'
]]

X = df[feature_cols]
y = df['sales']

# 時系列を考慮した分割（最後の20%をテストデータに）
split_date = df['date'].quantile(0.8)
train_mask = df['date'] < split_date
test_mask = df['date'] >= split_date

X_train, X_test = X[train_mask], X[test_mask]
y_train, y_test = y[train_mask], y[test_mask]

print(f"訓練データ: {X_train.shape}")
print(f"テストデータ: {X_test.shape}")
print(f"特徴量の数: {len(feature_cols)}")

LightGBMモデルの構築と特徴量重要度

# モデルの構築
params = {
    'objective': 'regression',
    'metric': 'rmse',
    'boosting_type': 'gbdt',
    'num_leaves': 100,
    'learning_rate': 0.05,
    'feature_fraction': 0.8,
    'bagging_fraction': 0.8,
    'bagging_freq': 5,
    'verbose': -1,
    'random_state': 42
}

train_data = lgb.Dataset(X_train, label=y_train)
valid_data = lgb.Dataset(X_test, label=y_test, reference=train_data)

model = lgb.train(
    params,
    train_data,
    valid_sets=[valid_data],
    num_boost_round=1000,
    callbacks=[lgb.early_stopping(50), lgb.log_evaluation(100)]
)

# 基本的な特徴量重要度
importance_df = pd.DataFrame({
    'feature': feature_cols,
    'importance': model.feature_importance(importance_type='gain')
}).sort_values('importance', ascending=False)

# 上位20個の特徴量を可視化
plt.figure(figsize=(10, 8))
top_features = importance_df.head(20)
plt.barh(range(len(top_features)), top_features['importance'])
plt.yticks(range(len(top_features)), top_features['feature'])
plt.xlabel('重要度（Gain）')
plt.title('特徴量重要度（上位20）')
plt.gca().invert_yaxis()
plt.tight_layout()
plt.show()

print("最も重要な特徴量TOP10:")
print(importance_df.head(10))

SHAP値による高度な解釈

SHAPのインストール

pip install shap

SHAP値の計算と可視化

import shap

# SHAP値の計算
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

# 1. Summary Plot（全体的な特徴量の影響）
plt.figure(figsize=(10, 8))
shap.summary_plot(shap_values, X_test, max_display=20)
plt.tight_layout()
plt.show()

# 2. 特徴量重要度（SHAP値の絶対値平均）
shap_importance = pd.DataFrame({
    'feature': feature_cols,
    'shap_importance': np.abs(shap_values).mean(axis=0)
}).sort_values('shap_importance', ascending=False)

# LightGBMの重要度と比較
comparison_df = pd.merge(
    importance_df[['feature', 'importance']],
    shap_importance[['feature', 'shap_importance']],
    on='feature'
)

# 正規化して比較
comparison_df['lgb_norm'] = comparison_df['importance'] / comparison_df['importance'].max()
comparison_df['shap_norm'] = comparison_df['shap_importance'] / comparison_df['shap_importance'].max()

# 比較プロット
fig, ax = plt.subplots(figsize=(10, 8))
x = np.arange(len(comparison_df.head(15)))
width = 0.35

ax.bar(x - width/2, comparison_df.head(15)['lgb_norm'], width, label='LightGBM重要度')
ax.bar(x + width/2, comparison_df.head(15)['shap_norm'], width, label='SHAP重要度')

ax.set_xlabel('特徴量')
ax.set_ylabel('正規化された重要度')
ax.set_title('LightGBM重要度 vs SHAP重要度')
ax.set_xticks(x)
ax.set_xticklabels(comparison_df.head(15)['feature'], rotation=45, ha='right')
ax.legend()

plt.tight_layout()
plt.show()

個別予測の解釈

# 特定のサンプルの予測を詳しく分析
sample_idx = 0
sample = X_test.iloc[sample_idx:sample_idx+1]

# Force plot（個別予測の説明）
shap.force_plot(
    explainer.expected_value, 
    shap_values[sample_idx], 
    sample,
    matplotlib=True
)
plt.show()

# Waterfall plot（より詳細な個別説明）
shap.waterfall_plot(
    shap.Explanation(
        values=shap_values[sample_idx],
        base_values=explainer.expected_value,
        data=sample.values[0],
        feature_names=feature_cols
    ),
    max_display=15
)
plt.tight_layout()
plt.show()

特徴量間の相互作用分析

# 依存性プロット（特徴量の値と影響の関係）
# 最も重要な特徴量について
top_feature = shap_importance.iloc[0]['feature']
second_feature = shap_importance.iloc[1]['feature']

plt.figure(figsize=(10, 6))
shap.dependence_plot(
    top_feature, 
    shap_values, 
    X_test,
    interaction_index=second_feature,
    alpha=0.5
)
plt.title(f'{top_feature}の影響（{second_feature}との相互作用）')
plt.tight_layout()
plt.show()

# 相互作用の強さを計算
interaction_values = []
for i in range(len(feature_cols)):
    for j in range(i+1, len(feature_cols)):
        # 簡易的な相互作用の推定
        interaction = np.abs(shap_values[:, i] * shap_values[:, j]).mean()
        interaction_values.append({
            'feature1': feature_cols[i],
            'feature2': feature_cols[j],
            'interaction_strength': interaction
        })

interaction_df = pd.DataFrame(interaction_values).sort_values(
    'interaction_strength', ascending=False
)

print("最も強い相互作用TOP10:")
print(interaction_df.head(10))

特徴量選択の自動化

def automated_feature_selection(X, y, initial_features, min_features=10, importance_threshold=0.01):
    """重要度に基づく自動特徴量選択"""
    selected_features = initial_features.copy()
    
    while len(selected_features) > min_features:
        # モデルの学習
        train_data = lgb.Dataset(X[selected_features], label=y)
        model = lgb.train(
            params,
            train_data,
            num_boost_round=100,
            callbacks=[lgb.log_evaluation(0)]
        )
        
        # 重要度の計算
        importance = pd.DataFrame({
            'feature': selected_features,
            'importance': model.feature_importance(importance_type='gain')
        })
        
        # 重要度の正規化
        importance['importance_ratio'] = (
            importance['importance'] / importance['importance'].sum()
        )
        
        # 閾値以下の特徴量を除外
        low_importance_features = importance[
            importance['importance_ratio'] < importance_threshold
        ]['feature'].tolist()
        
        if not low_importance_features:
            break
            
        # 最も重要度の低い特徴量を削除
        feature_to_remove = importance.nsmallest(1, 'importance')['feature'].values[0]
        selected_features.remove(feature_to_remove)
        
        print(f"削除: {feature_to_remove} (重要度: {importance[importance['feature']==feature_to_remove]['importance_ratio'].values[0]:.4f})")
    
    return selected_features

# 特徴量選択の実行
selected_features = automated_feature_selection(
    X_train, y_train, feature_cols.copy(), min_features=20
)

print(f"\n選択された特徴量数: {len(selected_features)}")
print(f"削除された特徴量数: {len(feature_cols) - len(selected_features)}")

実務的なTips

1. 特徴量の品質チェック

def feature_quality_check(df, features):
    """特徴量の品質をチェック"""
    quality_report = []
    
    for feature in features:
        report = {
            'feature': feature,
            'missing_ratio': df[feature].isna().mean(),
            'unique_values': df[feature].nunique(),
            'zero_ratio': (df[feature] == 0).mean(),
            'std': df[feature].std(),
            'skewness': df[feature].skew()
        }
        quality_report.append(report)
    
    quality_df = pd.DataFrame(quality_report)
    
    # 問題のある特徴量を特定
    print("欠損値が多い特徴量（>50%）:")
    print(quality_df[quality_df['missing_ratio'] > 0.5]['feature'].tolist())
    
    print("\n分散が非常に小さい特徴量:")
    print(quality_df[quality_df['std'] < 0.01]['feature'].tolist())
    
    print("\n歪度が大きい特徴量（|skewness| > 3）:")
    print(quality_df[abs(quality_df['skewness']) > 3]['feature'].tolist())
    
    return quality_df

2. 特徴量の保存と管理

# 特徴量エンジニアリングのパイプライン保存
feature_engineering_config = {
    'time_features': ['year', 'month', 'day', 'dayofweek', 'weekofyear'],
    'lag_features': [1, 7, 14, 30],
    'categorical_encoding': label_encoders,
    'selected_features': selected_features,
    'feature_importance': importance_df.to_dict()
}

# JSONやPickleで保存
import json
with open('feature_config.json', 'w') as f:
    json.dump({k: v for k, v in feature_engineering_config.items() 
               if k not in ['categorical_encoding', 'feature_importance']}, f)

まとめ

今回は、LightGBMで高精度なモデルを構築するための特徴量エンジニアリングとSHAPによる解釈について学びました。重要なポイントは：

1. 時系列特徴量: ラグ、移動平均、トレンドなど
2. 集約特徴量: カテゴリ別の統計量や相対的な位置
3. 交互作用特徴量: ビジネス知識に基づく特徴量の組み合わせ
4. SHAP値: モデルの解釈性を高める最新技術
5. 自動特徴量選択: 効率的なモデル構築

次回は「交差検証とアンサンブル学習」について、より堅牢なモデル構築手法を解説します。

演習問題

1. 他の時系列特徴量（曜日×時間帯など）を追加してみましょう
2. カテゴリ変数に対して異なるエンコーディング手法を試してみましょう
3. SHAP値を使って、特定のカテゴリの商品の売上要因を分析してみましょう

---

ご質問やフィードバックがございましたら、ぜひコメント欄でお知らせください。次回もお楽しみに！