PFTC3 Wayqecha (PERÚ)

Exploring how the data are organized

Trait data

We suggest starting exploring the data sets of traits, asking for the names of the columns, the levels of each factor, and dimensions.

colnames(traits)

Explore a bit the levels of each factor, and the variable names.

unique(traits$Elevation)
#>  [1] 3890.0 3714.8 3882.9 3893.1 3714.6 3714.7 3714.3 3444.1 3489.0 3446.3
#> [11] 3447.7 3490.0 3889.8 3078.8 3071.7 3073.1 3884.3     NA 3077.1 3687.0
#> [21] 3444.2 3487.7 3118.0 3125.0 3492.4 3123.3 3121.0 3487.4 3120.3 3683.5
#> [31] 3685.5 3686.8 3686.9 3626.7 3673.3 3665.4 3694.3 3687.4 3665.5 3671.2
#> [41] 3669.6 3447.3
unique(traits$Treatment)
#> [1] "B"  "C"  "BB" NA
unique(traits$Site)## just change here the name of the column you want to explore
#> [1] "QUE" "TRE" "ACJ" "WAY" "PIL"

Multivariate exploratory

PCA

This document explains how to perform a PCA using vegan::rda() and ggvegan::autoplot.rda()

First, subset the big “traits” tibble into a smaller one. Log-transform some of the traits, select a few traits that are quite complete and the columns Taxon and Site.

traits.sel <- traits %>% # trait selection
    mutate(Leaf_Thickness_Ave_mm.log = log(Leaf_Thickness_Ave_mm), 
           Leaf_Area_cm2.log = log(Leaf_Area_cm2)) %>%
    select(Taxon, Site, Plant_Height_cm,
           Leaf_Thickness_Ave_mm.log, Leaf_Area_cm2.log)

Select only those data from the site Wayquecha. For doing the PCA we only need the traits values. Eliminate the NAs too.

traits.WAY <- traits.sel %>%
  filter(Site == "WAY") %>%
  drop_na()

Now, store the result of the PCA into an object pca.WAY. Eliminate the columns Site and Taxon, we don’t need them for the PCA. Then, do the PCA using vegan::rda(). We use the argument scale = TRUE so that each variable is transformed to have a mean of zero and unit variance. This is necessary because the variables have different scales.

pca.WAY <- traits.WAY %>% # store result as "pca.WAY"
  select(-Site, -Taxon) %>%   # remove Site & Taxon columns
  rda(scale = TRUE) #do the PCA

And plot it using the autoplot function (this is coming from the ggvegan package, and is using autoplot.rda).

autoplot(pca.WAY) +
  theme(legend.position = "none")

Every leaf appears as a red dot. The traits are shown by arrows. Unfortunately the defaults have made the arrows too long, but we can control that with the const argument to autoplot.rda (You will need the latest version of ggvegan for this to work - install with remotes::install_github("richardjtelford/ggvegan")).

autoplot(pca.WAY, const = c(0.2, 1)) + 
  xlim(-.22, NA) + # make more space for the labels
  theme(legend.position = "none")

It would be nice to show PCA with different colours for the different species. We cannot do that directly with autoplot.rda(), but we can do it manually, using fortify() to extract the scores from the PCA. This gives us lots of control over the plot.

#extract scores from ordination
fpca.WAY <- fortify(pca.WAY, const = c(0.2, 1), axes = 1:2)

head(fpca.WAY)#inspect the object to see the differences between "sites" and "species" in vegan terminology
#>     Score                     Label         PC1          PC2
#> 1 species           Plant_Height_cm -0.09070637  0.060360514
#> 2 species Leaf_Thickness_Ave_mm.log -0.08924325 -0.067130140
#> 3 species         Leaf_Area_cm2.log -0.09659404  0.005340169
#> 4   sites                      sit1 -0.02044190 -0.007783237
#> 5   sites                      sit2  0.00772506 -0.048062284
#> 6   sites                      sit3  0.07404017 -0.039406178

#get the "site" scores (i.e. each leaf)
fpca.WAY_leaves <- fpca.WAY %>% 
  filter(Score == "sites") %>% 
  bind_cols(traits.WAY)

# get the trait/variable scores ("species" in vegan terminology)
fpca.WAY_traits <- fpca.WAY %>% 
  filter(Score == "species") 

#plot the leaf scores         
g <- ggplot(fpca.WAY_leaves, aes(x = PC1, y = PC2)) +
  geom_point(aes(colour = Taxon), show.legend = FALSE) +
  coord_equal() # force scale to be correct

g

Add trait arrows to plot:

g <- g +
  geom_segment(data = fpca.WAY_traits,
               aes(x = 0, y = 0, xend = PC1, yend = PC2),
               arrow = arrow(length = unit(0.2, "cm")),
               colour = "navy") +
  geom_text(data = fpca.WAY_traits, # add labels
            aes(x = PC1, y = PC2, label = Label), nudge_x = -0.005, size = 3, hjust = 1)

g <- g + 
  xlim(-.20, NA) + #space for labels
  theme(legend.position = "none")# no legend as too large
g

We can finish off the plot by adding the proportion of variance explained by each axis to the plot. We can calculate this with the function eigenvals. We can make the labels with glue::glue() or with paste().

explained <- eigenvals(pca.WAY)[1:2] / sum(eigenvals(pca.WAY)) * 100
explained <- round(explained, 1)

g + 
  labs(x = glue("PC1 [{explained[1]}%]"), y = glue("PC1 [{explained[2]}%]"))