doc/orca_tutorial_3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. [Threshold models](#threshold-models)\n",
    "\t1. [Proportional odds model (POM)](#proportional-odds-model-pom)\n",
    "\t2. [Neural network based on POM (NNPOM)](#neural-network-based-on-pom-nnpom)\n",
    "\t3. [Support vector for ordinal regression (SVOREX and SVORIM)](#support-vector-for-ordinal-regression-svorex-and-svorim)\n",
    "\t4. [Reduction from ordinal regression to binary SVM classifiers (REDSVM)](#reduction-from-ordinal-regression-to-binary-svm-classifiers-redsvm)\n",
    "\t5. [Kernel discriminant learning for ordinal regression (KDLOR)](#kernel-discriminant-learning-for-ordinal-regression-kdlor)\n",
    "\t6. [Ordinal regression boosting (ORBoost)](#ordinal-regression-boosting-orboost)\n",
    "\t7. [Custom Ensemble based on several projections](#custom-ensemble-based-on-several-projections)\n",
    "2. [References](#references)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Threshold models\n",
    "\n",
    "This tutorial will cover how to apply threshold models in the framework ORCA. It is highly recommended to have previously completed the [how to tutorial](orca_tutorial_1.md).\n",
    "\n",
    "Moreover, we are going to work again with melanoma diagnosis dataset. You should complete at least [the second section of the second tutorial](orca_tutorial_2.md#loading-the-dataset-and-performing-some-preliminary-experiments) in order to follow this third tutorial.\n",
    "\n",
    "All threshold models are designed with a very reasonable idea: the categories to be predicted in ordinal classification comes from the discretization of an underlying latent variable, so that we can try to model the latent variables and use a total of *Q-1* thresholds (for *Q* classes) to divide this variable in categories. In this way, the order of categories will be considered, because the intervals defined for each will be arranged in the same order, and a lot of flexibility will be given to the model by simply moving these thresholds.\n",
    "\n",
    "Because of this, there are many threshold model proposals in the literature, and ORCA includes some of the most popular ones:\n",
    "- One linear model (POM) [1].\n",
    "- One neural network model (NNPOM) [1,2].\n",
    "- Two support vector machine proposals (SVOREX and SVORIM) [3].\n",
    "- One reduction from ordinal regression to binary SVM (REDSVM) [4].\n",
    "- One discriminant analysis proposal (KDLOR) [5].\n",
    "- One ensemble model (ORBoost) [6].\n",
    "\n",
    "The corresponding script for this tutorial, ([exampleMelanomaTM.m](../src/code-examples/exampleMelanomaTM.m)), can be found and run in the [code example folder](../src/code-examples). Octave code is not shown here to simplify code, but it is included in ([exampleMelanomaTM.m](../src/code-examples/exampleMelanomaTM.m)).\n",
    "\n",
    "## Proportional odds model (POM)\n",
    "\n",
    "The POM  [1] arose from a statistical background, and it is based on an extension of binary logistic regression. Instead of using one single threshold to discriminate negative and positive class, the model includes as many thresholds as the number of classes minus one. The model uses one logistic function for each class, where the thresholds are arranged in ascending order and the linear part (projection vector) is common for all of them.\n",
    "\n",
    "A strong probabilistic basis is inherent to this approach: each logistic function is modelling the cumulative probability that a pattern has of belonging to the corresponding class or any of the previous classes. Cumulative probabilities can be easily transformed to standard probabilities by simply subtracting them.\n",
    "\n",
    "We can try the POM model using the melanoma dataset to check its performance. This method does not involve any parameter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "addpath('../src/Measures')\n",
    "addpath('../src/Algorithms')\n",
    "addpath('../src/Utils')\n",
    "\n",
    "% Disable warnings\n",
    "% MATLAB warnings\n",
    "%warning('off','MATLAB:nearlySingularMatrix')\n",
    "%warning('off','stats:mnrfit:IterOrEvalLimit')\n",
    "% Octave warnings\n",
    "warning('off','all');\n",
    "\n",
    "% Loading the data\n",
    "trainMelanoma = load('../exampledata/10-fold/melanoma-5classes-abcd-100/matlab/train_melanoma-5classes-abcd-100.5');\n",
    "testMelanoma = load('../exampledata/10-fold/melanoma-5classes-abcd-100/matlab/test_melanoma-5classes-abcd-100.5');\n",
    "% Form structures for training/test\n",
    "train.patterns = trainMelanoma(:,1:(end-1));\n",
    "train.targets = trainMelanoma(:,end);\n",
    "test.patterns = testMelanoma(:,1:(end-1));\n",
    "test.targets = testMelanoma(:,end);\n",
    "% Preprocess the data\n",
    "[train,test] = DataSet.deleteConstantAtributes(train,test);\n",
    "[train,test] = DataSet.standarizeData(train,test);\n",
    "[train,test] = DataSet.deleteNonNumericValues(train,test);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Apply POM model\n",
    "% Create the POM object\n",
    "algorithmObj = POM();\n",
    "\n",
    "% Train POM\n",
    "info = algorithmObj.fitpredict(train,test);\n",
    "\n",
    "% Evaluate the model\n",
    "fprintf('POM method\\n---------------\\n');\n",
    "fprintf('POM Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('POM MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An accuracy of ~67% is obtained in this fold using POM. Now, we can check the projection the model is doing and the set of thresholds."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Visualize the projection\n",
    "figure; hold on;\n",
    "\n",
    "hist(info.projectedTest', 30);\n",
    "y1=get(gca,'ylim');\n",
    "\n",
    "for i=1:size(info.model.thresholds,2)\n",
    "    line([+info.model.thresholds(i) +info.model.thresholds(i)],y1,'Color',[1 0 0]);\n",
    "end\n",
    "hold off;\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be checked no pattern is projected beyond the last threshold, so that the last class is ignored. Note that POM is a linear model and this can limit its accuracy. We can check this in the confusion matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "confusionmat(test.targets,info.predictedTest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Moreover, we can recode the projections, so that the patterns of each class appear in different colours:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Palette of colors. More at https://www.color-hex.com/\n",
    "palette=num2cell([255,98,0;255,192,0;160,255,0;0,234,255;66,0,116]/255,2);\n",
    "\n",
    "figure; hold on;\n",
    "Q = size(info.model.thresholds,2)+1;\n",
    "for i=1:Q\n",
    "    hist(info.projectedTest(test.targets==i), 30, 'FaceColor', palette{i});\n",
    "end\n",
    "y1=get(gca,'ylim');\n",
    "for i=1:size(info.model.thresholds,2)\n",
    "    line([info.model.thresholds(i) info.model.thresholds(i)],y1,'Color',[1 0 0]);\n",
    "end\n",
    "%legend('C1','C2','C3','C4','C5');\n",
    "legend(arrayfun(@(num) sprintf('C%d', num), 1:Q, 'UniformOutput', false))\n",
    "hold off;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be observed the three patterns from the last class are never correctly classified.\n",
    "\n",
    "We can also plot the probability estimations generated by POM, which are based on the `logit` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Visualize the cummulative probabilities\n",
    "figure; hold on;\n",
    "numPoints=300;\n",
    "x = linspace(min(info.model.thresholds)-3,max(info.model.thresholds)+3,numPoints);\n",
    "f = repmat(info.model.thresholds,numPoints,1) - repmat(x',1,Q-1);\n",
    "cumProb = [1./(1+exp(-f)) ones(numPoints,1)]; %logit function\n",
    "plot(x,cumProb,'-');\n",
    "y1=get(gca,'ylim');\n",
    "for i=1:size(info.model.thresholds,2)\n",
    "    line([info.model.thresholds(i) info.model.thresholds(i)],y1,'Color',[1 0 0]);    \n",
    "end\n",
    "hold off;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Visualize the individual probabilities\n",
    "figure; hold on;\n",
    "prob = cumProb;\n",
    "prob(:,2:end) = prob(:,2:end) - prob(:,1:(end-1));\n",
    "plot(x,prob,'-');\n",
    "y1=get(gca,'ylim');\n",
    "for i=1:size(info.model.thresholds,2)\n",
    "    line([info.model.thresholds(i) info.model.thresholds(i)],y1,'Color',[1 0 0]);    \n",
    "end\n",
    "hold off;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "As can be seen, those projections close to the thresholds can be classified in different classes according to the probability distribution. However, following the spirit of threshold models, the implementation of POM included in ORCA classify the patterns according to their position with respect to the thresholds.\n",
    "\n",
    "---\n",
    "\n",
    "***Exercise 1***: POM is a member of the Cumulative Link Models (CLM) family, where a `logit` function is considered. Now, you can try other probability distributions used in this models, such as the following ones (table from the [documentation](https://mran.microsoft.com/snapshot/2014-11-11/web/packages/ordinal/vignettes/clm_intro.pdf) of this [R package](https://cran.r-project.org/web/packages/ordinal/)):\n",
    "\n",
    "![Different link functions](tutorial/images/links.png)\n",
    "\n",
    "Compare graphically the different results.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural network based on POM (NNPOM)\n",
    "\n",
    "The idea of NNPOM is to extend POM [1,2] by considering a nonlinear projection, instead of a linear one. To do so, the projection function is a linear combination of nonlinear basis function (i.e. a neural network of one output node). As in NNOP, only one hidden layer is considered. Given that POM estimate a proper probability distribution, the cross entropy is used for gradient descent. The algorithm used for gradient descent is the iRProp+ algorithm.\n",
    "\n",
    "Three parameters must be specified in this case:\n",
    "- Parameter `hiddenN`, number of hidden neurons of the model.\n",
    "- Parameter `iter`, number of iterations for gradient descent.\n",
    "- Parameter `lambda`, regularization parameter in the error function (L2 regularizer), to avoid overfitting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Create the NNPOM object\n",
    "algorithmObj = NNPOM();\n",
    "% Train NNPOM\n",
    "info = algorithmObj.fitpredict(train,test,struct('hiddenN',20,'iter',500,'lambda',0.1));\n",
    "info = algorithmObj.fitpredict(train,test,struct('hiddenN',20,'iter',500,'lambda',0.1));\n",
    "% Evaluate the model\n",
    "fprintf('NNPOM method\\n---------------\\n');\n",
    "fprintf('NNPOM Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('NNPOM MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The source code of NNPOM clearly shows how the prediction is performed, which clarify a lot the model:\n",
    "\n",
    "```MATLAB\n",
    "function [projected, predicted]= privpredict(obj,test)\n",
    "    %PRIVPREDICT predicts labels of TEST patterns labels using fitted MODEL.\n",
    "    m = size(test,1);\n",
    "    a1 = [ones(m, 1) test];\n",
    "    z2 = a1*obj.model.Theta1';\n",
    "    a2 =  1.0 ./ (1.0 + exp(-z2));\n",
    "    projected=a2*obj.model.Theta2';\n",
    "\n",
    "    z3=repmat(obj.model.thresholds,m,1)-repmat(projected,1,obj.model.num_labels-1);\n",
    "    a3T =  1.0 ./ (1.0 + exp(-z3));\n",
    "    a3 = [a3T ones(m,1)];\n",
    "    a3(:,2:end) = a3(:,2:end) - a3(:,1:(end-1));\n",
    "    [M,predicted] = max(a3,[],2);            \n",
    "end\n",
    "```\n",
    "\n",
    "In this case, the results are not very good and even worse than POM. However, the configuration of the parameters is critical.\n",
    "\n",
    "You can also examine the threshold obtained and the projections:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info.model.thresholds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info.projectedTest(1:10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "***Exercise 2***: try to improve the results in this dataset by applying nested crossvalidation with an `ini` file.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Support vector for ordinal regression (SVOREX and SVORIM)\n",
    "\n",
    "Now we examine support vector approaches based on the threshold model structure [3]. The main methods are:\n",
    "- SVOREX, which is an ordinal formulation of the SVM paradigm, which computes discriminant parallel hyperplanes for the data and a set of thresholds by imposing explicit constraints in the optimization problem.\n",
    "- SVORIM, which is similar but computes discriminant parallel hyperplanes for the data and a set of thresholds by imposing implicit constraints in the optimization problem.\n",
    "\n",
    "Both are threshold models so that they obtain a projection together with a set of thresholds. In ORCA, they are implemented considering the RBF kernel. Consequently, the main parameters of both methods are:\n",
    "- Parameter `C`, importance given to errors.\n",
    "- Parameter `k`, inverse of the width of the RBF kernel.\n",
    "\n",
    "We are going to test both methods:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Apply the SVORIM model\n",
    "% Create the SVORIM object\n",
    "algorithmObj = SVORIM();\n",
    "\n",
    "% Train SVORIM\n",
    "info = algorithmObj.fitpredict(train,test,struct('C',10,'k',0.001));\n",
    "\n",
    "% Evaluate the model\n",
    "fprintf('SVORIM method\\n---------------\\n');\n",
    "fprintf('SVORIM Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('SVORIM MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest));\n",
    "\n",
    "% Store projections and thresholds (we will used them later)\n",
    "svorimProjections = info.projectedTest;\n",
    "svorimProjectionsTrain = info.projectedTrain;\n",
    "svorimThresholds = info.model.thresholds;\n",
    "\n",
    "%% Apply the SVOREX model\n",
    "% Create the SVOREX object\n",
    "algorithmObj = SVOREX();\n",
    "\n",
    "% Train SVOREX\n",
    "info = algorithmObj.fitpredict(train,test,struct('C',10,'k',0.001));\n",
    "\n",
    "% Evaluate the model\n",
    "fprintf('---------------\\n');\n",
    "fprintf('SVOREX method\\n---------------\\n');\n",
    "fprintf('SVOREX Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('SVOREX MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest));\n",
    "\n",
    "% Store projections and thresholds (we will used them later)\n",
    "svorexProjections = info.projectedTest;\n",
    "svorexProjectionsTrain = info.projectedTrain;\n",
    "svorexThresholds = info.model.thresholds;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot both projections by:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Represent both projections and thresholds\n",
    "figure; hold on; \n",
    "subplot(2,1,1)\n",
    "plot(svorimProjections,test.targets, 'o');\n",
    "y1=get(gca,'ylim');\n",
    "for i=1:size(svorimThresholds,2)\n",
    "    line([svorimThresholds(i) svorimThresholds(i)],y1,'Color',[1 0 0]);    \n",
    "end\n",
    "legend('SVORIM');\n",
    "subplot(2,1,2)\n",
    "plot(svorexProjections,test.targets, 'o');\n",
    "y1=get(gca,'ylim');\n",
    "for i=1:size(svorexThresholds,2)\n",
    "    line([svorexThresholds(i) svorexThresholds(i)],y1,'Color',[1 0 0]);    \n",
    "end\n",
    "legend('SVOREX');\n",
    "hold off;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fine tuning a bit the parameters, we can improve the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Apply the SVORIM model improved\n",
    "% Create the SVORIM object\n",
    "algorithmObj = SVORIM();\n",
    "\n",
    "% Train SVORIM\n",
    "info = algorithmObj.fitpredict(train,test,struct('C',200,'k',0.001));\n",
    "\n",
    "% Evaluate the model\n",
    "fprintf('SVORIM method improved\\n---------------\\n');\n",
    "fprintf('SVORIM Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('SVORIM MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Reduction from ordinal regression to binary SVM classifiers (REDSVM)\n",
    "\n",
    "The reduction from ordinal regression to binary SVM classifiers (REDSVM) [4] is a method that can be categorized both as threshold method or as decomposition method. The hyper-parameters are the well-known `k` and `C` of SVM variants."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Apply the REDSVM model\n",
    "% Create the REDSVM object\n",
    "algorithmObj = REDSVM();\n",
    "\n",
    "% Train REDSVM\n",
    "info = algorithmObj.fitpredict(train,test,struct('C',10,'k',0.001));\n",
    "\n",
    "% Evaluate the model\n",
    "fprintf('REDSVM method\\n---------------\\n');\n",
    "fprintf('REDSVM Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('REDSVM MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To better understand the relevance of parameters selection process, the following code optimizes parameters `k` and `C` using a 3Fold for each combination. Only in Matlab: the second cell plots corresponding validation results for `Acc` and `AMAE`. Note that the optimal combination may differ depending of the selected performance metric. Depending on your version of Matlab/Octave, a `contourf` or a `heatmap` is used for each metric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% REDSVM optimization\n",
    "clear T Ts;\n",
    "\n",
    "Metrics = {@MZE,@AMAE};\n",
    "setC = 10.^(-3:1:3);\n",
    "setk = 10.^(-3:1:3);\n",
    "% TODO: fix for Octave since table() is not supported\n",
    "Ts = cell(size(Metrics,2),1);\n",
    "nFolds = 3;\n",
    "\n",
    "if (exist ('OCTAVE_VERSION', 'builtin') > 0)\n",
    "    pkg load statistics;\n",
    "    CVO = cvpartition(train.targets,'KFold',nFolds);\n",
    "else\n",
    "    CVO = cvpartition(train.targets,'k',nFolds);\n",
    "end\n",
    "\n",
    "for m = 1:size(Metrics,2)\n",
    "    mObj = Metrics{m}();\n",
    "    fprintf('Grid search to optimize %s for REDSVM\\n', mObj.name);\n",
    "    bestError=Inf;\n",
    "    if (~exist ('OCTAVE_VERSION', 'builtin') > 0)\n",
    "      T = table();\n",
    "    end\n",
    "    for C=10.^(-3:1:3)\n",
    "        for k=10.^(-3:1:3)\n",
    "            error=0;\n",
    "            for ff = 1:nFolds\n",
    "                param = struct('C',C,'k',k);\n",
    "                info = algorithmObj.fitpredict(train,test,param);\n",
    "                error = error + mObj.calculateMetric(test.targets,info.predictedTest);\n",
    "\n",
    "            end\n",
    "            error = error / nFolds;\n",
    "            if error < bestError\n",
    "                bestError = error;\n",
    "                bestParam = param;\n",
    "            end\n",
    "            param.error = error;\n",
    "            if (~exist ('OCTAVE_VERSION', 'builtin') > 0)\n",
    "              T = [T; struct2table(param)];\n",
    "            end\n",
    "            fprintf('.');\n",
    "        end\n",
    "    end\n",
    "    if (~exist ('OCTAVE_VERSION', 'builtin') > 0)\n",
    "      Ts{m} = T;\n",
    "    end\n",
    "    fprintf('\\nBest Results REDSVM C %f, k %f --> %s: %f\\n', bestParam.C, bestParam.k, mObj.name, bestError);\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (exist ('OCTAVE_VERSION', 'builtin') > 0)\n",
    "  fprintf('This type of graphic is not supported in Octave\\n');\n",
    "else\n",
    "if verLessThan('matlab', '9.2')\n",
    "    % Use contours\n",
    "    figure;\n",
    "    hold on;\n",
    "    for m = 1:size(Metrics,2)\n",
    "        mObj = Metrics{m}();\n",
    "        subplot(size(Metrics,2),1,m)\n",
    "        x = Ts{m}{:,1};\n",
    "        y = Ts{m}{:,2};\n",
    "        z = Ts{m}{:,3};\n",
    "        numPoints=100;\n",
    "        [xi, yi] = meshgrid(linspace(min(x),max(x),numPoints),linspace(min(y),max(y),numPoints));\n",
    "        zi = griddata(x,y,z, xi,yi);\n",
    "        contourf(xi,yi,zi,15);\n",
    "        set(gca, 'XScale', 'log');\n",
    "        set(gca, 'YScale', 'log');\n",
    "        colorbar;\n",
    "        title([mObj.name ' optimization for REDSVM']);\n",
    "    end\n",
    "    hold off;\n",
    "else\n",
    "    % Use heatmaps\n",
    "    fprintf('Generating heat maps\\n');\n",
    "    figure;\n",
    "    subplot(2,1,1)\n",
    "    heatmap(Ts{1},'C','k','ColorVariable','error');\n",
    "    title('MZE optimization for REDSVM');\n",
    "\n",
    "    subplot(2,1,2)\n",
    "    heatmap(Ts{2},'C','k','ColorVariable','error');\n",
    "    title('AMAE optimization for REDSVM');\n",
    "end\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Kernel discriminant learning for ordinal regression (KDLOR)\n",
    "\n",
    "This method [5] adapts discriminant learning to the context of ordinal classification. The original discriminant learning problem is transformed by considering the minimum difference between the averages of two consecutive classes (in the ordinal scale). If this minimum difference is positive, the classes are correctly ranked according to the projection.\n",
    "\n",
    "After obtaining the projection, the thresholds are positioned in the mean point of the average projection of the classes they are separating. This makes the method quite sensitive to classes with low frequencies (being a good method for imbalanced ordinal regression problems).\n",
    "\n",
    "You can use KDLOR with three different kernels (`rbf`, `sigmoid` or `linear`) and with different optimization engines for the problem to be solved (`quadprog`, `qp` or `cvx`). Please, take into account that the `cvx` engine is not included by default in ORCA, so you will have to [install it](http://cvxr.com/cvx/). These two parameters are specified in the constructor (with the parameters `optimizationMethod` and `kernelType`, respectively), as they are not intended to be fine tuned by cross validation. During the optimization, KDLOR needs the following three parameters:\n",
    "- Parameter `C`, importance given to the minimum distance between average projections of classes.\n",
    "- Parameter `k`, inverse of the width of the RBF kernel.\n",
    "- Parameter `u`, constant summed to the main diagonal of the kernel matrix, aimed at avoiding singularities.\n",
    "\n",
    "Let try KDLOR in the melanoma dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Apply the KDLOR model\n",
    "% Create the KDLOR object\n",
    "algorithmObj = KDLOR('kernelType','rbf');\n",
    "\n",
    "% Train KDLOR\n",
    "info = algorithmObj.fitpredict(train,test,struct('C',1,'k',0.001,'u',0.01));\n",
    "\n",
    "% Evaluate the model\n",
    "fprintf('KDLOR method\\n---------------\\n');\n",
    "fprintf('KDLOR Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('KDLOR MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "MS.calculateMetric(test.targets,info.predictedTest)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "confusionmat(test.targets,info.predictedTest)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "as previously discussed, because of the way the thresholds are derived, this method is quite robust for minority classes.\n",
    "\n",
    "Let check it in the projections:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "% Visualize the projection with colors\n",
    "figure; hold on;\n",
    "Q = size(info.model.thresholds,2)+1;\n",
    "for i=1:Q\n",
    "    hist(info.projectedTest(test.targets==i),30, 'FaceColor', palette{i});\n",
    "end\n",
    "y1=get(gca,'ylim');\n",
    "for i=1:size(info.model.thresholds,2)\n",
    "    line([info.model.thresholds(i) info.model.thresholds(i)],y1,'Color',[1 0 0]);    \n",
    "end\n",
    "%legend('C1','C2','C3','C4','C5');\n",
    "legend(arrayfun(@(num) sprintf('C%d', num), 1:Q, 'UniformOutput', false))\n",
    "hold off;"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "\n",
    "***Exercise 3***: Compare the results obtained in KDLOR by using different kernel functions.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Ordinal regression boosting (ORBoost)\n",
    "\n",
    "Ordinal regression boosting (ORBoost) is a thresholded-ensemble model, which is composed of confidence functions, and their weighted linear combination is used as the one-dimensional mapping. A set of thresholds for this mapping is also included in the model and iteratively updated with the rest of parameters.\n",
    "\n",
    "As proposed by the authors, the total number of ensemble members is set to `T=2000`, and normalised sigmoid functions are used as the base classifier, where the smoothness parameter is `gamma=4`. Large margin bounds of the classification error and the absolute error are derived, from which two algorithms are presented: ORBoost with all margins and ORBoost with left-right margins. The `weights` parameter in the constructor configures whether the All margins versions is used (`weights=true`) or the Left-Right margin is used (`weights=false`).\n",
    "\n",
    "This is the code for running ORBoost with the melanoma diagnosis dataset (this will take some time):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Apply the ORBoost model\n",
    "% Create the ORBoost object\n",
    "algorithmObj = ORBoost('weights',true);\n",
    "\n",
    "% Train ORBoost\n",
    "info = algorithmObj.fitpredict(train,test);\n",
    "\n",
    "% Evaluate the model\n",
    "fprintf('ORBoost method\\n---------------\\n');\n",
    "fprintf('ORBoost Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('ORBoost MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Custom Ensemble based on several projections\n",
    "\n",
    "We can construct an ensemble by using the projection given by different models. In the following example, we combine the projections given by SVORIM and SVOREX method to construct a new dataset. With this dataset, we apply the POM algorithm:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%% Make an ensemble with SVORIM and SVOREX. The final decission by POM\n",
    "% Construct a new dataset with the projections\n",
    "newTrain.patterns = [svorexProjectionsTrain svorimProjectionsTrain];\n",
    "newTrain.targets = train.targets;\n",
    "newTest.patterns = [svorexProjections svorimProjections];\n",
    "newTest.targets = train.targets;\n",
    "% Preprocess the dataset\n",
    "[newTrain,newTest] = DataSet.deleteConstantAtributes(newTrain,newTest);\n",
    "[newTrain,newTest] = DataSet.standarizeData(newTrain,newTest);\n",
    "[newTrain,newTest] = DataSet.deleteNonNumericValues(newTrain,newTest);\n",
    "% Train the final POM model\n",
    "algorithmObj = POM();\n",
    "info = algorithmObj.fitpredict(newTrain,newTest);\n",
    "% Evaluate the ensemble\n",
    "fprintf('SVORIM+SVOREX+POM method\\n---------------\\n');\n",
    "fprintf('SVORIM+SVOREX+POM Accuracy: %f\\n', CCR.calculateMetric(test.targets,info.predictedTest));\n",
    "fprintf('SVORIM+SVOREX+POM MAE: %f\\n', MAE.calculateMetric(test.targets,info.predictedTest))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we check the dataset used for POM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "scatter(newTrain.patterns(:,1),newTrain.patterns(:,2),150,newTrain.targets, \"filled\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that, although the correlation of both projections is quite high, some patterns can be refined by considering both projections.\n",
    "\n",
    "---\n",
    "\n",
    "***Exercise 4***: construct a similar ensemble but using different SVORIM projections with different subsets of input variables (a 40% of randomly chosen variables). The number of members of the ensemble should be as a parameter (try 50).\n",
    "\n",
    "----\n",
    "\n",
    "***Exercise 5***: construct a similar ensemble but using different SVORIM projections with different parameters for the `C` value.\n",
    "\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "1. P. McCullagh, \"Regression models for ordinal data\",  Journal of the Royal Statistical Society. Series B (Methodological), vol. 42, no. 2, pp. 109–142, 1980.\n",
    "1. M. J. Mathieson, \"Ordinal models for neural networks\", in Proc. 3rd Int. Conf. Neural Netw. Capital Markets, 1996, pp. 523-536.\n",
    "1. W. Chu and S. S. Keerthi, \"Support Vector Ordinal Regression\", Neural Computation, vol. 19, no. 3, pp. 792–815, 2007. http://10.1162/neco.2007.19.3.792\n",
    "1. H.-T. Lin and L. Li, \"Reduction from cost-sensitive ordinal ranking to weighted binary classification\" Neural Computation, vol. 24, no. 5, pp. 1329-1367, 2012. http://10.1162/NECO_a_00265\n",
    "1. B.-Y. Sun, J. Li, D. D. Wu, X.-M. Zhang, and W.-B. Li, \"Kernel discriminant learning for ordinal regression\", IEEE Transactions on Knowledge and Data Engineering, vol. 22, no. 6, pp. 906-910, 2010. https://doi.org/10.1109/TKDE.2009.170\n",
    "1. H.-T. Lin and L. Li, \"Large-margin thresholded ensembles for ordinal regression: Theory and practice\", in Proc. of the 17th Algorithmic Learning Theory International Conference, ser. Lecture Notes in Artificial Intelligence (LNAI), J. L. Balcazar, P. M. Long, and F. Stephan, Eds., vol. 4264. Springer-Verlag, October 2006, pp. 319-333."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Octave",
   "language": "octave",
   "name": "octave"
  },
  "language_info": {
   "file_extension": ".m",
   "help_links": [
    {
     "text": "GNU Octave",
     "url": "https://www.gnu.org/software/octave/support.html"
    },
    {
     "text": "Octave Kernel",
     "url": "https://github.com/Calysto/octave_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-octave",
   "name": "octave",
   "version": "4.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}