Goal

Today we are going to practice linear discriminant analysis, or LDA.

We are going to work with our Federalist data from previous labs. As a reminder, in this data set we have three authors: Jay, Madison, and Hamilton. To load the data, put the following into a chunk in Markdown and press play.

train <- read.csv("https://www.dropbox.com/scl/fi/9q980le0xizsvsmzeknmt/train_Federalist.csv?rlkey=vem4v5yft4d5ufbkvbwd6kn3g&st=bycfifdu&dl=1")
train <- train[,1:80]

Today, the data set has 85 rows and 80 columns. There is a response variable, author, and 79 features from LIWC. Our goal is to describe how these 79 features can be used to tell the papers of the different authors apart.

Starting Small

Before we use all 79 features, let’s start with just a few. This will help us dig a little deeper into what LDA is doing. We are going to start with just 5 features.

To create a data set called trainSmall with just 5 features, run the following code:

trainSmall <- train[,c(1,2,8,46,47,68)]

The goal of LDA is to take these 5 features and use them to create 2 new features called discriminants. Each discriminant will be a weighted mixture of the 5 original features. For example, linear discriminant 1 (LD1) will be:

\[LD1_i = \beta_1 x_{1i} + \beta_2 x_{2i} + \beta_3 x_{3i} + \beta_4 x_{4i} + \beta_5 x_{5i}\]

Here, \(x_{1i}\) is the value of the first feature for row \(i\), and so on.

Using the same set up as LD1, linear discriminant 2 (LD2) will be:

\[LD2_i = \gamma_1 x_{1i} + \gamma_2 x_{2i} + \gamma_3 x_{3i} + \gamma_4 x_{4i} + \gamma_5 x_{5i}\]

As we can see from Questions 2 and 3, the value of LD1 and LD2 can be different even with the same values of the features! It all depends on the values of the coefficients that are used.

This means that in order to use LDA, we have to estimate the coefficients, meaning the \(\beta\) and \(\gamma\) terms. If this sounds a lot like what we do in regression models, it should!! The difference in LDA versus other regression models is how we choose the coefficients.

In LDA, we choose the coefficients to:

• 1. Make the LD1 and LD2 scores of the papers by the same author be close to one another
• 2. Make the average LD1 and LD2 score of the papers of different authors be far from one another.

In other words, when we convert from the features to LD1 and LD2, we want to do this in such a way that papers by the same authors are clustered together, and papers by different authors are kept far apart.

The math is to find the coefficients with these properties is beyond the scope of this course, but luckily R can help us do it!

# Load the libraries we need 
require("tidyverse")
require("caret")
require("MASS")
require("ggthemes")
require("rsample")
library(dplyr)

# We skip column 1 because it is Y, which is categorical!
# This means the features are in columns 2 -6, so we use 2:6 
for(i in 2:6){
  trainSmall[,i] <- scale(trainSmall[,i])
  trainSmall[,i] <- as.numeric(trainSmall[,i])
}

set.seed(279)
lda_mod <- lda(author ~ ., data = trainSmall)

lda_mod$scaling

trainSmall|>
  cbind(predict(lda_mod)$x) |>
  ggplot(aes(x = LD1, y = LD2, color = author, shape = author)) +
  geom_point(size = 1, alpha = 0.8) +
  labs(x = "LD1 (NUMBER%)", y = "LD2 (NUMER%)", color = "author") +
  scale_color_colorblind() +
  guides(color = guide_legend(override.aes = list(alpha = 1, size = 4))) +
  theme_bw() +
  theme(legend.title = element_text(face = "bold"),
        axis.title = element_text(face = "bold"),
        axis.text = element_text(color = "black"))

tibble(discriminant = c("LD1", "LD2"),
                               svd = lda_mod$svd) |>
  mutate(sigma_squared = svd^2,
         pct_rounded = round(sigma_squared/sum(sigma_squared)*100, 1)) |>
         dplyr::select(discriminant, pct_rounded)

LDA with all the features

To use all the features in LDA, we just switch from trainSmall to train.

# We skip column 1 because it is Y, which is categorical!
# This means the features are in columns 2 -80, so we use 2:80
for(i in 2:80){
  train[,i] <- scale(train[,i])
  train[,i] <- as.numeric(train[,i])
}

set.seed(279)
lda_mod <- lda(author ~ ., data = train)

At this point, we are ready to start interpreting. Remember that in LDA, we generally look at the top 10 largest (in absolute value) coefficients in LDA1 and LD2. To see those, we use the code below, but we will need to adapt it!!

features <- list(colnames(train)[-1])

top_features <- data.frame(lda_mod$scaling) |>
  rownames_to_column("features") |>
  pivot_longer(starts_with("LD"),
               names_to = "discriminant",
               values_to = "coefficient") |>
  group_by(discriminant) |>
  arrange(desc(abs(coefficient))) |>
  mutate(rank = row_number()) |>
  pivot_wider(id_cols = features,
              names_from = discriminant,
              values_from = c(coefficient, rank)) |>
  mutate(association = case_when(
    # ADAPT THIS
    coefficient_LD1 & coefficient_LD2  ~ "Hamilton",
    # ADAPT THIS
    coefficient_LD1 & coefficient_LD2 ~ "Madison",
    # ADAPT THIS
    coefficient_LD1 & coefficient_LD2 ~ "Jay",
    TRUE ~ NA
  )) |>
  filter(rank_LD2 <= 10 | rank_LD1 <= 10) |>
  arrange(rank_LD1)

knitr::kable(top_features)

At this point, we are ready to interpret! It takes more work to get to this point then it does in linear regression, but we can work with a lot more features at once.

Prediction

The last thing we can do with LDA is prediction. This is especially important in LDA because there is really no other way to see how the model is performing. We don’t have an R-squared or anything like that, so getting something like a confusion matrix can help us see how well LDA is reflecting the patterns in the data!

We do this for all the papers in the training data, and we can create the following confusion matrix using this code (no adaptation needed).

preds <- factor(predict(lda_mod, newdata = train)$class)
actual <- factor(train$author) 

overall_accuracy <- tibble(preds, actual) |>
  mutate(match = if_else(preds == actual, TRUE, FALSE)) |>
  summarize(overall_accuracy = mean(match) * 100)

conf_matrix <- confusionMatrix(data = preds, reference = actual)

as_tibble(conf_matrix$table) |>
  group_by(Reference) |>
  mutate(freq = n/sum(n)) |>
  ungroup() |>
  ggplot(aes(y = Prediction, x = Reference, fill = freq)) +
  geom_tile(color = "white") +
  geom_text(aes(label = paste0(round(freq * 100, 2),"%\n(",n,")")),
            color = "white", size = 3) +
  scale_fill_gradient(low  = "#9ecae1",high = "#08306b",
                      labels = scales::percent, limits = c(0, 1)) +
  labs(fill = "Percent\n", x = "True Author",
       y = "Predicted Author") +
  coord_equal() +
  theme_linedraw() +
  theme(panel.grid.major = element_blank(),
        axis.ticks = element_line(color = "black"),
        axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text = element_text(color = "black"),
        axis.title = element_text(color = "black", face = "bold"),
        legend.title = element_text(color = "black", face = "bold"),
        legend.text = element_text(color = "black"))