{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Machine Learning pipeline: Shu *et al*, 2020\n", "\n", "This notebook attempts to replicate the Machine Learning pipeline of the following paper:\n", "\n", "
\n", "Shu, Zhen‐Yu, et al. Predicting the progression of Parkinson's disease using conventional MRI and machine learning: An application of radiomic biomarkers in whole‐brain white matter. Magnetic Resonance in Medicine 85.3 (2021): 1611-1624.
" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import sys, os\n", "sys.path.append('code/scripts')\n", "import os\n", "import glob\n", "import pandas as pd\n", "import numpy as np\n", "from sklearn.model_selection import GridSearchCV, StratifiedKFold, cross_val_score\n", "from sklearn.base import BaseEstimator, TransformerMixin\n", "from sklearn.decomposition import PCA\n", "from sklearn.pipeline import Pipeline\n", "from sklearn.tree import DecisionTreeClassifier\n", "from sklearn.neighbors import KNeighborsClassifier\n", "from sklearn.naive_bayes import GaussianNB\n", "from sklearn.svm import SVC\n", "from sklearn.preprocessing import StandardScaler\n", "from mrmr.pandas import mrmr_classif\n", "from verify import check_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Usage\n", "\n", "Assign the name of the CSV file containing your PPMI data to the `data` variable. Additionaly, enter cohort information in the `cohort` list:\n", "\n", "Specify the directory with cohort details. Ensure the following for each cohort:\n", "\n", "- **cohort.csv**: Must have columns - PATNO, EVENT_ID, Description.\n", "\n", "- **demographics**.csv: Should include columns - PATNO, Age, Stage, SEX, output, UPDRS_TOT (optional)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "data = \"radiomics_results.csv\"\n", "\n", "cohorts = [\n", " \"cohorts/max_cohorts/cohort_max_off_2_35_matched\",\n", " \"cohorts/max_cohorts/cohort_max_off_off_65_75_unmatched\",\n", " \"cohorts/max_cohorts/cohort_max_both_2_59_matched\",\n", " \"cohorts/max_cohorts/cohort_max_both_102_92_unmatched\",\n", " \"cohorts/max_cohorts/cohort_max_on_2_43_matched\",\n", " \"cohorts/max_cohorts/cohort_max_on_on_66_81_unmatched\",\n", " \"cohorts/max_cohorts/cohort_max_2_72_matched\"\n", "]\n", "\n", "# Verify data file\n", "assert os.path.exists(data), f\"{data} does not exist.\"\n", "\n", "# Verify data\n", "assert all(check_data(cohort) for cohort in cohorts), \"Not all cohorts have valid data files.\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Helper functions" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def get_df(data, cohort, demographics, useUPDRS=False):\n", " \"\"\"\n", " Merge and preprocess data from multiple CSV files to create a consolidated dataframe.\n", "\n", " Parameters:\n", " - data (str): The path to the primary data CSV file.\n", " - cohort (str): The path to the 'cohort.csv' file containing cohort information.\n", " - demographics (str): The path to the 'demographics*.csv' file containing demographic information.\n", " - useUPDRS (bool, optional): Whether to include the 'UPDRS_TOT' column in the output dataframe. Defaults to False.\n", "\n", " Returns:\n", " - pandas.DataFrame: A consolidated dataframe containing merged data from the specified files.\n", " \"\"\"\n", " # Collect data files\n", " data_df = pd.read_csv(data)\n", " cohort_df = pd.read_csv(cohort)\n", " demographics_df = pd.read_csv(demographics)\n", "\n", " # Create output dataframe\n", " columnsToDrop = [\"EVENT_ID\", \"Description\", \"Age\", \"Stage\", \"SEX\", \"UPDRS_TOT\"] if useUPDRS else [\"EVENT_ID\", \"Description\", \"Age\", \"Stage\", \"SEX\"]\n", " output_df = cohort_df \\\n", " .merge(demographics_df, on=\"PATNO\") \\\n", " .drop(columnsToDrop, axis=1)\n", "\n", " # Merge together\n", " df = data_df \\\n", " .merge(output_df, on=\"PATNO\") \\\n", " .drop([\"Unnamed: 0\",], axis=1)\n", "\n", " return df\n", "\n", "\n", "def getModel(model):\n", " \"\"\"\n", " Helper function that returns a machine learning model and its parameters for tuning based on the top parameters from Shu et al. models.\n", "\n", " Parameters:\n", " - model (str): The type of machine learning model to be returned. Options: \"SVM\", \"DecisionTree\", \"kNN\", \"GNB\".\n", "\n", " Returns:\n", " - tuple: A tuple containing the machine learning model (classifier) and a dictionary of hyperparameter grid for tuning.\n", " \"\"\"\n", " prefix = \"train__\"\n", " \n", " if model == \"SVM\":\n", " param_grid = {f'{prefix}C': [10, 100, 1000],\n", " f'{prefix}gamma': [0.1, 0.01, 0.001],\n", " f'{prefix}kernel': ['rbf'],\n", " f'{prefix}class_weight': [None, 'balanced']}\n", " clf = SVC(probability=True)\n", " elif model == \"DecisionTree\":\n", " param_grid = {f'{prefix}max_depth': [4, 5, 6, 8],\n", " f'{prefix}max_leaf_nodes': list(range(5, 20, 1)),\n", " f'{prefix}min_samples_split': [5, 8, 16],\n", " f'{prefix}class_weight': [None, 'balanced']}\n", " clf = DecisionTreeClassifier()\n", " elif model == \"kNN\":\n", " param_grid = {f'{prefix}n_neighbors': list(range(5, 20)),\n", " f'{prefix}p': [1, 2],\n", " f'{prefix}weights': [\"uniform\", \"distance\"]}\n", " clf = KNeighborsClassifier()\n", " elif model == \"GNB\":\n", " param_grid = {f'{prefix}var_smoothing': np.logspace(0, -9, num=100)}\n", " clf = GaussianNB()\n", "\n", " param_grid[f'pca__n_components'] = [2, 3, 5, 7, 10]\n", "\n", " return clf, param_grid\n", "\n", "\n", "class MRMRFeatureSelector(BaseEstimator, TransformerMixin):\n", " \"\"\"\n", " Custom scikit-learn transformer for feature selection using MRMR (Maximum Relevance, Minimum Redundancy) algorithm.\n", "\n", " Parameters:\n", " - K (int, optional): The number of top features to select. Defaults to 7.\n", " - n_jobs (int, optional): The number of CPU cores to use for parallel computation. Defaults to -1 (using all available cores).\n", " \"\"\"\n", " def __init__(self, K=7, n_jobs=-1):\n", " self.K = K\n", " self.n_jobs = n_jobs\n", "\n", " def fit(self, X, y):\n", " \"\"\"\n", " Fit the MRMR feature selector to the input data.\n", "\n", " Parameters:\n", " - X (pandas.DataFrame): The feature matrix.\n", " - y (pandas.Series): The target variable.\n", "\n", " Returns:\n", " - self (MRMRFeatureSelector): The fitted transformer object.\n", " \"\"\"\n", " self.selected_features_ = mrmr_classif(X, y, K=self.K, n_jobs=self.n_jobs)\n", " return self\n", "\n", " def transform(self, X):\n", " \"\"\"\n", " Transform the input data by selecting the relevant features.\n", "\n", " Parameters:\n", " - X (pandas.DataFrame): The feature matrix.\n", "\n", " Returns:\n", " - pandas.DataFrame: The transformed dataframe containing only the selected features.\n", " \"\"\"\n", " return X[self.selected_features_]\n", "\n", "\n", "def nested_cross_validation(df, clf, param_grid, usePCA=True):\n", " \"\"\"\n", " Perform nested cross-validation for model evaluation.\n", "\n", " Parameters:\n", " - df (pandas.DataFrame): The input dataframe containing features and target variable.\n", " - clf (BaseEstimator): The machine learning classifier to be evaluated.\n", " - param_grid (dict): The hyperparameter grid for the grid search.\n", " - usePCA (bool, optional): Flag indicating whether to include PCA in the pipeline. Defaults to True.\n", "\n", " Returns:\n", " - tuple: A tuple containing the fitted grid search object and an array of ROC AUC scores from outer cross-validation.\n", " \"\"\"\n", " # Get data\n", " features = df.loc[:, ~df.columns.isin(['PATNO', 'output'])].columns\n", " X = df[features]\n", " y = df[\"output\"]\n", "\n", " # Define outer and inner cross-validation\n", " cv_outer = StratifiedKFold(n_splits=7, shuffle=True, random_state=42)\n", " cv_inner = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)\n", "\n", " # Define pipeline with or without PCA\n", " if usePCA:\n", " pipeline = Pipeline([\n", " ('scale', StandardScaler()),\n", " ('pca', PCA()),\n", " ('train', clf)\n", " ])\n", " else:\n", " del param_grid[\"pca__n_components\"]\n", " pipeline = Pipeline([\n", " ('scale', StandardScaler().set_output(transform=\"pandas\")),\n", " ('mrmr', MRMRFeatureSelector(K=7)),\n", " ('train', clf)\n", " ])\n", "\n", " # Train model using grid search\n", " grid = GridSearchCV(pipeline, param_grid, scoring='roc_auc', cv=cv_inner, n_jobs=1, verbose=0, refit=True)\n", " scores = cross_val_score(grid, X, y, scoring='roc_auc', cv=cv_outer, n_jobs=-1, verbose=0)\n", "\n", " return grid, scores\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Training loop" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def run():\n", " \"\"\"\n", " Perform nested cross-validation for multiple cohorts and machine learning models, recording results in a dataframe.\n", "\n", " Returns:\n", " - pandas.DataFrame: A dataframe containing the mean and standard deviation of ROC AUC scores for each cohort,\n", " model type, and feature selection method.\n", " \"\"\"\n", " results_df = pd.DataFrame(columns=[\"SVM\", \"Decision Tree\", \"kNN\", \"GNB\"])\n", "\n", " for cohort in cohorts:\n", " # Skip missing files\n", " if not check_data(cohort):\n", " print(f\"Skipping {cohort}\")\n", " continue\n", "\n", " # Get data\n", " df = get_df(data, f\"{cohort}/cohort.csv\", glob.glob(f\"{cohort}/demographics*\")[0])\n", " cohortName = cohort.split(\"/\")[-1]\n", "\n", " # Define features\n", " features = df.loc[:, ~df.columns.isin(['PATNO'])].columns\n", " df_features = df[features]\n", " \n", " print(f\"COHORT: {cohortName}\\n========================================\")\n", " for MODEL_TYPE in [\"SVM\", \"DecisionTree\", \"kNN\", \"GNB\"]:\n", " print(f\"\\tTraining the following model: {MODEL_TYPE}\")\n", "\n", " # Define classifier\n", " clf, param_grid = getModel(MODEL_TYPE)\n", "\n", " # Train\n", " _, scores_pca = nested_cross_validation(df_features, clf, param_grid, usePCA=True)\n", " _, scores_mrmr = nested_cross_validation(df_features, clf, param_grid, usePCA=False)\n", " \n", " if MODEL_TYPE == \"SVM\":\n", " results_df.loc[f\"PCA_{cohortName}\", \"SVM\"] = f\"{round(scores_pca.mean(), 3)} +/- {round(scores_pca.std(), 3)}\"\n", " if MODEL_TYPE == \"DecisionTree\":\n", " results_df.loc[f\"PCA_{cohortName}\", \"Decision Tree\"] = f\"{round(scores_pca.mean(), 3)} +/- {round(scores_pca.std(), 3)}\"\n", " if MODEL_TYPE == \"kNN\":\n", " results_df.loc[f\"PCA_{cohortName}\", \"kNN\"] = f\"{round(scores_pca.mean(), 3)} +/- {round(scores_pca.std(), 3)}\"\n", " if MODEL_TYPE == \"GNB\":\n", " results_df.loc[f\"PCA_{cohortName}\", \"GNB\"] = f\"{round(scores_pca.mean(), 3)} +/- {round(scores_pca.std(), 3)}\"\n", "\n", " if MODEL_TYPE == \"SVM\":\n", " results_df.loc[f\"MRMR_{cohortName}\", \"SVM\"] = f\"{round(scores_mrmr.mean(), 3)} +/- {round(scores_mrmr.std(), 3)}\"\n", " if MODEL_TYPE == \"DecisionTree\":\n", " results_df.loc[f\"MRMR_{cohortName}\", \"Decision Tree\"] = f\"{round(scores_mrmr.mean(), 3)} +/- {round(scores_mrmr.std(), 3)}\"\n", " if MODEL_TYPE == \"kNN\":\n", " results_df.loc[f\"MRMR_{cohortName}\", \"kNN\"] = f\"{round(scores_mrmr.mean(), 3)} +/- {round(scores_mrmr.std(), 3)}\"\n", " if MODEL_TYPE == \"GNB\":\n", " results_df.loc[f\"MRMR_{cohortName}\", \"GNB\"] = f\"{round(scores_mrmr.mean(), 3)} +/- {round(scores_mrmr.std(), 3)}\"\n", " \n", " return results_df" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "COHORT: cohort_max_off_2_35_matched\n", "========================================\n", "\tTraining the following model: SVM\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "\tTraining the following model: DecisionTree\n", "\tTraining the following model: kNN\n", "\tTraining the following model: GNB\n", "COHORT: cohort_max_off_off_65_75_unmatched\n", "========================================\n", "\tTraining the following model: SVM\n", "\tTraining the following model: DecisionTree\n", "\tTraining the following model: kNN\n", "\tTraining the following model: GNB\n", "COHORT: cohort_max_both_2_59_matched\n", "========================================\n", "\tTraining the following model: SVM\n", "\tTraining the following model: DecisionTree\n", "\tTraining the following model: kNN\n", "\tTraining the following model: GNB\n", "COHORT: cohort_max_both_102_92_unmatched\n", "========================================\n", "\tTraining the following model: SVM\n", "\tTraining the following model: DecisionTree\n", "\tTraining the following model: kNN\n", "\tTraining the following model: GNB\n", "COHORT: cohort_max_on_2_43_matched\n", "========================================\n", "\tTraining the following model: SVM\n", "\tTraining the following model: DecisionTree\n", "\tTraining the following model: kNN\n", "\tTraining the following model: GNB\n", "COHORT: cohort_max_on_on_66_81_unmatched\n", "========================================\n", "\tTraining the following model: SVM\n", "\tTraining the following model: DecisionTree\n", "\tTraining the following model: kNN\n", "\tTraining the following model: GNB\n", "COHORT: cohort_max_2_72_matched\n", "========================================\n", "\tTraining the following model: SVM\n", "\tTraining the following model: DecisionTree\n", "\tTraining the following model: kNN\n", "\tTraining the following model: GNB\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
SVMDecision TreekNNGNB
PCA_cohort_max_off_2_35_matched0.617 +/- 0.2440.597 +/- 0.1980.529 +/- 0.1840.566 +/- 0.181
PCA_cohort_max_off_off_65_75_unmatched0.541 +/- 0.1370.6 +/- 0.1420.526 +/- 0.1860.527 +/- 0.098
PCA_cohort_max_both_2_59_matched0.586 +/- 0.1890.515 +/- 0.0890.502 +/- 0.1520.415 +/- 0.109
PCA_cohort_max_both_102_92_unmatched0.458 +/- 0.1120.553 +/- 0.0910.504 +/- 0.060.553 +/- 0.134
PCA_cohort_max_on_2_43_matched0.506 +/- 0.1730.495 +/- 0.1460.545 +/- 0.140.486 +/- 0.192
PCA_cohort_max_on_on_66_81_unmatched0.531 +/- 0.0890.519 +/- 0.070.518 +/- 0.0660.48 +/- 0.134
PCA_cohort_max_2_72_matched0.405 +/- 0.0990.461 +/- 0.1050.467 +/- 0.1050.509 +/- 0.082
\n", "
" ], "text/plain": [ " SVM Decision Tree \\\n", "PCA_cohort_max_off_2_35_matched 0.617 +/- 0.244 0.597 +/- 0.198 \n", "PCA_cohort_max_off_off_65_75_unmatched 0.541 +/- 0.137 0.6 +/- 0.142 \n", "PCA_cohort_max_both_2_59_matched 0.586 +/- 0.189 0.515 +/- 0.089 \n", "PCA_cohort_max_both_102_92_unmatched 0.458 +/- 0.112 0.553 +/- 0.091 \n", "PCA_cohort_max_on_2_43_matched 0.506 +/- 0.173 0.495 +/- 0.146 \n", "PCA_cohort_max_on_on_66_81_unmatched 0.531 +/- 0.089 0.519 +/- 0.07 \n", "PCA_cohort_max_2_72_matched 0.405 +/- 0.099 0.461 +/- 0.105 \n", "\n", " kNN GNB \n", "PCA_cohort_max_off_2_35_matched 0.529 +/- 0.184 0.566 +/- 0.181 \n", "PCA_cohort_max_off_off_65_75_unmatched 0.526 +/- 0.186 0.527 +/- 0.098 \n", "PCA_cohort_max_both_2_59_matched 0.502 +/- 0.152 0.415 +/- 0.109 \n", "PCA_cohort_max_both_102_92_unmatched 0.504 +/- 0.06 0.553 +/- 0.134 \n", "PCA_cohort_max_on_2_43_matched 0.545 +/- 0.14 0.486 +/- 0.192 \n", "PCA_cohort_max_on_on_66_81_unmatched 0.518 +/- 0.066 0.48 +/- 0.134 \n", "PCA_cohort_max_2_72_matched 0.467 +/- 0.105 0.509 +/- 0.082 " ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "results_df = run()\n", "results_df" ] } ], "metadata": { "kernelspec": { "display_name": "research", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.4" }, "orig_nbformat": 4 }, "nbformat": 4, "nbformat_minor": 2 }