Python(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'

# pandasで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'})

# SMILES列からRDKitのMolオブジェクトを生成し、'ROMol'という新しい列を追加
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)
    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"
]

X_desc = np.array(df["ROMol"].apply(calculate_all_descriptors).tolist())

y = df["logS"].values
X_desc_df = pd.DataFrame(X_desc, columns=new_feature_names)


# 特徴量同士の相関行列を計算
correlation_matrix = X_desc_df.corr()

# ヒートマップで可視化
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, cmap='coolwarm', fmt='.2f')
plt.title('Feature Correlation Matrix')
plt.show()

import pandas as pd
import numpy as np
from rdkit import Chem
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
import matplotlib.pyplot as plt
import seaborn as sns

url = 'https://raw.githubusercontent.com/deepchem/deepchem/master/datasets/delaney-processed.csv'
df = pd.read_csv(url)
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')

# 9つの特徴量をすべて計算する関数
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]

# 特徴量のリスト
all_feature_names = [
    "MolLogP", "MolWt", "NumRotatableBonds", "AromaticProportion",
    "TPSA", "NumHDonors", "NumHAcceptors", "NumRings", "NumHeteroatoms"
]

# 全特徴量と目的変数を計算し、欠損値のある行を削除
X_all_list = df["ROMol"].apply(calculate_all_descriptors).tolist()
X_all_df = pd.DataFrame(X_all_list, columns=all_feature_names)
X_all_df['logS'] = df['logS']

y = X_all_df['logS'].values
X_all_df = X_all_df.drop('logS', axis=1)

# 'TPSA' を残し、相関の強い 'NumHDonors', 'NumHAcceptors', 'NumHeteroatoms' を除外
features_strategy1 = [
    "MolLogP", "MolWt", "NumRotatableBonds", "AromaticProportion",
    "TPSA", "NumRings"
]
X_strategy1_df = X_all_df[features_strategy1]

X_train1, X_test1, y_train1, y_test1 = train_test_split(X_strategy1_df, y, test_size=0.2, random_state=42)
model1 = LinearRegression().fit(X_train1, y_train1)
y_pred1 = model1.predict(X_test1)
rmse1 = np.sqrt(mean_squared_error(y_test1, y_pred1))
r2_1 = r2_score(y_test1, y_pred1)

print(f"選択した特徴量: {features_strategy1}")
print(f"RMSE: {rmse1:.4f}")
print(f"R² (決定係数): {r2_1:.4f}")

print("\n--- 戦略2：目的変数との相関の強さで選択 ---")

# 相関の高いグループとlogSの相関を計算
correlated_group = ["TPSA", "NumHDonors", "NumHAcceptors", "NumHeteroatoms"]
temp_df = X_all_df[correlated_group].copy()
temp_df['logS'] = y
logS_correlations = temp_df.corr()['logS'].abs().drop('logS')
best_feature = logS_correlations.idxmax()

print(f"相関の高いグループ: {correlated_group}")
print("各特徴量と logS との相関係数（絶対値）:")
print(logS_correlations)
print(f"\nlogSと最も相関が強い特徴量: '{best_feature}'")

# 'best_feature' を残し、相関グループの他の特徴量を除外
base_features = ["MolLogP", "MolWt", "NumRotatableBonds", "AromaticProportion", "NumRings"]
features_strategy2 = base_features + [best_feature]
features_strategy2 = list(dict.fromkeys(features_strategy2)) 
X_strategy2_df = X_all_df[features_strategy2]

X_train2, X_test2, y_train2, y_test2 = train_test_split(X_strategy2_df, y, test_size=0.2, random_state=42)
model2 = LinearRegression().fit(X_train2, y_train2)
y_pred2 = model2.predict(X_test2)
rmse2 = np.sqrt(mean_squared_error(y_test2, y_pred2))
r2_2 = r2_score(y_test2, y_pred2)

print(f"\n選択した特徴量: {features_strategy2}")
print(f"RMSE: {rmse2:.4f}")
print(f"R² (決定係数): {r2_2:.4f}")

# statsmodelsライブラリのインポート
# ※もし未インストールなら、お使いの環境で !pip install statsmodels を実行してください
from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant

# X_all_dfは9つの特徴量すべてを含むDataFrame
X_const = add_constant(X_all_df) 

# 各特徴量のVIFを計算して、見やすいようにDataFrameに格納
vif_df = pd.DataFrame()
vif_df["feature"] = X_const.columns
vif_df["VIF"] = [variance_inflation_factor(X_const.values, i) for i in range(X_const.shape[1])]

print(vif_df[vif_df["feature"] != 'const'])

X_df_round1 = X_all_df.drop('MolWt', axis=1)
X_const = add_constant(X_df_round1) 
vif_df = pd.DataFrame()
vif_df["feature"] = X_const.columns
vif_df["VIF"] = [variance_inflation_factor(X_const.values, i) for i in range(X_const.shape[1])]
print(vif_df[vif_df["feature"] != 'const'])

# MolWtを削除済みのデータフレームから、さらに 'TPSA' を削除
X_df_round2 = X_df_round1.drop('TPSA', axis=1)

# VIFを再々計算する
X_const_round2 = add_constant(X_df_round2)
vif_df_round2 = pd.DataFrame()
vif_df_round2["feature"] = X_const_round2.columns
vif_df_round2["VIF"] = [variance_inflation_factor(X_const_round2.values, i) for i in range(X_const_round2.shape[1])]

print("--- TPSAを削除後のVIFスコア ---")
print(vif_df_round2[vif_df_round2["feature"] != 'const'])

import pandas as pd
import numpy as np
from rdkit import Chem
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
import matplotlib.pyplot as plt
import seaborn as sns

url = 'https://raw.githubusercontent.com/deepchem/deepchem/master/datasets/delaney-processed.csv'
df = pd.read_csv(url)
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')

# 9つの特徴量をすべて計算する関数
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]

# 特徴量のリスト
all_feature_names = [
    "MolLogP", "MolWt", "NumRotatableBonds", "AromaticProportion",
    "TPSA", "NumHDonors", "NumHAcceptors", "NumRings", "NumHeteroatoms"
]

X_all_list = df["ROMol"].apply(calculate_all_descriptors).tolist()
X_all_df = pd.DataFrame(X_all_list, columns=all_feature_names)
X_all_df['logS'] = df['logS']

y = X_all_df['logS'].values
X_all_df = X_all_df.drop('logS', axis=1)

# 選択する特徴量のリスト
features_strategy_VIF = [
    "MolLogP", "NumRotatableBonds", "AromaticProportion",
    "NumHDonors", "NumHAcceptors", "NumRings", "NumHeteroatoms"
]
X_strategy1_df = X_all_df[features_strategy_VIF]

X_train, X_test, y_train, y_test = train_test_split(X_strategy1_df, y, test_size=0.2, random_state=42)
model_VIF = LinearRegression().fit(X_train, y_train)
y_pred_VIF = model_VIF.predict(X_test)
rmse_VIF = np.sqrt(mean_squared_error(y_test, y_pred_VIF))
r2_VIF = r2_score(y_test, y_pred_VIF)

print(f"選択した特徴量: {features_strategy_VIF}")
print(f"RMSE: {rmse_VIF:.4f}")
print(f"R² (決定係数): {r2_VIF:.4f}")

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