LDA

• Go to the sta-363-s20 GitHub organization and search for appex-02-lda
• Clone this repository into RStudio Cloud

• ${\mu }_1=-1.5$
• ${\mu }_2=1.5$
• ${\pi }_1={\pi }_2=0.5$
• ${\sigma }^2=1$
• typically we don't know the true parameters, we just use our training data to estimate them

Estimating parameters

$${\hat{\pi }}_k=\frac{n_k}{n}$$ $${\hat{\mu }}_k=\frac{1}{n_k}\sum _{i:y_i=k}{x}_i$$

$$\begin{array}{cc}{\hat{\sigma }}^2& =\frac{1}{n-K}\sum _{k=1}^K\sum _{i:y_i=k}(x_i-\widehat{{\mu }_k}{)}^2\\ & =\sum _{k=1}^K\frac{n_k-1}{n-K}{\hat{\sigma }}_k^2\end{array}$$

$${\hat{\sigma }}_k^2=\frac{1}{n_k-1}\sum _{i:y_i=k}(x_i-{\hat{\mu }}_k{)}^2$$

Estimating parameters (in R!)

$${\hat{\pi }}_k=\frac{n_k}{n}$$

df %>%
  group_by(y) %>%
  summarise(n = n()) %>%
  mutate(pi = n / sum(n))

## # A tibble: 3 x 3
##       y     n    pi
##   <dbl> <int> <dbl>
## 1     1     5 0.333
## 2     2     5 0.333
## 3     3     5 0.333

Estimating parameters (in R!)

$${\hat{\pi }}_k=\frac{n_k}{n}$$

df %>%
  group_by(y) %>%
  summarise(n = n()) %>%
  mutate(pi = n / sum(n))

## # A tibble: 3 x 3
##       y     n    pi
##   <dbl> <int> <dbl>
## 1     1     5 0.333
## 2     2     5 0.333
## 3     3     5 0.333

Estimating parameters (in R!)

$${\hat{\pi }}_k=\frac{n_k}{n}$$

df %>%
  group_by(y) %>%
  summarise(n = n()) %>%
  mutate(pi = n / sum(n))

## # A tibble: 3 x 3
##       y     n    pi
##   <dbl> <int> <dbl>
## 1     1     5 0.333
## 2     2     5 0.333
## 3     3     5 0.333

Estimating parameters (in R!)

$${\hat{\pi }}_k=\frac{n_k}{n}$$

df %>%
  group_by(y) %>%
  summarise(n = n()) %>% 
  mutate(pi = n / sum(n))

## # A tibble: 3 x 3
##       y     n    pi
##   <dbl> <int> <dbl>
## 1     1     5 0.333
## 2     2     5 0.333
## 3     3     5 0.333

Estimating parameters (in R!)

$${\hat{\pi }}_k=\frac{n_k}{n}$$

df %>%
  group_by(y) %>%
  summarise(n = n()) %>% 
  mutate(pi = n / sum(n))

Estimating parameters (in R!)

$${\hat{\pi }}_k=\frac{n_k}{n}$$

df %>%
  group_by(y) %>%
  summarise(n = n()) %>% 
  mutate(pi = n / sum(n)) %>%
  pull(pi) -> pi

Estimating parameters (in R!)

$${\hat{\pi }}_k=\frac{n_k}{n}$$

pi

## [1] 0.3333333 0.3333333 0.3333333

Estimating parameters (in R!)

$${\hat{\mu }}_k=\frac{1}{n_k}\sum _{i:y_i=k}{x}_i$$

df %>%
  group_by(y) %>%
  summarise(mu = mean(x))

## # A tibble: 3 x 2
##       y    mu
##   <dbl> <dbl>
## 1     1 -1.46
## 2     2  1.5 
## 3     3  3.54

Estimating parameters (in R!)

$${\hat{\mu }}_k=\frac{1}{n_k}\sum _{i:y_i=k}{x}_i$$

df %>%
  group_by(y) %>%
  summarise(mu = mean(x)) %>%
  pull(mu) -> mu

Estimating parameters (in R!)

$$\begin{array}{c}{\hat{\sigma }}^2=\sum _{k=1}^K\frac{n_k-1}{n-K}{\hat{\sigma }}_k^2\end{array}$$

df %>%
  group_by(y) %>%
  summarise(var_k = var(x),
            n = n()) %>%
  mutate(v = ((n - 1) / (sum(n) - 3)) * var_k) %>%
  summarise(sigma_sq = sum(v))

## # A tibble: 1 x 1
##   sigma_sq
##      <dbl>
## 1     1.47

Estimating parameters (in R!)

$$\begin{array}{c}{\hat{\sigma }}^2=\sum _{k=1}^K\frac{n_k-1}{n-K}{\hat{\sigma }}_k^2\end{array}$$

df %>%
  group_by(y) %>%
  summarise(var_k = var(x),
            n = n()) %>%
  mutate(v = ((n - 1) / (sum(n) - 3)) * var_k) %>%
  summarise(sigma_sq = sum(v)) %>%
  pull(sigma_sq) -> sigma_sq

Estimating parameters (in R!)

$${\delta }_k\left(x\right)=x\frac{{\mu }_k}{{\sigma }^2}-\frac{{\mu }_k^2}{2{\sigma }^2}+\log \left({\pi }_k\right)$$

• Let's predict the class for $x=2$

x <- 2
x * (mu / sigma_sq) - mu^2 / (2 * sigma_sq) + log(pi)

## [1] -3.8155857  0.1795063 -0.5436021

Estimating parameters (in R!)

$${\delta }_k\left(x\right)=x\frac{{\mu }_k}{{\sigma }^2}-\frac{{\mu }_k^2}{2{\sigma }^2}+\log \left({\pi }_k\right)$$

• Let's predict the class for $x=6$

x <- 6
x * (mu / sigma_sq) - mu^2 / (2 * sigma_sq) + log(pi)

## [1] -7.796499  4.269486  9.108750

From the discriminant score to probabilities

We can turn ${\hat{\delta }}_k\left(x\right)$ into estimates for class probabilities

$$\hat{P}(Y=k|X=x)=\frac{e^{{\hat{\delta }}_k\left(x\right)}}{{\sum }_{l=1}^K{e}^{{\hat{\delta }}_l\left(x\right)}}$$

• Classifying the largest ${\hat{\delta }}_k\left(x\right)$ is the same as classifying to the class with the largest $\hat{P}(Y=k|X=x)$
• For $K=2$:
  • classify to 2 if $\hat{P}(Y=2|X=x)\ge 0.5$
  • classify to 1 otherwise

Estimating parameters (in R!)

$$\hat{P}(Y=k|X=x)=\frac{e^{{\hat{\delta }}_k\left(x\right)}}{{\sum }_{l=1}^K{e}^{{\hat{\delta }}_l\left(x\right)}}$$

• Let's get the posterior probability of each class for $x=6$

x <- 6
d <- x * (mu / sigma_sq) - mu^2 / (2 * sigma_sq) + log(pi)
exp(d) / sum(exp(d))

## [1] 4.515655e-08 7.850755e-03 9.921492e-01

Estimating parameters (in R!)

• There is a function to do this in R called lda() in the MASS package

library(MASS)
model <- lda(y ~ x, data = df)

Estimating parameters (in R!)

• There is a function to do this in R called lda() in the MASS package

library(MASS) 
model <- lda(y ~ x, data = df)
predict(model, newdata = data.frame(x = 6))

## $class
## [1] 3
## Levels: 1 2 3
## 
## $posterior
##              1           2         3
## 1 4.515655e-08 0.007850755 0.9921492
## 
## $x
##        LD1
## 1 3.968523

LDA

• Go to the sta-363-s20 GitHub organization and search for appex-02-lda
• Clone this repository into RStudio Cloud
• Complete the exercises
• Knit, Commit, Push

Linear discriminant analysis $p>1$

• When $p>1$ the density takes on the multivariate normal density

$$f\left(x\right)=\frac{1}{(2\pi {)}^{p/2}|\Sigma {|}^{1/2}}{e}^{-\frac{1}{2}(x-\mu {)}^T{\Sigma }^{-1}(x-\mu )}$$

Linear discriminant analysis $p>1$

• The discriminant function is now

$${\delta }_k\left(x\right)=x^T{\Sigma }^{-1}{\mu }_k-\frac{1}{2}{\mu }_k^T{\Sigma }^{-1}{\mu }_k+\log {\pi }_k$$

• This is still a linear function!

Example $p=2$, $K=3$

• Here ${\pi }_1={\pi }_2={\pi }_3=1/3$
• The dashed lines the Bayes decision boundaries
  • If they were known, they would yield the fewest misclassification errors, among all possible classifiers.

LDA on Credit Data

• $\frac{23+252}{10000}$ errors - $2.75\%$ misclassification

  Since this is training error what is a possible concern?

• This could be overfit

LDA on Credit Data

• $\frac{23+252}{10000}$ errors - $2.75\%$ misclassification
• This could be overfit
• Since we have a large n and small p ( $n=10,000$, $p=4$ ) we aren't too worried about overfitting

LDA on Credit Data

• $\frac{23+252}{10000}$ errors - $2.75\%$ misclassification

  What would the error rate be if we classified to the prior, No default?

• $333/10000$ - $3.33\%$

LDA on Credit Data

• $\frac{23+252}{10000}$ errors - $2.75\%$ misclassification
• Since we have a large n and small p ( $n=10,000$, $p=4$ ) we aren't too worried about overfitting
• Of the true No's, we make $23/9667=0.2\%$ errors; of the true Yes's, we make $252/333=75.7\%$ errors!

Types of errors

• False positive rate: The fraction of truly negative that are classified as positive
• False negative rate: The fraction of truly positive that are classified as negative

• 0.2%

• 75.7%

Types of errors

• False positive rate: The fraction of truly negative that are classified as positive
• False negative rate: The fraction of truly positive that are classified as negative
• The Credit Default table was created by predicting the Yes class if

$$\hat{P}\left(\text{Default}|\text{Balance, Student}\right)\ge 0.5$$ We can change the two error rates by changing the *threshold from 0.5 to some other number between 0 and 1

$$\hat{P}\left(\text{Default}|\text{Balance, Student}\right)\ge threshold$$

Varying the threshold

• To reduce the false negative rate we may want the threshold to be 0.1 or less

ROC

Let's see it in R

library(MASS)
model <- lda(default ~ balance + student + income, data = Default)

• Use the lda() function in R from the MASS package

Let's see it in R

library(MASS)
model <- lda(default ~ balance + student + income, data = Default)
predictions <- predict(model)

• Use the lda() function in R from the MASS package
• Get the predicted classes along with posterior probabilities using the predict() function

Let's see it in R

library(MASS)
model <- lda(default ~ balance + student + income, data = Default)
predictions <- predict(model)
Default %>%
  mutate(predicted_class = predictions$class)

• Use the lda() function in R from the MASS package
• Get the predicted classes along with posterior probabilities using the predict() function
• Add the predicted class using the mutate() function

Let's see it in R

library(MASS)
model <- lda(default ~ balance + student + income, data = Default)
predictions <- predict(model)
Default %>%
  mutate(predicted_class = predictions$class) %>%
  summarise(fpr = 
              sum(default == "No" & predicted_class == "Yes") /
              sum(default == "No"),
            fnr = 
              sum(default == "Yes" & predicted_class == "No") /
              sum(default == "Yes"))

##           fpr       fnr
## 1 0.002275784 0.7627628

• Use the summarise() function to add the false positive and false negative rates

Let's see it in R

library(MASS)
library(tidymodels)
model <- lda(default ~ balance + student + income, data = Default)
predictions <- predict(model)
Default %>%
  mutate(predicted_class = predictions$class) %>%
  conf_mat(default, predicted_class) %>%
  autoplot(type = "heatmap")

Let's see it in R

• conf_mat() expects your outcome to be a factor variable

library(MASS)
library(tidymodels)
model <- lda(default ~ balance + student + income, data = Default)
predictions <- predict(model)
Default %>%
  mutate(predicted_class = predictions$class,
         default = as.factor(default)) %>%
  conf_mat(default, predicted_class) %>% 
  autoplot(type = "heatmap")