I am comparing several attributes between two types of things. My goals are (1) hypothesis testing and (2) to summarize the information from multiple attributes.

PCA is an option, but this method is agnostic to the groups which I would like to directly compare. Discriminant analysis also seems not appropriate as this is geared towards classification. I think that constrained ordination, e.g. RDA or CCA should be appropriate, and so I am specifying the matrix of attributes as the response variables and a dummy variable for type as the explanatory variable.

My question is: how can there be a positive correlation between the constrained axis and the first residual axis? The code I am using is this:

> library(vegan)

> sample data

   type   trait1   trait2    trait3
1     1 1.758808 2.435267 1.0193438
2     2 1.856269 2.531002 1.0926777
3     2 1.792604 2.470863 1.0203008
4     1 1.842596 2.419813 1.0078421
5     2 1.867928 2.540830 1.0499309
6     2 1.899811 2.530313 1.1007924
7     2 1.943430 2.432563 1.0344759
8     1 1.808498 2.465469 1.1102114
9     1 1.773386 2.349574 0.9659113
10    1 1.800390 2.444265 1.0380732
11    2 1.826511 2.533739 1.1153933
12    1 1.895856 2.539948 1.0701774
13    1 1.873971 2.543685 1.0622794
14    1 1.819450 2.508187 0.9820533
15    1 1.848569 2.512621 0.9966571
16    2 1.855092 2.490199 1.0482279
17    2 1.907723 2.500366 1.0839503

ord <- rda(sample.data[, -1] ~ sample.data[,1])

enter image description here

There is a clear positive association between the constrained axis (called RDA by the output) and the first residual axis (called PCA1). I thought that a feature of constrained ordination was that residual variance should be uncorrelated with the constrained variance. Is this somehow an artifact of having a single explanatory variable and multiple response variables?


2 Answers 2


My understanding of constrained ordination is that, rather than performing a singular value decomposition as is the case in PCA and using a sparse representation of the Sigma matrix, "constrained ordination" displays the residuals from a regression analysis adjusting for the other $p-2$ variables in the covariance matrix. These factors appear as a point cloud if they are in fact conditionally independent. If there are other unmeasured variables that you have not accounted for (mediators or moderators), you will see a trend in the point cloud like what you have shown above.

  • $\begingroup$ Is there any reason why this would affect the interpretation of % variance explained by the first axis, and the results of a permutation test to test the significance of the first axis? $\endgroup$
    – user243404
    Apr 4, 2019 at 3:01

This is an old and stale issue, but I happened to stumble on this first now.

There are two kind of constrained row scores in vegan::rda: the WA scores which are found from the response data, and the LC (linear combination scores) that are the linear combinations of your constraints. The LC scores are the PCA of (scaled) fitted values of your response variables, and WA scores are the values that you get when you use response variables to predict the values of constraints. You had only one constraint ("type" in your data), and the LC score is only that one constraint rescaled. You can plot the ordination with LC scores by requesting them, or with your example:

plot(ord, display="lc")

You can also see that the LC scores of RDA1 are uncorrelated with residual unconstrained axes PC1 to PC3:

 zapsmall(cor(scores(ord, display="lc", choices=1:4)))
     RDA1 PC1 PC2 PC3
RDA1    1   0   0   0
PC1     0   1   0   0
PC2     0   0   1   0
PC3     0   0   0   1

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.