LDA

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

• μ1=−1.5
• μ2=1.5
• π1=π2=0.5
• σ2=1
• typically we don't know the true parameters, we just use our training data to estimate them

Estimating parameters

ˆπk=nkn ˆμk=1nk∑i:yi=kxi

ˆσ2=1n−KK∑k=1∑i:yi=k(xi−^μk)2=K∑k=1nk−1n−Kˆσ2k

ˆσ2k=1nk−1∑i:yi=k(xi−ˆμk)2

Estimating parameters (in R!)

ˆπk=nkn

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!)

ˆπk=nkn

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!)

ˆπk=nkn

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!)

ˆπk=nkn

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!)

ˆπk=nkn

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

Estimating parameters (in R!)

ˆπk=nkn

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

Estimating parameters (in R!)

ˆπk=nkn

pi

## [1] 0.3333333 0.3333333 0.3333333

Estimating parameters (in R!)

ˆμk=1nk∑i:yi=kxi

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!)

ˆμk=1nk∑i:yi=kxi

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

Estimating parameters (in R!)

ˆσ2=K∑k=1nk−1n−Kˆσ2k

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!)

ˆσ2=K∑k=1nk−1n−Kˆσ2k

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!)

δk(x)=xμkσ2−μ2k2σ2+log(πk)

• 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!)

δk(x)=xμkσ2−μ2k2σ2+log(πk)

• 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 ˆδk(x) into estimates for class probabilities

ˆP(Y=k|X=x)=eˆδk(x)∑Kl=1eˆδl(x)

• Classifying the largest ˆδk(x) is the same as classifying to the class with the largest ˆP(Y=k|X=x)
• For K=2:
  • classify to 2 if ˆP(Y=2|X=x)≥0.5
  • classify to 1 otherwise

Estimating parameters (in R!)

ˆP(Y=k|X=x)=eˆδk(x)∑Kl=1eˆδl(x)

• 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(x)=1(2π)p/2|Σ|1/2e−12(x−μ)TΣ−1(x−μ)

Linear discriminant analysis p>1

• The discriminant function is now

δk(x)=xTΣ−1μk−12μTkΣ−1μk+logπk

• This is still a linear function!

Example p=2, K=3

• Here π1=π2=π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

• 23+25210000 errors - 2.75% misclassification

  Since this is training error what is a possible concern?

• This could be overfit

LDA on Credit Data

• 23+25210000 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

• 23+25210000 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

• 23+25210000 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

ˆP(Default|Balance, Student)≥0.5 We can change the two error rates by changing the *threshold from 0.5 to some other number between 0 and 1

ˆP(Default|Balance, Student)≥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")