Function to return percent inertia explained for cca.objects

[MikkoVihtakari]

First, I want to thank for making such a great package for multivariate and ecological community analyses! I use vegan a lot in my work. One thing I constantly miss is the possibility to quickly glance the inertia explained by axes in CA and CCA plots. Therefore, I normally have to modify the plots. Something that takes time, if you want to mass produce biplots to look at your data from several different angles. I thought of sharing my solution here in the hope that something similar would appear in the package in the future :)

#' @title Get percent of total inertia explained by principal or constrained axes
#' @param mod cca.object
#' @param axes A numeric vector indicating the axes for which percent explained inertia should be returned for
#' @return Returns a named vector containing explained inertia as percentages. Returns constrained inertia fo CCA and unconstrained for CA

axis.expl <- function(mod, axes = 1:2) {
  
  if(is.null(mod$CCA)) {
    sapply(axes, function(i) {
    100*mod$CA$eig[i]/mod$tot.chi
    })
  } else {
    sapply(axes, function(i) {
    100*mod$CCA$eig[i]/mod$tot.chi
    })
  }
    
}

## Lets try it:

data("dune")

x <- cca(dune)

(labs <- axis.expl(x))

#      CA1      CA2 
# 25.33987 18.91696 

## Add to plot

plot(x, xlab = "", ylab = "")

title(xlab = paste0(names(labs[1]), " (", sprintf("%.1f", labs[1]), "%)"))
title(ylab = paste0(names(labs[2]), " (", sprintf("%.1f", labs[2]), "%)"))

If there already is such a feature, please correct me :)

[jarioksa]

There are two (mutually contradictory) points here:

1. I hate this. I detest the practice of showing percentages explained in plots. That is the reason why it is not shown in vegan plots. The whole concept of "percentage explained" is against my statistical thinking. That said, if somebody supplies a function to do this, and vegan developers want to merge that new feature, I won't resist -- but that is up to @gavinsimpson and @psolymos.

2. If such a feature is added, it should be easier to use. Adding titles with sprintf is something no normal vegan user should need to know about.

[MikkoVihtakari]

1. Well, you do not need to call it "percentage explained", but name it as some other measure and make clear to the user that it really is not the same than R2 for regressions, for instance. A fact is that this percentage does tell stories about data and that in my view is the point of statistics.

2. The example above was just an example. Wrapping it into the plot.cca should not be a huge challenge.

[jarioksa]

Function summary.eigenvals() does this. The problem (if it is a problem) is easy use in labelling axes.

I do think that in many cases the signal strength (eigenvalues of shown axes) is more informative measure than the proportion of noise not shown (= sum of eigenvalues of axes not used) that you have in R². Just yesterday I saw scary example of case where R² was used with detrended correspondence analysis (DCA) where the measure is absolutely meaningless (and the DCA axis was also labeled with p-value which has no justification at all) -- and this was in a recently published paper. I sill dislike the idea of R², but if a code for easy labelling of axes is suggested, it may be merged. Depends quite a lot on @gavinsimpson who recently has developed summary.eigenvals() capabilities that are central here.

[MikkoVihtakari]

OK. I see. Yes, it is always possible to use tools wrong. I guess the philosophical question is whether it is the toolmaker's or user's fault (I have a clear opinion on this, but anyhow). You seem to ask for a suggestion how to implement the axis titles.

Looking at code here: https://github.com/vegandevs/vegan/blob/master/R/plot.cca.R

One could modify line 72 to:

plot(g[[1]], xlim = xlim, ylim = ylim, type = "n", asp = 1,
         xlab = "", ylab = "", ...)

Then add an argument to plot.cca. For example name it "axis.eigenvals".

Then one could add somewhere under line 75 the following:

evs <- eigenvals(x)
if(is.null(axis.eigenvals)) {
    title(xlab = names(evs[choices[1]]), ylab = names(evs[choices[2]]))
  } else {
    if(axis.eigenvals == "abs") {
      title(xlab = paste0(names(evs[choices[1]]), " (", sprintf("%.2f", evs[choices[1]]), ")"), ylab = paste0(names(evs[choices[2]]), " (", sprintf("%.2f", evs[choices[2]]), ")"))
  } else {
    if(axis.eigenvals == "rel") {
      title(xlab = paste0(names(evs[choices[1]]), " (", sprintf("%.1f", 100*evs[choices[1]]/x$tot.chi), "%)"), ylab = paste0(names(evs[choices[2]]), " (", sprintf("%.1f", 100*evs[choices[2]]/x$tot.chi), "%)"))
  } else {
    stop("axis.eigenvals has to be either NULL, 'abs', or 'rel'")
  }}
}

You can test it with following parameters:

x <- cca(dune)
choices = c(1,2)
plot(x, type = "n", xlab = "", ylab = "")
axis.eigenvals = NULL # change this between NULL, "abs" or "rel" and run the code above. 

Just an idea. Not a big deal, though.