Discriminant Analysis for Group Separation in R

November 17, 2016
By

(This article was first published on R – Aaron Schlegel, and kindly contributed to R-bloggers)

The term ‘discriminant analysis’ is often used interchangeably to represent two different objectives. These objectives of discriminant analysis are:

  • Description of group separation. Linear combinations of variables, known as discriminant functions, of the dependent variables that maximize the separation between the groups are used to identify the relative contribution of the p variables that best predict group membership.

  • Prediction of observations to groups using either linear or quadratic discriminant functions, known as LDA and QDA, respectively.

This post will explore the first objective of discriminant analysis with two groups. Future posts will examine the classification and prediction objective of discriminant analysis.

Discriminant Analysis for Two Groups

Discriminant analysis assumes the two samples or populations being compared have the same covariance matrix \Sigma but distinct mean vectors \mu_1 and \mu_2 with p variables. The discriminant function that maximizes the separation of the groups is the linear combination of the p variables. The linear combination denoted z = a'y transforms the observation vectors to a scalar. The discriminant functions thus take the form:

 z_{1i} = a'y_{1i} = a_1 y_{1i1} + a_2 y_{1i2} + \cdots + a_p y_{1ip} \qquad i = 1, 2, \cdots, n_1
 z_{2i} = a'y_{2i} = a_2 y_{2i1} + a_2 y_{2i2} + \cdots + a_p y_{2ip} \qquad i = 1, 2, \cdots, n_2

To compute the discriminant function coefficients, first find the sample means \bar{z}_1 and \bar{z}_2. The mean can be found by averaging the n values or as a linear combination of the sample mean vector y_1, \bar{y}.

 \bar{z}_1=\frac{1}{n_1}\sum^{n_1}_{i=1}z_{1i}=a'\bar{y}_1
 \bar{z}_2 = \frac{1}{n_2} \sum^n_{i=1} z_{2i} = a'\bar{y}_2

Where,

\bar{y}_1 = \sum^{n_1}_{i=1} \frac{y_{1i}}{n_1}
\bar{y}_{1} = \sum^{n_2}_{i=1} \frac{y_{2i}}{n_2}

The goal is to then find a vector a that maximizes the standardized squared difference (\bar{z}_1 - \bar{z}_2)^2 / s^2_z. The sample variance s^2_z is the sample variance of z_1, z_2, \cdots, z_n or from the vector a and the sample covariance matrix of the mean vectors y_1, y_2, \cdots, y_n, denoted by S.

 s^2_z = \frac{\sum^n_{i=1}(z_i - \bar{z})^2}{n - 1} = a'Sa

Thus the standardized squared distance (\bar{z}_1 - \bar{z}_2)^2 / s^2_z can also be written as the following:

 \frac{(\bar{z}_1 - \bar{z}_2)^2}{s^2_z} = \frac{[a'(\bar{y}_1 - \bar{y}_2)]^2}{a'S_{p1}a}

Where S_{p1} is an unbiased estimator of the covariance matrix \Sigma. S_{p1} is defined as:

 S_{p1} = \frac{1}{n_1 + n_2 - 2}(W_1 + W_2)

Where W_1 and W_2 are defined as matrices of the sample sum of squares and cross products.

 W_1 = \sum^{n_1}_{i=1}(y_{1i} - \bar{y}_1)(y_{1i} - \bar{y}_1)' = (n_1 - 1)S_1
 W_2 = \sum^{n_2}_{i=1}(y_{2i} - \bar{y}_2)(y_{2i} - \bar{y}_2)' = (n_2 - 1)S_2

For S_{p1} to exist, n_1 + n_2 - 2 > p must be satisified.

The maximum of the above function is found when a is equivalent or a multiple of the following:

 a = S_{p1}^{-1}(\bar{y}_1 - \bar{y}_2)

Since a can be a multiple of the above, it is not unique; however, its direction is unique. By ‘direction’, it is implied the relative values of the vector a, a_1, a_2, \cdots, a_p are unique.

Discriminant Analysis in R

The data we are interested in is four measurements of two different species of flea beetles. All measurements are in micrometers (\mu m) except for the elytra length which is in units of .01 mm. The data were obtained from the companion FTP site of the book Methods of Multivariate Analysis by Alvin Rencher.

Read the data and give the columns names for reference.

beetles <- read.table('BEETLES.DAT', col.names = c('Measurement.Number', 'Species', 'transverse.groove.dist', 'elytra.length', 'second.antennal.joint.length', 'third.antennal.joint.length'))

The dplyr package will be used for simple data manipulation.

library(dplyr)

Separate the two groups into different data frames.

beetle1 <- filter(beetles, Species == 1)[,3:6]
beetle2 <- filter(beetles, Species == 2)[,3:6]

Store the sample size and means of the two groups for later.

n1 <- nrow(beetle1)
n2 <- nrow(beetle2)

beetle1.means <- apply(beetle1, 2, mean)
beetle2.means <- apply(beetle2, 2, mean)

First, S_{p1} must be calculated.

w1 <- (n1 - 1) * var(beetle1)
w2 <- (n2 - 1) * var(beetle2)

sp1 <- 1 / (n1 + n2 - 2) * (w1 + w2)

As mentioned above, the groups are maximally separated when a = S_{p1}^{-1}(\bar{y}_1 - \bar{y}_2).

a <- solve(sp1) %*% (beetle1.means - beetle2.means)
a
##                                    [,1]
## transverse.groove.dist        0.3452490
## elytra.length                -0.1303878
## second.antennal.joint.length -0.1064338
## third.antennal.joint.length  -0.1433533

The output of which gives us the linear discriminant function coefficients. However, as noted earlier, the data is not commensurate and therefore needs to be scaled to provide any meaningful interpretation. The linear discriminant analysis coefficients can be standardized by diag(S_{p1})^{1/2}a.

diag(sp1)^(1/2) * a
##                                   [,1]
## transverse.groove.dist        4.136640
## elytra.length                -2.500550
## second.antennal.joint.length -1.157705
## third.antennal.joint.length  -2.067833

Which gives us the following discriminant function:

 z = 4.137y_1 - 2.501y_2 - 1.158y_3 - 2.068y_4

The interpretation of the discriminant function can be made in several ways. The most simple is to rank the absolute value of the coefficients and determine contribution based on the order of the coefficients. Another method is to perform a partial F-test to find the significance of the variables.

Judging from our discriminant function, it appears the first measurement is the most significant while the second and third measurements have similar contribution to group separation.

The MASS package contains the function lda() for performing linear discriminant analysis.

library(MASS)

The lda() function takes a formula argument.

beet.lda <- lda(Species ~ .-Measurement.Number, data = beetles)
beet.lda$scaling
##                                      LD1
## transverse.groove.dist       -0.09327642
## elytra.length                 0.03522706
## second.antennal.joint.length  0.02875538
## third.antennal.joint.length   0.03872998

Note the discriminant function coefficients are different than what we computed earlier. This difference is due to another scaling method employed by the lda() function. Since any multiple of a can be taken as the maximum vector, either vector would suffice as the solution. We can see the coefficients are scaled differently than the a vector we found earlier. Despite this difference in scaling, output of the lda() function would still provide the same interpretation of the coefficients if they were ordered by their absolute values.

beet.lda$scaling / a
##                                     LD1
## transverse.groove.dist       -0.2701715
## elytra.length                -0.2701715
## second.antennal.joint.length -0.2701715
## third.antennal.joint.length  -0.2701715

The group separation can be plotted by using the plot() function from the MASS package. Note since we are only concerned with one discriminant function the plot will be a histogram rather than a scatterplot.

plot(beet.lda)

histogram of group separation by discriminant analysis

Summary

This post explored the objective of group separation using discriminant function analysis. By performing and interpreting a discriminant analysis function, one can get a better sense of what contributes the most distinction between the sample groups. As we will see in future posts, the discriminant function can also be used to classify and predict future observations.

References

Rencher, A. (n.d.). Methods of Multivariate Analysis (2nd ed.). Brigham Young University: John Wiley & Sons, Inc.

The post Discriminant Analysis for Group Separation in R appeared first on Aaron Schlegel.

To leave a comment for the author, please follow the link and comment on their blog: R – Aaron Schlegel.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.

Search R-bloggers


Sponsors

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)