[Colab]CatBoost.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/ByungjunKim/EnergyTransitionKorea/blob/main/%5BColab%5DCatBoost.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "8mSeVBe8vXZB"
      },
      "source": [
        "# CatBoost\n",
        "This code contains the CatBoost modeling and plotting visualization process used in the paper."
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "!pip install -q catboost ipywidgets shap\n",
        "!jupyter nbextension enable --py widgetsnbextension"
      ],
      "metadata": {
        "id": "U2J7bSpEvhd7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7hW0YT48vXZD"
      },
      "outputs": [],
      "source": [
        "from catboost import CatBoostRegressor, Pool\n",
        "from sklearn.model_selection import train_test_split\n",
        "import pandas as pd\n",
        "import numpy as np\n",
        "from sklearn.metrics import r2_score\n",
        "# from sklearn.metrics import mean_absolute_error\n",
        "from sklearn.metrics import mean_absolute_percentage_error\n",
        "from sklearn.metrics import mean_squared_error\n",
        "import shap\n",
        "import matplotlib.pyplot as plt\n",
        "plt.rcParams[\"figure.dpi\"] = 200 # Increase to high DPI resolution.\n",
        "# plt.rcParams['font.family'] = 'NanumGothic'\n",
        "plt.rcParams['font.family'] = 'DejaVu Sans'\n",
        "plt.rcParams['font.sans-serif'] = ['sans-serif', 'DejaVu Sans', 'sans']\n",
        "plt.rcParams['axes.unicode_minus'] = False\n",
        "# plt.rcParams['text.usetex'] = True\n",
        "shap.initjs()"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Clone\n",
        "!git clone https://github.com/ByungjunKim/EnergyTransitionKorea.git"
      ],
      "metadata": {
        "id": "x-sMP_zCyPsS"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "R9elKm9ovXZH"
      },
      "outputs": [],
      "source": [
        "# Load the comment data from news articles related to energy transition.\n",
        "df = pd.read_csv('./EnergyTransitionKorea/data/reply_df_catboost.csv')\n",
        "df.head()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eAQkfX-vvXZH"
      },
      "outputs": [],
      "source": [
        "df.columns"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "yr0TwpgcvXZI"
      },
      "outputs": [],
      "source": [
        "# dropna\n",
        "df = df.dropna().reset_index(drop=True)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iYqxBe0zvXZI"
      },
      "outputs": [],
      "source": [
        "# Categorize news article titles.\n",
        "df.loc[df['title_p/n']==-1.0,'title_p/n'] = 'negative'\n",
        "df.loc[df['title_p/n']==0.0,'title_p/n'] = 'neutral'\n",
        "df.loc[df['title_p/n']==1.0,'title_p/n'] = 'positive'\n",
        "df['title_p/n'] = df['title_p/n'].astype('category')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "kqm-Uaz8vXZI"
      },
      "outputs": [],
      "source": [
        "# Categorize news article policy.\n",
        "df.loc[df['e_policy_p/n']==-1.0,'e_policy_p/n'] = 'negative'\n",
        "df.loc[df['e_policy_p/n']==0.0,'e_policy_p/n'] = 'neutral'\n",
        "df.loc[df['e_policy_p/n']==1.0,'e_policy_p/n'] = 'positive'\n",
        "df['e_policy_p/n'] = df['e_policy_p/n'].astype('category')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "fhnWxD6NvXZI"
      },
      "outputs": [],
      "source": [
        "# Categorize user type (heavy or regular) as a nominal variable.\n",
        "df['user_type'] = df['user_type'].astype('category')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zEbxXB8ivXZI"
      },
      "source": [
        "### Model 1: Dependent variable \"Condemning\""
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "J25v9ZkPvXZI"
      },
      "outputs": [],
      "source": [
        "# Creating data features.\n",
        "X_1 = df[['user_type', 'title_p/n', 'e_policy_p/n',\n",
        "         'politics', 'economy','society', 'culture', 'international', 'sports', 'IT_science', \\\n",
        "        'Condemning_past','Praising_past', 'Suffering_past', 'Self-Conscious_past',\n",
        "        'tokens_len_past', 'ttr_past']]\n",
        "y_1 = df['Condemning'].tolist()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HrspaBZGvXZJ"
      },
      "outputs": [],
      "source": [
        "# Split into training and test sets\n",
        "X_train_1, X_test_1, y_train_1, y_test_1 = train_test_split(X_1, y_1, test_size=0.2, random_state=2023)\n",
        "\n",
        "# Specify indices for categorical variables (0,1,2)\n",
        "cat_features = [0,1,2]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "jiHEiO94vXZJ"
      },
      "outputs": [],
      "source": [
        "# Convert to CatBoost Pool format\n",
        "train_pool_1 = Pool(X_train_1, y_train_1, cat_features=cat_features)\n",
        "test_pool_1 = Pool(X_test_1, y_test_1, cat_features=cat_features)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "PXa3sXpEvXZK"
      },
      "outputs": [],
      "source": [
        "# Model initialization\n",
        "model_1 = CatBoostRegressor(iterations=1000, learning_rate=0.1, verbose=200)"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# Train\n",
        "model_1.fit(train_pool_1, eval_set=test_pool_1)"
      ],
      "metadata": {
        "id": "NehIyrV_yixr"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "A95Z_tyXvXZK"
      },
      "outputs": [],
      "source": [
        "# save model\n",
        "# model_1.save_model('./catboost_model/condemning.model')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3tMTpgxxvXZM"
      },
      "outputs": [],
      "source": [
        "# Predict\n",
        "y_pred_1 = model_1.predict(X_test_1)\n",
        "\n",
        "# R-squared\n",
        "print(\"R-squared: {:.3f}\".format(r2_score(y_test_1, y_pred_1)))\n",
        "\n",
        "# Adjusted R2\n",
        "n = X_test_1.shape[0]\n",
        "k = X_test_1.shape[1]\n",
        "adjusted_r2 = 1 - (1 - r2_score(y_test_1, y_pred_1)) * (n - 1) / (n - k - 1)\n",
        "print(\"Adj R-squared: {:.3f}\".format(adjusted_r2))\n",
        "\n",
        "# MAPE\n",
        "print(\"MAPE: {:.3f}\".format(mean_absolute_percentage_error(y_test_1, y_pred_1)))\n",
        "\n",
        "# Normalized MAE\n",
        "# Calculate MAE\n",
        "mae = np.mean(np.abs(y_test_1 - y_pred_1))\n",
        "# Normalize MAE by the range of the dependent variable (Max - Min)\n",
        "normalized_mae_range = mae / (np.max(y_test_1) - np.min(y_test_1))\n",
        "\n",
        "# Normalize MAE by the standard deviation of the dependent variable\n",
        "normalized_mae_std = mae / np.std(y_test_1)\n",
        "print(\"Normalized MAE: {:.3f}\".format(normalized_mae_range))\n",
        "print(\"Normalized_std MAE: {:.3f}\".format(normalized_mae_std))\n",
        "\n",
        "# Calculate Coefficient of Variation of the RMSE (CV-RMSE)\n",
        "print(\"CV-RMSE: {:.3f}\".format(mean_squared_error(y_test_1, y_pred_1, squared=False) / np.mean(y_test_1)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "3R6JbNW4vXZM"
      },
      "outputs": [],
      "source": [
        "# Get and plot SHAP values\n",
        "explainer_1 = shap.Explainer(model_1)\n",
        "shap_values_1 = explainer_1(X_1)\n",
        "shap.summary_plot(shap_values_1, X_1,show=False)\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"a)\", fontsize=15)\n",
        "title.set_position([-0.5,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "dMt-e9qHvXZN"
      },
      "outputs": [],
      "source": [
        "# https://github.com/conorosully/medium-articles/blob/master/src/interpretable%20ml/SHAP/SHAP_catboost.ipynb\n",
        "\n",
        "#Create for placeholder SHAP values\n",
        "shap_values_cate_1 = explainer_1(X_1)\n",
        "\n",
        "#Get shaply values and feature values for odor\n",
        "odor_values = np.array(shap_values_1[:,'title_p/n'].values)\n",
        "odor_data = X_1['title_p/n']\n",
        "\n",
        "#Create new SHAP values array\n",
        "\n",
        "#Split odor SHAP values by unique odor categories\n",
        "unique_odor = list(X_1['title_p/n'].unique())\n",
        "new_shap_values = [np.array(pd.Series(odor_values)[odor_data==odor]) for odor in unique_odor]\n",
        "\n",
        "#Each sublist needs to be the same length\n",
        "max_len = max([len(v) for v in new_shap_values])\n",
        "new_shap_values = [np.append(vs,[np.nan]*(max_len - len(vs))) for vs in new_shap_values]\n",
        "new_shap_values = np.array(new_shap_values)\n",
        "\n",
        "#transpost matrix so categories are columns and SHAP values are rows\n",
        "new_shap_values = new_shap_values.transpose()\n",
        "\n",
        "#replace shap values\n",
        "shap_values_cate_1.values = np.array(new_shap_values)\n",
        "\n",
        "#replace data with placeholder array\n",
        "shap_values_cate_1.data = np.array([[0]*len(unique_odor)]*max_len)\n",
        "\n",
        "#replace base data with placeholder array\n",
        "shap_values_cate_1.base = np.array([0]*max_len)\n",
        "\n",
        "#replace feature names with category labels\n",
        "# odor_labels = {'a':'almond',\n",
        "#                'l':'anise',\n",
        "#                'c':'creosote',\n",
        "#                'y':'fishy',\n",
        "#                'f':'foul',\n",
        "#                'm':'musty',\n",
        "#                'n':'none',\n",
        "#                'p':'pungent',\n",
        "#                's':'spicy'}\n",
        "# labels = [\"{} ({})\".format(odor_labels[u],u) for u in unique_odor]\n",
        "shap_values_cate_1.feature_names = unique_odor\n",
        "\n",
        "#Use besswarm as before\n",
        "shap.plots.beeswarm(shap_values_cate_1, color_bar=False,show=False, color='#808080')\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"b)\", fontsize=15)\n",
        "title.set_position([-0.18,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UtmpyZgMvXZN"
      },
      "outputs": [],
      "source": [
        "#get shaply values and data\n",
        "# odor_values = shap_values_1[:,'title_p/n'].values\n",
        "# odor_data = X_1['title_p/n']\n",
        "# unique_odor = set(X['title_p/n'])\n",
        "\n",
        "#split odor shap values based on odor category\n",
        "odor_categories = sorted(list(set(odor_data)))\n",
        "\n",
        "odor_groups = []\n",
        "for o in odor_categories:\n",
        "    relevant_values = odor_values[odor_data == o]\n",
        "    odor_groups.append(relevant_values)\n",
        "\n",
        "# #replace categories with labels\n",
        "# odor_labels = {'a':'almond',\n",
        "#                'l':'anise',\n",
        "#                'c':'creosote',\n",
        "#                'y':'fishy',\n",
        "#                'f':'foul',\n",
        "#                'm':'musty',\n",
        "#                'n':'none',\n",
        "#                'p':'pungent',\n",
        "#                's':'spicy'}\n",
        "\n",
        "# labels = [odor_labels[u] for u in unique_odor]\n",
        "\n",
        "#plot boxplot\n",
        "plt.figure(figsize=(8, 5))\n",
        "\n",
        "plt.boxplot(odor_groups,labels=odor_categories)\n",
        "\n",
        "plt.ylabel('SHAP values (Condemning)',size=15)\n",
        "plt.xlabel('title_p/n',size=15)\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"c)\", fontsize=15)\n",
        "title.set_position([-0.13,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fVigkGCOvXZO"
      },
      "source": [
        "### Model 2: Dependent variable \"token_len\""
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "bPLc3caOvXZO"
      },
      "outputs": [],
      "source": [
        "# Creating data features.\n",
        "X_2 = df[['user_type', 'title_p/n', 'e_policy_p/n',\n",
        "        'politics', 'economy','society', 'culture', 'international', 'sports', 'IT_science', \\\n",
        "        # 'Condemning',\n",
        "        # 'Praising', 'Suffering','Self-Conscious',\\\n",
        "        'Condemning_past','Praising_past', 'Suffering_past', 'Self-Conscious_past',\n",
        "        'tokens_len_past', 'ttr_past']]\n",
        "y_2 = df['tokens_len']"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OlnrFt2RvXZO"
      },
      "outputs": [],
      "source": [
        "# Splitting data into training and testing sets\n",
        "X_train_2, X_test_2, y_train_2, y_test_2 = train_test_split(X_2, y_2, test_size=0.2, random_state=2023)\n",
        "\n",
        "# Specifying indices of categorical features\n",
        "cat_features = [0,1,2]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5YF27gdGvXZO"
      },
      "outputs": [],
      "source": [
        "# Convert to CatBoost Pool format\n",
        "train_pool_2 = Pool(X_train_2, y_train_2, cat_features=cat_features)\n",
        "test_pool_2 = Pool(X_test_2, y_test_2, cat_features=cat_features)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "MS37sn4GvXZO"
      },
      "outputs": [],
      "source": [
        "# Model initialization\n",
        "model_2 = CatBoostRegressor(iterations=1000, learning_rate=0.1, verbose=200)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "ntLjkmZKvXZO"
      },
      "outputs": [],
      "source": [
        "# Train\n",
        "model_2.fit(train_pool_2, eval_set=test_pool_2)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "WMtigC10vXZP"
      },
      "outputs": [],
      "source": [
        "# model save\n",
        "# model_2.save_model('./catboost_model/token_len.model')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "iGYqExAsvXZP"
      },
      "outputs": [],
      "source": [
        "# Predict\n",
        "y_pred_2 = model_2.predict(X_test_2)\n",
        "\n",
        "# R-squared\n",
        "print(\"R-squared: {:.3f}\".format(r2_score(y_test_2, y_pred_2)))\n",
        "\n",
        "# Adjusted R2\n",
        "n = X_test_2.shape[0]\n",
        "k = X_test_2.shape[1]\n",
        "adjusted_r2 = 1 - (1 - r2_score(y_test_2, y_pred_2)) * (n - 1) / (n - k - 1)\n",
        "print(\"Adj R-squared: {:.3f}\".format(adjusted_r2))\n",
        "\n",
        "# MAPE\n",
        "print(\"MAPE: {:.3f}\".format(mean_absolute_percentage_error(y_test_2, y_pred_2)))\n",
        "\n",
        "# Normalized MAE\n",
        "# Calculate MAE\n",
        "mae = np.mean(np.abs(y_test_2 - y_pred_2))\n",
        "# Normalize MAE by the range of the dependent variable (Max - Min)\n",
        "normalized_mae_range = mae / (np.max(y_test_2) - np.min(y_test_2))\n",
        "\n",
        "# Normalize MAE by the standard deviation of the dependent variable\n",
        "normalized_mae_std = mae / np.std(y_test_2)\n",
        "print(\"Normalized MAE: {:.3f}\".format(normalized_mae_range))\n",
        "print(\"Normalized_std MAE: {:.3f}\".format(normalized_mae_std))\n",
        "\n",
        "# Calculate Coefficient of Variation of the RMSE (CV-RMSE)\n",
        "print(\"CV-RMSE: {:.3f}\".format(mean_squared_error(y_test_2, y_pred_2, squared=False) / np.mean(y_test_2)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xdsFARRtvXZP"
      },
      "outputs": [],
      "source": [
        "# Get and plot SHAP values\n",
        "explainer_2 = shap.Explainer(model_2)\n",
        "shap_values_2 = explainer_2(X_2)\n",
        "shap.summary_plot(shap_values_2, X_2,show=False)\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"a)\", fontsize=15)\n",
        "title.set_position([-0.5,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "mesFFkGSvXZQ"
      },
      "outputs": [],
      "source": [
        "#Create for placeholder SHAP values\n",
        "shap_values_cate_2 = explainer_2(X_2)\n",
        "\n",
        "#Get shaply values and feature values for odor\n",
        "odor_values = np.array(shap_values_2[:,'e_policy_p/n'].values)\n",
        "odor_data = X_2['e_policy_p/n']\n",
        "\n",
        "#Create new SHAP values array\n",
        "\n",
        "#Split odor SHAP values by unique odor categories\n",
        "unique_odor = list(X_2['e_policy_p/n'].unique())\n",
        "new_shap_values = [np.array(pd.Series(odor_values)[odor_data==odor]) for odor in unique_odor]\n",
        "\n",
        "#Each sublist needs to be the same length\n",
        "max_len = max([len(v) for v in new_shap_values])\n",
        "new_shap_values = [np.append(vs,[np.nan]*(max_len - len(vs))) for vs in new_shap_values]\n",
        "new_shap_values = np.array(new_shap_values)\n",
        "\n",
        "#transpost matrix so categories are columns and SHAP values are rows\n",
        "new_shap_values = new_shap_values.transpose()\n",
        "\n",
        "#replace shap values\n",
        "shap_values_cate_2.values = np.array(new_shap_values)\n",
        "\n",
        "#replace data with placeholder array\n",
        "shap_values_cate_2.data = np.array([[0]*len(unique_odor)]*max_len)\n",
        "\n",
        "#replace base data with placeholder array\n",
        "shap_values_cate_2.base = np.array([0]*max_len)\n",
        "\n",
        "#replace feature names with category labels\n",
        "# odor_labels = {'a':'almond',\n",
        "#                'l':'anise',\n",
        "#                'c':'creosote',\n",
        "#                'y':'fishy',\n",
        "#                'f':'foul',\n",
        "#                'm':'musty',\n",
        "#                'n':'none',\n",
        "#                'p':'pungent',\n",
        "#                's':'spicy'}\n",
        "# labels = [\"{} ({})\".format(odor_labels[u],u) for u in unique_odor]\n",
        "shap_values_cate_2.feature_names = unique_odor\n",
        "\n",
        "#Use besswarm as before\n",
        "shap.plots.beeswarm(shap_values_cate_2, color_bar=False,show=False,color='#808080')\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"b)\", fontsize=15)\n",
        "title.set_position([-0.18,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "HeU-v9mJvXZQ"
      },
      "outputs": [],
      "source": [
        "#get shaply values and data\n",
        "# odor_values = shap_values[:,'e_policy_p/n'].values\n",
        "# odor_data = X['e_policy_p/n']\n",
        "# unique_odor = set(X['e_policy_p/n'])\n",
        "\n",
        "#split odor shap values based on odor category\n",
        "odor_categories = sorted(list(set(odor_data)))\n",
        "\n",
        "odor_groups = []\n",
        "for o in odor_categories:\n",
        "    relevant_values = odor_values[odor_data == o]\n",
        "    odor_groups.append(relevant_values)\n",
        "\n",
        "# #replace categories with labels\n",
        "# odor_labels = {'a':'almond',\n",
        "#                'l':'anise',\n",
        "#                'c':'creosote',\n",
        "#                'y':'fishy',\n",
        "#                'f':'foul',\n",
        "#                'm':'musty',\n",
        "#                'n':'none',\n",
        "#                'p':'pungent',\n",
        "#                's':'spicy'}\n",
        "\n",
        "# labels = [odor_labels[u] for u in unique_odor]\n",
        "\n",
        "#plot boxplot\n",
        "plt.figure(figsize=(8, 5))\n",
        "\n",
        "plt.boxplot(odor_groups,labels=odor_categories)\n",
        "\n",
        "plt.ylabel('SHAP values (Token Len)',size=15)\n",
        "plt.xlabel('e_policy_p/n',size=15)\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"c)\", fontsize=15)\n",
        "title.set_position([-0.13,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5MzkBKuhvXZQ"
      },
      "source": [
        "### Model 3: Dependent variable \"ttr\""
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "OnzCK9ijvXZQ"
      },
      "outputs": [],
      "source": [
        "# Creating data features.\n",
        "X_3 = df[['user_type', 'title_p/n', 'e_policy_p/n', 'politics', 'economy',\\\n",
        "       'society', 'culture', 'international', 'sports', 'IT_science', \\\n",
        "        # 'Praising', 'Suffering','Self-Conscious',\\\n",
        "        'Condemning_past','Praising_past', 'Suffering_past', 'Self-Conscious_past',\n",
        "        'tokens_len_past', 'ttr_past']]\n",
        "y_3 = df['ttr']"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "w8Gh21-zvXZR"
      },
      "outputs": [],
      "source": [
        "# Splitting data into training and testing sets\n",
        "X_train_3, X_test_3, y_train_3, y_test_3 = train_test_split(X_3, y_3, test_size=0.2, random_state=2023)\n",
        "\n",
        "# Specifying indices of categorical features\n",
        "cat_features = [0,1,2]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "5rcBr3CpvXZR"
      },
      "outputs": [],
      "source": [
        "# Convert to CatBoost Pool format.\n",
        "train_pool_3 = Pool(X_train_3, y_train_3, cat_features=cat_features)\n",
        "test_pool_3 = Pool(X_test_3, y_test_3, cat_features=cat_features)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "K_2QusicvXZR"
      },
      "outputs": [],
      "source": [
        "# Model initialization\n",
        "model_3 = CatBoostRegressor(iterations=1000, learning_rate=0.1, verbose=200)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uD5BVjuivXZR"
      },
      "outputs": [],
      "source": [
        "# Train\n",
        "model_3.fit(train_pool_3, eval_set=test_pool_3)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4BYFHKSHvXZU"
      },
      "outputs": [],
      "source": [
        "# model save\n",
        "# model_3.save_model('./catboost_model/token_ttr.model')"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "lXDrcxUdvXZU"
      },
      "outputs": [],
      "source": [
        "# Predict\n",
        "y_pred_3 = model_3.predict(X_test_3)\n",
        "\n",
        "# R-squared\n",
        "print(\"R-squared: {:.3f}\".format(r2_score(y_test_3, y_pred_3)))\n",
        "\n",
        "# Adjusted R2\n",
        "n = X_test_3.shape[0]\n",
        "k = X_test_3.shape[1]\n",
        "adjusted_r2 = 1 - (1 - r2_score(y_test_3, y_pred_3)) * (n - 1) / (n - k - 1)\n",
        "print(\"Adj R-squared: {:.3f}\".format(adjusted_r2))\n",
        "\n",
        "# MAPE\n",
        "print(\"MAPE: {:.3f}\".format(mean_absolute_percentage_error(y_test_3, y_pred_3)))\n",
        "\n",
        "# Normalized MAE\n",
        "# Calculate MAE\n",
        "mae = np.mean(np.abs(y_test_3 - y_pred_3))\n",
        "# Normalize MAE by the range of the dependent variable (Max - Min)\n",
        "normalized_mae_range = mae / (np.max(y_test_3) - np.min(y_test_3))\n",
        "\n",
        "# Normalize MAE by the standard deviation of the dependent variable\n",
        "normalized_mae_std = mae / np.std(y_test_3)\n",
        "print(\"Normalized MAE: {:.3f}\".format(normalized_mae_range))\n",
        "print(\"Normalized_std MAE: {:.3f}\".format(normalized_mae_std))\n",
        "\n",
        "# Calculate Coefficient of Variation of the RMSE (CV-RMSE)\n",
        "print(\"CV-RMSE: {:.3f}\".format(mean_squared_error(y_test_3, y_pred_3, squared=False) / np.mean(y_test_3)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "r0p7tTiYvXZU"
      },
      "outputs": [],
      "source": [
        "# Get and plot SHAP values\n",
        "explainer_3 = shap.Explainer(model_3)\n",
        "shap_values_3 = explainer_3(X_3)\n",
        "shap.summary_plot(shap_values_3, X_3, show=False)\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"a)\", fontsize=15)\n",
        "title.set_position([-0.5,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "leN6f6LKvXZV"
      },
      "outputs": [],
      "source": [
        "#Create for placeholder SHAP values\n",
        "shap_values_cate_3 = explainer_3(X_3)\n",
        "\n",
        "#Get shaply values and feature values for odor\n",
        "odor_values = np.array(shap_values_3[:,'e_policy_p/n'].values)\n",
        "odor_data = X_3['e_policy_p/n']\n",
        "\n",
        "#Create new SHAP values array\n",
        "\n",
        "#Split odor SHAP values by unique odor categories\n",
        "unique_odor = list(X_3['e_policy_p/n'].unique())\n",
        "new_shap_values = [np.array(pd.Series(odor_values)[odor_data==odor]) for odor in unique_odor]\n",
        "\n",
        "#Each sublist needs to be the same length\n",
        "max_len = max([len(v) for v in new_shap_values])\n",
        "new_shap_values = [np.append(vs,[np.nan]*(max_len - len(vs))) for vs in new_shap_values]\n",
        "new_shap_values = np.array(new_shap_values)\n",
        "\n",
        "#transpost matrix so categories are columns and SHAP values are rows\n",
        "new_shap_values = new_shap_values.transpose()\n",
        "\n",
        "#replace shap values\n",
        "shap_values_cate_3.values = np.array(new_shap_values)\n",
        "\n",
        "#replace data with placeholder array\n",
        "shap_values_cate_3.data = np.array([[0]*len(unique_odor)]*max_len)\n",
        "\n",
        "#replace base data with placeholder array\n",
        "shap_values_cate_3.base = np.array([0]*max_len)\n",
        "\n",
        "#replace feature names with category labels\n",
        "# odor_labels = {'a':'almond',\n",
        "#                'l':'anise',\n",
        "#                'c':'creosote',\n",
        "#                'y':'fishy',\n",
        "#                'f':'foul',\n",
        "#                'm':'musty',\n",
        "#                'n':'none',\n",
        "#                'p':'pungent',\n",
        "#                's':'spicy'}\n",
        "# labels = [\"{} ({})\".format(odor_labels[u],u) for u in unique_odor]\n",
        "shap_values_cate_3.feature_names = unique_odor\n",
        "\n",
        "#Use besswarm as before\n",
        "shap.plots.beeswarm(shap_values_cate_3, color_bar=False,show=False,color='#808080')\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"b)\", fontsize=15)\n",
        "title.set_position([-0.18,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "4vbgfVlvvXZV"
      },
      "outputs": [],
      "source": [
        "#get shaply values and data\n",
        "# odor_values = shap_values[:,'e_policy_p/n'].values\n",
        "# odor_data = X['e_policy_p/n']\n",
        "# unique_odor = set(X['e_policy_p/n'])\n",
        "\n",
        "#split odor shap values based on odor category\n",
        "odor_categories = sorted(list(set(odor_data)))\n",
        "\n",
        "odor_groups = []\n",
        "for o in odor_categories:\n",
        "    relevant_values = odor_values[odor_data == o]\n",
        "    odor_groups.append(relevant_values)\n",
        "\n",
        "# #replace categories with labels\n",
        "# odor_labels = {'a':'almond',\n",
        "#                'l':'anise',\n",
        "#                'c':'creosote',\n",
        "#                'y':'fishy',\n",
        "#                'f':'foul',\n",
        "#                'm':'musty',\n",
        "#                'n':'none',\n",
        "#                'p':'pungent',\n",
        "#                's':'spicy'}\n",
        "\n",
        "# labels = [odor_labels[u] for u in unique_odor]\n",
        "\n",
        "#plot boxplot\n",
        "plt.figure(figsize=(8, 5))\n",
        "\n",
        "plt.boxplot(odor_groups,labels=odor_categories)\n",
        "\n",
        "plt.ylabel('SHAP values (TTR)',size=15)\n",
        "plt.xlabel('e_policy_p/n',size=15)\n",
        "\n",
        "# Setting the title with custom alignment\n",
        "title = plt.title(\"c)\", fontsize=15)\n",
        "title.set_position([-0.13,1]) # You can adjust the [0, 1.02] values as needed\n",
        "title.set_ha('left')\n",
        "\n",
        "plt.show()"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "python311",
      "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.11.6"
    },
    "orig_nbformat": 4,
    "colab": {
      "provenance": [],
      "include_colab_link": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}