LDA

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

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

Estimating parameters

ˆπk=nkn$${\hat{\pi }}_k=\frac{n_k}{n}$$ ˆμk=1nk∑i:yi=kxi$${\hat{\mu }}_k=\frac{1}{n_k}\sum _{i:y_i=k}{x}_i$$

ˆσ2=1n−KK∑k=1∑i:yi=k(xi−^μk)2=K∑k=1nk−1n−Kˆσ2k$$\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}$$

ˆσ2k=1nk−1∑i:yi=k(xi−ˆμk)2$${\hat{\sigma }}_k^2=\frac{1}{n_k-1}\sum _{i:y_i=k}(x_i-{\hat{\mu }}_k{)}^2$$

Estimating parameters (in R!)

ˆπk=nkn$${\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!)

ˆπk=nkn$${\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!)

ˆπk=nkn$${\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!)

ˆπk=nkn$${\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!)

ˆπk=nkn$${\hat{\pi }}_k=\frac{n_k}{n}$$

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

Estimating parameters (in R!)

ˆπk=nkn$${\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!)

ˆπk=nkn$${\hat{\pi }}_k=\frac{n_k}{n}$$

pi

## [1] 0.3333333 0.3333333 0.3333333

Estimating parameters (in R!)

ˆμk=1nk∑i:yi=kxi$${\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!)

ˆμk=1nk∑i:yi=kxi$${\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!)

ˆσ2=K∑k=1nk−1n−Kˆσ2k$$\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!)

ˆσ2=K∑k=1nk−1n−Kˆσ2k$$\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!)

δk(x)=xμkσ2−μ2k2σ2+log(πk)$${\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 <- 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)$${\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 <- 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)${\hat{\delta }}_k\left(x\right)$ into estimates for class probabilities

ˆP(Y=k|X=x)=eˆδk(x)∑Kl=1eˆδl(x)$$\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 ˆδk(x)${\hat{\delta }}_k\left(x\right)$ is the same as classifying to the class with the largest ˆP(Y=k|X=x)$\hat{P}(Y=k|X=x)$
• For K=2$K=2$:
  • classify to 2 if ˆP(Y=2|X=x)≥0.5$\hat{P}(Y=2|X=x)\ge 0.5$
  • classify to 1 otherwise

Estimating parameters (in R!)

ˆP(Y=k|X=x)=eˆδk(x)∑Kl=1eˆδl(x)$$\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$

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$p>1$

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

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

• The discriminant function is now

δk(x)=xTΣ−1μk−12μTkΣ−1μk+logπk$${\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$p=2$, K=3$K=3$

• Here π1=π2=π3=1/3${\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

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

  Since this is training error what is a possible concern?

• This could be overfit

LDA on Credit Data

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

LDA on Credit Data

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

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

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

LDA on Credit Data

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