第12回 PLS回帰（Partial Least Squares Regression）¶

科目：ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕 担当：小野原彩香 対象：3年生（後期）

学習目標¶

この回を終えると、以下のことができるようになります。

PCAとPLSの目的の違いを明確に説明できる
潜在変数（Latent Variable）の概念を理解する
scikit-learnのPLSRegressionを使ってモデルを実装できる
VIPスコアで重要な特徴量を特定できる
交差検証でモデルの汎化性能を評価できる
発展：PCAとPLSが最大化している量の違いと、VIPの閾値が1である理由を、自分で導ける

1. PLSとは（PCAとの違い）¶

1.1 PCAとPLSの目的比較¶

PCAとPLSは、どちらも「高次元データを潜在変数に圧縮する」手法ですが、目的が根本的に異なります。

比較軸	PCA	PLS（PLS回帰）
学習の種類	教師なし学習	教師あり学習
最大化する情報	Xの分散	XとYの共分散
目的	データの構造把握・可視化	Yの予測
潜在変数の方向	Xが最もばらつく方向	XとYが最も共変する方向
ラベルの使用	不要	必要

直感的な説明：

PCAは「データをできるだけ広く広げる方向」を探します。 PLSは「Yを予測するために最も役立つ方向」を探します。

これらは一致しないことが多いです。Xの分散が大きくても、Yとの相関が低い方向は予測に役立ちません。

導出①：PCAとPLSが最大化している量を、式で並べる

表にある「最大化する情報」の行を式にします。どちらも長さ1の重みベクトル \(w\) でスコア \(t = Xw\) を作るところまでは同じで、何を大きくする \(w\) を選ぶかだけが違います。

\[\text{PCA}: \ \max_{\|w\|=1} \mathrm{Var}[Xw] = \max_{\|w\|=1} w^{T}\Sigma w\]

\[\text{PLS}: \ \max_{\|w\|=1} \mathrm{Cov}[Xw,\ y] = \max_{\|w\|=1} w^{T}c, \qquad c = \frac{X^{T}y}{n-1}\]

PLS側は \(w\) について一次式なので、答えがすぐ書けます。長さ1のベクトルと、決まったベクトル \(c\) との内積が最大になるのは、\(w\) が \(c\) と同じ向きのときです。つまり

\[w \propto X^{T}y\]

第1潜在変数の向きは「各変数と \(y\) の共分散を並べたベクトル」そのものです。\(y\) を使う（教師あり）とは、この意味です。PCA側の式には \(y\) がどこにも出てこないので、\(y\) が何であっても同じ向きが出ます。

数値で確かめる：\(x_1\) はばらつきが大きいが \(y\) と無関係、\(x_2\) はばらつきが小さいが \(y\) を決める、という200点のデータを作ると、

	\(x_1\)	\(x_2\)
分散	19.1187	0.9460
\(y\) との共分散	0.1557	1.8806

2つの向きは

\[w_{\text{PCA}} = (-0.999992,\ -0.004082), \qquad w_{\text{PLS}} = (0.082483,\ 0.996592)\]

となり、なす角は 85.03°。ほぼ直交します。それぞれのスコアと \(y\) の関係を見ると、

	スコアの分散	\(y\) との相関
PCA方向	19.1190	−0.0193
PLS方向	1.0819	0.9373

PCAが選んだ向きは、ばらつきは20倍近く大きいのに \(y\) とはほとんど無相関です。「Xの分散が大きくても、Yとの相関が低い方向は予測に役立たない」を、数値がそのまま示しています。

第2以降の潜在変数は、抽出済みの成分を差し引いてから同じことを繰り返して求めます（次節のNIPALS）。その収束の性質までは、この授業で扱う範囲外です。

使った数学：分散と共分散の定義。内積が同じ向きのとき最大になること。

1.2 地理空間LLMへの応用文脈¶

本授業での具体的な問い：

「どの空間特徴量が、LLMのエラー発生率と最も関連しているか？」

X：空間特徴量（クエリ密度、POI密度、駅距離、人口密度など）
Y：エラー発生率（グリッドセルごとの平均エラー率）

PLS回帰を使うことで、エラー発生率を説明・予測する「空間的な構造」を発見できます。

1.3 多重共線性への強さ¶

空間特徴量は互いに相関していることが多いです（例：人口密度とPOI密度は都市部で共に高い）。

通常の重回帰は多重共線性に弱いですが、PLSは潜在変数に圧縮してからYに回帰するため、多重共線性に強いという利点があります。

2. PLS回帰の原理¶

2.1 潜在変数（Latent Variable）¶

PLSは、XとYの両方を少数の潜在変数（LV）に圧縮します。

X行列（n×p）  →  T（スコア行列, n×h）
Y行列（n×q）  →  U（スコア行列, n×h）

ただし h << p（h: 潜在変数数, p: 元の特徴量数）

第1潜在変数は「XとYの共分散が最大になる方向」を抽出します。第2潜在変数はそれに直交する方向で、次に共分散が大きい方向です。

2.2 NIPALS（非線形反復最小二乗法）¶

PLSは通常、NIPALS（Nonlinear Iterative Partial Least Squares）アルゴリズムで解きます。

主なステップ（簡略版）： 1. YのスコアベクトルUを初期化 2. XのウェイトW = X'U / ||X'U||を計算 3. XのスコアT = XWを計算 4. YのローディングQ = Y'T / ||T||²を計算 5. Uを更新 6. 収束するまで2〜5を繰り返す 7. 残差（X, Yを潜在変数で説明した残り）を計算し、次の潜在変数へ

2.3 VIP（Variable Importance in Projection）スコア¶

PLSでは各特徴量の重要度をVIPスコアで定量化できます。

\[\text{VIP}_j = \sqrt{p \cdot \frac{\sum_{h=1}^{H} (w_{jh}^2 \cdot \text{SS}_h)}{\sum_{h=1}^{H} \text{SS}_h}}\]

\(w_{jh}\)：第h潜在変数における特徴量jのウェイト
\(\text{SS}_h\)：第h潜在変数が説明するYの分散
\(p\)：特徴量数

判断基準： - VIP > 1.0 → 重要な特徴量 - VIP > 1.5 → 特に重要な特徴量 - VIP < 0.8 → 重要度が低い特徴量

導出②：なぜ VIP の閾値が「1」なのか

VIPの判断基準が 1.0 なのは慣習ではなく、式から出てきます。定義の両辺を二乗します。

\[\text{VIP}_j^2 = p \cdot \frac{\sum_{h} w_{jh}^2\,\text{SS}_h}{\sum_{h} \text{SS}_h}\]

これを全特徴量について足します。\(j\) についての和を中に入れると、

\[\sum_{j=1}^{p} \text{VIP}_j^2 = p \cdot \frac{\sum_{h} \left(\text{SS}_h \sum_{j} w_{jh}^2\right)}{\sum_{h} \text{SS}_h}\]

ここで各潜在変数のウェイトベクトル \(w_h\) は長さ1に正規化されています（\(\sum_j w_{jh}^2 = 1\)）。したがって括弧の中は \(\text{SS}_h\) そのものになり、

\[\sum_{j=1}^{p} \text{VIP}_j^2 = p \cdot \frac{\sum_{h}\text{SS}_h}{\sum_{h}\text{SS}_h} = p\]

両辺を \(p\) で割ると

\[\frac{1}{p}\sum_{j=1}^{p} \text{VIP}_j^2 = 1\]

VIPの二乗の平均は、どんなデータでも必ず1になります。つまりVIPは「平均的な特徴量を1としたときの相対的な重要度」であり、1を超えるとは「平均より貢献が大きい」という意味です。閾値1は、この性質から自然に出てくる基準です。

同時に注意も読み取れます。合計が \(p\) で固定されている以上、VIP>1の特徴量は必ず存在します。仮に全部の変数が \(y\) と無関係でも、その中で相対的に大きいものは1を超えます。VIPが1を超えたことは「この変数は重要だ」ではなく「この変数はほかより相対的に効いている」を意味するだけです。まずモデル全体が \(y\) を説明できているか（\(R^2\) や交差検証のRMSE）を確認してから読む必要があります。

使った数学：総和の順序の入れ替え。単位ベクトルの成分の二乗和が1であること。

3. 空間特徴量×エラー発生率への適用¶

3.1 データ準備¶

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import cross_val_score, KFold
from sklearn.metrics import r2_score, mean_squared_error

plt.rcParams['font.family'] = 'Hiragino Sans'
plt.rcParams['axes.unicode_minus'] = False

# データ読み込み
df = pd.read_parquet('data/spatial_features.parquet')

# 特徴量とターゲットの定義
feature_cols = [
    'query_density',
    'poi_density',
    'road_density',
    'nearest_station_dist',
    'green_area_ratio',
    'building_coverage',
    'elevation_mean',
    'population_density',
    'weekday_query_ratio',     # 平日クエリ比率
    'night_query_ratio',       # 夜間クエリ比率
    'query_text_avg_length',   # クエリテキスト平均文字数
]
target_col = 'error_rate'

# 欠損値を除外
data = df[feature_cols + [target_col]].dropna()
X = data[feature_cols].values
y = data[target_col].values.reshape(-1, 1)

print(f"サンプル数: {len(X)}")
print(f"特徴量数: {len(feature_cols)}")
print(f"エラー率 統計量:")
print(f"  平均: {y.mean():.4f}, 標準偏差: {y.std():.4f}")
print(f"  最小: {y.min():.4f}, 最大: {y.max():.4f}")

# 標準化（X, Yともに）
scaler_X = StandardScaler()
scaler_y = StandardScaler()
X_scaled = scaler_X.fit_transform(X)
y_scaled = scaler_y.fit_transform(y)

3.2 事前分析：相関ヒートマップ¶

# 特徴量とターゲットの相関を確認
df_analysis = pd.DataFrame(X_scaled, columns=feature_cols)
df_analysis['error_rate'] = y_scaled.ravel()

corr_matrix = df_analysis.corr()

plt.figure(figsize=(12, 10))
mask = np.triu(np.ones_like(corr_matrix, dtype=bool))
sns.heatmap(
    corr_matrix, mask=mask, annot=True, fmt='.2f',
    cmap='coolwarm', center=0, vmin=-1, vmax=1,
    linewidths=0.5, cbar_kws={'label': 'ピアソン相関係数'}
)
plt.title('特徴量間の相関行列（error_rate含む）')
plt.tight_layout()
plt.savefig('output/correlation_matrix.png', dpi=150, bbox_inches='tight')
plt.show()

# エラー率との相関が高い特徴量
corr_with_error = corr_matrix['error_rate'].drop('error_rate').abs().sort_values(ascending=False)
print("\nerror_rateとの相関係数（絶対値順）:")
print(corr_with_error)

4. モデルの評価と解釈¶

4.1 潜在変数数の選択（交差検証）¶

# 潜在変数数（LV数）をCVで選択
kf = KFold(n_splits=5, shuffle=True, random_state=42)
lv_range = range(1, min(len(feature_cols), 10) + 1)
cv_rmse = []

for n_lv in lv_range:
    pls = PLSRegression(n_components=n_lv, scale=False)
    # neg_MSEを返すため、-1倍してsqrtでRMSEに変換
    neg_mse_scores = cross_val_score(
        pls, X_scaled, y_scaled.ravel(),
        scoring='neg_mean_squared_error',
        cv=kf
    )
    rmse = np.sqrt(-neg_mse_scores.mean())
    cv_rmse.append(rmse)
    print(f"LV={n_lv}: CV-RMSE={rmse:.4f}")

# 最適LV数
best_n_lv = lv_range[np.argmin(cv_rmse)]
print(f"\n最適潜在変数数: {best_n_lv}")

# グラフ
plt.figure(figsize=(8, 5))
plt.plot(list(lv_range), cv_rmse, 'bo-', markersize=8)
plt.axvline(x=best_n_lv, color='red', linestyle='--', label=f'最適LV={best_n_lv}')
plt.xlabel('潜在変数数（LV数）')
plt.ylabel('CV-RMSE')
plt.title('LV数の選択（5-Fold CV）')
plt.legend()
plt.grid(alpha=0.3)
plt.savefig('output/pls_cv_lv.png', dpi=150, bbox_inches='tight')
plt.show()

4.2 最終モデルの学習と評価¶

# 最適LV数でモデル学習
pls_final = PLSRegression(n_components=best_n_lv, scale=False)
pls_final.fit(X_scaled, y_scaled)

# 予測
y_pred_scaled = pls_final.predict(X_scaled)
y_pred = scaler_y.inverse_transform(y_pred_scaled)
y_true = scaler_y.inverse_transform(y_scaled)

# 評価指標
r2 = r2_score(y_true, y_pred)
rmse = np.sqrt(mean_squared_error(y_true, y_pred))
print(f"=== モデル評価 ===")
print(f"R²: {r2:.4f}")
print(f"RMSE: {rmse:.6f}")

# 実測 vs 予測 散布図
plt.figure(figsize=(8, 8))
plt.scatter(y_true, y_pred, alpha=0.3, s=20, c='steelblue')
min_val = min(y_true.min(), y_pred.min())
max_val = max(y_true.max(), y_pred.max())
plt.plot([min_val, max_val], [min_val, max_val], 'r--', lw=2, label='完全予測線')
plt.xlabel('実測エラー率')
plt.ylabel('予測エラー率')
plt.title(f'PLS回帰：実測 vs 予測（R²={r2:.3f}）')
plt.legend()
plt.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('output/pls_actual_vs_predicted.png', dpi=150, bbox_inches='tight')
plt.show()

4.3 VIPスコアの計算と可視化¶

def calculate_vip(pls_model: PLSRegression) -> np.ndarray:
    """
    PLSモデルからVIPスコアを計算する

    参考: Wold et al. (2001)
    """
    t = pls_model.x_scores_       # スコア行列 (n, h)
    w = pls_model.x_weights_      # ウェイト行列 (p, h)
    q = pls_model.y_loadings_     # Yローディング (q, h)

    p = w.shape[0]  # 特徴量数
    h = t.shape[1]  # 潜在変数数

    # 各潜在変数がYを説明する分散（SS）
    ss = np.diag(t.T @ t @ q.T @ q)

    # VIPスコアの計算
    vip = np.sqrt(p * np.sum(w**2 * ss, axis=1) / np.sum(ss))
    return vip

vip_scores = calculate_vip(pls_final)
vip_df = pd.DataFrame({
    'feature': feature_cols,
    'vip': vip_scores
}).sort_values('vip', ascending=False)

print("=== VIPスコア（重要度順）===")
print(vip_df.to_string(index=False))

# VIPスコアの棒グラフ
plt.figure(figsize=(10, 6))
colors = ['crimson' if v >= 1.0 else 'steelblue' for v in vip_df['vip']]
plt.barh(vip_df['feature'], vip_df['vip'], color=colors)
plt.axvline(x=1.0, color='black', linestyle='--', lw=2, label='VIP=1.0（重要性の閾値）')
plt.axvline(x=1.5, color='gray', linestyle=':', lw=2, label='VIP=1.5（高重要性）')
plt.xlabel('VIPスコア')
plt.title('PLS回帰：特徴量重要度（VIPスコア）\nエラー発生率の予測に対する貢献度')
plt.legend()
plt.gca().invert_yaxis()  # 上位が上に来るよう反転
plt.tight_layout()
plt.savefig('output/pls_vip_scores.png', dpi=150, bbox_inches='tight')
plt.show()

# VIP > 1.0の特徴量を抽出
important_features = vip_df[vip_df['vip'] >= 1.0]['feature'].tolist()
print(f"\n重要特徴量（VIP≥1.0）: {important_features}")

4.4 PLS潜在変数スコアの可視化¶

# スコアプロット（LV1 vs LV2）
T = pls_final.x_scores_

plt.figure(figsize=(10, 8))
scatter = plt.scatter(T[:, 0], T[:, 1],
                      c=y_true.ravel(),
                      cmap='coolwarm', alpha=0.5, s=30)
plt.colorbar(scatter, label='エラー率')
plt.xlabel(f'LV1スコア')
plt.ylabel(f'LV2スコア')
plt.title('PLS スコアプロット（LV1 vs LV2）\n色: エラー率')
plt.axhline(y=0, color='gray', linestyle='-', alpha=0.3)
plt.axvline(x=0, color='gray', linestyle='-', alpha=0.3)
plt.tight_layout()
plt.savefig('output/pls_score_plot.png', dpi=150, bbox_inches='tight')
plt.show()

4.5 空間的な残差の分析¶

import geopandas as gpd
from shapely.geometry import Point

# 残差を計算
residuals = y_true.ravel() - y_pred.ravel()

# GeoDataFrameに統合
df_valid = data.copy().reset_index(drop=True)
df_valid['predicted_error_rate'] = y_pred.ravel()
df_valid['residual'] = residuals

gdf_result = gpd.GeoDataFrame(
    df_valid,
    geometry=[Point(lon, lat) for lon, lat in
              zip(df['longitude'].loc[data.index],
                  df['latitude'].loc[data.index])],
    crs='EPSG:4326'
)

# 残差の地図表示
fig, axes = plt.subplots(1, 2, figsize=(16, 7))

# 左：予測エラー率
gdf_result.plot(column='predicted_error_rate', ax=axes[0],
                cmap='Reds', markersize=20, alpha=0.7, legend=True,
                legend_kwds={'label': '予測エラー率'})
axes[0].set_title('予測エラー率（PLS）')

# 右：残差（正: 過小予測、負: 過大予測）
gdf_result.plot(column='residual', ax=axes[1],
                cmap='RdBu_r', markersize=20, alpha=0.7, legend=True,
                legend_kwds={'label': '残差（実測-予測）'})
axes[1].set_title('予測残差の空間分布')

plt.tight_layout()
plt.savefig('output/pls_spatial_residual.png', dpi=150, bbox_inches='tight')
plt.show()

5. ハンズオン手順¶

5.1 実行ステップの確認¶

# 出力ディレクトリ確認
ls output/

# データ確認
python -c "
import pandas as pd
df = pd.read_parquet('data/spatial_features.parquet')
print(df.shape)
print(df['error_rate'].describe())
"

5.2 分析の流れまとめ¶

1. データ読み込み（data/spatial_features.parquet）
   ↓
2. 欠損値処理・標準化
   ↓
3. 相関行列でXとYの関係を事前確認
   ↓
4. 5-Fold CVで最適LV数を決定
   ↓
5. 最終モデル学習（PLSRegression）
   ↓
6. R²・RMSEで性能評価
   ↓
7. VIPスコアで重要特徴量を特定
   ↓
8. 残差の空間分布を可視化
   ↓
9. 結果の解釈・考察

6. 演習課題¶

課題1（基本）¶

data/spatial_features.parquetを読み込み、以下を実行してください。

error_rateと各特徴量のピアソン相関係数を計算し、散布図行列（pairplot）を描く
5-Fold CVでLV数を選択し、CV-RMSEのグラフを描く

課題2（応用）¶

最適LV数でPLS回帰を学習し、R²とRMSEを求める
VIPスコアを計算し、VIP≥1.0の特徴量を特定する
VIPスコアが最も高い特徴量とerror_rateの散布図を描き、関係を考察する

課題3（発展）¶

PLSで選択した重要特徴量のみを使った通常の重回帰（LinearRegression）と、PLS全特徴量の予測精度を比較する
VIPスコアが高かった特徴量について「なぜエラー発生率と関係があるのか」を地理的・技術的に考察する（200字以上）
学習データと検証データを空間的に分割（緯度で分ける）して、空間的汎化性能を評価する

提出形式：Jupyterノートブック、ファイル名：12_pls_regression_[学籍番号].ipynb

まとめ¶

項目	ポイント
PCA vs PLS	PCA：Xの分散最大化（教師なし）、PLS：XYの共分散最大化（教師あり）
潜在変数	XとYが最も共変する低次元方向
多重共線性	空間特徴量は互いに相関しているが、PLSはそれに強い
VIPスコア	VIP≥1.0が重要特徴量の目安、エラー発生の原因特定に活用
交差検証	LV数の選択にCV-RMSEを使用、過学習を防ぐ
空間残差	残差の空間パターンから「モデルが説明できない要因」を探索
次回への接続	分析基盤となるDuckDB/Polarsパイプラインを第13回で深掘りする

次回予告： 第13回では、これまでの分析で使ってきたDuckDBとPolarsについて、より体系的に理解し、パイプラインの高速化・効率化を学びます。

作成：小野原彩香 / ビッグデータ応用Ⅱ〔特徴量エンジニアリング〕