はじめに

ここまで，化合物の水溶解度を予測する機械学習モデルを作成してきました．

単純な6記述子から作成したベースラインモデル（①）から，RDKit記述子の追加（②），非線形モデルへの切り替え（③），フィンガープリントの導入（④）と段階的にモデルを改善してきました．

最後にSHAP解析を使って「どの記述子がどう予測に効いているか」を可視化し，モデルの中身を解釈します．
また，Delaney の元論文の予測式（ESOL）との性能比較も行い全体をまとめます．
J. S. Delaney, "ESOL: Estimating Aqueous Solubility Directly from Molecular Structure",
J. Chem. Inf. Comput. Sci., 2004, 44, 1000–1005.
https://pubs.acs.org/doi/10.1021/ci034243x

SHAP解析とは

SHAP（SHapley Additive exPlanations）は，「各記述子が、個々の予測にどれだけ貢献したか」を数値化する手法です．

tree系の重要度（.feature_importances_）や線形モデルの係数（.coef_）は「全体として、どの特徴量が重要か」しか分かりません．

SHAPを使うと，さらに踏み込んで：

この分子の予測が低いのはLogP が高いからだ
この分子の予測が高いのは分子量が小さいからだ
というように1分子ごとに「なぜその予測になったか」の理由が見えます．

summary plot では、全分子の SHAP 値を重ねて表示することで，「記述子の値が高い（赤）とき予測を上げるか下げるか」という記述子の効き方の方向まで可視化できます．

SHAP計算の実行

import shap
# 最良モデル(LightGBM チューニング済み)で全訓練データを学習
best_lgb = LGBMRegressor(
    n_estimators=200, max_depth=-1, learning_rate=0.1,
    num_leaves=15, random_state=42, verbose=-1
)
best_lgb.fit(X_train2, y_train2)

# SHAP 値を計算
explainer = shap.TreeExplainer(best_lgb) # LightGBMの木構造を解析してSHAP値計算の準備
shap_values = explainer.shap_values(X_test2) # テストデータの各サンプル×各特徴量について計算

# SHAP summary plot(全体像)
plt.figure()
shap.summary_plot(shap_values, X_test2, feature_names=desc_cols_v2, show=False)
plt.title("SHAP Summary Plot")
plt.tight_layout()
plt.show()

このプロットは「各記述子が予測にどう影響しているか」を示します．

横軸:SHAP 値(正=溶解度を高く予測する方向．負=低く予測する方向)
縦軸:記述子(上から影響が大きい順)
色:赤 = その記述子の値が高い、青 = 低い
点の分布:その記述子の影響のばらつき

MolLogPが圧倒的に重要

横方向の広がりが最も大きい：他の記述子と比べて圧倒的に影響力が大きい．

色と方向の関係
赤い点(LogP が高い = 疎水的) → 負の SHAP(溶解度を低く予測)
青い点(LogP が低い = 親水的) → 正の SHAP(溶解度を高く予測)

Polar Surface Area が2位

極性表面積が大きい(赤) → 正の SHAP(溶解度を上げる方向)
極性表面が大きい分子は水と相互作用しやすいので溶けやすい，という．

また，相関係数では r = +0.12 と弱かったのにSHAP では2位でした．
これは他の記述子との組み合わせで効いている(非線形な効果)ことを意味します．

Molecular Weight が3位

分子量が大きい(赤) → 負の SHAP(溶解度を下げる)
「大きい分子ほど溶けにくい」ということです．
相関(r = -0.64)で2位だったのが，SHAP では3位になりました．
MolLogP と情報が重複しているため(大きい分子は疎水的になりやすい)，LogP がある分，分子量の独自の貢献が減ったのだと思われます．

LightGBMのfeature_importancesとの比較

LightGBMのfeature_importances_（gainベース）との比較も行います．

LightGBM gain: その特徴量で分岐したとき、誤差がどれだけ減ったかの累積
→ 「モデルの学習にどれだけ使われたか」

SHAP: その特徴量が、各予測値をどれだけ動かしたか
→ 「最終的な予測にどれだけ影響したか」
ということで，SHAPの方が信頼性が高いとされているそうです．

# LightGBM 組み込みの重要度(gain ベース)
lgb_gain = LGBMRegressor(
    n_estimators=200, max_depth=-1, learning_rate=0.1,
    num_leaves=15, random_state=42, verbose=-1,
    importance_type="gain"
)
lgb_gain.fit(X_train2, y_train2)

print("=== 特徴量重要度の比較 ===")
print(f"{'記述子':25s} {'LightGBM gain':>15s} {'SHAP':>10s}")
print("-" * 52)

# SHAP の平均絶対値
shap_importance = np.abs(shap_values).mean(axis=0)

for i, col in enumerate(desc_cols_v2):
    print(f"  {col:23s} {lgb_gain.feature_importances_[i]:>13.1f} {shap_importance[i]:>10.3f}")

1位は一致して，2〜5位で少し順位が入れ替わるりますが，大枠は同じだとわかりました．
両者で大きく矛盾していないということも，モデルの信頼性をサポートしています．

文献モデルとの比較

これで一通りモデル作成は終了しました．
このデータセットの元論文である Delaney (2004) の予測式と比較したいと思います．

esol_col = "ESOL predicted log solubility in mols per litre"
# ESOL の R²
esol_r2 = r2_score(y, df[esol_col])

# 自分のモデルの CV R²
best_lgb = LGBMRegressor(
    n_estimators=200, max_depth=-1, learning_rate=0.1,
    num_leaves=15, random_state=42, verbose=-1
)
cv_scores = cross_val_score(best_lgb, X_v2, y, cv=5, scoring="r2")

print(f"ESOL(Delaney):  R² = {esol_r2:.4f}")
print(f"LightGBM(自分): CV R² = {cv_scores.mean():.4f}")

自分のモデルが Delaney の元論文の式を上回りました．
文献を上回った！と喜んでもいいのかもしれませんが，
記述子の数もモデルも違いますし（元論文：4記述子+線形回帰，今回：9記述子+LightGBM），2004年の論文ですしね．．．

また，ESOL の R² は全1128分子で計算した値であるのに対し，自分のモデルの CV R² は5分割交差検証の平均値なので，厳密に同じ条件での比較ではない点に注意が必要です．

記述子の増加と非線形モデルの使用は大きくモデルの改善につながったということがわかります．

比較プロットも作成します．

# 最良モデルで予測を再計算
best_lgb = LGBMRegressor(
    n_estimators=200, max_depth=-1, learning_rate=0.1,
    num_leaves=15, random_state=42, verbose=-1
)
best_lgb.fit(X_train2, y_train2)
y_pred_tuned = best_lgb.predict(X_test2)

# ESOL vs 自分のモデル 比較プロット
fig, axes = plt.subplots(1, 2, figsize=(14, 6))

# 左:ESOL
esol_col = "ESOL predicted log solubility in mols per litre"
axes[0].scatter(y, df[esol_col], alpha=0.5, s=20)
axes[0].plot([-12, 2], [-12, 2], "r--", label="ideal")
axes[0].set_xlabel("Measured logS")
axes[0].set_ylabel("Predicted logS")
esol_r2 = r2_score(y, df[esol_col])
axes[0].set_title(f"ESOL (Delaney's equation, R²={esol_r2:.3f})")
axes[0].legend()

# 右:自分のモデル
axes[1].scatter(y_test2, y_pred_tuned, alpha=0.5, s=20)
axes[1].plot([-12, 2], [-12, 2], "r--", label="ideal")
axes[1].set_xlabel("Measured logS")
axes[1].set_ylabel("Predicted logS")
axes[1].set_title(f"LightGBM (our model, test R²={r2_score(y_test2, y_pred_tuned):.3f})")
axes[1].legend()

plt.suptitle("ESOL vs Our Model", fontsize=14)
plt.tight_layout()
plt.show()

右側のモデルの方が全体的に赤いラインに乗っていることがわかります．（左（ESOL）は全1128分子，右（LightGBM）はテスト226分子でプロットの点の数が異なる点に注意が必要です．）

まとめ

全体をまとめます．

ステップ	内容	CV R²	改善
1	6記述子 + 線形回帰	0.676	―
2	9記述子 + 線形回帰（+LogP）	0.783	+0.107
3	9記述子 + LightGBM	0.888	+0.105
3-2	9記述子 + LightGBM（チューニング）	0.889	+0.001
4	9記述子 + FP + LightGBM	0.906	+0.018

ESOL（Delaney の式）のR² = 0.811を上回りました．

今回のシリーズを通じて実感した，精度改善の優先順位は以下の通りです:

良い特徴量を作る（最大の効果 +0.107）
適切なモデルを選ぶ（大きな効果 +0.105）
フィンガープリントの追加（中程度の効果 +0.018）
ハイパーパラメータの調整（小さな効果 +0.001）

今回の学習を通して，
EDA → 記述子計算 → ベースライン → モデル改善 → 解釈
という一連の流れを学ぶことができました．

モデルの検討や最適化は時間をかければ正直誰でもできうることなので，やっぱり良い特徴量をいかに作成・設計するということが重要だと思いました．そのためにはドメイン知識が必要で，自分なりの観点を持つのが他の人との差別化にも繋がるし，研究としての面白さかなと思います．

AIのおかげで新しい分野について効率よく学べるのは楽しいです．
引き続き色々試していきながら，仕事にも活かしていけたらなと思います．