1

I am looking for the exact formula for the coefficients of linear discriminants in lda() from the package MASS. Consider the following data set.

set.seed(1)
sample<-data.frame(X=c(rnorm(40,-1),rnorm(40,1)),Class=c(rep("First",40),rep("Second",40)))
library(MASS)
lda.fit<-lda(Class~X,data=sample)

lda.fit contains the following data.

Call:
lda(Class ~ X, data = sample)

Prior probabilities of groups:
 First Second 
   0.5    0.5 

Group means:
                X
First  -0.9079738
Second  1.1202668

Coefficients of linear discriminants:
       LD1
X 1.103193

What is the exact formula used to calculate the coefficients of linear discriminants? How is 1.103193 calculated?

Any help is much appreciated!

P.S. I saw many answers on this website about LDA, but I am not able to find the exact formula used in lda() from the MASS package.

Cm7F7Bb
  • 247
  • 2
  • 9
  • I had exactly the same problem and found a very satisfactory answer at https://stats.stackexchange.com/questions/308872/difference-in-scaling-of-linear-discriminant-analysis-coefficients-between-manua – mindconnect.cc Mar 04 '20 at 00:31

1 Answers1

1

Instead of the formula for computing the coefficient of LDA, I'll show you an R session.

> set.seed(1)
> sample<-data.frame(X=c(rnorm(40,-1),rnorm(40,1)),Class=c(rep("First",40),rep("Second",40)))
> library(MASS)
> lda.fit<-lda(Class~X,data=sample)
> 
> group1 = sample[sample$Class=="First",]
> group2 = sample[sample$Class=="Second",]
> var1 = var(group1[,c(1)])
> var2 = var(group2[,c(1)])
> 
> n1 = nrow(group1)
> n2 = nrow(group2)
> n = n1+n2
> K = 2
> 
> var <- 1/(n - K) * (var1*(n1-1) + var2*(n2-1))
> 
> 1/sqrt(var) # same as below
[1] 1.103193
> 
> lda.fit$scaling 
       LD1
X 1.103193

Actually, this code does not reveal the whole story because $p$, the number of predictors is 1. I made an example for the case $p=2$ below.

set.seed(1)
sample <- data.frame(X=c(rnorm(40,-1),rnorm(40,1)), Y=c(rnorm(40,-1),rnorm(40,1)), Class=c(rep("First",40),rep("Second",40)));
library(MASS)
lda.fit<-lda(Class~X+Y,data=sample); lda.fit

Call:
lda(Class ~ X + Y, data = sample)

Prior probabilities of groups:
 First Second 
   0.5    0.5 

Group means:
                X          Y
First  -0.9079738 -0.8831604
Second  1.1202668  0.6896082

Coefficients of linear discriminants:
        LD1
X 0.9892123
Y 0.8650886

group1 = sample[sample$Class=="First",]
group2 = sample[sample$Class=="Second",]
cov1 = cov(group1[,c(1,2)])
cov2 = cov(group2[,c(1,2)])
n1 = nrow(group1)
n2 = nrow(group2)
n = n1+n2
K = 2
Sigma <- 1/(n - K) * (cov1*(n1-1) + cov2*(n2-1))
SigmaInv = solve(Sigma);
mu2_mu1 = matrix(lda.fit$means[2,] - lda.fit$means[1,]);
coeff.vec = SigmaInv %*% mu2_mu1;
# coeff.vec and lda.fit$scaling are linearly dependent
v.scalar = sqrt(t(coeff.vec) %*% Sigma %*% coeff.vec);
coeff.vec/drop(v.scalar) 
       [,1]
X 0.9892123
Y 0.8650886