【RDKit機械学習入門】特徴量エンジニアリングで予測精度を上げてみる

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from rdkit.Chem import Descriptors
from rdkit.Chem.PandasTools import AddMoleculeColumnToFrame
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score

%matplotlib inline
sns.set_style('whitegrid')


# データセットのURL
url = 'https://raw.githubusercontent.com/deepchem/deepchem/master/datasets/delaney-processed.csv'
df = pd.read_csv(url)

# 目的変数とsmilesデータ以外を削除
cols_to_keep = ["measured log solubility in mols per litre", "smiles"]
df = df[cols_to_keep].copy() 
# カラム名が長いため、簡潔な名前に変更
df = df.rename(columns={'measured log solubility in mols per litre': 'logS'})

AddMoleculeColumnToFrame(df, smilesCol='smiles', molCol='ROMol')

def calculate_all_descriptors(mol):
    # --- 前回の記事の特徴量 ---
    mw   = Descriptors.MolWt(mol)
    logp = Descriptors.MolLogP(mol)
    rb   = Descriptors.NumRotatableBonds(mol)
    aromatic_flags   = [atom.GetIsAromatic() for atom in mol.GetAtoms()]
    num_heavy_atoms  = mol.GetNumHeavyAtoms()
    ap = aromatic_flags.count(True) / num_heavy_atoms if num_heavy_atoms else 0

    # --- 新しい特徴量 ---
    tpsa = Descriptors.TPSA(mol)
    h_donors = Descriptors.NumHDonors(mol)
    h_acceptors = Descriptors.NumHAcceptors(mol)
    
    return [logp, mw, rb, ap, tpsa, h_donors, h_acceptors]

# 新しい特徴量リスト
new_feature_names = [
    "MolLogP", "MolWt", "NumRotatableBonds", "AromaticProportion",
    "TPSA", "NumHDonors", "NumHAcceptors"
]

# すべての分子について記述子を計算
X_desc = np.array(df["ROMol"].apply(calculate_all_descriptors).tolist())

# 目的変数 logS
y = df["logS"].values

X_desc_df = pd.DataFrame(X_desc, columns=new_feature_names)
X_desc_df.head()
X_train, X_test, y_train, y_test = train_test_split(
    X_desc_df, y, test_size=0.2, random_state=42
)

lr_model = LinearRegression()
lr_model.fit(X_train, y_train)
y_pred_lr = lr_model.predict(X_test)
rmse_lr = np.sqrt(mean_squared_error(y_test, y_pred_lr))
r2_lr   = r2_score(y_test, y_pred_lr)

# 実測値と予測値を散布図で比較
plt.figure(figsize=(6, 6))
plt.scatter(y_test, y_pred_lr, alpha=0.5)
plt.plot([-8, 2], [-8, 2], "r--")  # y = x
plt.title("Actual vs. Predicted (Linear Regression)")
plt.xlabel("Actual log S")
plt.ylabel("Predicted log S")
plt.axis("square")
plt.show()

print("--- モデル精度---")
print(f"RMSE  : {rmse_lr:.3f}")
print(f"R²    : {r2_lr:.3f}")

print("\n--- 特徴量の係数 ---")
for name, coef in zip(new_feature_names, lr_model.coef_):
    print(f"{name}: {coef:+.3f}")

def calculate_all_descriptors(mol):
    # --- 元の特徴量 ---
    mw   = Descriptors.MolWt(mol)
    logp = Descriptors.MolLogP(mol)
    rb   = Descriptors.NumRotatableBonds(mol)
    aromatic_flags   = [atom.GetIsAromatic() for atom in mol.GetAtoms()]
    num_heavy_atoms  = mol.GetNumHeavyAtoms()
    ap = aromatic_flags.count(True) / num_heavy_atoms if num_heavy_atoms else 0

    # --- 新しい特徴量 ---
    tpsa = Descriptors.TPSA(mol)
    h_donors = Descriptors.NumHDonors(mol)
    h_acceptors = Descriptors.NumHAcceptors(mol)
    num_rings = Descriptors.RingCount(mol)
    num_hetero = Descriptors.NumHeteroatoms(mol)
    
    return [logp, mw, rb, ap, tpsa, h_donors, h_acceptors, num_rings, num_hetero]

# 新しい特徴量の名前リストも更新する
new_feature_names = [
    "MolLogP", "MolWt", "NumRotatableBonds", "AromaticProportion",
    "TPSA", "NumHDonors", "NumHAcceptors", "NumRings", "NumHeteroatoms"
]