\documentclass[ignorenonframetext,]{beamer}
\setbeamertemplate{caption}[numbered]
\setbeamertemplate{caption label separator}{: }
\setbeamercolor{caption name}{fg=normal text.fg}
\beamertemplatenavigationsymbolsempty
\usepackage{lmodern}
\usepackage{amssymb,amsmath}
\usepackage{ifxetex,ifluatex}
\usepackage{fixltx2e} % provides \textsubscript
\ifnum 0\ifxetex 1\fi\ifluatex 1\fi=0 % if pdftex
  \usepackage[T1]{fontenc}
  \usepackage[utf8]{inputenc}
\else % if luatex or xelatex
  \ifxetex
    \usepackage{mathspec}
  \else
    \usepackage{fontspec}
  \fi
  \defaultfontfeatures{Ligatures=TeX,Scale=MatchLowercase}
\fi
% use upquote if available, for straight quotes in verbatim environments
\IfFileExists{upquote.sty}{\usepackage{upquote}}{}
% use microtype if available
\IfFileExists{microtype.sty}{%
\usepackage{microtype}
\UseMicrotypeSet[protrusion]{basicmath} % disable protrusion for tt fonts
}{}
\newif\ifbibliography
\hypersetup{
            pdftitle={Module 4: CLASSIFICATION},
            pdfauthor={Mette Langaas and Julia Debik, Department of Mathematical Sciences, NTNU},
            pdfborder={0 0 0},
            breaklinks=true}
\urlstyle{same}  % don't use monospace font for urls
\usepackage{color}
\usepackage{fancyvrb}
\newcommand{\VerbBar}{|}
\newcommand{\VERB}{\Verb[commandchars=\\\{\}]}
\DefineVerbatimEnvironment{Highlighting}{Verbatim}{commandchars=\\\{\}}
% Add ',fontsize=\small' for more characters per line
\usepackage{framed}
\definecolor{shadecolor}{RGB}{248,248,248}
\newenvironment{Shaded}{\begin{snugshade}}{\end{snugshade}}
\newcommand{\KeywordTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\DataTypeTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{#1}}
\newcommand{\DecValTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\BaseNTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\FloatTok}[1]{\textcolor[rgb]{0.00,0.00,0.81}{#1}}
\newcommand{\ConstantTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\CharTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\SpecialCharTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\StringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\VerbatimStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\SpecialStringTok}[1]{\textcolor[rgb]{0.31,0.60,0.02}{#1}}
\newcommand{\ImportTok}[1]{#1}
\newcommand{\CommentTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\DocumentationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AnnotationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\CommentVarTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\OtherTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{#1}}
\newcommand{\FunctionTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\VariableTok}[1]{\textcolor[rgb]{0.00,0.00,0.00}{#1}}
\newcommand{\ControlFlowTok}[1]{\textcolor[rgb]{0.13,0.29,0.53}{\textbf{#1}}}
\newcommand{\OperatorTok}[1]{\textcolor[rgb]{0.81,0.36,0.00}{\textbf{#1}}}
\newcommand{\BuiltInTok}[1]{#1}
\newcommand{\ExtensionTok}[1]{#1}
\newcommand{\PreprocessorTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textit{#1}}}
\newcommand{\AttributeTok}[1]{\textcolor[rgb]{0.77,0.63,0.00}{#1}}
\newcommand{\RegionMarkerTok}[1]{#1}
\newcommand{\InformationTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\WarningTok}[1]{\textcolor[rgb]{0.56,0.35,0.01}{\textbf{\textit{#1}}}}
\newcommand{\AlertTok}[1]{\textcolor[rgb]{0.94,0.16,0.16}{#1}}
\newcommand{\ErrorTok}[1]{\textcolor[rgb]{0.64,0.00,0.00}{\textbf{#1}}}
\newcommand{\NormalTok}[1]{#1}
\usepackage{longtable,booktabs}
\usepackage{caption}
% These lines are needed to make table captions work with longtable:
\makeatletter
\def\fnum@table{\tablename~\thetable}
\makeatother
\usepackage{graphicx,grffile}
\makeatletter
\def\maxwidth{\ifdim\Gin@nat@width>\linewidth\linewidth\else\Gin@nat@width\fi}
\def\maxheight{\ifdim\Gin@nat@height>\textheight0.8\textheight\else\Gin@nat@height\fi}
\makeatother
% Scale images if necessary, so that they will not overflow the page
% margins by default, and it is still possible to overwrite the defaults
% using explicit options in \includegraphics[width, height, ...]{}
\setkeys{Gin}{width=\maxwidth,height=\maxheight,keepaspectratio}

% Prevent slide breaks in the middle of a paragraph:
\widowpenalties 1 10000
\raggedbottom

\AtBeginPart{
  \let\insertpartnumber\relax
  \let\partname\relax
  \frame{\partpage}
}
\AtBeginSection{
  \ifbibliography
  \else
    \let\insertsectionnumber\relax
    \let\sectionname\relax
    \frame{\sectionpage}
  \fi
}
\AtBeginSubsection{
  \let\insertsubsectionnumber\relax
  \let\subsectionname\relax
  \frame{\subsectionpage}
}

\setlength{\parindent}{0pt}
\setlength{\parskip}{6pt plus 2pt minus 1pt}
\setlength{\emergencystretch}{3em}  % prevent overfull lines
\providecommand{\tightlist}{%
  \setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}
\setcounter{secnumdepth}{0}

\title{Module 4: CLASSIFICATION}
\subtitle{TMA4268 Statistical Learning V2018}
\author{Mette Langaas and Julia Debik, Department of Mathematical Sciences, NTNU}
\date{week 5 2018 (Version 28.01.2018)}

\begin{document}
\frame{\titlepage}

\begin{frame}

\textbf{Learning material for this module:}

\begin{itemize}
\tightlist
\item
  James et al (2013): An Introduction to Statistical Learning. Chapter
  2.2.3, 4.
\end{itemize}

\begin{block}{Move to:}

\begin{itemize}
\tightlist
\item
  \protect\hyperlink{PartA}{Part A: Introduction to classification, and
  modelling class densities}
\item
  \protect\hyperlink{PartB}{Part B: Modelling posterior probabilites}
\item
  \protect\hyperlink{RecEx}{Recommended exercises}
\item
  \protect\hyperlink{further}{Further reading}
\item
  \protect\hyperlink{Rpackages}{Packages to install before knitting this
  R Markdown file}
\end{itemize}

\end{block}

\end{frame}

\begin{frame}{Topics in this module}

\textbf{Part A: Introduction to classification, and modelling class
densities }

\begin{itemize}
\tightlist
\item
  Aim of this module
\item
  What is classification and what is discrimination?
\item
  Loss function and the Bayes classifier
\item
  Modelling class densities

  \begin{itemize}
  \tightlist
  \item
    Linear discriminant analysis
  \item
    Quadratic discriminant analysis
  \item
    Naive Bayes (optional)
  \end{itemize}
\item
  Modelling posterior probabilities

  \begin{itemize}
  \tightlist
  \item
    KNN-classifier
  \end{itemize}
\end{itemize}

\end{frame}

\begin{frame}

\textbf{Part B: Modelling posterior probabilites }

\begin{itemize}
\tightlist
\item
  Modelling posterior probabilities (cont.)

  \begin{itemize}
  \tightlist
  \item
    Linear regression for classification problems
  \item
    Logistic regression
  \end{itemize}
\item
  Sensitivity, specificity, ROC and AUC
\item
  Comparisons
\item
  Extensions
\item
  Quiz (team mode)
\end{itemize}

\textbf{Recommended exercises}

\textbf{Further reading}

\textbf{Rpackages to install}

\end{frame}

\begin{frame}{ Part A}

\begin{block}{What will you learn?}

\begin{itemize}
\tightlist
\item
  What is classification and discrimination?
\item
  What is the Bayes classifier and the Bayes risk?
\item
  What is the sampling paradigm and what is modelled then?
\item
  Linear discriminant analysis: model, method, results.
\item
  Quadratic discriminant analysis: model, method, results.
\item
  Naive Bayes - when and why?
\item
  That is the diagnostic paradigm and what is modelled then?
\item
  KNN - majority vote or estimate posterior class probability?
\end{itemize}

\end{block}

\end{frame}

\begin{frame}[fragile]{Example: Which type of iris species?}

The \texttt{iris} flower data set was introduced by the British
statistician and biologist Ronald Fisher in 1936.

\begin{itemize}
\tightlist
\item
  Three plant species \{setosa, virginica, versicolor\} (50 observation
  of each), and
\item
  four features: \texttt{Sepal.Length}, \texttt{Sepal.Width},
  \texttt{Petal.Length} and \texttt{Petal.Width}.
\end{itemize}

\includegraphics[width=100px]{iris}

Image taken from
\url{http://blog.kaggle.com/2015/04/22/scikit-learn-video-3-machine-learning-first-steps-with-the-iris-dataset/}

\end{frame}

\begin{frame}{Discrimination and classification}

\begin{itemize}
\tightlist
\item
  Discrimination is concerned with the description and the separation of
  the classes and
\item
  classification with the allocation and the assignment of observations
  to predefined classes\\
\item
  The two tasks are closely related and often overlap.
\end{itemize}

Suppose we have a qualititative response value that can be a member in
one of \(K\) classes \(\mathcal{C} = \{c_1, c_2, ..., c_k, ..., c_K\}\).
Further, suppose these elements are numbered \(1, 2, ..., K\).

In classification we build a function \(f(X)\) that takes a vector of
input variables \(X\) and predicts its class membership, such that
\(Y \in \mathcal{C}\).

We would also assess the \emph{uncertainty} in this classification.
Sometimes the role of the different predictors may be of main interest.
Often we are interested in estimating the \emph{probability} that \(X\)
belongs to a class rather than the classification.

In discrimination we might also focus on describing the class boundaries
using \emph{discriminant functions}.

\end{frame}

\begin{frame}{Set-up}

\textbf{Training set:} observations (independent pairs)
\(\{(x_1, y_1), ..., (x_n, y_n)\}\) where the response variable \(Y\) is
\emph{qualitative} and labelled \(1, 2, ..., K\).

The training set is used to construct the classification rule (or
similar).

\textbf{Test set:} observations (independent pairs), same format as the
training set.

The test set is used to evaluate the classification rule.

\textbf{Loss function:}: The misclassifications are given the loss 1 and
the correct classifications loss 0 - this is called \emph{0/1-loss}.
(Quadratic loss is not used for classification.)

But, what if the \emph{cost} of misclassification is different for the
classes? In a more general set-up we could also take that into account.
Read more about this: Decision theoretic framework (sample space, action
space, loss function, risk function, Bayes risk, cost function).

\end{frame}

\begin{frame}{Bayes classifier}

Assume that we know or can estimate the probability that a new
observation \(x_0\) belongs to class \(k\):
\[p_k(x_0) = \text{Pr}(Y=k | X=x_0), \quad k = 1, 2, ... K.\] This is
the probability that \(Y=k\) given the observation \(x_0\). The
\emph{Bayes classifier assigns an observation to the most likely class},
given its predictor values.

It can be proven that the Bayes classifier is the classifier minimizing
the expected 0/1-loss. This classifier is also called the maximum
posterior probability (MAP) classifier.

\end{frame}

\begin{frame}

\textbf{The Bayes classifier}

\begin{itemize}
\tightlist
\item
  has the \emph{smallest test error rate}.
\item
  However, we do not know the conditional distribution of \(Y\) given
  \(X\) for real data. Computing the Bayes classifier is thus impossible
  for real situations.
\item
  The class boundaries using the Bayes classifier is called the
  \emph{Bayes decision boundary}.
\item
  The error rate at \(X=x_0\) is
  \(1-\text{max}\text{Pr}(Y=j\mid X=x_0)\).
\item
  The overall Bayes error rate therefore is given as
  \[1-\text{E}(\text{max} \text{Pr}(Y=j\mid X))\] where the expectation
  is over \(X\).
\item
  The Bayes error rate is comparable to the \emph{irreducible error} in
  the regression setting - we can't go lower!
\end{itemize}

\end{frame}

\begin{frame}

\begin{block}{Synthetic example: what is the Bayes error?}

Suppose we have observations coming from two classes: \{{green},
{orange}\}, where
\[X_{\text{green}}\sim \mathcal{N}(-2, 1.5) \text{ and }
X_{\text{orange}}\sim \mathcal{N}(2, 1.5) \] and that the probability of
observing each class is equal.

\begin{itemize}
\tightlist
\item
  Where is the Bayes decision boundary here?
\item
  How can we calculate the Bayes error rate (intuitively- see the
  graph)?
\item
  What would you estimate the Bayes error rate to be?
\item
  What if someone constructs a classifier and it has a lower error rate
  than the Bayes error rate - is this possible?
\end{itemize}

\end{block}

\end{frame}

\begin{frame}

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-2-1.pdf}

\end{frame}

\begin{frame}[fragile]

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-3-1.pdf}

Bayes error: \texttt{round(0.5*2*pnorm(0,mean=2,sd=1.5),2)}=0.09

\end{frame}

\begin{frame}{Bayes theorem - for discrimination and classification}

Bayes theorem

\begin{align*}
\text{Pr}(Y=k | X=x) &= 
\frac{\text{Pr}(X=x \cap Y=k)}{\text{Pr}(X=x)}\\
&=\frac{\text{Pr}(X=x | Y=k) \text{Pr}(Y=k)}{\text{Pr}(X=x)}\\
&= \frac{\pi_k f_k(x)}{\sum_{l=1}^K \pi_l f_l(x)}
\end{align*}

Here \(f_k(x) = \text{Pr}(X=x | Y=k)\) is the \emph{density} for \(X\)
in class \(k\) and \(\pi_k = \text{Pr}(Y=k)\) is the prior probability
for class \(k\).

\end{frame}

\begin{frame}{Bayes decision rule - two paradigms}

The most popular approach to classication is to use the Bayes decision
rule with the 0/1 loss function= classify to the class with the largest
\(P(Y=k\mid X=x_0)\).

Two approaches:

\begin{itemize}
\item
  The \textbf{diagnostic paradigm}: We focus on \emph{directly}
  estimating the posterior distribution for the classes
  \(\text{Pr}(Y=k \mid X=x)\).
\item
  The \textbf{sampling paradigm}: There focus is on estimating the prior
  probabilities \(\pi_k\) for the classes and the class conditional
  distributions \(f_k(x)\). We classify to the class with the maximal
  product \(\pi_k f_k(x)\).
\end{itemize}

We \emph{first} look at the sampling paradigm - and then we need to
model the pdf for each class. Popular: the multivariate normal
distribution!

\end{frame}

\begin{frame}

\begin{block}{Univariate normal class distributions - the role of the
prior}

Suppose (again) we have observations coming from two classes: \{{green},
{orange}\}, where
\[y_{\text{green}}\sim \mathcal{N}(-2, 1.5) \text{ and }
y_{\text{orange}}\sim \mathcal{N}(2, 1.5) \]

In the figure below we have specified the prior probabilities to be
equal, \(\pi_1 = \pi_2 = 0.5\) and have plotted \(\pi_k f_k(x)\) for the
two classes.

The decision boundary is where the point of intersection of the two
lines is, because here \(\pi_1 f_1(x)=\pi_2 f_2(x)\). Thus all points to
left of the decision boundary will be classified as green and similarly,
all points to the right of the decision boundary will be classified as
orange.

\end{block}

\end{frame}

\begin{frame}

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-4-1.pdf}

\end{frame}

\begin{frame}

We now specify different priors, such that \(\pi_1 = 0.3\) and
\(\pi_2 = 0.7\). We see that this has affected the decision boundary,
which now has shifted to the left.

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-5-1.pdf}

\end{frame}

\begin{frame}{Linear Discriminant Analysis (LDA) (p=1)}

Using Linear discriminant analysis we assume the class conditional
distributions are normal (Gaussian).

The univariate normal pdf
\[f(x) = \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{1}{2}\big(\frac{x-\mu}{\sigma}\big)^2},\]
has parameters \(\mu\) (mean) and \(\sigma\) (standard deviation).

Assume that we observe \(K\) classes, where each class conditional
distribution is normal, where the classes have mean \(\mu_k\) and
standard deviation \(\sigma_k\):
\[f_k(x) = \frac{1}{\sqrt{2\pi}\sigma_k} e^{-\frac{1}{2}\big(\frac{x-\mu_k}{\sigma_k}\big)^2}\]
With LDA we assume that all of the classes have the \emph{same standard
deviation} \(\sigma_k = \sigma\).

\end{frame}

\begin{frame}

We can insert the expression for each class distribution into Bayes
formula to obtain the posterior probability
\(p_k(x) = \text{Pr}(Y = k | X = x)\):
\[p_k(x) = \frac{\text{Pr}(X=x | Y=k) \text{Pr}(Y=k)}{\text{Pr}(X=x)}=\frac{\pi_k \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{1}{2}\big(\frac{x-\mu_k}{\sigma}\big)^2}}{\sum_{l=1}^K \pi_l \frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{1}{2}\big(\frac{x-\mu_l}{\sigma}\big)^2}} \]

Our rule is to classify to the class for which \(p_k(x)\) is largest. It
can be shown that this is equivalent to assigning \(x\) to the class
with the largest \emph{discriminant score} \(\delta_k(x)\):
\[\delta_k(x) = x\cdot \frac{\mu_k}{\sigma^2} - \frac{\mu_k^2}{2 \sigma^2}+\log(\pi_k).\]
This decision boundaries between the classes are \emph{linear} in \(x\).

\textbf{Q:} So, what do we need to use this result in practice?

\end{frame}

\begin{frame}

\begin{block}{Parameter estimators}

In real life situations we will not know the prior probabilities, the
class mean or standard deviation=need parameter estimators.

\begin{itemize}
\item
  Prior probability for class \(k\) is (often) estimated by taking the
  fraction of observations \(n_k\) (out of \(n\)) coming from class
  \(k\): \(\hat{\pi}_k = \frac{n_k}{n}.\)
\item
  The mean value for class \(k\) is simply the sample mean of all
  observations from class \(k\):
  \[\hat{\mu}_k = \frac{1}{n_k}\sum_{i:y_i=k} x_i.\]
\item
  The standard deviation: sample standard deviation across all classes:
  \[\hat{\sigma}^2=\frac{1}{n-K}\sum_{k=1}^K \sum_{i: y_i=k} (x_i-\hat{\mu}_k)^2 = \sum_{k=1}^K \frac{n_k - 1}{n - K} \cdot \hat{\sigma}_k^2.\]
  \(\hat{\sigma}_k\): estimated standard deviation of all observations
  from class \(k\).
\end{itemize}

\end{block}

\end{frame}

\begin{frame}[fragile]

\begin{block}{How would the estimation affect our misclassification
rate?}

If \(\mu_1=-2\), \(\mu_2=2\), \(\sigma=1.5\), \(\pi_1=\pi_2=0.5\) we
found that the class boundary is at 0 and the Bayes error is
\texttt{round(2*0.5*pnorm(0,2,1.5))}=0.09.

But, in a real life situation

\begin{itemize}
\tightlist
\item
  we estimate the class boundary
\item
  and we do not know the true distribution for the classes.
\end{itemize}

How can we then estimate the goodness of our estimator?

\end{block}

\end{frame}

\begin{frame}[fragile]

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  Use the training set to estimate parameters and class boundary:
\item
  Use the test set to estimate misclassification rate.
\end{enumerate}

\footnotesize

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{n=}\DecValTok{1000}\NormalTok{;pi1=pi2=}\FloatTok{0.5}\NormalTok{;mu1=}\OperatorTok{-}\DecValTok{2}\NormalTok{;mu2=}\DecValTok{2}\NormalTok{;sigma=}\FloatTok{1.5}\NormalTok{;}\KeywordTok{set.seed}\NormalTok{(}\DecValTok{1}\NormalTok{)}
\NormalTok{n1train=}\KeywordTok{rbinom}\NormalTok{(}\DecValTok{1}\NormalTok{,n,pi1);n2train=n}\OperatorTok{-}\NormalTok{n1train}
\NormalTok{n1test=}\KeywordTok{rbinom}\NormalTok{(}\DecValTok{1}\NormalTok{,n,pi1);n2test=n}\OperatorTok{-}\NormalTok{n1test}
\NormalTok{train1=}\KeywordTok{rnorm}\NormalTok{(n1train,mu1,sigma);train2=}\KeywordTok{rnorm}\NormalTok{(n2train,mu2,sigma)}
\NormalTok{test1=}\KeywordTok{rnorm}\NormalTok{(n1test,mu1,sigma);test2=}\KeywordTok{rnorm}\NormalTok{(n2test,mu2,sigma)}
\NormalTok{sigma2.}\DecValTok{1}\NormalTok{=}\KeywordTok{var}\NormalTok{(train1);sigma2.}\DecValTok{2}\NormalTok{=}\KeywordTok{var}\NormalTok{(train2)}
\NormalTok{estsigma2=((n1train}\OperatorTok{-}\DecValTok{1}\NormalTok{)}\OperatorTok{*}\NormalTok{sigma2.}\DecValTok{1}\OperatorTok{+}\NormalTok{(n2train}\OperatorTok{-}\DecValTok{1}\NormalTok{)}\OperatorTok{*}\NormalTok{sigma2.}\DecValTok{2}\NormalTok{)}\OperatorTok{/}\NormalTok{(n}\OperatorTok{-}\DecValTok{2}\NormalTok{)}

\NormalTok{rule=}\FloatTok{0.5}\OperatorTok{*}\NormalTok{(}\KeywordTok{mean}\NormalTok{(train1)}\OperatorTok{+}\KeywordTok{mean}\NormalTok{(train2))}\OperatorTok{+}
\StringTok{  }\NormalTok{estsigma2}\OperatorTok{*}\NormalTok{(}\KeywordTok{log}\NormalTok{(n2train}\OperatorTok{/}\NormalTok{n)}\OperatorTok{-}\KeywordTok{log}\NormalTok{(n1train}\OperatorTok{/}\NormalTok{n))}\OperatorTok{/}\NormalTok{(}\KeywordTok{mean}\NormalTok{(train1)}\OperatorTok{-}\KeywordTok{mean}\NormalTok{(train2))}

\KeywordTok{c}\NormalTok{((}\KeywordTok{sum}\NormalTok{(train1}\OperatorTok{>}\NormalTok{rule)}\OperatorTok{+}\KeywordTok{sum}\NormalTok{(train2}\OperatorTok{<}\NormalTok{rule))}\OperatorTok{/}\NormalTok{n, (}\KeywordTok{sum}\NormalTok{(test1}\OperatorTok{>}\NormalTok{rule)}\OperatorTok{+}\KeywordTok{sum}\NormalTok{(test2}\OperatorTok{<}\NormalTok{rule))}\OperatorTok{/}\NormalTok{n)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] 0.105 0.115
\end{verbatim}

\normalsize

\end{frame}

\begin{frame}

\begin{block}{Training error rate}

the proportion of mistakes that are made if we apply classifier
\(\hat{f}\) to the training observations, i.e.
\(\hat{y}_i=\hat{f}(x_i)\).
\[\frac{1}{n}\sum_{i=1}^n \text{I}(y_i \neq \hat{y}_i).\] Here I is the
indicator function (to give our 0/1 loss) which is defined as:
\[\text{I}(a=\hat{a}) = \begin{cases} 1 \text{ if } a = \hat{a} \\ 0 \text{ else } \end{cases}\]
The indicator function counts the number of times our model has made a
wrong classification. The training error rate is the fraction of
misclassifications made on our training set.

A very low training error rate may imply overfitting.

\end{block}

\end{frame}

\begin{frame}

\begin{block}{Test error rate}

Here the fraction of misclassifications is calculated when our model is
applied on a test set. From what we have learned about regression we can
deduct that this gives a better indication of the true performance of
the classifier (than the traning error).
\[\text{Ave}(I(y_0-\hat{y}_0))\] where the average is over all the test
observations \((x_0,y_0)\).

We assume that a \emph{good} classifier is a classifier that has a
\emph{low} test error.

\end{block}

\end{frame}

\begin{frame}

\begin{block}{The confusion matrix}

The confusion matrix is a table that can show the performance of
classifier, given that the true values are known.

We can make a confusion matrix from the training or test set, but will
in most cases do so for the test set only.

The rows represent the true classes, while the columns represent the
predicted classes. Inside the table we have counts (just labelled
``correct'' and ``wrong'' below - but should be numbers). The sum of the
diagonal is the total number of correct classifications. The sum of all
elements off the diagonal is the total number of misclassifications.

\begin{longtable}[]{@{}lllll@{}}
\toprule
& Predicted 1 & Predicted 2 & \ldots{} & Predicted K\tabularnewline
True 1 & correct & wrong & \ldots{} & wrong\tabularnewline
True 2 & wrong & correct & \ldots{} & wrong\tabularnewline
\ldots{} & \ldots{} & \ldots{} & \ldots{} & \ldots{}\tabularnewline
True K & wrong & wrong & \ldots{} & correct\tabularnewline
\bottomrule
\end{longtable}

\end{block}

\end{frame}

\begin{frame}[fragile]

The confusion matrix can be obtained in \texttt{R} by using the
\texttt{table} function with the predicted classes and the true classes
need to be given as function arguments.

We will soon come back to the special case of two classes - where we may
think of this as ``-'' (non disease) and ``+'' (disease).

\end{frame}

\begin{frame}

\begin{block}{LDA when p\textgreater{}1}

Linear discriminant analysis can be generalized to situations when
\(p>1\) covariates are used. The decision boundaries are still linear.

The multivariate normal distribution function:
\[f(x) = \frac{1}{(2 \pi)^{p/2}|\boldsymbol{\Sigma}|^{1/2}}\exp({-\frac{1}{2}({\bf x}-\boldsymbol\mu)^T \boldsymbol{\Sigma}^{-1}({\bf x}-\boldsymbol\mu)})\]
This gives the following expression for the discriminant function:
\[\delta_k(x) = {\bf x}^T \boldsymbol{\Sigma}^{-1}\boldsymbol\mu_k - \frac{1}{2}\boldsymbol\mu_k^T \boldsymbol{\Sigma}^{-1}\boldsymbol\mu_k + \log \pi_k.\]

Compulsory Exercise 1, Problem 3: you are asked to derive this
expression for K=2!

\end{block}

\end{frame}

\begin{frame}

\begin{block}{Posterior probabilites}

Sometimes the probability that an observation comes from a class \(k\)
is more interesting than the actual classification itself. These class
probabilities can be estimated from the priors and class conditional
distributions, or from the discriminant functions:

\begin{align*}\hat{\text{Pr}}(Y=k | X=x)&=
\frac{\hat{\pi}_k \cdot \frac{1}{(2 \pi)^{p/2}|\hat{\boldsymbol{\Sigma}}|^{1/2}} \exp(-\frac{1}{2}
({\bf x}-\hat{\boldsymbol\mu}_k)^T \hat{\boldsymbol{\Sigma}}^{-1}
({\bf x}-\hat{\boldsymbol\mu}_k))}
{\sum_{l=1}^K \hat{\pi}_l 
\frac{1}{(2 \pi)^{p/2}|\hat{\boldsymbol{\Sigma}}|^{1/2}}
\exp(-\frac{1}{2}
({\bf x}-\hat{\boldsymbol\mu}_l)^T 
\hat{\boldsymbol{\Sigma}}^{-1}
({\bf x}-\hat{\boldsymbol\mu}_l))}\\
&=
\frac{e^{\hat{\delta}_k(x)}}{\sum_{l=1}^K e^{\hat{\delta}_l(x)}}.\end{align*}

\end{block}

\end{frame}

\begin{frame}

\begin{block}{Quadratic Discriminant Analysis (QDA)}

In quadratic discriminant analysis we assume that the distributions of
the classes is Gaussian, but with covariance matrix
\(\boldsymbol{\Sigma}_k\) for each class.

The discriminant functions are now given by:

\begin{align*} \delta_k(x) &= -\frac{1}{2}(x-\mu_k)^T \boldsymbol{\Sigma}_k^{-1}(x-\mu_k)-\frac{1}{2}\log |\boldsymbol{\Sigma}_k| + \log \pi_k \\ &= -\frac{1}{2} x^T \boldsymbol{\Sigma}_k^{-1}x + x^T \boldsymbol{\Sigma}_k^{-1}\mu_k - \frac{1}{2} \mu_k^T \boldsymbol{\Sigma}_k^{-1}\mu_k - \frac{1}{2}\log |\boldsymbol{\Sigma}_k | + \log \pi_k.\end{align*}

These decision boundaries are \emph{quadratic} functions of \(x\).

\end{block}

\end{frame}

\begin{frame}

\begin{block}{LDA vs QDA}

In QDA we are allowed for the covariance matrices to be different for
the classes, but for LDA they are assumed to be the same, so QDA is more
flexible than LDA.

\textbf{Q:}

\begin{itemize}
\tightlist
\item
  But, if the covariance matrices in theory are equal - will they not be
  estimated equal?
\item
  Should we not always prefer QDA to LDA?
\end{itemize}

\end{block}

\end{frame}

\begin{frame}

\textbf{A:}

Bias-variance trade-off:

If the assumption of equal covariance matrices

\begin{itemize}
\tightlist
\item
  is wrong, then LDA may suffer from high bias for the parameter
  estimators.
\item
  is correct, then QDA is better off. But, for small sample sizes the
  covariance matrices might be poorly estimated (high variance of
  estimators).
\end{itemize}

If the number of covariates is high:

\begin{itemize}
\tightlist
\item
  then QDA requires estimating \(K*p*(p+1)/2\) parameters,
\item
  while LDA only requires \(p*(p+1)/2\).
\end{itemize}

Thus, LDA is less flexible than QDA and might therefore have much less
variance.

\end{frame}

\begin{frame}{Naive Bayes - when we have many covariates (optional)}

In Naive Bayes (Idiot's Bayes) we assume that each class density is the
product of marginal densities - i.e.~inputs are conditionally
independent in each class. \[f_k(x)=\prod_{j=1}^p f_{kj}(x_j)\] This is
generally not true, but it simplifies the estimation dramatically.

The original naive Bayes used univariate normal marginal distributions,
but generalizations can be made.

Instead of estimating a full covariance matrix, only the diagonal
elements are estimated.

This method often produces good results, even though the joint pdf is
not the product of the marginal pdf. This might be because we are not
focussing on estimation of class pdfs, but class boundaries.

\end{frame}

\begin{frame}[fragile]

\begin{block}{Example: Classification of iris plants}

We will use \texttt{sepal\ width} and \texttt{sepal\ length} to build a
classificator, here 50 observations from each class available.

\begin{longtable}[]{@{}rrrrl@{}}
\toprule
Sepal.Length & Sepal.Width & Petal.Length & Petal.Width &
Species\tabularnewline
\midrule
\endhead
5.1 & 3.5 & 1.4 & 0.2 & setosa\tabularnewline
4.9 & 3.0 & 1.4 & 0.2 & setosa\tabularnewline
4.7 & 3.2 & 1.3 & 0.2 & setosa\tabularnewline
4.6 & 3.1 & 1.5 & 0.2 & setosa\tabularnewline
5.0 & 3.6 & 1.4 & 0.2 & setosa\tabularnewline
5.4 & 3.9 & 1.7 & 0.4 & setosa\tabularnewline
\bottomrule
\end{longtable}

\end{block}

\end{frame}

\begin{frame}

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-8-1.pdf}

\end{frame}

\begin{frame}

\begin{block}{Iris: LDA}

\includegraphics{4Classif_files/figure-beamer/irislda-1.pdf}

\end{block}

\end{frame}

\begin{frame}

\begin{block}{Iris: QDA}

\includegraphics{4Classif_files/figure-beamer/irisqda-1.pdf}

\end{block}

\end{frame}

\begin{frame}[fragile]

\begin{block}{Iris: compare LDA and QDA}

We now want to compare the predictive performance of our two
classifiers. We do this by dividing the original \texttt{iris} data set,
randomly, into train and test samples of equal size:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{set.seed}\NormalTok{(}\DecValTok{1}\NormalTok{)}
\NormalTok{train =}\StringTok{ }\KeywordTok{sample}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\DecValTok{150}\NormalTok{, }\DecValTok{75}\NormalTok{)}

\NormalTok{iris_train =}\StringTok{ }\NormalTok{iris[train, ]}
\NormalTok{iris_test =}\StringTok{ }\NormalTok{iris[}\OperatorTok{-}\NormalTok{train, ]}

\NormalTok{iris_lda2 =}\StringTok{ }\KeywordTok{lda}\NormalTok{(Species}\OperatorTok{~}\NormalTok{Sepal.Length }\OperatorTok{+}\StringTok{ }\NormalTok{Sepal.Width, }
                \DataTypeTok{data=}\NormalTok{iris_train, }
                \DataTypeTok{prior=}\KeywordTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{,}\DecValTok{1}\NormalTok{,}\DecValTok{1}\NormalTok{)}\OperatorTok{/}\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\end{block}

\end{frame}

\begin{frame}[fragile]

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# Training error}
\KeywordTok{table}\NormalTok{(}\KeywordTok{predict}\NormalTok{(iris_lda2, }\DataTypeTok{newdata=}\NormalTok{iris_train)}\OperatorTok{$}\NormalTok{class, }
\NormalTok{      iris_train}\OperatorTok{$}\NormalTok{Species)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##             
##              setosa versicolor virginica
##   setosa         25          0         0
##   versicolor      1         21         7
##   virginica       0          7        14
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# Test error}
\NormalTok{iris_lda2_predict =}\StringTok{ }\KeywordTok{predict}\NormalTok{(iris_lda2, }\DataTypeTok{newdata=}\NormalTok{iris_test)}
\KeywordTok{table}\NormalTok{(iris_lda2_predict}\OperatorTok{$}\NormalTok{class, iris}\OperatorTok{$}\NormalTok{Species[}\OperatorTok{-}\NormalTok{train])}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##             
##              setosa versicolor virginica
##   setosa         24          0         0
##   versicolor      0         13         6
##   virginica       0          9        23
\end{verbatim}

\end{frame}

\begin{frame}

The LDA classifier has a training error rate of 15/75, that is 20 \%.
When tested on a new set it correctly classified 24+13+23 times and
misclassified 15 times. This gives a misclassification rate of
\[\text{Test error}_\text{LDA}  = \frac{15}{75} =0.2.\]

Using a different division into training and test set will give (small)
changes to these numbers.

\end{frame}

\begin{frame}[fragile]

\begin{block}{Iris: training and test error for QDA}

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{iris_qda2 =}\StringTok{ }\KeywordTok{qda}\NormalTok{(Species}\OperatorTok{~}\NormalTok{Sepal.Length }\OperatorTok{+}\StringTok{ }\NormalTok{Sepal.Width, }\DataTypeTok{data=}\NormalTok{iris_train, }
                \DataTypeTok{prior=}\KeywordTok{c}\NormalTok{(}\DecValTok{1}\NormalTok{,}\DecValTok{1}\NormalTok{,}\DecValTok{1}\NormalTok{)}\OperatorTok{/}\DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

Important: use the same division into training and test set for methods
we want to compare.

\end{block}

\end{frame}

\begin{frame}[fragile]

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# Training error}
\KeywordTok{table}\NormalTok{(}\KeywordTok{predict}\NormalTok{(iris_qda2, }\DataTypeTok{newdata=}\NormalTok{iris_train)}\OperatorTok{$}\NormalTok{class, iris_train}\OperatorTok{$}\NormalTok{Species)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##             
##              setosa versicolor virginica
##   setosa         26          0         0
##   versicolor      0         22         6
##   virginica       0          6        15
\end{verbatim}

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# Test error}
\NormalTok{iris_qda2_predict =}\StringTok{ }\KeywordTok{predict}\NormalTok{(iris_qda2, }\DataTypeTok{newdata=}\NormalTok{iris_test)}
\KeywordTok{table}\NormalTok{(iris_qda2_predict}\OperatorTok{$}\NormalTok{class, iris}\OperatorTok{$}\NormalTok{Species[}\OperatorTok{-}\NormalTok{train])}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##             
##              setosa versicolor virginica
##   setosa         24          0         0
##   versicolor      0         14        10
##   virginica       0          8        19
\end{verbatim}

\end{frame}

\begin{frame}

The QDA classifier has a training error rate of \(16\%\). When tested on
a new set, the misclassification error rate was
\[\text{Test error}_\text{QDA}  = \frac{18}{75}=.24\]

The LDA classifier has given the smallest test error for classifying
iris plants based on sepal width and sepal length for our test set and
should be preferred in this case.

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  Would another division of the data into training and test set give the
  same conclusion?
\end{enumerate}

We will look into other choice than dividing into one training and one
test set in Module 5 (crossvalidation).

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\setcounter{enumi}{1}
\tightlist
\item
  What about the other two covariates? Would adding them to the model (4
  covariates) give a better classification rule?
\end{enumerate}

\end{frame}

\begin{frame}

\begin{block}{Optional: Fishers initial idea}

In 1936 R. A. Fisher developed LDA.

\begin{itemize}
\tightlist
\item
  His aim was to find a linear combination of the explanatory variables
  which \emph{maximized the ratio of its between group to within group
  variance}.
\item
  In this way the observations are transformed so that they are
  separated as much as possible.
\item
  His approach has the advantage that it is suited for visual inspection
  and graphical description , it ``separates'' the populations.
\end{itemize}

Let the between-class variance be denoted
\(\boldsymbol{B}=\sum_{m=1}^{M}(\boldsymbol{\mu}_{m}- \bar{\boldsymbol{\mu}})(\boldsymbol{\mu}_{m}-\bar{\boldsymbol{\mu}})^{T}\),
where \(\boldsymbol{\mu}_{m}\) denotes the mean of class \(\omega_{m}\)
and \(\bar{\boldsymbol{\mu}}\) the overall mean.

The within-class variance is assumed equal for all classes and is
denoted \(\boldsymbol{\Sigma}\) (\(\boldsymbol{\Sigma}\) is assumed to
have full rank).

The linear combination \(\boldsymbol{l}^{T}\boldsymbol{X}\) that
maximize
\({\boldsymbol{l}^{T}\boldsymbol{B}\boldsymbol{l}}/{\boldsymbol{l}^{T}\boldsymbol{\Sigma}\boldsymbol{l}}\)
under the constraint that
\(\boldsymbol{l}^{T}\boldsymbol{\Sigma}\boldsymbol{l}=1\) is found to be
the scaled eigenvectors of \(\boldsymbol{\Sigma}^{-1}\boldsymbol{B}\)
corresponding to the nonzero eigenvalues of
\(\boldsymbol{\Sigma}^{-1}\boldsymbol{B}\). The eigenvector
corresponding to the largest eigenvalue defines the first discriminant
\(\boldsymbol{l}_{1}^{T}\boldsymbol{X}\). The second linear discriminant
\(\boldsymbol{l}_{2}^{T}\boldsymbol{X}\) is constructed from the
eigenvector corresponding to the second largest eigenvalue and so on.
(We also have
\(Cov(\boldsymbol{l}_{j}^{T}\boldsymbol{X},\boldsymbol{l}_{i}^{T}\boldsymbol{X})=0\)
for \(i \ne j\) and \(Var(\boldsymbol{l}_{j}^{T}\boldsymbol{X})=1\).)

The number of linear discriminants equals the number of nonzero
eigenvalues, which is smaller than or equal to \(\min(D,M-1)\).
Observations are assigned to the class of the nearest (Euclidean
distance) class mean in the discriminant space.

This equals classification to the nearest Mahalanobis distance
population mean in the input space. Again, the parameters
\(\boldsymbol{\mu}_{i}\) and \(\boldsymbol{\Sigma}\) and the
between-class variance \(\boldsymbol{B}\) are usually unavailable.
Replacing the parameters by estimates from the training set leads to
Fisher's sample linear discriminants. \cite{bdr:natoasi} advocates the
use of robust estimates.

Classification using Fisher's linear discriminant (LDA) is equivalent to
using (\ref{eq:lda}) with equal prior probabilities
\(\pi(\omega_{i})\)=\(\pi\). To see this, \cite{johnswich}, as above
\(\boldsymbol{l}_{j}^{T}\) is the \(j\)th linear discriminant and
eigenvector of \(\boldsymbol{\Sigma}^{-1}\boldsymbol{B}\). Let
\(\boldsymbol{e}_{j}=\boldsymbol{\Sigma}^{\frac{1}{2}}\boldsymbol{l}_{j}\)
be a scaled eigenvector of
\(\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{B}\boldsymbol{\Sigma}^{-\frac{1}{2}}\)
and \(\boldsymbol{E}\) be the orthogonal matrix of eigenvectors
\(\boldsymbol{e}_{j}\). Using LDA an observation \(\boldsymbol{x}\) is
assigned to the class of the nearest class mean in the discriminant
space.
\[ \sum_{j=1}^{D} (\boldsymbol{l}_{j}^{T}\boldsymbol{x}-\boldsymbol{l}_{j}^{T}\boldsymbol{\mu}_{i})^{2} =
\sum_{j=1}^{D} \boldsymbol{l}_{j}^{T}(\boldsymbol{x}-\boldsymbol{\mu}_{i})^{2} \]
\[
=(\boldsymbol{x}-\boldsymbol{\mu}_{i})^{T}\boldsymbol{\Sigma}^{- \frac{1}{2}} 
\boldsymbol{E}\boldsymbol{E}^{T}\boldsymbol{\Sigma}^{- \frac{1}{2}}(\boldsymbol{x}-\boldsymbol{\mu}_{i})
=(\boldsymbol{x}-\boldsymbol{\mu}_{i})^{T}\boldsymbol{\Sigma}^{-1}(\boldsymbol{x}-\boldsymbol{\mu}_{i})\]
Removing terms independent of \(i\) leads to the linear discriminant
function \textbackslash{}
\(\boldsymbol{\mu}_{i}^{T}\boldsymbol{\Sigma}^{-1}\boldsymbol{x}- \frac{1}{2}\boldsymbol{\mu}_{i}^{T}\boldsymbol{\Sigma}^{-1}\boldsymbol{\mu}_{i}\)
which equals equation (\ref{eq:lda}) when the prior probabilities of the
classes are equal.

\end{block}

\end{frame}

\begin{frame}{Bayes decision rule - over to diagnostic paradigm}

Remember:

\begin{itemize}
\item
  The \textbf{diagnostic paradigm}: We focus on \emph{directly}
  estimating the posterior distribution for the classes
  \(\text{Pr}(Y=k \mid X=x)\).
\item
  The \textbf{sampling paradigm}: There focus is on estimating the prior
  probabilities for the classes and the class conditional distributions.
  We classify to the class with the maximal product \(\pi_k f_k(x)\).
\end{itemize}

Now we move to the \emph{diagnostic paradigm} and the \(K\)-nearest
neighbor classifier.

The \(K\)-nearest neighbour classifier estimates
\(\text{Pr}(Y=k \mid X=x)\) and classifies a new observation based on
this estimated probability

\end{frame}

\begin{frame}

\begin{block}{Synthetic example for KNN classification}

Consider a two-class example, and equal class prior probabilites.

A new observation \(x_0\) will be classified to \(A\) if
\(\text{Pr}(Y=A | X=x_0) > 0.5\) and to class \(B\) otherwise.

\begin{itemize}
\tightlist
\item
  The figure below shows a plot of 100 observations from two classes
  \(A\) (red dots) and \(B\) (turquoise dots),
\item
  simulated from a bivariate normal distribution with mean vectors
  \(\mu_A = (1, 1)^T\) and \(\mu_B = (3, 3)^T\) and a covariance matrix
  \(\Sigma_A = \Sigma_B = \begin{pmatrix} 2\hspace{2mm} 0 \\ 0 \hspace{2mm} 2 \end{pmatrix}\).
\item
  We want to find a rule to classify a new observation to class \(A\) or
  \(B\).
\end{itemize}

\end{block}

\end{frame}

\begin{frame}

\includegraphics{4Classif_files/figure-beamer/synthAB-1.pdf}

\end{frame}

\begin{frame}

\begin{block}{K-nearest neighbour classifier}

(warning: K is not the number of classes, but neighbours\ldots{})

The \(K\)-nearest neighbour classifier (KNN) works in the following way:

\begin{itemize}
\tightlist
\item
  Given a new observation \(x_0\) it searches for the \(K\) points in
  our training data that are closest to it.
\item
  These points make up the neighborhood of \(x_0\), \(\mathcal{N}_0\).
\item
  The point \(x_0\) is classified by taking a majority vote of the
  neighbors.
\item
  This means that KNN estimate the posterior class probability as:
  \[\hat{\text{Pr}}(Y=j | X = x_0) = \frac{1}{K} \sum_{i \in \mathcal{N}_0} I(y_i = j).\]
\end{itemize}

\end{block}

\end{frame}

\begin{frame}

\begin{block}{Synthetic data for KNN - continued}

\begin{itemize}
\tightlist
\item
  Assume we have a new observation \(X_0 = (x_{01}, x_{02})^T\) which we
  want to classify as belonging to the class \(A\) or \(B\).
\item
  To illustrate this problem we fit the \(K\)-nearest neighbor
  classifier to our simulated data set with \(K = 1, 3, 10\) and \(150\)
  and observe what happens.
\end{itemize}

In our plots, the small colored dots show the predicted classes for an
evenly-spaced grid. The lines show the decision boundaries. If our new
observation falls into the region within the red decision boundary, it
will be classified as \(A\). If it falls into the region within the
green decision boundary, it will be classified as \(B\).

\end{block}

\end{frame}

\begin{frame}

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-11-1.pdf}

\end{frame}

\begin{frame}

We see that the choice of \(K\) has a big influence on the result of our
classification. By choosing \(K=1\) the classification is made to the
same class as the one nearest neighbor. When \(K=3\) a majority vote is
taken among the three nearest neighbors, and so on. We see that as \(K\)
gets very large, the decision boundary tends towards a straight line
(which is the Bayes boundary in this set-up).

To find the optimal value of \(K\) the typical procedure is to try
different values of \(K\) and then test the predictive power of the
different classifiers, for example by cross-validation, which will be
discussed in Module 5.

We see that after trying all choices for \(K\) between 1 and 50, we see
that a few choices of \(K\) gave the smallest misclassification error
rate, estimating by leave-one out cross-validation (Leave-one-out
cross-validation will be discussed in Module 5). The smallest error rate
is equal to 0.165. This means that the classifier makes a
misclassification 16.5\% of the time and a correct classification 83.5\%
of the time.

\end{frame}

\begin{frame}

\includegraphics{4Classif_files/figure-beamer/knnerror-1.pdf}

\end{frame}

\begin{frame}

This above example showed the bias-variance trade-off in a
classification setting. Choosing a value of \(K\) amounts to choosing
the correct level of flexibility of the classifier. This again is
critical to the success of the classifier. A too low value of \(K\) will
give a very flexible classifier (with high variance and low bias) which
will fit the training set too well (it will overfit) and make poor
predictions for new observations. Choosing a high value for \(K\) makes
the classifier loose its flexibility and the classifier will have low
variance but high bias.

\end{frame}

\begin{frame}

\begin{block}{The curse of dimensionality}

The nearest neighbor classifier can be quite good if the number of
predictor \(p\) is small and the number of observations \(n\) is large.
We need enough close neighbors to make a good classification.

The effectiveness of the KNN classifier falls quickly when the dimension
of the preditor space is high. This is because the nearest neighbors
tend to be far away in high dimensions and the method no longer is
local. This is referred to as the \emph{curse of dimensionality}.

\end{block}

\end{frame}

\begin{frame}{What did we cover today - and what is left?}

\begin{itemize}
\tightlist
\item
  What is classification and discrimination?
\item
  What is the Bayes classifier and the Bayes risk?
\item
  What is the sampling paradigm and what is modelled then?
\item
  Linear discriminant analysis: model, method, results.
\item
  Quadratic discriminant analysis: model, method, results.
\item
  Naive Bayes - when and why?
\item
  That is the diagnostic paradigm and what is modelled then?
\item
  KNN - majority vote or estimate posterior class probability?
\end{itemize}

Part B: more on methods that model and estimate the posterior class
probabilites.

\end{frame}

\begin{frame}{ Part B: Modelling posterior probabilites}

What to remember from Part A?

\end{frame}

\begin{frame}

\begin{block}{Can linear regression be used for classification?}

There are two main reasons for why linear regression is badly suited for
classification problems:

\begin{enumerate}
\def\labelenumi{\arabic{enumi}.}
\tightlist
\item
  Linear regression can produce probabilities smaller than zero or
  larger than one.
\item
  Linear regression cannot be used when the response variable has more
  than two possible classes.
\end{enumerate}

We illustrate the problems of using linear regression for classification
by examples.

\end{block}

\end{frame}

\begin{frame}[fragile]

\textbf{Example 1:} This example uses the \texttt{Default} data set from
the \texttt{ISLR} package. Suppose we want to predict if a new customer
will \texttt{default} or not based on his/her \texttt{balance} or
\texttt{income}. We try to model this using a simple linear regression
and a dummy variable:
\[Y = \begin{cases} 1 \quad \text{ if default = "Yes"} \\ 0 \quad \text{ if default = "No"} \end{cases}.\]
It would be tempting to do the classification according to the rule:
classify as \emph{yes} if \(\hat{Y}>0.5\), else as \emph{no}.

\end{frame}

\begin{frame}

\begin{longtable}[]{@{}llrr@{}}
\toprule
default & student & balance & income\tabularnewline
\midrule
\endhead
No & No & 729.5265 & 44361.625\tabularnewline
No & Yes & 817.1804 & 12106.135\tabularnewline
No & No & 1073.5492 & 31767.139\tabularnewline
No & No & 529.2506 & 35704.494\tabularnewline
No & No & 785.6559 & 38463.496\tabularnewline
No & Yes & 919.5885 & 7491.559\tabularnewline
\bottomrule
\end{longtable}

\end{frame}

\begin{frame}

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-12-1.pdf}

\end{frame}

\begin{frame}[fragile]

The above plots shows \texttt{default} as a function of \texttt{balance}
and \texttt{default} as a function of \texttt{income} with corresponding
fitted linear regression lines (red for \texttt{x=balance} and orange
for \texttt{x=income}). Notice that linear regression in this case
produces probabilities that are zero or bigger than one. It does not
estimate the probabilities well. This can be avoided by using logistic
regression.

\end{frame}

\begin{frame}[fragile]

Logistic regression ensures that the estimated probabilities lie in the
interval between 0 and 1. This is illustrated in the figure below. The
blue line shows the fitted line when performing logistic regression on
\texttt{default} as a function of \texttt{balance}.

\includegraphics{4Classif_files/figure-beamer/defaultglm-1.pdf}

\end{frame}

\begin{frame}

\textbf{Example 2: } Suppose we want to classify the a film. We have
defined three classes: \{ drama, comedy, science-fiction\}. We could try
to model this using linear regression and the following coding:
\[Y = \begin{cases} 1 \quad \text{if drama}, \\ 2 \quad \text{if comedy}, \\ 3 \quad \text{if science-fiction}.\end{cases}\]
However, this coding implies an ordering of the variables and that the
difference between the classes is equal. There is in general no natural
way to code a quantitative variable with more than two classes such that
it can be used with linear regression.

\end{frame}

\begin{frame}{Logistic regression}

Logistic regression models the \emph{probability} that \(Y\) belongs to
a particular class (\(\mathcal{C} = \{1, 0\}\) or \{success, failure\})
given one or many covariates \(X\).
\[Y_i = \begin{cases} 1 \text{ with probability } p(X), \\ 0 \text{ with probability } 1-p(X). \end{cases}\]
In order to ensure that the probabilities lie in the interval between 0
and 1, for all values of \(X\), the \emph{logistic function} is used. In
the case of one covariate, the logistic function has the form:
\[ \text{Pr}(Y=1 | X=x) = p(X) = \frac{e^{\beta_0+\beta_1 X}}{1 + e^{\beta_0 + \beta_1 X}}.\]
The form of this function is illustrated in the figure in example 1 (the
blue line): It is S-shaped, and ranges between 0 and 1. The parameter
\(\beta_1\) determines the rate of increase or decrease of the S-shaped
curve, and the sign indicates whether the curve ascends or descends.

\end{frame}

\begin{frame}

\begin{block}{The odds ratio}

By manipulating the logistic function, we can write it in a form which
is referred to as the \emph{odds ratio}:
\[\frac{p(X)}{1-p(X)} = e^{\beta_0 + \beta_1 X}.\] The above equation
gives the odds of the case \(Y=1\). It is the probability of success
(\(Y=1\)) divided by the probability of failure (\(Y=0\)). The odds
increase multiplicatively by \(e^{\beta_1}\) for every one-unit increase
in \(X\).

\end{block}

\begin{block}{The logit link function}

By taking logarithms on both sides of the odds ratio we obtain the logit
link function:
\[\eta(X) = \log \Big (\frac{p(X)}{1-p(X)} \Big) = \beta_0 + \beta_1 X.\]
This function is linear in \(X\).

\end{block}

\end{frame}

\begin{frame}

\begin{block}{Maximum Likelihood}

We assume that pairs of covariates and responses \(\{x_i, y_i\}\) are
measured independently of each other. Given \(n\) such observation
pairs, the likelihood function of a logistic regression model can be
written as:
\[L(\boldsymbol{\beta}) = \prod_{i=1}^n L_i(\boldsymbol{\beta}) = \prod_{i=1}^n f(y_i; \boldsymbol{\beta}) = \prod_{i=1}^n p(x_i)^{y_i}(1-p(x_i))^{1-y_i},\]
where \(\boldsymbol{\beta} = (\beta_0, \beta_1)^T\).

To find estimates of the parameters we need to maximize this function.
We start by taking the logarithm of the likelihood function:

\begin{align*} l(\boldsymbol{\beta}) =& \sum_{i=1}^n \Big ( y_i \log p(x_i) + (1-y_i) \log(1 - p(x_i) )\Big ) \\ &= \sum_{i=1}^n \Big ( y_i \log \Big (\frac{p(x_i)}{1-p(x_i)} \Big) + \log(1-p(x_i)) \Big ) \\
&= \sum_{i=1}^n \Big (y_i (\beta_0 + \beta_1 x_i) - \log(1 + e^{\beta_0 + \beta_1 x_i} \Big ).\end{align*}

To maximize the log-likelihood function we set the derivatives equal to
zero. This gives us nonlinear equations in \(\boldsymbol{\beta}\). To
solve these equations the Newton-Raphson algorithm is usually used.

\end{block}

\end{frame}

\begin{frame}

\begin{block}{Predictions}

Assuming we have fit a logistic regression model to our data set, and
have obtained parameter estimates \(\hat{\beta}_0\) and
\(\hat{\beta}_1\). We want to use this model to make a prediction when
given a new observation of \(X\). This is done by inserting the
parameter estimates and the \(X\) value into the logistic function:
\[\hat{p}(X) = \frac{e^{\hat{\beta}_0 + \hat{\beta}_1 X}}{1+e^{\hat{\beta}_0 + \hat{\beta}_1 X}}. = \frac{e^{\hat{\eta}(X)}}{1+e^{\hat{\eta}(X)}}\]
This \(\hat{p}(X)\) is the estimated probability that the new
observation \(X\) belongs to the class defined by \(Y=1\).

In the case of qualitative covariates, a dummy variable needs to be
introduced. This is done in a similar fashion as for linear regression.
Assume we have a qualitative covariate with two classes \{A, B\}. We
create a new variable:
\[X_i = \begin{cases}1, \quad \text{if obs. }i \text{ from class A} \\ 0, \quad   \text{if obs. } i \text{ from class B}.\end{cases}\]
And insert it into the link function:
\[\hat{\eta}_i = \hat\beta_0 + \hat\beta_1 X_i = \begin{cases} \hat\beta_0 + \hat\beta_1, \quad &\text{if obs. }i \text{ from class A}  \\ \hat\beta_0, \quad &\text{if obs. }i \text{ from class B} \end{cases}\]

\end{block}

\end{frame}

\begin{frame}

\begin{block}{Multiple logistic regression}

The simple logistic regression model easily generalizes to the case of
predicting a binary response with multiple covariates. Assume we have an
observation vector \(\mathbf{X}\), which consists of observations from
\(p\) covariates: \[\mathbf{X} = (X_1, X_2, \dots, X_p)^T.\] A multiple
logistic regression is defined through the following link function
\[\eta(\mathbf{X}) = \beta_0 + \beta_1 X_1 + \beta_2 X_2, \dots, \beta_p X_p.\]
We insert this link function into the logistic function, and obtain the
probabilities as before:
\[p (\mathbf{X}) = \frac{e^{\beta_0 + \beta_1 X + \dots + \beta_p X_p}}{1 +e^{\beta_0 + \beta_1 X + \dots + \beta_p X_p}}. \]
As in the case of simple logistic regression, we use maximum likelihood
to obtain parameter estimates.

\end{block}

\end{frame}

\begin{frame}[fragile]

\begin{block}{Example: South African heart disease data set}

In this example we use the \texttt{SAhert} data set from the
\texttt{ElemStatLearn} package. This is retrospective sample of males in
a heart-disease high-risk region in South Africa. It consists of 462
observations on the 10 variables. All subjects are male in the age range
15-64. There are 160 cases (individuals who have suffered from a
conorary heart disease) and 302 controls (individuals who have not
suffered from a conorary heart disease).

The response value:

\begin{itemize}
\tightlist
\item
  \texttt{chd} : conorary heart disease \{yes, no\} coded by the numbers
  \{1, 0\}
\end{itemize}

Covariates:

\begin{itemize}
\tightlist
\item
  \texttt{sbp} : systolic blood pressure\\
\item
  \texttt{tobacco} : cumulative tobacco (kg)\\
\item
  \texttt{ldl} : low density lipoprotein cholesterol
\item
  \texttt{adiposity} : a numeric vector
\item
  \texttt{famhist} : family history of heart disease. Categorical
  variable with two levels: \{Absent, Present\}.
\item
  \texttt{typea} : type-A behavior
\item
  \texttt{obesity} : a numerical value
\item
  \texttt{alcohol} : current alcohol consumption
\item
  \texttt{age} : age at onset
\end{itemize}

\end{block}

\end{frame}

\begin{frame}[fragile]

The goal is to identify important risk factors. We start by loading the
data and looking at the data:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ElemStatLearn)}
\KeywordTok{library}\NormalTok{(knitr)}
\KeywordTok{library}\NormalTok{(kableExtra)}
\NormalTok{heartds =}\StringTok{ }\NormalTok{SAheart}
\NormalTok{heartds}\OperatorTok{$}\NormalTok{chd =}\StringTok{ }\KeywordTok{as.factor}\NormalTok{(heartds}\OperatorTok{$}\NormalTok{chd)}
\KeywordTok{kable}\NormalTok{(}\KeywordTok{head}\NormalTok{(heartds))}
\end{Highlighting}
\end{Shaded}

\begin{longtable}[]{@{}rrrrlrrrrl@{}}
\toprule
sbp & tobacco & ldl & adiposity & famhist & typea & obesity & alcohol &
age & chd\tabularnewline
\midrule
\endhead
160 & 12.00 & 5.73 & 23.11 & Present & 49 & 25.30 & 97.20 & 52 &
1\tabularnewline
144 & 0.01 & 4.41 & 28.61 & Absent & 55 & 28.87 & 2.06 & 63 &
1\tabularnewline
118 & 0.08 & 3.48 & 32.28 & Present & 52 & 29.14 & 3.81 & 46 &
0\tabularnewline
170 & 7.50 & 6.41 & 38.03 & Present & 51 & 31.99 & 24.26 & 58 &
1\tabularnewline
134 & 13.60 & 3.50 & 27.78 & Present & 60 & 25.99 & 57.34 & 49 &
1\tabularnewline
132 & 6.20 & 6.47 & 36.21 & Present & 62 & 30.77 & 14.14 & 45 &
0\tabularnewline
\bottomrule
\end{longtable}

\end{frame}

\begin{frame}[fragile]

In order to investigate the data further, we use the \texttt{ggpairs}
function from the \texttt{GGally} library, to make scatter plots of the
covariates. The coloring is done according to the response variable,
where green represents a case \(Y=1\) and orange represents a control
\(Y=0\).

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(ggplot2)}
\KeywordTok{library}\NormalTok{(GGally)}
\KeywordTok{ggpairs}\NormalTok{(heartds, ggplot2}\OperatorTok{::}\KeywordTok{aes}\NormalTok{(}\DataTypeTok{color=}\NormalTok{chd), }\DataTypeTok{upper=}\StringTok{"blank"}\NormalTok{,  }
        \DataTypeTok{lower =} \KeywordTok{list}\NormalTok{(}\DataTypeTok{continuous =} \KeywordTok{wrap}\NormalTok{(}\StringTok{"points"}\NormalTok{, }\DataTypeTok{alpha =} \FloatTok{0.3}\NormalTok{, }\DataTypeTok{size=}\FloatTok{0.2}\NormalTok{)))}
\end{Highlighting}
\end{Shaded}

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-14-1.pdf}

\end{frame}

\begin{frame}[fragile]

We now fit a multiple logistic regression model using the \texttt{glm}
function and the full data set. In order to fit a logistic model, the
\texttt{family} argument must be set equal to \texttt{="binomial"}. The
\texttt{summary} function prints out the estimates of the coefficients,
their standard errors and z-values. As for a linear regression model,
the significant coefficients are indicated by stars where the
significant codes are included in the \texttt{R} outprint.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{glm_heart =}\StringTok{ }\KeywordTok{glm}\NormalTok{(chd}\OperatorTok{~}\NormalTok{., }\DataTypeTok{data=}\NormalTok{heartds, }\DataTypeTok{family=}\StringTok{"binomial"}\NormalTok{)}
\KeywordTok{summary}\NormalTok{(glm_heart)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = chd ~ ., family = "binomial", data = heartds)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.7781  -0.8213  -0.4387   0.8889   2.5435  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -6.1507209  1.3082600  -4.701 2.58e-06 ***
## sbp             0.0065040  0.0057304   1.135 0.256374    
## tobacco         0.0793764  0.0266028   2.984 0.002847 ** 
## ldl             0.1739239  0.0596617   2.915 0.003555 ** 
## adiposity       0.0185866  0.0292894   0.635 0.525700    
## famhistPresent  0.9253704  0.2278940   4.061 4.90e-05 ***
## typea           0.0395950  0.0123202   3.214 0.001310 ** 
## obesity        -0.0629099  0.0442477  -1.422 0.155095    
## alcohol         0.0001217  0.0044832   0.027 0.978350    
## age             0.0452253  0.0121298   3.728 0.000193 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 596.11  on 461  degrees of freedom
## Residual deviance: 472.14  on 452  degrees of freedom
## AIC: 492.14
## 
## Number of Fisher Scoring iterations: 5
\end{verbatim}

\end{frame}

\begin{frame}[fragile]

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(knitr)}
\KeywordTok{library}\NormalTok{(kableExtra)}
\KeywordTok{kable}\NormalTok{(}\KeywordTok{coef}\NormalTok{(glm_heart), }\DataTypeTok{col.names =} \KeywordTok{paste}\NormalTok{(}\StringTok{"Estimated coeffs"}\NormalTok{), }
      \DataTypeTok{digits =} \DecValTok{3}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\begin{longtable}[]{@{}lr@{}}
\toprule
& Estimated coeffs\tabularnewline
\midrule
\endhead
(Intercept) & -6.151\tabularnewline
sbp & 0.007\tabularnewline
tobacco & 0.079\tabularnewline
ldl & 0.174\tabularnewline
adiposity & 0.019\tabularnewline
famhistPresent & 0.925\tabularnewline
typea & 0.040\tabularnewline
obesity & -0.063\tabularnewline
alcohol & 0.000\tabularnewline
age & 0.045\tabularnewline
\bottomrule
\end{longtable}

\end{frame}

\begin{frame}[fragile]

\begin{block}{Multinomial logistic regression}

The logistic regression model can be generalized for a response variable
with more than two classes. Assume we have a response variable with
\(K\) possible classes and \(p\) covariates. The probability that \(Y\)
belongs to class \(k\), given an observation vector
\(\mathbf{X} = (X_1, X_2, \dots, X_p)^T\) is modeled by:
\[\text{Pr}(Y = k | \mathbf{X}) = \frac{e^{\beta_{0k} + \beta_{1k} X_1 + \dots + \beta_{pk} X_p}}{\sum_{l=1}^K e^{\beta_{0l}+\beta_{1l}X_1 + \dots + \beta_{pl}X_p}}.\]

We have a separate link function for each of the \(K\) classes. The
multinomial logistic regression model is implemented in the
\texttt{glmnet} package in \texttt{R}.

\end{block}

\end{frame}

\begin{frame}

\begin{block}{The Akaike Information Criterion (AIC)}

The AIC score is given by: \[AIC = 2 \cdot p -2 \cdot \text{loglik},\]
where \(p\) is the number of model parameters. The loglik is the
maximized log-likelihood and is given by
\[\text{loglik} = \sum_{i}^n \log \text{Pr}_{\hat{\boldsymbol{\beta}}} (y_i),\]
where \(\hat{\boldsymbol{\beta}}\) is the maximum-likelihood estimate of
the parameter-vector
\(\boldsymbol{\beta} = (\beta_0, \beta_1, ..., \beta_{p-1})^T\). The
role of \(p\) is to penalize models with many parameters as a high
number of parameters may lead to overfitting.

The AIC value can be used to choose between candidate logistic
regression models, where the model with the lowest AIC value is the one
expected to give the best fit.

\end{block}

\end{frame}

\begin{frame}{Confusion - sensitivity, specificity}

In a two class problem - assume the classes are labelled ``-'' (non
disease) and ``+'' (disease).

\begin{longtable}[]{@{}llll@{}}
\toprule
& Predicted - & Predicted + & Total\tabularnewline
True - & True Negative TN & False Positive FP & N\tabularnewline
True + & False Negative FN & True Positive TP & P\tabularnewline
Total & N* & P* &\tabularnewline
\bottomrule
\end{longtable}

\end{frame}

\begin{frame}

For this situation the following quantities are defined:

\begin{longtable}[]{@{}lll@{}}
\toprule
Name & Definition & Synonyms\tabularnewline
False positive rate & FP/N & Type I error, 1-specificity\tabularnewline
True positive rate & TP/P & 1-Type II error, power, sensitivity,
recall\tabularnewline
Positive predictive value (PPV) & TP/P* & Precision, 1-false discovery
proportion\tabularnewline
Negative predictive value (NPV) & TN/N* &\tabularnewline
\bottomrule
\end{longtable}

These two tables are tables 4.6 and 4.7 in our ISL-textbook.

\end{frame}

\begin{frame}

\textbf{Sensitivity} is the proportion of correctly classified positive
observations:
\(\frac{\# \text{True Positive}}{\# \text{Condition Positive}}\).

\textbf{Specificity} is the proportion of correctly classified negative
observations:
\(\frac{\# \text{True Negative}}{\# \text{Condition Negative}}\).

\end{frame}

\begin{frame}

\begin{block}{ROC curves}

The receiver operating characteristics (ROC) curve gives a graphical
display of the sensitivity against specificity, as the threshold value
is moved over the range of all possible values. An ideal classifier will
give a ROC curve which hugs the top left corner, while a straight line
represents a classifier with a random guess of the outcome. We will see
an ROC curve, in the continuation of the South African heart data
example, below.

The \textbf{AUC} score is the area under the AUC curve. It ranges
between the values 0 and 1, where a higher value indicates a better
classifier. An AUC score equal to 1 would imply that all observations
are correctly classified. The AUC score is useful for comparing the
performance of different classifiers, as all possible threshold values
are taken into account.

\end{block}

\end{frame}

\begin{frame}[fragile]

\begin{block}{Example Continued: South African heart disease}

We want to evaluate our multiple logistic model for the \texttt{SAheart}
data set. In order to investigate the training error and the test error,
we divide the original data set, randomly, into two samples of equal
size.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{set.seed}\NormalTok{(}\DecValTok{20}\NormalTok{)}
\NormalTok{train_ID =}\StringTok{ }\KeywordTok{sample}\NormalTok{(}\DecValTok{1}\OperatorTok{:}\KeywordTok{nrow}\NormalTok{(heartds), }\KeywordTok{nrow}\NormalTok{(heartds)}\OperatorTok{/}\DecValTok{2}\NormalTok{)}
\NormalTok{train_SA =}\StringTok{ }\NormalTok{heartds[train_ID, ]}
\NormalTok{test_SA =}\StringTok{ }\NormalTok{heartds[}\OperatorTok{-}\NormalTok{train_ID, ]}
\end{Highlighting}
\end{Shaded}

\end{block}

\end{frame}

\begin{frame}[fragile]

We now fit a logistic regression model, using the training set only:

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{glm_SA =}\StringTok{ }\KeywordTok{glm}\NormalTok{(chd}\OperatorTok{~}\NormalTok{. , }\DataTypeTok{data=}\NormalTok{train_SA, }\DataTypeTok{family=}\StringTok{"binomial"}\NormalTok{)}
\KeywordTok{summary}\NormalTok{(glm_SA)}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## 
## Call:
## glm(formula = chd ~ ., family = "binomial", data = train_SA)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.9715  -0.7993  -0.4098   0.8780   2.2163  
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    -7.425033   1.919850  -3.868  0.00011 ***
## sbp             0.013101   0.008822   1.485  0.13755    
## tobacco         0.088854   0.037542   2.367  0.01794 *  
## ldl             0.160858   0.082623   1.947  0.05155 .  
## adiposity       0.010770   0.038713   0.278  0.78086    
## famhistPresent  1.039578   0.335824   3.096  0.00196 ** 
## typea           0.042366   0.018254   2.321  0.02029 *  
## obesity        -0.044412   0.058290  -0.762  0.44611    
## alcohol        -0.004820   0.006672  -0.722  0.47000    
## age             0.045777   0.016873   2.713  0.00667 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 301.69  on 230  degrees of freedom
## Residual deviance: 232.54  on 221  degrees of freedom
## AIC: 252.54
## 
## Number of Fisher Scoring iterations: 5
\end{verbatim}

\end{frame}

\begin{frame}[fragile]

By comparing this outprint with the corresponding outprint above, we see
that the estimated coefficients slightly differ. This is because a
different data set has been used to fit the model. We previously used
the full data set.

We want to estimate the probability of a \texttt{chd} event for the
observations in the test set. To do this we can insert the estimated
coefficient into the logistic equation, remembering that
\texttt{famhist} is a categorical covariate, modeled by a dummy
variable:
\[x_\text{famhist} = \begin{cases} 1, \quad \text{ if Present}, \\ 0, \quad \text{ if Absent}.\end{cases}\]
The estimated probability of \(Y=1\) if \texttt{famhist\ =\ "Present"},
given a vector \(\mathbf{X}\) of covariate observations is:
\[\hat{p}(\mathbf{X}) =\frac{e^{-7.43 + 0.01 x_{\text{sbp}} + 0.09 x_{\text{tobacco}} + 0.16 x_\text{ldh} + 0.01 x_\text{adiposity}+1.04 \cdot 1 + 0.04 x_\text{typea}-0.04 x_\text{obesity}-0.00 x_\text{alcohol} +0.05 x_{\text{age}}}}{1+e^{-7.43 + 0.01 x_{\text{sbp}} + 0.09 x_{\text{tobacco}} + 0.16 x_\text{ldh} + 0.01 x_\text{adiposity}+1.04 \cdot 1 + 0.04 x_\text{typea}-0.04 x_\text{obesity}-0.00 x_\text{alcohol} +0.05 x_{\text{age}}}} .\]
Whereas, if \texttt{famhist\ =\ "Absent"} the estimated probability is:
\[\hat{p}(\mathbf{X}) =\frac{e^{-7.43 + 0.01 x_{\text{sbp}} + 0.09 x_{\text{tobacco}} + 0.16 x_\text{ldh} + 0.01 x_\text{adiposity}+1.04 \cdot 0 + 0.04 x_\text{typea}-0.04 x_\text{obesity}-0.00 x_\text{alcohol} +0.05 x_{\text{age}}}}{1+e^{-7.43 + 0.01 x_{\text{sbp}} + 0.09 x_{\text{tobacco}} + 0.16 x_\text{ldh} + 0.01 x_\text{adiposity}+1.04 \cdot 0 + 0.04 x_\text{typea}-0.04 x_\text{obesity}-0.00 x_\text{alcohol} +0.05 x_{\text{age}}}} .\]

\end{frame}

\begin{frame}[fragile]

The \texttt{predict} function does these calculations for us. When
specifying \texttt{type="response"} the function returns the
probabilities for \(Y=1\).

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{probs_SA =}\StringTok{ }\KeywordTok{predict}\NormalTok{(glm_SA, }\DataTypeTok{newdata=}\NormalTok{test_SA, }\DataTypeTok{type=}\StringTok{"response"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

From these probabilities we can obtain classifications, by specifying a
threshold value. We have here chosen a threshold value of 0.5. By using
the \texttt{ifelse} function we specify that all probabilities larger
than 0.5 are to be classified as 1, while the remaining probabilities
are to be classified as 0.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{pred_SA =}\StringTok{ }\KeywordTok{ifelse}\NormalTok{(probs_SA }\OperatorTok{>}\StringTok{ }\FloatTok{0.5}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{0}\NormalTok{)}

\NormalTok{predictions_SA =}\StringTok{ }\KeywordTok{data.frame}\NormalTok{(probs_SA, pred_SA, test_SA[,}\DecValTok{10}\NormalTok{])}
\KeywordTok{colnames}\NormalTok{(predictions_SA) =}\StringTok{ }\KeywordTok{c}\NormalTok{(}\StringTok{"Estim. prob. of Y=1"}\NormalTok{,}\StringTok{"Predicted class"}\NormalTok{,}\StringTok{"True class"}\NormalTok{)}
\KeywordTok{kable}\NormalTok{(}\KeywordTok{head}\NormalTok{(predictions_SA))}
\end{Highlighting}
\end{Shaded}

\begin{longtable}[]{@{}lrrl@{}}
\toprule
& Estim. prob. of Y=1 & Predicted class & True class\tabularnewline
\midrule
\endhead
2 & 0.3547764 & 0 & 1\tabularnewline
4 & 0.7732669 & 1 & 1\tabularnewline
5 & 0.6889170 & 1 & 1\tabularnewline
6 & 0.6404794 & 1 & 0\tabularnewline
10 & 0.6507839 & 1 & 1\tabularnewline
11 & 0.7241305 & 1 & 1\tabularnewline
\bottomrule
\end{longtable}

\end{frame}

\begin{frame}[fragile]

We can now use the confusion matrix to count the number of
misclassifications. The below confusion matrix is calculated using the
test set and comparing the predicted classes with the true classes.

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{table}\NormalTok{(pred_SA, SAheart[}\OperatorTok{-}\NormalTok{train_ID,}\DecValTok{10}\NormalTok{])}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
##        
## pred_SA   0   1
##       0 130  37
##       1  24  40
\end{verbatim}

The logistic model has correctly classified 130+40 times, and
misclassified 24+37 times. The misclassification test error rate is
thus: \[\text{Test error} = \frac{24+37}{231} \approx 0.264\]

\end{frame}

\begin{frame}[fragile]

The training error can be calculated in a similar fashion, but now we
use the fitted model to make prediction for the training set.

\begin{Shaded}
\begin{Highlighting}[]
\NormalTok{SA_train_prob =}\StringTok{ }\NormalTok{glm_SA}\OperatorTok{$}\NormalTok{fitted.values}
\NormalTok{SA_train_pred =}\StringTok{ }\KeywordTok{ifelse}\NormalTok{(SA_train_prob}\OperatorTok{>}\FloatTok{0.5}\NormalTok{, }\DecValTok{1}\NormalTok{, }\DecValTok{0}\NormalTok{)}
\NormalTok{conf_train =}\StringTok{ }\KeywordTok{table}\NormalTok{(SA_train_pred, SAheart[train_ID, }\DecValTok{10}\NormalTok{])}
\NormalTok{misclas_train =}\StringTok{ }\NormalTok{(}\DecValTok{231}\OperatorTok{-}\KeywordTok{sum}\NormalTok{(}\KeywordTok{diag}\NormalTok{(conf_train)))}\OperatorTok{/}\DecValTok{231}
\NormalTok{misclas_train}
\end{Highlighting}
\end{Shaded}

\begin{verbatim}
## [1] 0.2510823
\end{verbatim}

The train misclassification error rate is \(\approx 25.1\%\). In order
to see how our model performs for different threshold values, we can
plot a ROC curve:

\begin{Shaded}
\begin{Highlighting}[]
\KeywordTok{library}\NormalTok{(pROC)}
\NormalTok{SA_roc =}\StringTok{ }\KeywordTok{roc}\NormalTok{(SAheart[}\OperatorTok{-}\NormalTok{train_ID, }\DecValTok{10}\NormalTok{], probs_SA, }\DataTypeTok{legacy.axes=}\OtherTok{TRUE}\NormalTok{)}
\KeywordTok{ggroc}\NormalTok{(SA_roc)}\OperatorTok{+}\KeywordTok{ggtitle}\NormalTok{(}\StringTok{"ROC curve"}\NormalTok{)}\OperatorTok{+}\StringTok{ }
\StringTok{  }\KeywordTok{annotate}\NormalTok{(}\StringTok{"text"}\NormalTok{, }\DataTypeTok{x =} \FloatTok{0.25}\NormalTok{, }\DataTypeTok{y =} \FloatTok{0.30}\NormalTok{, }\DataTypeTok{label =} \StringTok{"AUC = 0.7762"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\includegraphics{4Classif_files/figure-beamer/unnamed-chunk-25-1.pdf}

\end{frame}

\begin{frame}{Which classification method is the best?}

\begin{block}{Advantages of discriminant analysis}

\begin{itemize}
\tightlist
\item
  Discriminant analysis is more stable than logistic regression when the
  classes are well-separated.
\item
  Discriminant analysis is more stable than logistic regression if the
  number of observations \(n\) is small and the distribution of the
  predictors \(X\) is approximately normal.
\end{itemize}

Assume a binary classification problem with one covariate. Recall that
multiple logistic regression is defined using the logit link function:
\[\log \Big ( \frac{p(X)}{1-p(X)}\Big ) = \beta_0 + \beta_1 X.\] For LDA
we have that \(p_0(X)\) is the probability that the observation \(X\)
belongs to class 0, while \(p_1(X)=1-p_0(X)\) is the probability that it
belongs to class 1.

\begin{align*}\log \frac{\text{Pr}(Y=0 |  X = x)}{\text{Pr}(Y=1 | X=x)} &= \log \frac{\pi_0}{\pi_1} + \log \frac{f_0(x)}{f_1(x)} \\ &= \log \frac{\pi_0}{\pi_1} -\frac{1}{2 \sigma^2}(x-\mu_0)^2 + \frac{1}{2 \sigma^2}(x-\mu_1)^2 \\ &= \log \frac{\pi_0}{\pi_1} - \frac{1}{2 \sigma^2} (x^2 - 2 x \mu_0 + \mu_0^2 - x^2 + 2 x \mu_1 - \mu_1^2) \\ &= \log \frac{\pi_0}{\pi_1} -\frac{1}{2 \sigma^2}(\mu_0^2 - \mu_1^2)+\frac{1}{\sigma^2}(\mu_0 - \mu_1) x \\ &= \alpha_0 + \alpha_1 x \end{align*}

\end{block}

\end{frame}

\begin{frame}

The two methods can thus be written in the same form. The difference is
in \emph{how} the parameters are estimated.

The \emph{joint} density of \(X\) and \(Y\) can be written as:
\[\text{Pr}(X, Y) = \text{Pr}(X)\text{Pr}(Y | X).\] For both logistic
regression and LDA the conditional probability of \(Y\) given \(X\) has
a logit-liner form.

\begin{itemize}
\tightlist
\item
  Logistic regression fits the parameters by maximizing the
  \emph{conditional} likelihood Pr(\(Y|X\)).
\item
  Linear discriminant analysis fits the parameters by maximizing the
  \emph{full} likelihood based on the joint density Pr(\(X, Y\)).
\item
  The results are usually quite similar.
\end{itemize}

\end{frame}

\begin{frame}[fragile]{ Recommended Exercises }

\begin{block}{From our ISL-textbook:}

\begin{block}{Exercise 4.7.6}

Suppose we collect data for a group of students in a statistics class
with variables \(X_1\) = hours studied, \(X_2\) = undergrad GPA, and
\(Y\) = receive an A. We fit a logistic regression and produce estimated
coefficient,
\(\hat{\beta}_0 = -6, \hat{\beta}_1 = 0.05, \hat{\beta}_2 = 1\).

\begin{enumerate}
\def\labelenumi{\alph{enumi}.}
\tightlist
\item
  Estimate the probability that a student who studies for 40 h and has
  an undergrad GPA of 3.5 gets an A in the class.\\
\item
  How many hours would the student in part a) need to study to have a
  50\% chance of getting an A in the class?
\end{enumerate}

\end{block}

\begin{block}{Exercise 4.7.9}

This problem has to do with \emph{odds}.

\begin{enumerate}
\def\labelenumi{\alph{enumi}.}
\tightlist
\item
  On average, what fraction of people with an odds of 0.37 of defaulting
  on their credit card payment will in fact default?\\
\item
  Suppose that an individual has a 16\% change of defaulting on her
  credit card payment. What are the odds that she will default?
\end{enumerate}

\end{block}

\begin{block}{Exercise 4.7.10}

This question should be answered using the \texttt{Weekly} data set,
which is part of the \texttt{ISLR}packages. This data is similar in
nature to the \texttt{Smarket} data from this chapter's lab, except that
it contains 1,089 weekly returns for 21 years, from the beginning of
1990 to the end of 2010.

\begin{enumerate}
\def\labelenumi{\alph{enumi}.}
\tightlist
\item
  Produce some numerical and graphical summaries of the \texttt{Weekly}
  data. Do there appear to be any patterns?\\
\item
  Use the full data set to perform a logistic regression with
  \texttt{Direction} as the response and the five lag variables plus
  \texttt{Volume}as predictors. Use the \texttt{summary} function to
  print the results. Do any of the predictors appear to be statistically
  significant? If so, which ones?\\
\item
  Compute the confusion matrix and overall fraction of correct
  predictions. Explain what the confusion matrix is telling you about
  the types of mistakes made by logistic regression.\\
\item
  Now fit the logistic regression model using a training data period
  from 1990 to 2008, with \texttt{Lag2} as the only predictor. Compute
  the confusion matrix and the overall fraction of correct predictions
  for the held out data (that is, the data from 2009 and 2010).\\
\item
  Repeat d) using LDA.
\item
  Repeat d) using QDA.
\item
  Repeat d) using KNN with \(K=1\).
\item
  Which of these methods appear to provide the best results on this
  data?
\item
  Experiment with different combinations of predictors, including
  possible transformations and interaction, for each of the methods.
  Report the variables, method, and associated confusion matrix that
  appears to provide the best results on the held out data. Note that
  you should also experiment with values for \(K\) in the KNN
  classifier.
\end{enumerate}

\end{block}

\begin{block}{Exercise 4.7.11}

In this problem, you will develop a model to predict whether a given car
gets high or low gas mileage based on the \texttt{Auto} data set.

\begin{enumerate}
\def\labelenumi{\alph{enumi}.}
\tightlist
\item
  Create a binary variable, \texttt{mpg01}, that contains a 1 if
  \texttt{mpg} contains a value above its median, and a 0 if
  \texttt{mpg} contains a value below its median. You can compute the
  median using the \texttt{median()} function. Note that you may find it
  helpful to use the \texttt{data.frame()} function to create a single
  data set containing both \texttt{mpg01}and the other \texttt{Auto}
  variables.
\item
  Explore the data graphically in order to investigate the association
  between \texttt{mpg01} and the other features. Which of the other
  features seems most likely to be useful in predicting \texttt{mpg01}?
  Scatter plots and boxplots may be useful tools to answer this
  question. Describe your findings.
\item
  Split the data into a training set and a test set.
\item
  Perform LDA on the training data in order to predict \texttt{mpg01}
  using the variables that seemed most associated with \texttt{mpg01} in
  b. What is the test error of the model obtained?
\item
  Repeat d) using LDA.
\item
  Repeat d) using logistic regression.
\item
  Repeat d) using KNN with different values of \(K\). Which value of
  \(K\) seems to perform the best on this data?
\end{enumerate}

\end{block}

\end{block}

\begin{block}{Modified version of Trial exam 2017 at University of Oslo}

\begin{enumerate}
\def\labelenumi{\alph{enumi}.}
\item
  Install the \texttt{DAAG} package and load the \texttt{frogs} data
  set. This data set consists of 212 observations of the following
  variables:

  \begin{itemize}
  \tightlist
  \item
    \texttt{pres.abs}: a binary variable (0/1) indicating the
    presence/absence of frogs at a particular location.\\
  \item
    \texttt{northing} : reference point\\
  \item
    \texttt{easting} : reference point\\
  \item
    \texttt{altitude} : altitude in meters
  \item
    \texttt{distance} : distance to nearest extant population, in meters
  \item
    \texttt{NoOfPools} : number of potential breeding pools
  \item
    \texttt{NoOfSites} : Number of potential breeding sites within a
    radius of 2 km
  \item
    \texttt{avrain} : mean rainfall during Spring
  \item
    \texttt{meanmin} : mean minimum temperature during Spring
  \item
    \texttt{meanmax} : mean maximum temperature during Spring
  \end{itemize}

  \begin{longtable}[]{@{}rllllllllll@{}}
  \toprule
  pre & s.abs nor & thing eas & ting alt & itude dis & tance NoO &
  fPools NoO & fSites avr & ain mea & nmin mea & nmax\tabularnewline
  \midrule
  \endhead
  2 & 1 & 115 & 1047 & 1500 & 500 & 232 & 3 & 155.00 & 3.57 &
  14.00\tabularnewline
  3 & 1 & 110 & 1042 & 1520 & 250 & 66 & 5 & 157.67 & 3.47 &
  13.80\tabularnewline
  4 & 1 & 112 & 1040 & 1540 & 250 & 32 & 5 & 159.67 & 3.40 &
  13.60\tabularnewline
  5 & 1 & 109 & 1033 & 1590 & 250 & 9 & 5 & 165.00 & 3.20 &
  13.17\tabularnewline
  6 & 1 & 109 & 1032 & 1590 & 250 & 67 & 5 & 165.00 & 3.20 &
  13.17\tabularnewline
  7 & 1 & 106 & 1018 & 1600 & 500 & 12 & 4 & 167.33 & 3.13 &
  13.07\tabularnewline
  \bottomrule
  \end{longtable}
\item
  Fit a full logistic model to the \texttt{frogs} data set, where
  \texttt{pres.abs} is the response variable. Call this model
  \texttt{glm.full}.

  \begin{enumerate}
  \def\labelenumii{\roman{enumii}.}
  \tightlist
  \item
    What is the log-likelihood for the \texttt{glm.full} model?
  \item
    Compute the confusion matrix for \texttt{glm.full}.
  \item
    What type of error can you find from this confusion matrix?
  \item
    Plot an ROC curve for \texttt{glm.full}. What is the AUC score?
  \item
    Should you be satisfied with this model? Argument why.
  \end{enumerate}
\item
  Find an alternative model \texttt{glm.alt} using backward stepwise
  selection, until only significant variables are included. Use the AIC
  values in the model-selection procedure.

  \begin{enumerate}
  \def\labelenumii{\roman{enumii}.}
  \tightlist
  \item
    What is the log-likelihood the \texttt{glm.alt} model?
  \item
    Why is this value lower than for the \texttt{glm.full} model? Based
    on the AIC values, which of the models \texttt{glm.full} and
    \texttt{glm.alt} should you prefer?
  \item
    Compute the confusion matrix for \texttt{glm.alt}.
  \item
    Plot an ROC curve for \texttt{glm.alt}. What is the AUC score?
  \item
    Based on the AUC scores, which of the models \texttt{glm.full} and
    \texttt{glm.alt} should you prefer?
  \end{enumerate}
\item
  Use the \texttt{summary} function print a summary of the
  \texttt{glm.alt} model.

  \begin{enumerate}
  \def\labelenumii{\roman{enumii}.}
  \tightlist
  \item
    What do the p-values of the last column tell you?
  \item
    Give an explanation of why the p-values for the variables which are
    present in both models may differ across the models.
  \end{enumerate}
\end{enumerate}

\end{block}

\end{frame}

\begin{frame}{ Further reading }

\begin{itemize}
\tightlist
\item
  Need details on the simple linear regression: From TMA4240/TMA4245
  Statistics we have the thematic page for
  \href{https://wiki.math.ntnu.no/tma4245/tema/begreper/regression}{Simple
  linear regression (in Norwegian)}.
\item
  Need more advanced thory: Theoretical version (no simple linear
  regression) from TMA4315 Generalized linear models H2017:
  \href{https://www.math.ntnu.no/emner/TMA4315/2017h/2MLR.html}{TMA4315M2:
  Multiple linear regression}
\item
  Slightly different presentation (more focus on multivariate normal
  theory):
  \href{https://www.math.ntnu.no/emner/TMA4267/2017v/TMA4267V2017Part2.pdf}{Slides
  and written material from TMA4267}
\item
  And, same source, but now
  \href{http://www.math.ntnu.no/emner/TMA4267/2017v/TMA4267V2017Part3.pdf}{Slides
  and written material from TMA4267 Linear Statistical Models in 2017,
  Part 3: Hypothesis testing and ANOVA}
\item
  \href{https://www.youtube.com/playlist?list=PL5-da3qGB5IBSSCPANhTgrw82ws7w_or9}{Videoes
  on YouTube by the authors of ISL, Chapter 2}
\end{itemize}

\end{frame}

\begin{frame}[fragile]{ R packages to install before knitting this R
Markdown file}

\footnotesize

\begin{Shaded}
\begin{Highlighting}[]
\CommentTok{# packages to install before knitting this R Markdown file}
\CommentTok{# to knit the Rmd}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"knitr"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"rmarkdown"}\NormalTok{)}

\CommentTok{# nice tables in Rmd}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"kableExtra"}\NormalTok{)}

\CommentTok{# cool layout for the Rmd}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"prettydoc"}\NormalTok{) }\CommentTok{# alternative to github}

\CommentTok{#plotting}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"ggplot2"}\NormalTok{) }\CommentTok{# cool plotting}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"ggpubr"}\NormalTok{) }\CommentTok{# for many ggplots}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"GGally"}\NormalTok{) }\CommentTok{# for ggpairs}
\CommentTok{#datasets}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"ElemStatLearn"}\NormalTok{) }\CommentTok{# for ozone data set}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"gamlss.data"}\NormalTok{)}\CommentTok{#rent index data set here}
\CommentTok{#methods}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"nortest"}\NormalTok{)}\CommentTok{#test for normality - e.g. Anderson-Darling}

\KeywordTok{install.packages}\NormalTok{(}\StringTok{"car"}\NormalTok{) }\CommentTok{# vif}
\KeywordTok{library}\NormalTok{(Matrix)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"reshape"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"corrplot"}\NormalTok{)}
\KeywordTok{install.packages}\NormalTok{(}\StringTok{"tidyverse"}\NormalTok{)}
\end{Highlighting}
\end{Shaded}

\end{frame}

\end{document}
