大数据文摘出品
赛题任务为根据从小分子结构中提取的3177个维度特征，预测小分子的六个化学性质。作者将赛题归纳为一个回归问题，直接训练六个模型来分别预测对应的六个性质，此baseline评测得分为8.30。
药物研发是一项成本极高的工作。著名的医学期刊JAMA的一篇调查论文显示，研发一款癌症药物的成本在6.48亿美元左右。其中，大量成本都会用于待选药物分子的测试实验上。虚拟筛选等计算技术可以缩小筛选对象集，降低制药成本，而机器学习技术的出现极大地辅助了计算机辅助药物设计的进步。因此，化学信息学或药物研发领域开始使用各种机器学习技术，包括SVM、随机森林和深度学习，以及图神经网络等。根据分子结构信息预测分子可能的化学性质，将在化学研究和制药领域产生重要应用。
围绕这一课题，今年2月，北京智源研究院联合晶泰科技举办了一场药物研发小分子性质预测赛，赛题任务即根据从小分子结构中提取的3177个维度特征，预测小分子的以下六个化学性质：
- Dipole Moment
- HOMO energy
- LUMO energy
- zero-point vibrational energy
- atomization energy at zero kelvin
- atomization energy at room temperature
这六个性质的准确预测对于药物的发现和开发将提供重要价值。
赛程过半，许多优秀的选手相继开源了比赛的baseline代码模型，赛事承办方biendata整理出几份优秀的baseline文档供各位同学参考学习。
下面的叶枫旭同学的baseline文档，最后评测得分为8.30。
简单分析
这次赛题是一个回归问题，而且要预测的变量有六个。所以比较简单直接的方法就是训练六个模型来分别预测六个性质。
读取数据
代码语言：javascript复制import pandas as pd
import os
import tqdm
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
# 注意改成自己的路径
train_data = pd.read_csv('molecule_open_data/candidate_train.csv')
test_data = pd.read_csv('molecule_open_data/candidate_val.csv')
sample_submission = pd.read_csv('molecule_open_data/sample_submission.csv')
train_answer = pd.read_csv('molecule_open_data/train_answer.csv')
def show_info(df):
print(df.shape)
print(df.columns)
print()
# 查看数据的 shape 和 columns
show_info(train_data)
show_info(test_data)
show_info(sample_submission)
show_info(train_answer)
# 将 p1 ~ p6 拼接到原来的特征
train_data = train_data.merge(train_answer, on='id', how='left')
show_info(train_data)
train_data.head()
数据探索及可视化
代码语言：javascript复制# 查看特征的数据类型。有 21 个特征是浮点数，3156 个特征是整数
print(train_data.info())
print()
print(test_data.info())
# 查看是否有缺失值，结果显示没有缺失值。
train_data.isna().sum().sum()
p1-p6相关统计
代码语言：javascript复制# 查看 p1 ~ p6 的统计指标，如 min, max, mean 等。
train_answer.describe()
# 查看 p1 ~ p6 的分布，可以看到，p5、p6 十分接近
# 所有 property 离群点都不多，p1,2,4,5,6 接近正态分布。
# p1 在右边有个小长尾，也许可以试试对 p1 做一下 log 变换
from scipy.stats import norm
f, ax = plt.subplots(3, 2, figsize=(10, 10))
sns.distplot(train_answer['p1'], kde=True, fit=norm, ax=ax[0, 0])
sns.distplot(train_answer['p2'], kde=True, fit=norm, ax=ax[0, 1])
sns.distplot(train_answer['p3'], kde=True, fit=norm, ax=ax[1, 0])
sns.distplot(train_answer['p4'], kde=True, fit=norm, ax=ax[1, 1])
sns.distplot(train_answer['p5'], kde=True, fit=norm, ax=ax[2, 0])
sns.distplot(train_answer['p6'], kde=True, fit=norm, ax=ax[2, 1])
# 查看 p1 ~ p6 之间的相关性。可以看到，p5 p6 相关性为 1，p3 p4 相关性 0.65。除此之外， property 之间的相关性都比较低。
# 也许可以用某个 property 作为另一个 property 的特征。（这个也许可以尝试）
sns.heatmap(train_answer.corr(), annot=True, linewidths=0.2)
fig=plt.gcf()
fig.set_size_inches(7, 6)
部分可视化特征
代码语言：javascript复制select_feats = ['34', '3157', '3158', '3175']
f, ax = plt.subplots(len(select_feats), 2, figsize=(10, 5 * len(select_feats)))
for i in range(len(select_feats)):
sns.distplot(train_data[select_feats[i]], kde=True, ax=ax[i, 0])
sns.distplot(test_data[select_feats[i]], kde=True, ax=ax[i, 1])
ax[i, 0].set_title('Train data')
ax[i, 1].set_title('Test data')
# 注：可以看到，在这四个特征上，训练集和测试集的分布还是十分接近的。
# 3157，3158，3175 为数值型特征（连续的）。34 从图上看似乎只有 0,1,2 三个取值。
# 3157 和 3158 不仅取值范围相近，而且密度曲线也很接近。
# （3157,1358） 和 3175 虽然取值范围不同，但密度曲线比较接近。
# 使用 value_counts() 方法， 可以看到，特征 34 的取值实际上有 0,1,2,3,4。只不过特征值为 3、4 的样本数占比太少。
print(train_data['34'].value_counts())
print()
print(test_data['34'].value_counts())
# 特征 34 的分布
f, ax=plt.subplots(1, 2, figsize=(10,5))
explode = [0.03] * train_data['34'].nunique()
train_data['34'].value_counts().plot(kind='bar', ax=ax[0])
train_data['34'].value_counts().plot.pie(explode=explode, autopct='%1.1f%%', ax=ax[1],shadow=True)
# 注：在后面的模型训练中，我将所有特征当成连续数值型特征来处理。
# 但其实也可以尝试把某些特征（比如这里的 34）当成类别型特征处理，做一些特征工程，看效果有没有提升
# 直接选取特征 0 到 20。（可以自己修改 show_num 的值，查看更多特征的分布）
show_num = 20
select_feats = [str(i) for i in range(show_num)]
f, ax = plt.subplots(len(select_feats), 2, figsize=(10, 5 * len(select_feats)))
for i in range(len(select_feats)):
sns.distplot(train_data[select_feats[i]], kde=True, ax=ax[i, 0])
sns.distplot(test_data[select_feats[i]], kde=True, ax=ax[i, 1])
ax[i, 0].set_title('Train data')
# 各个特征在训练集和测试集上的分布依旧一致。
# 一个奇怪的现象是，某些特征（0,1,2,3,4,5,19）似乎只有一个取值。我们在下个 cell 验证一下
single_value_cols = []
columns = [str(i) for i in range(3177)]
for col in columns:
if train_data[col].nunique() == 1:
print(col, train_data[col].unique(), test_data[col].unique())
single_value_cols.append(col)
print('只有一个取值的 col 数：', len(single_value_cols))
# 可以看到，有 51 个特征，只有一个取值 0。可以删掉这些特征。（注：本文并没有删除这些特征）
特征工程
这一部分是比较重要的一步。这里我就举个简单粗暴的例子——使用几个强特进行多项式特征生成。
代码语言：javascript复制# 拼接训练集和测试集，做特征的时候就不用分开做了
data = pd.concat([train_data, test_data], axis=0, sort=False).reset_index(drop=True)
from sklearn.preprocessing import PolynomialFeatures
# 不做特征工程的时候跑过一次模型，筛选出来的每个 property 重要性排前 4 的特征
high_imps_list = [
['25', '3175', '34', '3158'],
['25', '16', '3175', '3157'],
['34', '10', '3157', '15'],
['3161', '3155', '3173', '3157'],
['3155', '3173', '3175', '3151'],
['3155', '3173', '3175', '3151']
]
# 每个 property 使用不同的多项式特征
pf_df_list = []
for i in range(len(high_imps_list)):
high_imps = high_imps_list[i]
pf = PolynomialFeatures(degree=5)
pf_array = pf.fit_transform(data[high_imps])
pf_df = pd.DataFrame(pf_array)
pf_df.columns = [f'p{i 1}_pf{col}' for col in pf_df.columns]
# 实际上，pf_0 取值全为 1，pf_1,2,3,4 其实就是 high_imps 的四个特征
drop_columns = [f'p{i 1}_pf{j}' for j in range(len(high_imps) 1)]
pf_df.drop(columns=drop_columns, inplace=True)
pf_df_list.append(pf_df)
# 可以看到，生成的特征比较多，但估计只有部分有用。可以先全部加上，后面慢慢删（根据特征重要性等）。
pf_df_list[0].head()
# 将生成的特征拼接到原特征中(注：后面分类时不同 property 使用不同的多项式特征)
data = pd.concat([data] pf_df_list, axis=1).reset_index(drop=True)
show_info(data)
模型训练（Catboost）和后续分析
代码语言：javascript复制import gc
import pickle
import datetime
import numpy as np
from catboost import CatBoostRegressor, Pool
from sklearn.model_selection import KFold模型训练
代码语言：javascript复制# 6 个 property 共用的初始 3177 特征
use_cols = [col for col in data.columns if col != 'id' and 'p' not in col]
print('Number of common used features:', len(use_cols))
%%time
# 交叉验证折数和待预测的 property 数
n_fold = 5
n_property = 6
train = data[data['p1'].notna()].reset_index(drop=True)
test = data[data['p1'].isna()].reset_index(drop=True)
print('Shape of train dataset:', train.shape)
print('Shape of test dataset:', test.shape)
# 测试集和训练集的预测结果
pred_test = np.zeros((len(test), n_property))
pred_train = np.zeros((len(train), n_property))
# 保存特征重要性的 DataFrame
imps = [pd.DataFrame() for _ in range(n_property)]
smape_score = np.zeros((n_fold, n_property))
all_index = range(len(train))
kfold = KFold(n_splits=n_fold, shuffle=True, random_state=1)
for fold_idx, (tr_idx, va_idx) in enumerate(kfold.split(all_index)):
time_stamp = datetime.datetime.now()
print('*' * 120)
print(f"Fold [{fold_idx}]: " time_stamp.strftime('%Y.%m.%d-%H:%M:%S'))
temp_train = train.iloc[tr_idx]
temp_valid = train.iloc[va_idx]
for p_idx in range(n_property):
print(f'Training for p{p_idx 1} ......')
# 当前 property 使用的特征
temp_cols = use_cols [col for col in data.columns if f'p{p_idx 1}_' in col]
print('Length of temp_cols:', len(temp_cols))
label = 'p' str(p_idx 1)
x_train, x_valid = temp_train[temp_cols], temp_valid[temp_cols]
y_train, y_valid = temp_train[label], temp_valid[label]
cate_features = []
train_pool = Pool(x_train, y_train, cat_features=cate_features)
eval_pool = Pool(x_valid, y_valid, cat_features=cate_features)
cbt_model = CatBoostRegressor(iterations=10000, # 注：baseline 提到的分数是用 iterations=60000 得到的，但运行时间有点久
learning_rate=0.1, # 注：事实上好几个 property 在 lr=0.1 时收敛巨慢。后面可以考虑调大
eval_metric='SMAPE',
use_best_model=True,
random_seed=42,
logging_level='Verbose',
task_type='GPU',
devices='0',
gpu_ram_part=0.5,
early_stopping_rounds=400,
)
cbt_model.fit(train_pool, eval_set=eval_pool, verbose=1000)
smape_score[fold_idx][p_idx] = cbt_model.best_score_['validation']['SMAPE']
# 模型保存
if not os.path.exists('models'):
os.mkdir('models')
with open(f'./models/cbt_f{len(use_cols)}_fold_{fold_idx}_p{p_idx 1}.pkl', 'wb') as f:
pickle.dump(cbt_model, f)
if fold_idx == 0:
imps[p_idx]['feat'] = temp_cols
imps[p_idx][f'score{fold_idx}'] = cbt_model.feature_importances_
# 注：pred_train 保存的是 oof 预测结果，可用于后面的 stacking 等操作。（baseline 里其实没有用到 pred_train）
pred_train[va_idx, p_idx] = cbt_model.predict(x_valid)
pred_test[:, p_idx] = cbt_model.predict(test[temp_cols]) / 5
print()
print('smape_score:', smape_score[fold_idx])
print(f'fold {fold_idx} round {cbt_model.best_iteration_} : smape_mean: {np.mean(smape_score[fold_idx]):.6f}')
del cbt_model, train_pool, eval_pool
del x_train, y_train, x_valid, y_valid
gc.collect()
print('*' * 120)
print('Mean smape in each fold:', smape_score.mean(1))
print()
print('Mean smape of each property:', smape_score.mean(0))
print()
print('Final mean smape:', np.mean(smape_score))特征的重要性
代码语言：javascript复制for i in range(len(imps)):
# 对五折的 importance score 进行平均
imps[i]['score_mean'] = imps[i].apply(lambda row: (row['score0'] row['score1'] row['score2'] row['score3']
row['score4']) / 5, axis=1)
imps[i] = imps[i].sort_values(by='score_mean', ascending=False).reset_index(drop=True)
for imp in imps:
print(imp.head())
print()
# 每个 property 重要性前四的特征
high_imps_list = []
for imp in imps:
high_imps_list.append(imp['feat'].values[:4])
'''
使用原始 3177 特征跑出来的 high_imps_list 如下
[array(['25', '3175', '34', '3158'], dtype=object),
array(['25', '16', '3175', '3157'], dtype=object),
array(['34', '10', '3157', '15'], dtype=object),
array(['3161', '3155', '3173', '3157'], dtype=object),
array(['3155', '3173', '3175', '3151'], dtype=object),
array(['3155', '3173', '3175', '3151'], dtype=object)]
'''
high_imps_list
生成提交文件
代码语言：javascript复制pred_df = pd.DataFrame(pred_test, columns=[f'p{i}' for i in range(1, 7)])
result = pd.concat([test_data[['id']], pred_df], axis=1).reset_index(drop=True)
result.to_csv('result.csv', index=False)
result.head()
总结和分析
实际上，在不做特征工程的情况下，线下得分 8.69，线上 8.39。做了文中的多项式特征生成之后，线下得分 8.62，线上 8.30，提升幅度有限，后续还有一些可提升的空间：
- 针对不同的 property 做不同的特征工程。
- 针对不同的 property 调参。比如，其实 lr=0.1 时收敛有点慢，有些 property 跑了 60000 iterations 还没收敛，这时候就可以调大 lr 看看效果。
- 可以试试将部分只有少数个整数取值的特征，当成类别型变量处理。
- p1 和 p3 的 smape 较大，可以重点分析下原因，做下特征工程。
- 多看看各个特征的分布，看能发现什么规律。
比赛介绍
2020年2月，北京智源人工智能研究院、专注于AI驱动药物研发的科技企业晶泰科技、数据评测平台biendata，共同发布了“智源抗疫-药物研发小分子性质预测赛”，并同步开放了评测竞赛，总奖金10万元。本次比赛要求选手根据从小分子结构中提取的3177个维度特征，预测对于药物发现和开发有重要意义的六个化学性质。本次比赛希望鼓励更多的程序员与算法爱好者发挥计算之所长，通过运用公开数据库或者文献数据与成果，探索如何以AI助力药物研发，为加速新药研究贡献一份力量。
