As mentioned in the post on classification with linear discriminant analysis, LDA assumes the groups in question have equal covariance matrices \((\Sigma_1 = \Sigma_2 = \cdots = \Sigma_k)\). Therefore, often when the groups do not have equal covariance matrices, observations are frequently assigned to groups with large variances on the diagonal of its corresponding covariance matrix (Rencher, n.d., pp. 321–321). Quadratic discriminant analysis is a modification of LDA that does not assume equal covariance matrices amongst the groups. In quadratic discriminant analysis, the group’s respective covariance matrix \(S_i\) is employed in predicting the group membership of an observation, rather than the pooled covariance matrix \(S_{p1}\) in linear discriminant analysis. The classification function in QDA is, therefore:

As in LDA, the observation \(y\) is assigned to the group for which \(D_i^2(y)\) is smallest.

One caveat to quadratic discriminant analysis is each group’s sample size \(n_i\) must be greater than the number of dependent variables \(p\).

Quadratic Discriminant Analysis in R

The beetles data, obtained from the companion FTP site of the book Methods of Multivariate Analysis by Alvin Rencher, will be analyzed by quadratic discriminant analysis.

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

head(beetles)

##   Measurement.Number Species transverse.groove.dist elytra.length
## 1                  1       1                    189           245
## 2                  2       1                    192           260
## 3                  3       1                    217           276
## 4                  4       1                    221           299
## 5                  5       1                    171           239
## 6                  6       1                    192           262
##   second.antennal.joint.length third.antennal.joint.length
## 1                          137                         163
## 2                          132                         217
## 3                          141                         192
## 4                          142                         213
## 5                          128                         158
## 6                          147                         173

The following function implements quadratic discriminant analysis to predict the group membership of the beetle observations. The function computes the group means and covariance matrices and then calculates \(D_i^2(y)\), as shown above, and outputs the predicted group, confusion matrix, and error rate.

two.group.quadratic.classification <- function(data, grouping, newdata) {
  dat.split <- split(data, grouping)
  g1 <- as.data.frame(dat.split[1])
  g2 <- as.data.frame(dat.split[2])
  g1.means <- apply(g1, 2, mean)
  g2.means <- apply(g2, 2, mean)
  g1.covar <- cov(g1)
  g2.covar <- cov(g2)
  
  prediction <- apply(newdata, 1, function(y) {
    d2.y1 <- (y - g1.means) %*% solve(g1.covar) %*% (y - g1.means)
    d2.y2 <- (y - g2.means) %*% solve(g2.covar) %*% (y - g2.means)
    ifelse(d2.y1^2 > d2.y2^2, 2, 1)
  })
  
  class.table <- table(grouping, prediction, dnn = c('Actual Group','Predicted Group'))
  pred.errors <- sum(diag(t(apply(class.table, 2, rev)))) / dim(data)[1]
  results <- list('Prediction'=prediction, 'Table of Predictions'=class.table, 'Error Rate'=pred.errors)
  
  return(results)
}

Run the function with the observed data as the newdata argument.

beetle.quad <- two.group.quadratic.classification(beetles[,3:6], beetles[,2], beetles[,3:6])
beetle.quad
## $Prediction
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2
## [36] 2 2 2 2
## 
## $`Table of Predictions`
##             Predicted Group
## Actual Group  1  2
##            1 19  0
##            2  1 19
## 
## $`Error Rate`
## [1] 0.02564103

The Mass package also supplies the qda() function to perform quadratic discriminant analysis.

library(MASS)

Similar to the lda() function in the MASS package, the qda() function takes a formula argument.

beetle.qda <- qda(Species ~.-Measurement.Number, data = beetles)

qda.pred <- predict(beetle.qda)$class
qda.pred
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 1 2 2 2 2 2 2 2 2
## [36] 2 2 2 2
## Levels: 1 2


table(beetles$Species, qda.pred, dnn = c('Actual Group','Predicted Group'))
##             Predicted Group
## Actual Group  1  2
##            1 19  0
##            2  1 19

Quadratic discriminant analysis predicted the same group membership as LDA. As noted in the previous post on linear discriminant analysis, predictions with small sample sizes, as in this case, tend to be rather optimistic and it is therefore recommended to perform some form of cross-validation on the predictions to yield a more realistic model to employ in practice.

Cross-Validation of Quadratic Discriminant Analysis Classifications

As before, we will use leave-one-out cross-validation to find a more realistic and less optimistic model for classifying observations in practice. Leave-one-out cross-validation is performed by using all but one of the sample observation vectors to determine the classification function and then using that classification function to predict the omitted observation’s group membership. The procedure is repeated for each observation so that each is classified by a function of the other observations. The approach to building the quadratic discriminant functions remains the same as before, with the exception that all but \(N-1\) observations are used.

cv.prediction <- c()

for (i in 1:dim(beetles)[1]) {
  holdout <- beetles[-i,]
  y <- as.numeric(beetles[i,3:6])
  
  holdout1 <- holdout[holdout$Species == 1,][,3:6]
  holdout2 <- holdout[holdout$Species == 2,][,3:6]
  
  holdout1.means <- apply(holdout1, 2, mean)
  holdout2.means <- apply(holdout2, 2, mean)
  
  holdout1.covar <- cov(holdout1)
  holdout2.covar <- cov(holdout2)

  d2.y1 <- (y - holdout1.means) %*% solve(holdout1.covar) %*% (y - holdout1.means)
  d2.y2 <- (y - holdout2.means) %*% solve(holdout2.covar) %*% (y - holdout2.means)
  
  group <- ifelse(d2.y1^2 > d2.y2^2, 2, 1)
  cv.prediction <- append(cv.prediction, group)
}

table(beetles$Species, cv.prediction, dnn = c('Actual Group','Predicted Group'))
##             Predicted Group
## Actual Group  1  2
##            1 17  2
##            2  2 18

The cross-validated results show two misclassified observations for each group, giving an error rate of approximately 10.3%. This error rate is slightly higher than the 7.7% error rate we found with cross-validated linear discriminant classifications. One could potentially test both models on new observations to determine their predictive power; however, as quadratic discriminant analysis makes fewer assumptions regarding the data and involves more parameters, it is likely that model would be more realistic in classifying observations.

Cross-validation can also be done with the qda() function. If the CV argument is set to TRUE, there is no need to call the predict() function on the qda object.

beetle.qda.cv <- qda(Species ~.-Measurement.Number, CV = TRUE, data = beetles)
table(beetles$Species, beetle.qda.cv$class, dnn = c('Actual Group','Predicted Group'))
##             Predicted Group
## Actual Group  1  2
##            1 17  2
##            2  2 18

The cross-validated results from the qda() function agree with our results.

References

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

The post Quadratic Discriminant Analysis of Two Groups appeared first on Aaron Schlegel.