{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "5b100c3b-c535-62d6-3404-0f815649c061"
      },
      "source": [
        "Matthews correlation coefficient\n",
        "--------------------------------\n",
        " \n",
        "Wiki: [https://en.wikipedia.org/wiki/Matthews_correlation_coefficient][1]\n",
        "\n",
        "Author: [Vopani][2]   \n",
        "R-implementation of [CPMP's fast mcc optimization][3] originally from [Anokas mcc optimization code][4].\n",
        "   \n",
        "This code takes true values and predicted probabilities, and selects the threshold/cutoff that yields the best MCC score. It is slower than the Python implementation, maybe someone can optimize this further, before trying to use it in xgboost or elsewhere.\n",
        "\n",
        "  [1]: https://en.wikipedia.org/wiki/Matthews_correlation_coefficient\n",
        "  [2]: https://www.kaggle.com/rohanrao\n",
        "  [3]: https://www.kaggle.com/cpmpml/bosch-production-line-performance/optimizing-probabilities-for-best-mcc\n",
        "  [4]: https://www.kaggle.com/c/bosch-production-line-performance/forums/t/22917/optimising-probabilities-binary-prediction-script"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "316f7880-2e25-f8b4-f8fb-de98a8d67828"
      },
      "outputs": [],
      "source": [
        "# This R environment comes with all of CRAN preinstalled, as well as many other helpful packages\n",
        "# The environment is defined by the kaggle/rstats docker image: https://github.com/kaggle/docker-rstats"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "7a3c32ba-0f06-3312-0d89-2dce728a8a31"
      },
      "source": [
        "**Evaluating MCC**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ec03e602-be87-753f-1c8e-6e9ab09eabb2"
      },
      "outputs": [],
      "source": [
        "mcc <- function(TP, FP, FN, TN)\n",
        "{\n",
        "    num <- (TP*TN) - (FP*FN)\n",
        "    den <- (TP+FP)*(TP+FN)*(TN+FP)*(TN+FN)\n",
        "\n",
        "    if (den == 0)\n",
        "    {\n",
        "      return(0)\n",
        "    }else\n",
        "    {\n",
        "      return(num / sqrt(den))\n",
        "    }\n",
        "}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "3a582b82-7f45-9aeb-eb51-404f3699c558"
      },
      "source": [
        "**Optimizing MCC**"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "e85d60bd-766b-b490-8180-66eaa60f7374"
      },
      "outputs": [],
      "source": [
        "eval_mcc <- function(y_true, y_prob, plot=T)\n",
        "{\n",
        "    y_true_order <- y_true[order(y_prob)]\n",
        "\n",
        "    n <- length(y_true)\n",
        "    best_mcc <- 0\n",
        "    best_cutoff <- -1\n",
        "    mccs <- c()\n",
        "\n",
        "    for (i in 1:(n-1))\n",
        "    {\n",
        "        tp <- as.numeric(sum(y_true_order[(i+1):n]))\n",
        "        fp <- n - i - tp\n",
        "        fn <- as.numeric(sum(y_true_order[1:i]))\n",
        "        tn <- i - fn\n",
        "\n",
        "        new_mcc <- mcc(tp, fp, fn, tn)\n",
        "        mccs[i] <- new_mcc\n",
        "\n",
        "        if (new_mcc >= best_mcc)\n",
        "        {\n",
        "            best_mcc <- new_mcc\n",
        "            best_cutoff <- y_prob[order(y_prob)[i]]\n",
        "        }\n",
        "    }\n",
        "\n",
        "    if (plot) { plot(mccs) }\n",
        "\n",
        "    cat(best_mcc, best_cutoff, \"\\n\")\n",
        "    \n",
        "    return(best_mcc)\n",
        "}"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "b0b38ff6-022c-fbc5-7f64-dc2e026629c2"
      },
      "source": [
        "Trying out with a sample"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "5107ebfe-ab41-e8d8-922e-e2a0e9f55fc1"
      },
      "outputs": [],
      "source": [
        "set.seed(235)\n",
        "y_prob0 <- runif(10000, 0, 1)\n",
        "y_prob <- y_prob0 + 0.5 * runif(10000, 0, 1) - 0.02\n",
        "\n",
        "y_true <- ifelse(y_prob0 > 0.6, 1, 0)\n",
        "\n",
        "eval_mcc(y_true, y_prob)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "46b64c96-88d0-73f5-b26e-36a37169e4ac"
      },
      "source": [
        "Use the following function in XGBoost's feval"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "a96d71b5-1004-fd0d-9aa3-642d16835184"
      },
      "outputs": [],
      "source": [
        "mcc_eval <- function(pred, dtrain)\n",
        "{    \n",
        "    y_true <- getinfo(dtrain, \"label\")\n",
        "    best_mcc <- eval_mcc(y_true, pred, plot=F)\n",
        "    \n",
        "    return(list(metric=\"mcc\", value=best_mcc))\n",
        "}"
      ]
    }
  ],
  "metadata": {
    "_change_revision": 0,
    "_is_fork": false,
    "kernelspec": {
      "display_name": "R",
      "language": "R",
      "name": "ir"
    },
    "language_info": {
      "codemirror_mode": "r",
      "file_extension": ".r",
      "mimetype": "text/x-r-source",
      "name": "R",
      "pygments_lexer": "r",
      "version": "3.3.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}