Regression scores in constrained ordination (CCA, RDA, db-RDA)?

[jarioksa]

Should we implement regression scores in constrained ordination? Currently we have biplot scores which are displayed as default for continuous variables and polynomial constraints of ordered factors. We do have coef function that returns the regression coefficients for original (centred and weighted) variables, but regression scores are similar coefficients for standardized variables. So the original scale (like the unit of measurement) does influence coef, but it does not influence regression scores.

For regression scores we should (or could: this can be leaved to scores()) add an argument to return coef of standardized variables. This is a trivial one-liner. To use these coefficients in ordination plots, we also need to add scaling of these scores similarly as biplot scores. Here we have a bit more work.

The interpretation of regression scores is up to user who requests them. Biplot scores are easier to interpret as they show for each variable independently the direction of its steepest change or the direction of the gradient (at least with scaling = "species"). Regression scores cannot be interpreted individually, but all arrows must be interpreted together: the LC score for sites is obtained by summing the effect of displayed arrows for the values of each environmental variable for each sampling unit. Moreover, in partial models the biplot arrows for Conditions (partialling variables) are zero and are not shown, but their regression coefficients are non-zero and should be shown (although Canoco does not show them). For Conditions, the arrows show what must be done to partial out these effects. As an example, here the graphs of LC scores (= linear combinations of environmental variables) and biplot or regression arrows for the model rda(dune ~ Condition(Management) + A1 + Manure, dune.env)

Here the function that extracts standardized regression scores from a vegan constrained ordination result, and additional scales them if scaling = 1 (== scaling = "sites") -- other scaling alternatives are not implemented. This shows how little we need to do to implement regression scores, and most things to be done are janitorial (how to handle these in plots):

Regr <-
    function(x, scaling = 2)
{
    X <- qr.X(x$CCA$QR)
    out <- sqrt(colSums(X^2)) * coef(x)
    if (abs(scaling) == 1) {
        eig <- x$CCA$eig
        if (inherits(x, "rda"))
            eig <- eig/x$tot.chi
        out <- sweep(out, 2, sqrt(eig), "*")
    }
    out
}

[jarioksa]

Implemented in merge 83389b4