16. Hyperparameter Search and Feature Selection with Scikit¶

# 取得訓練集與測試集
dsWdbc = WdbcDataset()
# NOTE: 回傳的維持是 DataFrame 與 Series 的類型
X_train, X_test, Y_train, Y_test = dsWdbc.get_xy()

# 顯示切割後的大小
print('X: train shape = {}, test shape = {}'.format(X_train.shape, X_test.shape))
print('Y: train shape = {}, test shape = {}'.format(Y_train.shape, Y_test.shape))

# 顯示切割後的類別分布
print('Train classes and counts:\n{}'.format(Y_train.value_counts()))
print('Test classes and counts:\n{}'.format(Y_test.value_counts()))

from sklearn.model_selection import ParameterGrid

# hyperparameters for LogisticRegression
hyparam_list = [{
    'solver': ['newton-cg', 'sag', 'lbfgs'],
    'penalty': ['l2'],
    'max_iter': [300, 200, 100, 90, 80],
    'n_jobs': [-1]
}, {
    'solver': ['saga'],
    'penalty': ['l1', 'l2', 'elasticnet'],
    'max_iter': [300, 200, 100, 90, 80],
    'n_jobs': [-1]
}]

hyparam_grid = ParameterGrid(hyparam_list)
print('Type of ParameterGrid object =', type(hyparam_grid))
print('Number of combinations =', len(list(hyparam_grid)))
for n,p in enumerate(hyparam_grid):
    print('#{}: {}'.format(n, p))

# use grid search
blr_hsgrid = bclass_hyparam_search(blrpipe, X_train, Y_train, hyparams=hyparam_list, randomize=False)

# 搜尋完的結果都記錄在 GridSearchCV 物件的屬性中
print('Hyperparameter grid search result:\n  * best index: {}\n  * best parameters: {}\n  * best score: {:.3f}'.format(
    blr_hsgrid.best_index_, blr_hsgrid.best_params_, blr_hsgrid.best_score_))

# 輸出最佳得分，把 "mean_test" 字串換成 "validate" ，目的是為了記錄成日誌時，避免跟實際與測試集測試的 test score 搞混。
for k in ['mean_test_roc_auc','mean_test_accuracy','mean_test_sensitivity','mean_test_specificity']:
    print('  * {}: {:.3f}'.format(k.replace('mean_test','validate'), blr_hsgrid.cv_results_[k][blr_hsgrid.best_index_]))

import pickle

# use random search
blr_hsrand = bclass_hyparam_search(blrpipe, X_train, Y_train, hyparams=hyparam_list, randomize=True)

print('Hyperparameter random search result:\n  * best index: {}\n  * best parameters: {}\n  * best score: {:.3f}'.format(
    blr_hsrand.best_index_, blr_hsrand.best_params_, blr_hsrand.best_score_))
for k in ['mean_test_roc_auc','mean_test_accuracy','mean_test_sensitivity','mean_test_specificity']:
    print('  * {}: {:.3f}'.format(k.replace('mean_test','validate'), blr_hsrand.cv_results_[k][blr_hsrand.best_index_]))

# save the result and tag with current datetime
jobid = time.strftime('%Y%m%d-%H%M%S')

# log result to file
dfresult = pd.DataFrame(blr_hsrand.cv_results_)
log_file = 'blr-hypsearch-random-{}.csv'.format(jobid)
dfresult.to_csv(log_file, index=False)
print('Hyperparameter random search result logged to file: {}'.format(log_file))

# save the best model to file
model_file = 'blr-hypsearch-random-{}.pickle'.format(jobid)
with open(model_file, 'wb') as f:
    pickle.dump(blr_hsrand.best_estimator_, f, pickle.HIGHEST_PROTOCOL)
print('The best model pickled to file: {}'.format(model_file))

from sklearn.metrics import confusion_matrix

# Load the best model back and test it
with open(model_file, 'rb') as f:
    blr_hsrand_model = pickle.load(f)

# predict with test set
Y_predict = blr_hsrand_model.predict(X_test)

# confusion matrix
TN, FP, FN, TP = confusion_matrix(Y_test, Y_predict).ravel()
print('[Confusion Matrix]:\n  | TP: {} | FP: {} |\n  | FN: {} | TN: {} |\n'.format(TP, FP, FN, TN))

# sensitivity and specificity
sensitivity = recall_score(Y_test, Y_predict, pos_label=1)
print('[Testing Scores]:\n  * Sensitivity: {:.3f}'.format(sensitivity))
specificity = recall_score(Y_test, Y_predict, pos_label=0)
print('  * Specificity: {:.3f}'.format(specificity))
# accuracy and ROC AUC
accuracy = accuracy_score(Y_test, Y_predict)
print('  * Accuracy: {:.3f}'.format(accuracy))
rocauc = roc_auc_score(Y_test, blr_hsrand_model.predict_proba(X_test)[:,1])
rocauc_ci = rocauc_ci_bootstrap(blr_hsrand_model, X_test, Y_test)
print('  * ROC AUC : {:.3f}, 95% CI = ({:.3f}, {:.3f})'.format(rocauc, rocauc_ci[0], rocauc_ci[1]))

# 畫出二元分類的 confusion matrix
plot_confusion_matrix(Y_test, Y_predict)

# 轉成 DataFrame 方便觀察
dfFeatureScore = pd.DataFrame({
    'χ²': wdbc_chi2,
    'χ² p-value': chi2_pvalue,
    'ANOVA f': wdbc_fstat,
    'ANOVA f p-value': fstat_pvalue,
    'Mutual Information': wdbc_mi
}, index=X_train01.columns)

dfFeatureScore

_, axs = plt.subplots(1, 3, figsize=(14,10))
sns.heatmap(dfFeatureScore[['χ²']], annot=True, yticklabels=True, ax=axs[0])
sns.heatmap(dfFeatureScore[['ANOVA f']], annot=True, yticklabels=False, ax=axs[1])
sns.heatmap(dfFeatureScore[['Mutual Information']], annot=True, yticklabels=False, ax=axs[2])

plt.tight_layout()
plt.show()

_, axcorr = plt.subplots(figsize=(16,12))
sns.heatmap(dsWdbc.df.corr(), annot=True, square=True, ax=axcorr)

plt.tight_layout()
plt.show()

chi2p_mask = chi2_pvalue < 1e-2

print('{} features selected by χ² p-value < 0.01:'.format(np.sum(chi2p_mask)))
list(X_train01.columns[chi2p_mask])

from sklearn.feature_selection import SelectFpr

chi2p_selectfpr = SelectFpr(chi2, alpha=0.01)
chi2p_selectfpr.fit(X_train01, Y_train)
chi2p_fpr_mask = chi2p_selectfpr.get_support()

print('{} features selected by χ² FPR test:'.format(np.sum(chi2p_fpr_mask)))
list(X_train01.columns[chi2p_fpr_mask])

from sklearn.feature_selection import SelectFdr

chi2p_selectfdr = SelectFdr(chi2, alpha=0.01)
chi2p_selectfdr.fit(X_train01, Y_train)
chi2p_fdr_mask = chi2p_selectfdr.get_support()

print('{} features selected by χ² FDR test:'.format(np.sum(chi2p_fdr_mask)))
list(X_train01.columns[chi2p_fdr_mask])

from sklearn.feature_selection import SelectFwe

chi2p_selectfwe = SelectFwe(chi2, alpha=0.01)
chi2p_selectfwe.fit(X_train01, Y_train)
chi2p_fwe_mask = chi2p_selectfwe.get_support()

print('{} features selected by χ² FWE test:'.format(np.sum(chi2p_fwe_mask)))
list(X_train01.columns[chi2p_fwe_mask])

k2select = 15
fstat_dec_index = np.flip(np.argsort(wdbc_fstat, kind='stable'))
print('Top {} ANOVA f-value:\n{}\n'.format(k2select, wdbc_fstat[fstat_dec_index[:k2select]]))
print('Top {} features selected from ANOVA f:'.format(k2select))
list(X_train01.columns[np.sort(fstat_dec_index[:k2select])])

from sklearn.feature_selection import SelectKBest

k2select = 15
fstat_selectk = SelectKBest(f_classif, k=k2select)
fstat_selectk.fit(X_train01, Y_train)
fstat_k_mask = fstat_selectk.get_support()

print('Top {} features selected from ANOVA f:'.format(k2select))
list(X_train01.columns[fstat_k_mask])

k2select = 15
mi_dec_index = np.flip(np.argsort(wdbc_mi, kind='stable'))
print('Top {} Mutual Information:\n{}\n'.format(k2select, wdbc_mi[mi_dec_index[:k2select]]))
print('Top {} features selected from Mutual Information:'.format(k2select))
list(X_train01.columns[np.sort(mi_dec_index[:k2select])])

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import StratifiedKFold
from sklearn.feature_selection import SequentialFeatureSelector

# starting from a smaller filtered subset
mi15_index = np.sort(mi_dec_index[:k2select])
X_train_mi15 = X_train.iloc[:, mi15_index]

blrpipe = Pipeline(
    steps=[
        ('zscaler', StandardScaler()),
        ('blr', LogisticRegression())
    ]
)

hyparams = {
    'blr__solver': 'saga',
    'blr__penalty': 'l2',
    'blr__tol': 1e-5,
    'blr__C': 1.0,
    'blr__max_iter': 100,
    'blr__n_jobs': None
}

blrpipe.set_params(**hyparams)

seqfwd_selector = SequentialFeatureSelector(
    blrpipe,
    n_features_to_select=8,
    direction='backward',
    scoring='roc_auc',
    cv=StratifiedKFold(n_splits=10, shuffle=True),
    n_jobs=-1
)

seqfwd_selector.fit(X_train_mi15, Y_train)

seqfwd_mask = seqfwd_selector.get_support()
print('Sequential Forward Selected mask:\n{}\n'.format(seqfwd_mask))
print('Sequential Forward Selected features:')
list(X_train_mi15.columns[seqfwd_mask])

from sklearn.feature_selection import SelectKBest, mutual_info_classif

# the classifier
blrfspipe = Pipeline(
    steps=[
        ('zscaler', StandardScaler()),
        ('fselectk', SelectKBest(mutual_info_classif)),
        ('blr', LogisticRegression())
    ]
)

# hyperparameter searching space for LogisticRegression
hypfs_list = [{
    'fselectk__k': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
    'blr__solver': ['newton-cg', 'sag', 'lbfgs'],
    'blr__penalty': ['l2', 'none'],
    'blr__tol': [1e-4],
    'blr__C': [1.0, 1e-1, 1e-2],
    'blr__max_iter': [120, 100, 80],
    'blr__n_jobs': [-1]
}, {
    'fselectk__k': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10],
    'blr__solver': ['saga'],
    'blr__penalty': ['l1', 'l2', 'elasticnet', 'none'],
    'blr__tol': [1e-4],
    'blr__C': [1.0, 1e-1, 1e-2],
    'blr__max_iter': [120, 100, 80],
    'blr__n_jobs': [-1]
}]

# use grid search
blr_hypfs = bclass_hyparam_search(blrfspipe, X_train, Y_train, hyparams=hypfs_list, randomize=False)

print('Hyperparameter feature search result:\n  * best index: {}\n  * best parameters: {}\n  * best score: {:.3f}'.format(
    blr_hypfs.best_index_, blr_hypfs.best_params_, blr_hypfs.best_score_))
for k in ['mean_test_roc_auc','mean_test_accuracy','mean_test_sensitivity','mean_test_specificity']:
    print('  * {}: {:.3f}'.format(k.replace('mean_test','validate'), blr_hypfs.cv_results_[k][blr_hypfs.best_index_]))

mi_mask = blr_hypfs.best_estimator_['fselectk'].get_support()

print('{} features selected from Mutual Information:'.format(np.sum(mi_mask)))
list(X_train.columns[mi_mask])