Search code examples
rpcavegan

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?


Solution

  • 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