Chapter_02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification\n",
    "\n",
    "linear Regression models assume a quantitative $Y$\n",
    "\n",
    "But: often $Y$ is qualitative/ categorical e.g. eye color \"blue, brown, green\"\n",
    "\n",
    "$\\Rightarrow$ predicting qualitative response ($Y$) is called classification\n",
    "\n",
    "classifier widely used:\n",
    " - logistic regression\n",
    " - linear discriminant analysis (LDA)\n",
    " - K-nearest neighbors (KNN)\n",
    " - quadratic discriminant analysis (QDA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Why not linear Regression?\n",
    "\n",
    "again linear regression models assume a quantitative response variable $Y$.\n",
    "\n",
    "But: qualitative response variables do not assume an ordering or the distance between classes in not equal\n",
    "\n",
    "Hence, quantitative $\\neq$ qualitative!\n",
    "\n",
    "Note: For a binary (two level) qualitative response we can implement a Dummy Variable approach!\n",
    "\n",
    "$$Y=\n",
    "\\begin{cases}\n",
    "0 \\quad False \\quad (A)\\\\\n",
    "1 \\quad True \\quad (B)\n",
    "\\end{cases}$$\n",
    "\n",
    "predict $\\hat{y}>0.5$ to class/ group B\n",
    "\n",
    "$$Pr(group=B|X)=X\\hat{\\beta}$$\n",
    "\n",
    "But: no guarantee that $X\\hat{\\beta}\\in [0,1]$\n",
    "\n",
    "When we would expect $X\\hat{\\beta}\\notin [0,1]$?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![alt text](cdf.svg \"cdf function\")\n",
    "\n",
    "If we would expect very high or very low probabilities, then a linear model would not be a good approximation.\n",
    "\n",
    "In contrast in the middle of the CDF function (Logistic Regression) the curve is very steep, thus the linear regression is able to approximate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Logistic Regression\n",
    "\n",
    "$Pr(group = B | X) = X\\hat{\\beta}$\n",
    "\n",
    "linear model: \n",
    " $$p(X)=\\beta_0 + \\beta_1 X$$\n",
    " \n",
    "logistic function:\n",
    " $$p(X)=\\frac{e^{\\beta_0 + \\beta_1 X}}{1+e^{\\beta_0 + \\beta_1 X}}$$\n",
    " \n",
    "Hence, the odds are:\n",
    " $$ \\frac{p(X)}{1-p(X)} = e^{\\beta_0 + \\beta_1 X} \\in \\{0,\\infty \\} \\\\\n",
    "     log \\left(\\frac{p(X)}{1-p(X)}\\right) =  \\beta_0 + \\beta_1 X $$\n",
    "is an S-shaped function\n",
    " - estimation using ML\n",
    " - logistic regression for $>2$ classes: \n",
    " $$ 1- Pr(group = A | X)- Pr(group = B | X)$$\n",
    " \n",
    "**Logistic Regression: conditional distribution of Y on X**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Discriminant Analysis (LDA)\n",
    "\n",
    " * model the distribution of the predictors $X$ spearately in each of the classes\n",
    " * use Bayes theorem to flip these around into estimes for $Pr(Y=k|X=x)$\n",
    " \n",
    "### Why not using logistic Regression?\n",
    " - if classes well-separated, logistic Regression estimates are unstable\n",
    " - if n is small and the distribution of the predictors X is approx. normal in each class, LDA is again more stable\n",
    " - LDA is popular if classes $>2$\n",
    " \n",
    "### Bayes Theorem for Classification Problems\n",
    "\n",
    " - $\\pi_k\\equiv Pr(Y=k)$  prior probability that random obs. belongs to class k\n",
    " - $f_k(X)\\equiv Pr(X=x|Y=k)$ is density \n",
    " \n",
    " $$\n",
    " p_k(X)=Pr(Y=k|X=x) =\\frac{Pr(X=x\\mid Y=k)Pr(Y=k)}{Pr(X=x)}\\\\\n",
    " = \\frac{\\pi_k f_k(X)}{\\sum_{l=1}^K \\pi_l f_l(X)}\n",
    " $$\n",
    " \n",
    "Thus, $p_k(X)$ is the posterior probability.\n",
    "\n",
    "Hence, for estimation we need to compute $\\pi_k$ using the sample and make assumptions about the functional form of $f_k$\n",
    "\n",
    "Deriving the LDA formula is done by setting up the Bayes theorem with the assumed functionalform, then taking the logarithm and rearrange these term!\n",
    "\n",
    "### Why linear?\n",
    "$$ \\delta_k(x)=x*\\frac{\\mu_k}{\\sigma^2}-\\frac{\\mu_k^2}{2\\sigma^2} + log(\\pi_k) $$\n",
    "\n",
    "LDA assignes $x$ to class $k$, where $\\delta_k(x)$ is maximized.\n",
    "\n",
    "Therefore the **decision boundaries are linear**!\n",
    "\n",
    "### Sensitivity, Specificity and total error\n",
    "\n",
    "* Sensitivity = e.g. percentage of true defaulters correctly predicted\n",
    "* Specificity = e.g. perentage of non-defaulters correctly identified \n",
    "\n",
    "**True-positive rate = sensitivity**\n",
    "\n",
    "**False-positive rate = 1-specificity**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quadratic Discriminant Analysis (QDA)\n",
    "\n",
    "in principle the same as LDA, but assumes that **each class** has it's **own covariate matrix**\n",
    "\n",
    "$$ X \\sim N(\\mu_k, \\Sigma_k) $$\n",
    "\n",
    "### Bias-variance Trade-off\n",
    "\n",
    "* $p$ predictors estimating $\\Sigma$ entails estimating \n",
    "  $$ \\frac{p(p+1)}{2} $$\n",
    "  \n",
    "* estimating $\\Sigma_k$ leads to estimating\n",
    "  $$ K*\\frac{p(p+1)}{2} $$\n",
    "  \n",
    "$\\Rightarrow$ LDA is much less flexible than QDA (assuming $\\Sigma$ instead of $\\Sigma_k$, thus LDA is linear in $X$)\n",
    "\n",
    "Note: $\\Sigma \\neq \\Sigma_k \\quad, \\forall k$ leads to high bias for LDA! \n",
    "\n",
    "### Why quadratic?\n",
    "$$ \\delta_k(x)= -\\frac{1}{2}(x-\\mu_k)'\\Sigma_k^{-1}(x-\\mu_k) -\\frac{1}{2}log \\mid \\Sigma_k \\mid + log(\\pi_k) $$\n",
    "\n",
    "QDA assignes $x$ to class $k$, where $\\delta_k(x)$ is maximized.\n",
    "\n",
    "Since $(x-\\mu_k)'\\Sigma_k^{-1}(x-\\mu_k)$ the **decision boundaries are quadratic**!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of Classifiers\n",
    "\n",
    "### Logistic Regression vs LDA\n",
    "\n",
    " * both linear decision boundaries\n",
    " * different estimation and assumptions\n",
    " \n",
    "if observations are Gausian with common covariance matrix $\\Sigma$:\n",
    "**LDA > Logistic Regression**\n",
    "\n",
    "### K-nearest Neighbors\n",
    "\n",
    "$$ Pr(Y=j|X=x_0)=\\frac{1}{K} \\sum_{i\\in N_0} I(y_i=j) $$\n",
    "\n",
    "is a non-parametric approach: no assumptions about the shape of decision boundaries\n",
    "\n",
    "is useful if decision boundary is nonlinear\n",
    "\n",
    "**KNN > LDA & Logistic Regression**\n",
    "\n",
    "But: no coefficient estimates!\n",
    "\n",
    "### QDA as compromise\n",
    "\n",
    "is a compromise between the non-parametric KNN & LDA/ Logistic Regression \n",
    " - not as flexbile as KNN\n",
    " - performs better with limited sample size (number of observations)\n",
    " \n",
    "### Comparison of Classifiers\n",
    " * Logistic regression - linear boundary, ML-estimation (no distributional assumption)\n",
    " \n",
    " * LDA - linear boundary, Gaussian distributional assumption\n",
    " \n",
    " * QDA - non-linear boundary, Gaussian distributional assumption\n",
    " \n",
    " * KNN - flexibility depending on $k$ (if $k=1$ most flexible), select $k$ by CV, no assumptions are made"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Classification Error assessment\n",
    "\n",
    "Suppose we have three different classes and an classification error rate of $5\\%$.\n",
    "\n",
    "Consider two examples, where this is a good/ bad result.\n",
    "\n",
    " * a bad result given these prior probabilities\n",
    " \n",
    "   $P(K_1)=0.99$\n",
    " \n",
    "   $P(K_2)=0.05$\n",
    "   \n",
    "   $P(K_3)=0.05$\n",
    "  \n",
    " * a good result given these prior probabilities\n",
    "   \n",
    "   $P(K_1)=0.33$\n",
    "   \n",
    "   $P(K_2)=0.33$\n",
    "   \n",
    "   $P(K_3)=0.33$\n",
    "   \n",
    "The second example is good, because the Classifier (only $5\\%$ error) performs better than the Null-Classifier!\n",
    "\n",
    "Hence, we have gained no more information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Derivation of LDA\n",
    "\n",
    "let's start with the Bayes Theorem:\n",
    "$$\n",
    "p_k(X)=Pr(Y=k|X=x)\\\\\n",
    "=\\frac{Pr(X=x\\mid Y=k)Pr(Y=k)}{Pr(X=x)}\n",
    "$$\n",
    "\n",
    "we now can define and insert some terms:\n",
    "$$\n",
    "P(Y=k)=\\pi_k\\\\\n",
    "P(X=x|Y=k)=f_k=\\frac{1}{\\sqrt{2\\pi\\sigma_k^2}} exp\\{-\\frac{1}{2\\sigma^2}(x-\\mu_k)^2\\}\\\\\n",
    "P(X=x)=C \\text{ does not depend on $k$}\n",
    "$$\n",
    "\n",
    "let's take the logarithm of both sides\n",
    "$$\n",
    "log(P(Y=k|X=x))=log(\\pi_k)+log(\\frac{1}{\\sqrt{2\\pi\\sigma_k^2}})-\\frac{1}{2\\sigma_k^2}(x-\\mu_k)^2 -log(C)\n",
    "$$\n",
    "\n",
    "since we are interested in maximizing the posterior probability of $x$ belonging to class $k$ and by assumption $\\sigma_k^2=\\sigma^2$\n",
    "$$\n",
    "log(P(Y=k|X=x))=\\delta_k(x)=log(\\pi_k)-\\frac{1}{2\\sigma_k^2}(x-\\mu_k)^2 \\\\\n",
    "\\delta_k(x)=x*\\frac{\\mu_k}{\\sigma^2}-\\frac{\\mu_k^2}{2\\sigma^2}+log(\\pi_k)\n",
    "$$\n",
    "\n",
    "Thus, the decision boundary is $\\delta_k=\\delta_l$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises for LDA and QDA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the following data generating process in which n observations belong to one of two classes. There are two covariates, drawn from normal distribution $x_1\\sim N$ and $x_2 \\sim N$ with class specific means. The class means are $\\mu_1 = (−3 \\quad 3)$ for class 1, and $\\mu_2 = (5 \\quad 5)$ for class 2 and $\\Sigma_1 = \\Sigma_2$. Initially, you may set $\\Sigma=(16 \\quad -2)(-2 \\quad 9)$ and $n_1 =300$ and $n_2 =500$.\n",
    "The goal of this exercise is to compare the performance of linear discriminant analysis and quadratic discriminant analysis when classifying observations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(MASS)       # to fit LDA and QDA analysis\n",
    "library(mvtnorm)\n",
    "# Note: packages above require the data to be saved as a data frame\n",
    "# class 1\n",
    "n1 = 300\n",
    "mu1= c(-3,3)\n",
    "# class 2\n",
    "n2 = 500\n",
    "mu2= c(5,5)\n",
    "# Variance-Covariance Matrix\n",
    "covmat=matrix(c(16,-2,-2,9), nrow=2, ncol=2)\n",
    "#total number of observations\n",
    "N=n1+n2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "a) Generate the covariates from a multivariate normal distribution using the $\\mu_k$ and $\\Sigma$ as described above and an indicator variable indicating class dependence for n observations and combine these in a data frame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "## a)\n",
    "# DGP - for both classes\n",
    "set.seed(123)\n",
    "X1=mvrnorm(n = n1, mu = mu1, Sigma = covmat)\n",
    "X2=mvrnorm(n = n2, mu = mu2, Sigma = covmat)\n",
    "df1 <- data.frame(X1)\n",
    "df1['class'] = 1\n",
    "df2 <- data.frame(X2)\n",
    "df2['class'] = 2\n",
    "df  <- merge(df1,df2, all=TRUE)\n",
    "rm(df1, df2, X1, X2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) Calculate the linear discriminant analysis and quadratic discriminant analysis, estimating all relevant quantities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "        Length Class  Mode     \n",
       "prior   2      -none- numeric  \n",
       "counts  2      -none- numeric  \n",
       "means   4      -none- numeric  \n",
       "scaling 2      -none- numeric  \n",
       "lev     2      -none- character\n",
       "svd     1      -none- numeric  \n",
       "N       1      -none- numeric  \n",
       "call    3      -none- call     \n",
       "terms   3      terms  call     \n",
       "xlevels 0      -none- list     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "        Length Class  Mode     \n",
       "prior   2      -none- numeric  \n",
       "counts  2      -none- numeric  \n",
       "means   4      -none- numeric  \n",
       "scaling 8      -none- numeric  \n",
       "ldet    2      -none- numeric  \n",
       "lev     2      -none- character\n",
       "N       1      -none- numeric  \n",
       "call    3      -none- call     \n",
       "terms   3      terms  call     \n",
       "xlevels 0      -none- list     "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## b)\n",
    "mod_lda  =  lda(class ~ X1 + X2, data=df)\n",
    "summary(mod_lda)\n",
    "class_lfit  <- as.numeric(predict(mod_lda)$class)\n",
    "mod_qda  =  qda(class ~ X1 + X2 , data=df)\n",
    "summary(mod_qda)\n",
    "class_qfit  <- as.numeric(predict(mod_qda)$class)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "c) Calculate the mean training error for both methods and compare."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 0.11875\n",
      "[1] 0.12\n",
      "[1] -0.00125\n"
     ]
    }
   ],
   "source": [
    "## c)\n",
    "# Note: since classes are ordinal scale we can not use MSE, due to the fact\n",
    "#       that the distance between class 1 and 3 are the same as between 1 and 2\n",
    "#       furthermore it is not appropriarte to use OLS\n",
    "MTE_LDA=sum(class_lfit!=df$class)/N\n",
    "MTE_QDA=sum(class_qfit!=df$class)/N\n",
    "print(MTE_LDA)\n",
    "print(MTE_QDA)\n",
    "print(MTE_LDA-MTE_QDA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Simulation Study**\n",
    "\n",
    "a) Evaluate the difference between the two methods through calculating classification training error in a simulation study for 100 different samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\f"
     ]
    },
    {
     "data": {
      "text/plain": [
       "   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. \n",
       " -1e-04  -1e-04  -1e-04  -1e-04  -1e-04  -1e-04 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# a)\n",
    "rm(list=ls())\n",
    "cat(\"\\014\")\n",
    "\n",
    "MCN=100\n",
    "MSE=matrix(NaN,MCN,2)\n",
    "\n",
    "\n",
    "n1 = 300\n",
    "mu1= c(-3,3)\n",
    "n2 = 500\n",
    "mu2= c(5,5)\n",
    "covmat=matrix(c(16,-2,-2,9), nrow=2, ncol=2)\n",
    "N=n1+n2\n",
    "\n",
    "set.seed(123)\n",
    "\n",
    "for (i in 1:MCN){\n",
    "  X1=mvrnorm(n = n1, mu = mu1, Sigma = covmat)\n",
    "  X2=mvrnorm(n = n2, mu = mu2, Sigma = covmat)\n",
    "  df1 <- data.frame(X1)\n",
    "  df1['class'] = 1\n",
    "  df2 <- data.frame(X2)\n",
    "  df2['class'] = 2\n",
    "  df  <- merge(df1,df2, all=TRUE)\n",
    "\n",
    "  mod_lda  =  lda(class ~ X1 + X2, data=df)\n",
    "  class_lfit  <- as.numeric(predict(mod_lda)$class)\n",
    "  mod_qda  =  qda(class ~ X1 + X2 , data=df)\n",
    "  class_qfit  <- as.numeric(predict(mod_qda)$class)\n",
    "\n",
    "  MSE[i,1]=sum(class_lfit!=df$class)/N\n",
    "  MSE[i,2]=sum(class_qfit!=df$class)/N\n",
    "}\n",
    "\n",
    "avg_MSE_LDA=mean(MSE[,1])\n",
    "avg_MSE_QDA=mean(MSE[,2])\n",
    "\n",
    "#par(mfrow=c(1,2))\n",
    "#plot(MSE[,1], ylab=\"LDA\")\n",
    "#abline(h=avg_MSE_LDA, col=\"red\")\n",
    "#plot(MSE[,2], ylab=\"QDA\")\n",
    "#abline(h=avg_MSE_QDA, col=\"red\")\n",
    "\n",
    "summary(avg_MSE_LDA-avg_MSE_QDA)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "b) Consider the theoretical properties of lda and qda that we discussed in the lecture: Which properties of the initial simulation set up could we manipulate in order to increase the difference between the classification error of lda and qda? Test your intuition by performing a suitable simulation study."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\f"
     ]
    },
    {
     "data": {
      "text/plain": [
       "    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. \n",
       "0.002275 0.002275 0.002275 0.002275 0.002275 0.002275 "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# b)\n",
    "# Note: since we have different covariate matrixes sigma 1 and sigma 2\n",
    "#       QDA is more precise than LDA\n",
    "#       From a theoretical perspective:\n",
    "#       * if LDAs assumption that the K classes share a common covariance matrix\n",
    "#          is badly off, then LDA --> high bias\n",
    "#       * LDA is a much less flexible classifier than QDA --> lower variance\n",
    "# Hence: Try with different covariance matrices\n",
    "\n",
    "rm(list=ls())\n",
    "cat(\"\\014\")\n",
    "\n",
    "MCN=100\n",
    "MSE=matrix(NaN,MCN,2)\n",
    "\n",
    "\n",
    "n1 = 300\n",
    "mu1= c(-3,3)\n",
    "n2 = 500\n",
    "mu2= c(5,5)\n",
    "covmat_1=matrix(c(16,-2,-2,9), nrow=2, ncol=2)\n",
    "covmat_2=matrix(c(10,-2,-2,5), nrow=2, ncol=2)\n",
    "N=n1+n2\n",
    "\n",
    "set.seed(123)\n",
    "\n",
    "for (i in 1:MCN){\n",
    "  X1=mvrnorm(n = n1, mu = mu1, Sigma = covmat_1)\n",
    "  X2=mvrnorm(n = n2, mu = mu2, Sigma = covmat_2)\n",
    "  df1 <- data.frame(X1)\n",
    "  df1['class'] = 1\n",
    "  df2 <- data.frame(X2)\n",
    "  df2['class'] = 2\n",
    "  df  <- merge(df1,df2, all=TRUE)\n",
    "\n",
    "  mod_lda  =  lda(class ~ X1 + X2, data=df)\n",
    "  class_lfit  <- as.numeric(predict(mod_lda)$class)\n",
    "  mod_qda  =  qda(class ~ X1 + X2 , data=df)\n",
    "  class_qfit  <- as.numeric(predict(mod_qda)$class)\n",
    "\n",
    "  MSE[i,1]=sum(class_lfit!=df$class)/N\n",
    "  MSE[i,2]=sum(class_qfit!=df$class)/N\n",
    "}\n",
    "\n",
    "avg_MSE_LDA=mean(MSE[,1])\n",
    "avg_MSE_QDA=mean(MSE[,2])\n",
    "\n",
    "#par(mfrow=c(1,2))\n",
    "#plot(MSE[,1], ylab=\"LDA\")\n",
    "#abline(h=avg_MSE_LDA, col=\"red\")\n",
    "#plot(MSE[,2], ylab=\"QDA\")\n",
    "#abline(h=avg_MSE_QDA, col=\"red\")\n",
    "\n",
    "summary(avg_MSE_LDA-avg_MSE_QDA)"
   ]
  }
 ],
 "metadata": {
  "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.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}