LDA: Biological Applications

This page provides worked examples of LDA applied to gene expression and microbiome data. The core method is described in the main LDA page.

A recurring theme in biological applications is that the number of variables (genes, OTUs) far exceeds the number of samples. LDA breaks down in this setting, so the standard approach is to first reduce dimensions with PCA and then run LDA on the PC scores.

Cancer Subtype Classification

Gene expression datasets typically have thousands of variables and tens to hundreds of samples. The PCA-then-LDA workflow handles this well and produces interpretable results.

library(MASS)
library(ggplot2)

# Simulate RNA-seq expression data: 60 samples, 100 genes, 3 subtypes
set.seed(42)
n_genes   <- 100
n_samples <- 60
subtypes  <- factor(rep(c("TypeA", "TypeB", "TypeC"), each = 20))

expr <- matrix(rnorm(n_genes * n_samples), nrow = n_samples)
expr[1:20,  1:15]  <- expr[1:20,  1:15]  + 2   # TypeA signature
expr[21:40, 16:30] <- expr[21:40, 16:30] + 2   # TypeB signature
expr[41:60, 31:45] <- expr[41:60, 31:45] + 2   # TypeC signature

# 1. Reduce dimensions with PCA first (p >> n)
pca <- prcomp(expr, scale. = TRUE)

# Retain enough PCs to explain 80% of variance
cumvar <- cumsum(pca$sdev^2 / sum(pca$sdev^2))
n_pcs  <- which(cumvar >= 0.80)[1]
pc_scores <- pca$x[, 1:n_pcs]

# 2. LDA on PC scores
lda_cancer <- lda(subtypes ~ pc_scores)

# 3. Cross-validate
lda_cv <- lda(subtypes ~ pc_scores, CV = TRUE)
cat("LOOCV accuracy:", mean(lda_cv$class == subtypes), "\n")
table(Predicted = lda_cv$class, Actual = subtypes)

# 4. Visualise
pred <- predict(lda_cancer)
lda_df <- data.frame(pred$x, Subtype = subtypes)

ggplot(lda_df, aes(LD1, LD2, colour = Subtype)) +
  geom_point(size = 3) +
  stat_ellipse() +
  labs(title = "LDA: Cancer Subtype Separation") +
  theme_minimal()

Identifying discriminant genes

To identify which genes drive the separation, map the LDA loadings back through the PCA rotation:

# LDA loadings are in PC space; back-project to gene space
lda_loadings_gene <- pca$rotation[, 1:n_pcs] %*% lda_cancer$scaling

# Top genes contributing to LD1
top_ld1 <- sort(abs(lda_loadings_gene[, 1]), decreasing = TRUE)
head(top_ld1, 10)

Microbiome Sample Classification

For microbiome data, apply CLR transformation before PCA to handle compositionality, then run LDA on PC scores.

library(phyloseq)
library(microbiome)
library(MASS)
library(ggplot2)

data(GlobalPatterns)
gp <- GlobalPatterns

# 1. Filter rare taxa
gp_filtered <- filter_taxa(gp,
  function(x) sum(x > 3) > (0.1 * length(x)), TRUE)

# 2. CLR transformation
gp_clr <- microbiome::transform(gp_filtered, "clr")
otu_clr <- t(as(otu_table(gp_clr), "matrix"))
sample_type <- get_variable(gp_clr, "SampleType")

# 3. PCA to reduce dimensions
pca <- prcomp(otu_clr, scale. = FALSE)  # CLR already standardised

cumvar <- cumsum(pca$sdev^2 / sum(pca$sdev^2))
n_pcs  <- which(cumvar >= 0.80)[1]
pc_scores <- pca$x[, 1:n_pcs]

# 4. LDA on PC scores
lda_micro <- lda(sample_type ~ pc_scores)

# 5. Cross-validate
lda_cv <- lda(sample_type ~ pc_scores, CV = TRUE)
cat("LOOCV accuracy:", mean(lda_cv$class == sample_type), "\n")
table(Predicted = lda_cv$class, Actual = sample_type)

# 6. Plot
pred <- predict(lda_micro)
lda_df <- data.frame(pred$x[, 1:2], SampleType = sample_type)
names(lda_df)[1:2] <- c("LD1", "LD2")

ggplot(lda_df, aes(LD1, LD2, colour = SampleType)) +
  geom_point(size = 3) +
  stat_ellipse() +
  labs(title = "LDA: Microbiome Sample Types") +
  theme_minimal()

Practical Considerations

Scaling. Always scale continuous predictors before LDA. CLR-transformed microbiome data and VST-transformed RNA-seq data do not require additional scaling; raw morphological measurements do.

Unbalanced groups. LDA defaults to prior probabilities proportional to group sizes. With very unequal groups, the larger group dominates predictions. Set equal priors explicitly if classification balance matters:

lda(groups ~ ., data = data, prior = rep(1/k, k))  # k = number of groups

QDA as an alternative. If Box's M test indicates unequal covariances between groups, fit QDA instead:

qda_model <- qda(Species ~ ., data = iris)

QDA allows each group to have its own covariance matrix but requires more data per group (roughly n > 10p per group).

Posterior probabilities for uncertain samples. Inspect posterior probabilities to identify samples the model is uncertain about. These are often biologically interesting: transitional states, mixed communities, or mislabelled samples:

pred <- predict(lda_model)
uncertain <- which(apply(pred$posterior, 1, max) < 0.80)
pred$posterior[uncertain, ]