{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "5b100c3b-c535-62d6-3404-0f815649c061"
      },
      "source": [
        "## Extreme benchmarking of MCC implementations in R.\n",
        "\n",
        "by Laurae\n",
        "\n",
        "## Quick results\n",
        "\n",
        "* Method 1: disconnects the kernel / too slow\n",
        "* Method 2: disconnects the kernel / too slow\n",
        "* Method 3: does not coerce integer to numeric, hence unscalable for large ints\n",
        "* Method 4: scales perfectly\n",
        "\n",
        "The original source (which this one has been forked from) shows all methods provide the same result.\n",
        "\n",
        "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_mcc1 <- 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": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "de013cf3-8c98-9ee4-3df3-256e2a52ec63"
      },
      "outputs": [],
      "source": [
        "eval_mcc2 <- 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",
        "    mccs[i] <- mcc(tp, fp, fn, tn)\n",
        "  }\n",
        "  \n",
        "  if (plot) { plot(mccs) }\n",
        "  \n",
        "  best = which.max(mccs)\n",
        "  best_mcc = mccs[best]\n",
        "  cat(best_mcc, y_prob[order(y_prob)[best]], \"\\n\")\n",
        "  \n",
        "  return(best_mcc)\n",
        "}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "ef8a6039-9353-96c8-db30-90d18d54e863"
      },
      "outputs": [],
      "source": [
        "require(ROCR)\n",
        "eval_mcc3 <- function(actual, predicted_prob, show=FALSE){\n",
        "  \n",
        "  pred2<- ROCR::prediction(predicted_prob, actual)\n",
        "  perf<- ROCR:: performance(pred2, \"mat\")\n",
        "  optim = which.max(perf@y.values[[1]])[1]\n",
        "  best_split <- perf@x.values[[1]][optim]\n",
        "  best_mcc = perf@y.values[[1]][optim]\n",
        "\n",
        "  if (show){\n",
        "    ROCR::plot(perf)\n",
        "  }\n",
        "  return(c(best_mcc))\n",
        "}"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "77444118-afeb-a1d0-03ba-0ae78b124634"
      },
      "outputs": [],
      "source": [
        "eval_mcc4 <- function(y_true, y_prob, show=F) {\n",
        "  idx <- order(y_prob)\n",
        "  y_prob_sort <- y_prob[idx]\n",
        "  y_true_sort <- y_true[idx]\n",
        "  n <- length(y_true)\n",
        "  nump <- sum(y_true)\n",
        "  numn <- n- nump\n",
        "\n",
        "  tn_v <- cumsum(y_true_sort == 0)\n",
        "  fp_v <- cumsum(y_true_sort == 1)\n",
        "  fn_v <- numn - tn_v\n",
        "  tp_v <- nump - fp_v\n",
        "  sup_v <- tp_v * tn_v - fp_v * fn_v\n",
        "  inf_v <- sqrt((tp_v + fp_v) * (tp_v + fn_v) * (tn_v + fp_v) * (tn_v + fn_v))\n",
        "  mcc_v <- sup_v/inf_v\n",
        "  mcc_v[!is.finite(mcc_v)] <- 0\n",
        "  best_id <- which.max(mcc_v)\n",
        "  best_mcc <- mcc_v[best_id]\n",
        "  best_proba <- y_prob_sort[best_id]\n",
        "\n",
        "  if (show) {\n",
        "    y_pred <- as.numeric(y_prob > best_proba)\n",
        "    score <- mcc(y_true, y_pred)\n",
        "    cat(\"\\n\",score, best_mcc)\n",
        "    return(list(best_proba, best_mcc, y_pred))\n",
        "  } else {\n",
        "    return(best_mcc)\n",
        "  }\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(11111)\n",
        "y_prob0 <- runif(1183748, 0, 1)\n",
        "y_prob <- y_prob0 + 0.5 * runif(1183748, 0, 1) - 0.02\n",
        "\n",
        "y_true <- ifelse(y_prob0 > 0.6, 1, 0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "_cell_guid": "284278b8-68ef-7140-0c09-782e4be3643c"
      },
      "source": [
        "Time them"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "9395e5bc-7f9e-b7a7-5452-748764be36a1"
      },
      "outputs": [],
      "source": [
        "# Disconnects the kernel, hence skipped\n",
        "#system.time(eval_mcc1(y_true, y_prob, plot = F))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "1dfc75ac-0e98-c856-8a35-dbc2b942431a"
      },
      "outputs": [],
      "source": [
        "# Disconnects the kernel, hence skipped\n",
        "#system.time(eval_mcc2(y_true, y_prob, plot = F)) # darragh"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "4b659943-bc11-50cc-5a46-87c75fca2fe7"
      },
      "outputs": [],
      "source": [
        "system.time(eval_mcc3(y_true, y_prob, show = F)) # mrooijer"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "bb867515-3ebd-a424-8da1-3587e8ce680a"
      },
      "outputs": [],
      "source": [
        "system.time(eval_mcc4(y_true, y_prob, show = F)) # daniel fg"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "65867a2f-692d-2579-ddeb-ae7bc2140d4a"
      },
      "outputs": [],
      "source": [
        "#all.equal(eval_mcc3(y_true, y_prob, plot = F), eval_mcc1(y_true, y_prob, plot = F))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "4dda754b-e92e-3d37-5086-b115f47a1ff2"
      },
      "outputs": [],
      "source": [
        "all.equal(eval_mcc3(y_true, y_prob, plot = F), eval_mcc2(y_true, y_prob, show = F))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "_cell_guid": "3063ed95-90cb-419e-adda-c09bf63b4aac"
      },
      "outputs": [],
      "source": [
        "all.equal(eval_mcc3(y_true, y_prob, plot = F), eval_mcc4(y_true, y_prob, show = F))"
      ]
    },
    {
      "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
}