CODEWITHSUNDEEP
rpcavegan
Mikko
Mikko

Reputation: 7755

How to calculate species contribution percentages for vegan rda/cca objects?

I am trying to reproduce the column ("variable" in FactoMineR::PCA, "species" in vegan::rda) contribution percentages to axes from the FactoMineR package in vegan. The contributions are coded in FactoMiner::PCA objects:

library(FactoMineR)
library(vegan)

data(dune)

fm <- FactoMineR::PCA(dune, scale.unit = FALSE, graph = FALSE)

head(round(sort(fm$var$contrib[,1], decreasing = TRUE), 3))
# Lolipere Agrostol Eleopalu Planlanc  Poaprat  Poatriv 
#  17.990   16.020   13.866    7.088    6.861    4.850

By looking at the code for FactoMiner::PCA, I found out that the contributions are calculated as squared axis coordinates divided by the axis eigenvalue and multiplied by 100%:

head(round(sort(100*fm$var$coord[,1]^2/fm$eig[1], decreasing = TRUE), 3))
# Lolipere Agrostol Eleopalu Planlanc  Poaprat  Poatriv 
#  17.990   16.020   13.866    7.088    6.861    4.850

I fail to replicate the calculations above using a vegan::rda object: 

vg <- rda(dune)

head(round(sort(100*scores(vg, choices = 1, display = "sp", 
scaling = 0)[,1]^2/vg$CA$eig[1], decreasing = TRUE), 3))
# Lolipere Agrostol Eleopalu Planlanc  Poaprat  Poatriv 
#   0.726    0.646    0.559    0.286    0.277    0.196

I am obviously doing something wrong, and the discrepancy is probably due to the difference how these two packages calculate the coordinates for columns as eigenvalues for axes are quite similar (identical for my actual dataset), but the coordinates are not:

# vegan eigenvalue for axis 1
vg$CA$eig[1]
#     PC1 
# 24.79532 

# FactoMineR eigenvalue for axis 1
fm$eig[1]
# [1] 23.55555

# vegan column coordinates for axis 1
head(round(scores(vg, choices = 1, display = "sp", scaling = 0)[,1], 3))
# Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere 
#  -0.176    0.400    0.007    0.155   -0.163   -0.097 

#FactoMineR, column coordinates for axis 1
head(round(fm$var$coord[,1], 3))
# Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere 
#   0.854   -1.943   -0.033   -0.751    0.791    0.472 

# Sum of column coordinates for vegan axis 1 to illustrate the difference
sum(scores(vg, choices = 1, display = "sp", scaling = 0)[,1])
# [1] -0.796912

# Sum of column coordinates for FactoMineR axis 1 to illustrate the difference
sum(fm$var$coord[,1])
# [1] 3.867738

How do I calculate column/species percentage contributions to ordination axes using a vegan rda object?

Upvotes: 4

Views: 1805

Answers (2)

Mikko
Mikko

Reputation: 7755

vegan's unscaled (scaling = 0) column coordinates have equal sums of squares (i.e. 1 for each axis). You can get "column contributions" by simply squaring the unscaled coordinates:

head(sort(round(100*scores(vg, display = "sp", scaling = 0)[,1]^2, 3), decreasing = TRUE))
# Lolipere Agrostol Eleopalu Planlanc  Poaprat  Poatriv 
#  17.990   16.020   13.866    7.088    6.861    4.850

Sum of these "contributions" equals 100% as said above:

sum(100*scores(vg, display = "sp", scaling = 0)[,1]^2)
# [1] 100

Upvotes: 3

Jari Oksanen
Jari Oksanen

Reputation: 3682

The unscaled scores in vegan are unscaled in the (normal) sense that their sum of squares is 1 -- independent of eigenvalues:

> colSums(scores(vg, choices=1:4,dis="sp", scaling=0)^2)
  PC1 PC2 PC3 PC4 
    1   1   1   1

I think this is documented. If you want to call these squared terms as contributions, that's OK with me. The same holds for cca, but there you need to study weighted sums of squares. Moreover, the unscaled scores for sites (dis = "si") will have the same unit sum of squares: that is the idea of being unscaled. If you scale either species or sites, then the same relationships no longer hold for the other set of scores. In general, being unscaled means that the scores are orthonormal so that their crossproduct is an identity matrix (diagonal or sum of squares 1 and non-diagonal elements 0). For scaled scores these sums of squares are proportional to eigenvalues (but it may be useful to read vegan vignette on design decisions for constant scaling of the scores).

vegan functions goodness and inertcomp may (or may not) give you the information you're looking for.

Upvotes: 6

Related Questions