Using a biplot of values obtained through principal component analysis, it is possible explore the explanatory variables that make up each principle component.  **Is this also possible with Linear Discriminant Analysis?**

Examples provided use the The data is "Edgar Anderson's Iris Data" (http://en.wikipedia.org/wiki/Iris_flower_data_set).

Example PCA biplot using the iris data set in R (code below):

![enter image description here][1]

This figure indicates that Petal length and Petal width are important in determining PC1 score and in discriminating between Species groups.  setosa has smaller petals and wider sepals.

Apparently, similar conclusions can be drawn from plotting linear discriminant predictor variable scaling scores on the linear discriminants (lda.fit$scaling: a matrix which transforms observations to discriminant functions, normalized so that within groups covariance matrix is spherical). I am not certain what the LDA plot presents, hence the question.  The axis are the two first linear discriminants (LD1 99% and LD2 1% of trace).  The coordinates of the red vectors are "Coefficients of linear discriminants" also described as "scaling".  Scaling is calculated as diag(1/f1, , p) and f1 is sqrt(diag(var(x - group.means[g, ]))).   Data can be projected data onto the linear discriminants (using predict) in conjuction with predictor variable scaling scores (code below, as demonstrated http://stackoverflow.com/a/17240647/742447. The data and the predictor variables are plotted together so that which species were defined by an increase in which predictor variables increased can be seen (as is done for usual PCA biplots).:

![Example LDA biplot using the iris data set in R][2]

From this plot, Sepal width, Petal Width and Petal Length all contribute to a similar level to LD1.  As expected, setosa appears to smaller petals and wider sepals.

There is no built-in way to plot such biplots from LDA in R and few discussions of this online, which makes me wary of this approach.

**Does this LDA plot (see code below) provide a statistically valid interpretation of predictor variable scaling scores ?**


Code for PCA:


    require(grid)
    
      iris.pca <- prcomp(iris[,-5])
      PC <- iris.pca
      x="PC1"
      y="PC2"
      PCdata <- data.frame(obsnames=iris[,5], PC$x)
      
      datapc <- data.frame(varnames=rownames(PC$rotation), PC$rotation)
      mult <- min(
        (max(PCdata[,y]) - min(PCdata[,y])/(max(datapc[,y])-min(datapc[,y]))),
        (max(PCdata[,x]) - min(PCdata[,x])/(max(datapc[,x])-min(datapc[,x])))
      )
      datapc <- transform(datapc,
                          v1 = 1.6 * mult * (get(x)),
                          v2 = 1.6 * mult * (get(y))
      )
      
      datapc$length <- with(datapc, sqrt(v1^2+v2^2))
      datapc <- datapc[order(-datapc$length),]
      
      p <- qplot(data=data.frame(iris.pca$x),
                 main="PCA",
                 x=PC1,
                 y=PC2,
                 shape=iris$Species)
      #p <- p + stat_ellipse(aes(group=iris$Species))
      p <- p + geom_hline(aes(0), size=.2) + geom_vline(aes(0), size=.2)
      p <- p + geom_text(data=datapc, 
                         aes(x=v1, y=v2,
                             label=varnames,
                             shape=NULL,
                             linetype=NULL,
                             alpha=length), 
                         size = 3, vjust=0.5,
                         hjust=0, color="red")
      p <- p + geom_segment(data=datapc, 
                            aes(x=0, y=0, xend=v1,
                                yend=v2, shape=NULL, 
                                linetype=NULL,
                                alpha=length),
                            arrow=arrow(length=unit(0.2,"cm")),
                            alpha=0.5, color="red")
      p <- p + coord_flip()
      
      
      print(p)


Code for LDA



    #Perform LDA analysis
    iris.lda <- lda(as.factor(Species)~.,
                     data=iris)
    
    #Project data on linear discriminants
    iris.lda.values <- predict(iris.lda, iris[,-5])
    
    #Extract scaling for each predictor and
    data.lda <- data.frame(varnames=rownames(coef(iris.lda)), coef(iris.lda))
    data.lda$length <- with(data.lda, sqrt(LD1^2+LD2^2))
    scale.para <- 0.75
    
    #Plot the results
    p <- qplot(data=data.frame(iris.lda.values$x),
               main="LDA",
               x=LD1,
               y=LD2,
               shape=iris$Species)#+stat_ellipse()
    p <- p + geom_hline(aes(0), size=.2) + geom_vline(aes(0), size=.2)
    p <- p + theme(legend.position="none")
    p <- p + geom_text(data=data.lda,
                       aes(x=LD1*scale.para, y=LD2*scale.para,
                           label=varnames, 
                           shape=NULL, linetype=NULL,
                           alpha=length),
                       size = 3, vjust=0.5,
                       hjust=0, color="red")
    p <- p + geom_segment(data=data.lda,
                          aes(x=0, y=0,
                              xend=LD1*scale.para, yend=LD2*scale.para,
                              shape=NULL, linetype=NULL,
                              alpha=length),
                          arrow=arrow(length=unit(0.2,"cm")),
                          color="red")
    p <- p + coord_flip()
    
    print(p)


  [1]: https://i.sstatic.net/urswP.jpg
  [2]: https://i.sstatic.net/5xWcW.jpg