基於機器學習模型預測甲狀腺癌之復發(python)

這是我在碩二上修中山資管系開設「巨量資料分析導論」課程的期末專題，取自真實世界中383名甲狀腺癌患者資料，用下採樣(undersampling)解決資料不平衡的問題。在python中以羅吉斯回歸、SVM、隨機森林、XGboost等模型預測患者是否復發，其中在k-fold交叉驗證下，以XGboost之0.97AUC表現最為優異。最後的Feature Importance中，所有模型皆將「癌症治療的反應」列為第一，因此是相當重要的特徵，然而這可能是過強的解釋變數(strong predictor)。

Gary Chen

26 min readApr 10, 2024

資料來源

資料來源為加州大學爾灣分校機器學習資料庫(UC Irvine Machine Learning Repository)當中的《分化型甲狀腺癌的復發》資料集。這個資料集是2023年10月才發佈的，非常新，沒什麼人用過。

共涵蓋383 名甲狀腺癌患者觀察值，有13 個臨床病理特徵，皆為類別變數，在此僅列出其中幾項：

「 Thyroid Function(甲狀腺功能)」：正常、亢進或低下；

「Adenopathy(淋巴結)」：無淋巴結腫大、右側腫大、左側腫大、雙側腫大、腫大位於後部等；

「Pathology(癌症病理)」：乳突狀癌細胞、濾泡狀癌細胞、Hurthel細胞等；

「T」：此為Tumor(腫瘤)的縮寫，根據癌腫瘤的大小與範圍分為7種等級

「N」：此為Node(淋巴結)的縮寫，根據淋巴結轉移的嚴重程度分為3種等級

「M」：此為Metastasis(轉移)的縮寫，分為「無遠處轉移」與「有遠處轉移」這2種情況

「Response(癌症治療的反應)」：接受癌症治療後患者會有「優異反應」、「不確定反應」等四種情況。

最後是outcome變數，也就是「有復發」與「無復發」兩種結果。

單變量分析(Univariate Analysis)

在拿到資料後，可先進行簡單的探索式資料分析(EDA)，快速判斷資料的分布以及變數間的相關性等等資訊。

年齡(Age): 30 到35 歲高峰的右尾分布

性別(gender) : 女性佔多數(81.5%

癌症分期(stage) : 癌症第一期的患者佔86%

癌症治療的反應(Response) : 接受癌症治療後，反應優良者佔54%

有無復發(Recurred) : 最終無復發者佔72%

雙變量分析(Bivariate Analysis)

上圖為相關係數heatmap，最下排即可觀察到與目標變數「有無復發(Recurred)」相關性前五高變數：

1.Response (癌症治療的反應) : 0.76

2. Risk(復發風險) : 0.73

3. N (Node, 淋巴結) 0.63

4. Adenopathy (淋巴結腫大) 0.61

5. T (Tumor, 腫瘤) 0.56

特徵選擇(Feature Selection)

並不是所有的特徵都要放進模型內，過多的特徵會引起過度擬和，或者引入噪音而干擾模型，故在此剔除掉不相關的變數。

以上將「 Thyroid Function(甲狀腺功能)」與「Physical Examination(身體檢查)」單獨對目標變數跑LPM與Logit模型後，可觀察到每個虛擬變數的P值皆未達統計顯著水準，因此判斷該變數對復發與否並沒有直接相關。

下圖中，特徵當中有一個變數是「Risk」，我判斷這個指的正是復發的風險。然而用「Risk(復發的風險)」去預測「Recurred(有無復發)」，結果必然會非常漂亮。這也是教授在課堂上一直很強調的，當模型太漂亮時要非常小心，你可能已經是拿Y預測Y了! 因此將其刪除。

另外圖中還有一個特徵是「Stage(癌症分期)」，有一到四期。然而癌症分期正是由前者的T、N、M的等級來決定的，這是醫界實務中常用的TNM分期系統(TNM Staging System)。可將「Stage」視為T、N、M的加權平均，兩者必須取一。該選擇誰呢?

「Stage」只是把癌症患者分為四種等級，而TNM分期系統的分類更加細膩。例如圖中皆為癌症一期的患者，但仔細看他們的T(癌腫瘤大小)、N(淋巴結轉移程度)、M(有無遠處轉移)都有所不同。因此剔除分類較粗糙的「Stage(癌症分期)」變數。

資料處理(data processing)

由於有復發者僅為108 人，佔了全樣本(383人)的28%，屬於類別不平衡資料。為避免模型都預測「無復發」，採用隨機下採樣(Undersampling)法，資料經隨機刪除後，「無復發」與「有復發」各為108 筆，全樣本為216 筆。

其中80% 作為訓練集(172 筆)，20% 作為測試集(44 筆) 。另外在特徵轉換上，將類別型變數one hot encoding。

模型預測

羅吉斯回歸(Logistics Regression)

邏輯斯回歸不僅可以拿來分類，也能用於變數上的解釋。結果如下：

年齡為連續變數，故每增加一單位年齡，復發的odd ratio上升e^0.078 =1.08 倍。Response(癌症治療的反應)變數中列出的是Excellent Response( 優異反應) 的結果，並以Indeterminate Response(不確定反應)為基準。即癌症治療反應為優異的患者，相較不確定反應的患者，其復發的機率更低。

再來進行預測，羅吉斯回歸python code與ROC Curve結果如下：

# 'Recurred' is the target column
selected_cols = ['Gender', 'Smoking', 'Hx Smoking', 'Hx Radiothreapy', 'Adenopathy', 
                 'Pathology', 'Focality', 'T', 'N', 'M','Response', 'Recurred']

# Extract the columns needed for one-hot encoding
data_selected = balanced_data[selected_cols]

# Split the data into features and target variable
X = data_selected.drop('Recurred', axis=1)
y = data_selected['Recurred']

# Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)

# Initialize and fit the logistic regression model
logreg = LogisticRegression(max_iter=1000) 
logreg.fit(X_train, y_train)

# Predict on the test set
y_pred = logreg.predict(X_test)

# Predict on training and testing data
y_pred_train = logreg.predict(X_train)
y_pred_test = logreg.predict(X_test)

# Evaluate the model
def evaluate_model(y_true, y_pred):
    precision = precision_score(y_true, y_pred)
    recall = recall_score(y_true, y_pred)
    f1 = f1_score(y_true, y_pred)
    accuracy = accuracy_score(y_true, y_pred)
    cm = confusion_matrix(y_true, y_pred)
    return precision, recall, f1, accuracy, cm

# Evaluate on training set
precision_train, recall_train, f1_train, accuracy_train, cm_train = evaluate_model(y_train, y_pred_train)

y_pred_proba_train = logreg.predict_proba(X_train)[:, 1]
auc_train = roc_auc_score(y_train, y_pred_proba_train)

# Evaluate on testing set
precision_test, recall_test, f1_test, accuracy_test, cm_test = evaluate_model(y_test, y_pred_test)

# Calculate AUC for testing set
y_pred_proba_test = logreg.predict_proba(X_test)[:, 1]
auc_test = roc_auc_score(y_test, y_pred_proba_test)

# Feature importance
feature_importance = pd.DataFrame({'Feature': X_train.columns, 'Importance': logreg.coef_[0]})
feature_importance.sort_values(by='Importance', ascending=False, inplace=True)

print("Training Set Metrics:")
print(f"AUC for Training Set: {auc_train:.4f}")
print(f"Precision: {precision_train}")
print(f"Recall: {recall_train}")
print(f"F1 Score: {f1_train}")
print(f"Accuracy: {accuracy_train}")
print(f"Confusion Matrix:\n{cm_train}")

print("\nTesting Set Metrics:")
print(f"Precision: {precision_test}")
print(f"AUC for Testing Set: {auc_test:.4f}")
print(f"Recall: {recall_test}")
print(f"F1 Score: {f1_test}")
print(f"Accuracy: {accuracy_test}")
print(f"Confusion Matrix:\n{cm_test}")
print(f"AUC for Testing Set: {auc_test:.4f}")

print("\nFeature Importance:")
print(feature_importance)

# K-fold cross-validation
# Define the metrics you want to evaluate
scoring = {
    'precision': make_scorer(precision_score),
    'recall': make_scorer(recall_score),
    'f1_score': make_scorer(f1_score),
    'accuracy': make_scorer(accuracy_score)
}

# Define the k-fold cross-validation
kf = KFold(n_splits=10, shuffle=True, random_state=42)

# Perform cross-validation and compute scores for each metric
cv_results = cross_validate(logreg, X, y, cv=kf, scoring=scoring)

# Calculate mean AUC across k-fold cross-validation
mean_auc_cv = np.mean(auc_scores)
std_auc_cv = np.std(auc_scores)

# Calculate mean scores for each metric
mean_precision = np.mean(cv_results['test_precision'])
mean_recall = np.mean(cv_results['test_recall'])
mean_f1_score = np.mean(cv_results['test_f1_score'])
mean_accuracy = np.mean(cv_results['test_accuracy'])

# Print the mean scores
print("\nK-fold cross validation:")
print(f"Mean AUC across k-fold CV: {mean_auc_cv:.4f} ± {std_auc_cv:.4f}")
print(f"Mean Precision: {mean_precision:.2f}")
print(f"Mean Recall: {mean_recall:.2f}")
print(f"Mean F1 Score: {mean_f1_score:.2f}")
print(f"Mean Accuracy: {mean_accuracy:.2f}")

K-fold(k=10)交叉驗證下，平均AUC為0.97。

支持向量機(SVM)

SVM是監督式的學習方法，原理是找到一個決策邊界(decision boundary)讓兩類之間的邊界(margins)最大化，使其可以完美區隔開來(SVM的推導詳見TommyHuang文章)。

SVM有兩個超參數(Hyperparameter)要調整 : Gamma和C

Gamma = How far the influence of a single training example reaches( 單一樣本所影響的範圍有多大)

圖中的 γ = gamma, 只有rbf核(SVM中最常用的核函數)才需要調整γ，線性核則不需要

若gamma小，資料點的影響力範圍比較遠，對超平面來說，較遠的資料點也有影響力。換句話說，低gamma值意味著模型會考慮更廣泛的點來定義決策邊界，因此能勾勒更加平滑的決策邊界。

若gamma大，資料點的影響力範圍比較近，要使得兩個點被視為屬於同一組，它們需要非常接近。高gamma值會導致模型更加關注每個資料點附近的細節，容易勾勒出更複雜、更不規則的決策邊界，也容易過度擬和。

C = How much you want to punish your model for each misclassified point (對於誤分類的點的懲罰有多大)

C愈大對誤分類的懲罰越大，即要求模型對訓練模型更契合，因此幾乎沒有分類錯誤。然而可能存在過度擬合(Overfitting)，使得模型泛化能力弱。

C愈小則邊界愈軟(The lower the C parameters, the softer the margin)，即允許更多的分類錯誤，此時可能存在欠擬合(Underfitting)，但可提高模型的泛化能力。

超參數Gamma和C的組合與調整都會影響模型的擬合程度，如下圖可見gamma與c若過大(見右下角)會使得SVM過度擬和。

如何選擇最佳的 Gamma和C : 網格搜尋(Grid Search)

# Train SVM model
svm_model = SVC(probability=True, kernel='rbf') 
svm_model.fit(X_train, y_train)

# Define the parameters for tuning
param_grid = {
    'C': [0.01, 0.1, 1, 10, 100],    # Range of regularization parameter C
    'gamma': [0.1, 1, 10]            # Range of gamma values
}

# Perform grid search with cross-validation 
grid_search = GridSearchCV(svm_model, param_grid, cv=10, scoring='accuracy')

# Fit the model for parameter tuning
grid_search.fit(X, y)

# Get the best parameters and the best score
best_params = grid_search.best_params_
best_score = grid_search.best_score_

print("Best Parameters:", best_params)
print("Best Score (Accuracy):", best_score)

交叉驗證結果中，c=1, gamma =0.1時accuracy的表現最好，如下：

Best Parameters: {'C': 1, 'gamma': 0.1}
Best Score (Accuracy): 0.9445887445887446

再來進行預測，SVM python code與ROC Curve結果如下：

# Train SVM model
svm_model = SVC(probability=True, C=1, kernel='rbf', gamma=0.1) 
svm_model.fit(X_train, y_train)

# Predict on the test set
y_pred_svm = svm_model.predict(X_test)
y_pred = svm_model.predict(X_test)

# Predict probabilities for AUC calculation
y_pred_proba_train_svm = svm_model.predict_proba(X_train)[:, 1]
y_pred_proba_test_svm = svm_model.predict_proba(X_test)[:, 1]

# Predict on training and testing data
y_pred_train_svm = svm_model.predict(X_train)
y_pred_test_svm = svm_model.predict(X_test)


# Evaluate the model
def evaluate_model_svm(y_true, y_pred):
    precision = precision_score(y_true, y_pred)
    recall = recall_score(y_true, y_pred)
    f1 = f1_score(y_true, y_pred)
    accuracy = accuracy_score(y_true, y_pred)
    cm = confusion_matrix(y_true, y_pred)
    return precision, recall, f1, accuracy, cm

# Evaluate on training set
precision_train_svm, recall_train_svm, f1_train_svm, accuracy_train_svm, cm_train_svm = evaluate_model_svm(y_train,y_pred_train_svm)

# Evaluate on testing set
precision_test_svm, recall_test_svm, f1_test_svm, accuracy_test_svm, cm_test_svm = evaluate_model_svm(y_test, y_pred_test_svm)

# Calculate AUC for training set
auc_train_svm = roc_auc_score(y_train, y_pred_proba_train_svm)

# Calculate AUC for testing set
auc_test_svm = roc_auc_score(y_test, y_pred_proba_test_svm)

print("Support Vector Machine Metrics:")
print(f"AUC for Training Set: {auc_train_svm:.4f}")
print(f"Precision: {precision_train_svm}")
print(f"Recall: {recall_train_svm}")
print(f"F1 Score: {f1_train_svm}")
print(f"Accuracy: {accuracy_train_svm}")
print(f"Confusion Matrix:\n{cm_train_svm}")

print("\nTesting Set Metrics:")
print(f"AUC for Testing Set: {auc_test_svm:.4f}")
print(f"Precision: {precision_test_svm}")
print(f"Recall: {recall_test_svm}")
print(f"F1 Score: {f1_test_svm}")
print(f"Accuracy: {accuracy_test_svm}")
print(f"Confusion Matrix:\n{cm_test_svm}")

# Feature importance
perm_importance = permutation_importance(svm_model, X_test, y_test, n_repeats=30, random_state=42)

# Get feature importance scores
importance_scores = perm_importance.importances_mean

# Create a DataFrame for feature importance
feature_importance = pd.DataFrame({'Feature': X_train.columns, 'Importance': importance_scores})
feature_importance.sort_values(by='Importance', ascending=False, inplace=True)

print("\nFeature Importance:")
print(feature_importance)

# K-fold cross-validation for SVM
cv_results_svm = cross_validate(svm_model, X, y, cv=kf, scoring=scoring)

# Calculate mean AUC across k-fold cross-validation for SVM
auc_scores_svm = cross_val_score(svm_model, X, y, cv=kf, scoring='roc_auc')
mean_auc_cv_svm = np.mean(auc_scores_svm)
std_auc_cv_svm = np.std(auc_scores_svm)

# Calculate mean scores for each metric
mean_precision_svm = np.mean(cv_results_svm['test_precision'])
mean_recall_svm = np.mean(cv_results_svm['test_recall'])
mean_f1_score_svm = np.mean(cv_results_svm['test_f1_score'])
mean_accuracy_svm = np.mean(cv_results_svm['test_accuracy'])

# Print the mean scores for SVM
print("\nSVM K-fold cross validation:")
print(f"Mean AUC: {mean_auc_cv_svm:.4f} ± {std_auc_cv_svm:.4f}")
print(f"Mean Precision: {mean_precision_svm:.2f}")
print(f"Mean Recall: {mean_recall_svm:.2f}")
print(f"Mean F1 Score: {mean_f1_score_svm:.2f}")
print(f"Mean Accuracy: {mean_accuracy_svm:.2f}")

k-fold交叉驗證下平均AUC為0.95。

隨機森林(Random Forest)

隨機森林(Random Forest)其實就是建構一棵棵各自獨立的決策樹，每棵樹會用到哪些訓練資料及特徵都是隨機決定。最後以投票方式(眾數/取平均)來決定最終的結果。隨機森林ROC Curve結果如下：

k-fold交叉驗證下平均AUC為0.98。

XGBoost(極限梯度提升)

XGBoost算法思路是不斷地添加樹，不斷地進行特征分裂來生長一棵樹，每次添加一個樹，其實是學習一個新函數，去擬合上次預測的殘差。XGBoost作者定義為：" a scalable machine learning system for tree boosting."

其中的"boosting"是關鍵，這和bagging都是集成式學習的方法。boosting的原理是「舊的分類器在訓練有些資料落在confusion area，如果再用全部的data下去訓練，錯的資料永遠都會判錯，因此我們需要針對錯誤的資料去學習，那這樣新訓練出來的分類器才能針對這些錯誤判讀的資料得到好的結果。」(詳見TommyHuang的文章)

XGBoost ROC Curve結果如下：

k-fold交叉驗證下平均AUC為0.97。

結論

本文使用383名甲狀腺癌患者之資料，將「甲狀腺功能」、「身體檢查」、「復發風險」、「癌症分期」等特徵排除後，以隨機下採樣方式將類別平衡後，最終樣本為216筆，即有復發與無復發各為108筆資料。

k=10的交叉驗證下，羅吉斯回歸與機器學習模型的AUC普遍都落在0.93至0.97之間，其中，以XGBoost表現最為優異。然而，太過精準的預測結果要歸咎過於強大的解釋變數(strong predictor)，說明如下。

特徵重要性(Feature Importance)

特徵重要性排名中，「癌症治療的反應(Response)」皆為第一。Why?

根據中國北京協和醫院之《分化型甲狀腺癌治療反應評估體系的解讀》文章，癌症治療的反應為優良者(Excellent)，該患者已達到臨床無瘤生存的狀態，癌症復發率僅為1%~4%。

而結構性不完全反應者(Structural Incomplete Response)，該患者將持續存在或出現新的頸部轉移病灶，若未予治療，50%~85%處於疾病持續狀態，並有50%發生遠處轉移的機率。

而在本文383筆樣本中，反應優良者為208人，其中有復發的僅有1人。結構性不完全反應者有91人，其中有復發的佔89人，這樣極端的分布都顯現該特徵是過於強大的解釋變數(strong predictor)，使得模型預測能力太強。

還記得期末報告時，教授馬上就指出這模型太漂亮「too good to be true」。也就是說這case不需要用到機器學習的方法，簡單跑邏輯斯回歸就有很好的結果了。另外教授也指出，當資料不平衡的時候用下採樣會喪失掉樣本，應考慮使用權重的方法，讓多數類別和少數類別設置不同的權重來實現。

我發現其他人做的結果也常會有模型太漂亮的問題，這也是大家在使用機器學習的一個誤區，認為一定要用fancy的模型才好，但其實不然。若發現有這種問題，就應該要想到有過於強大的解釋變數(strong predictor)，那其實寫個三層左右的if else分類器就有很好的結果了，這時候機器學習就不適合應用在這個資料集上，再換別的資料集吧!

參考文獻

陳欣慧(2023) ,「運用機器學習進行透析中低血壓之偵測」, 國立台灣大學資訊工程學系研究所碩士論文。
Nan Miles Xi(2022), “Improving the diagnosis of thyroid cancer by machine learning and clinical data”, scientific reports.
scikit learn, RBF SVM parameters
Visually Explained, The Kernel Trick in Support Vector Machine (SVM)

基於機器學習模型預測甲狀腺癌之復發(python)

Written by Gary Chen