Python_code_1_XGBOOST-ctliu2023.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "87d3e580",
   "metadata": {},
   "outputs": [],
   "source": [
    "##python codes for the XGBoost machine learning model\n",
    "#This model can calssify basaltic sample's erupted environment, subaerial or submerged eruption\n",
    "#Here, we extract the major and trace elements from \"Parental training dataset\" to train the XGBoost model\n",
    "\n",
    "\n",
    "\n",
    "from __future__ import division, print_function, unicode_literals\n",
    "\n",
    "# Check if the version of python is 3.5 and above\n",
    "import sys\n",
    "assert sys.version_info >= (3, 5)\n",
    "\n",
    "# Check to see if sklearn is version 0.20 and above\n",
    "import sklearn #\n",
    "assert sklearn.__version__ >= \"0.20\"\n",
    "from sklearn.decomposition import PCA\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt\n",
    "import os \n",
    "from xgboost import XGBClassifier#\n",
    "from sklearn.metrics import accuracy_score, r2_score, make_scorer, f1_score, recall_score, precision_score\n",
    "from sklearn.model_selection import GridSearchCV#\n",
    "from sklearn.model_selection import train_test_split#\n",
    "from sklearn.metrics import roc_auc_score#\n",
    "from sklearn.metrics import roc_curve, auc#\n",
    "from sklearn.metrics import confusion_matrix#\n",
    "import seaborn as sns\n",
    "\n",
    "np.random.seed(2022) \n",
    "\n",
    "# Make matplotlib diagrams work better\n",
    "# matplotlib inline\n",
    "import matplotlib as mpl\n",
    "mpl.rc('axes', labelsize=14)\n",
    "mpl.rc('xtick', labelsize=12)\n",
    "mpl.rc('ytick', labelsize=12)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "487b6b91",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ignoring Unnecessary Warnings\n",
    "import warnings\n",
    "warnings.filterwarnings(action=\"ignore\", message=\"^internal gelsd\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f28446c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# read Excel\n",
    "#Notes\n",
    "#1.put the \"Table-S2.Training dataset for XGBoost modeling\" and this jupyter python code into a same path\n",
    "\n",
    "#2.Please remove the first row (the head) before applying the \"Table-S2.Training dataset for XGBoost modeling\".\n",
    "\n",
    "orig_data = pd.read_excel('Table-S2.Training dataset for XGBoost modeling.xlsx')#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a59e94bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Separate features and labels\n",
    "X = orig_data.drop([\"TRUE_VALUE\"], axis=1).copy()\n",
    "y = orig_data[\"TRUE_VALUE\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "772a18b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(16, 6))\n",
    "sns.countplot(orig_data.TRUE_VALUE, palette=\"Set2\")\n",
    "plt.xticks(rotation=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c8c0e93",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_clf = XGBClassifier()\n",
    "xgb_clf.fit(X_train, y_train)\n",
    "y_pred_xgb = xgb_clf.predict(X_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a5d52e4",
   "metadata": {},
   "outputs": [],
   "source": [
    "#before running our model, check the imbalance problem of our trining dataset at the first\n",
    "from sklearn.model_selection import cross_val_score#\n",
    "\n",
    "# Ten-fold cross validation #\n",
    "scores = cross_val_score(xgb_clf, X, y,\n",
    "                        scoring = \"accuracy\", cv=10,#cv\n",
    "                        n_jobs=-1)\n",
    "def display_scores(scores):\n",
    "\n",
    "    print(\"Scores:\", scores)\n",
    "    print(\"Mean:\", scores.mean())\n",
    "    print(\"Standard deviation:\", scores.std())\n",
    "    \n",
    "display_scores(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4580ad6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Randomly seperate the traning dataset\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)\n",
    "\n",
    "X_train.shape, X_test.shape, y_train.shape, y_test.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efaabf11",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_clf.feature_importances_#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf75fbac",
   "metadata": {},
   "outputs": [],
   "source": [
    "# show feature importance\n",
    "for feature_name, score in zip(list(X.columns), xgb_clf.feature_importances_):\n",
    "    print(feature_name, \":\", score)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e541e7d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "learning_rate = [ 0.1, 0.2, 0.5,0.6,0.7]\n",
    "depth = [3, 4, 5, 6, 7]\n",
    "min_split = [0.1,0.2,0.3,1]#\n",
    "alpha1 = [0.1,0.3,0.5,0.7,0.9, 1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d4cf5a66",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb = XGBClassifier(objective='binary:logistic',\n",
    "                    eval_metric = 'auc', tree_method='hist', seed=2022,importance_type = 'cover')#"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e09f3245",
   "metadata": {},
   "outputs": [],
   "source": [
    "xgb_cv = GridSearchCV(xgb, param_grid = {'eta': learning_rate, 'gamma': min_split, 'max_depth': depth, 'alpha':alpha1}, \n",
    "                      cv=10, scoring='f1') #cv=10: tenfold cross-validation\n",
    "xgb_cv.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49277b28",
   "metadata": {},
   "outputs": [],
   "source": [
    "#calculate a mean value and standard deviation of the tenfold cross-validation scores of a best trained model\n",
    "from sklearn.model_selection import cross_val_score\n",
    "\n",
    "scores = cross_val_score(xgb_cv.best_estimator_, X_train, y_train,\n",
    "                        scoring = \"f1\", cv=10,\n",
    "                        n_jobs=-1)\n",
    "def display_scores(scores):\n",
    "    \n",
    "    print(\"Scores:\", scores)\n",
    "    print(\"Mean:\", scores.mean())\n",
    "    print(\"Standard deviation:\", scores.std())\n",
    "    \n",
    "display_scores(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "cc3a95e7",
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'xgb_cv' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32m~\\AppData\\Local\\Temp/ipykernel_8300/4038436740.py\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[1;31m#show the best score (mean value)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mxgb_cv\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mbest_score_\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m: name 'xgb_cv' is not defined"
     ]
    }
   ],
   "source": [
    "#show the best score (mean value)\n",
    "xgb_cv.best_score_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "18427128",
   "metadata": {},
   "outputs": [],
   "source": [
    "#show parameters of the best model\n",
    "xgb_cv.best_estimator_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be05f05e",
   "metadata": {},
   "outputs": [],
   "source": [
    "features = list(X.columns)\n",
    "importances = xgb_cv.best_estimator_.feature_importances_\n",
    "indices = np.argsort(importances)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d26166d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "df4 = pd.DataFrame({'features':features,'importances':importances})\n",
    "# df4.to_excel('ML_model_Feature_importance.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "427927db",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.pyplot import MultipleLocator\n",
    "plt.barh(range(len(indices)), importances[indices], color='c', align='center')\n",
    "plt.yticks(range(len(indices)), [features[i] for i in indices], fontsize=12)\n",
    "plt.xticks(fontsize=20)\n",
    "plt.xlabel('Relative Importance',fontsize=25)\n",
    "x_major_locator=MultipleLocator(0.02)#\n",
    "ax = plt.gca()\n",
    "ax.xaxis.set_major_locator(x_major_locator)\n",
    "plt.xlim((0,0.12))\n",
    "# plt.savefig('ML_model_Feature_importance.pdf', dpi=500)\n",
    "plt.show()\n",
    "plt.rcParams[\"figure.figsize\"] = (20, 10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78338196",
   "metadata": {},
   "outputs": [],
   "source": [
    "# predict the test data set\n",
    "xgb_test = xgb_cv.best_estimator_\n",
    "xgb_test.fit(X_train,y_train)\n",
    "y_test_pred = xgb_test.predict(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ace042f",
   "metadata": {},
   "outputs": [],
   "source": [
    "print('Accuracy: %.4f' % accuracy_score(y_test, y_test_pred))\n",
    "print('ROC AUC: %.4f' % roc_auc_score(y_test, y_test_pred))\n",
    "print('Precision: %.4f' % precision_score(y_true=y_test, y_pred=y_test_pred))\n",
    "print('Recall: %.4f' % recall_score(y_true=y_test, y_pred=y_test_pred))\n",
    "print('F1 Score: %.4f' % f1_score(y_true=y_test, y_pred=y_test_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8f21e48",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import classification_report\n",
    "\n",
    "print(classification_report(y_true=y_test, y_pred=y_test_pred))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ca94855",
   "metadata": {},
   "outputs": [],
   "source": [
    "confmat = confusion_matrix(y_true=y_test, y_pred=y_test_pred)\n",
    "\n",
    "print(confmat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b2ac52c",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.matshow(confmat, cmap=plt.cm.Blues, alpha=0.3)\n",
    "for i in range(confmat.shape[0]):\n",
    "    for j in range(confmat.shape[1]):\n",
    "        ax.text(x=j, y=i, s=confmat[i,j], va='center', ha='center')\n",
    "plt.xlabel('predicted label',fontsize=20)\n",
    "plt.ylabel('true label',fontsize=20)\n",
    "# plt.savefig('ML_model_confusion_matrix.pdf', dpi=800)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53679079",
   "metadata": {},
   "outputs": [],
   "source": [
    "#predict the application dataset that is non label data\n",
    "\n",
    "#We extract \"Table-S4.Application dataset for prediction\" from parental application dataset\n",
    "\n",
    "#Please remove the first row (the head) before applying the \"Table-S4.Application dataset for prediction\".\n",
    "predict_data = pd.read_excel ('Table-S4.Application dataset for prediction.xlsx')\n",
    "x_predict_data = predict_data # process.fit_transform(predict_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "923e9243",
   "metadata": {},
   "outputs": [],
   "source": [
    "predict_results = xgb_cv.best_estimator_.predict_proba(x_predict_data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1bcfd46c",
   "metadata": {},
   "outputs": [],
   "source": [
    "predict_results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6438063b",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame(data=predict_results)\n",
    "#df.to_excel('predict_results.xlsx')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bde6ffe6",
   "metadata": {},
   "outputs": [],
   "source": [
    "###Xgboost end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}