Three functions in this package have to do with displaying relationships between site ordinations (based on species presences or abundances, i.e. community composition) and environmental variables. vegan’s envfit function makes linear fits between the variables and the ordination axes. These can then be displayed as vectors originating at the center of the plot. Projections of sites on these vectors should order the sites with respect to the variable. But the fit may not be linear. For individual variables, linearity can be checked graphically by making contour (or surface) plots, and by making bubble plots. In the first case, contours of the value of the variable are fit to the ordination plot. In the second, the site symbols are sized in proportion to the value of the variable.
gg_envfit
Usage for this function is:
gg_envfit(ord, env, groups = NA, scaling = 1, choices = c(1, 2), perm = 999, alpha = 0.05,
angle = 20, len = 0.5, unit = "cm", arrow.col = "red", pt.size = 3, plot = TRUE)
And a simple example copied from the documentation is:
data(varespec)
data(varechem)
vare.dist <- vegdist(varespec, method = "bray")
set.seed(123)
vare.mds <- monoMDS(vare.dist)
set.seed(123)
gg_envfit(ord = vare.mds, env = varechem, perm = 9999, pt.size = 2, alpha = 0.2)
As with gg_ordiplot, the function silently returns a list of the data frames (df_ord and df_arrows) used for the plot and the plot itself. They can be captured, examined and used to make modified plots as explained in the other vignette (“Modifying ggordiplots Plots”) included in this package.
gg_ordisurf and gg_ordibubble
The print method for envfit gives a table including the r2 and p.value (Pr>r) for each variable.
set.seed(123)
envfit(vare.mds, env = varechem, permutations = 9999)
##
## ***VECTORS
##
## MDS1 MDS2 r2 Pr(>r)
## N -0.08438 0.99643 0.2244 0.0628 .
## P -0.79760 -0.60318 0.1893 0.1128
## K -0.91829 -0.39591 0.1549 0.1674
## Ca -0.87347 -0.48688 0.2836 0.0287 *
## Mg -0.95770 -0.28777 0.2134 0.0810 .
## S -0.46045 -0.88769 0.0502 0.5873
## Al 1.00000 0.00057 0.4819 0.0012 **
## Fe 0.99076 0.13563 0.4072 0.0033 **
## Mn -1.00000 0.00224 0.4183 0.0026 **
## Zn -0.99761 -0.06912 0.1391 0.2076
## Mo 0.77127 0.63651 0.0482 0.5980
## Baresoil -0.99749 0.07085 0.1640 0.1564
## Humdepth -0.91408 -0.40554 0.4758 0.0019 **
## pH 0.99413 -0.10823 0.2062 0.0926 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Permutation: free
## Number of permutations: 9999
I used a fairly large value of alpha in the gg_envfit plot above so that vectors for variables with weak correlations would still be included in the plot. And of course they all plotted as vectors, i.e. straight lines, because that is the way the model works. Let’s compare variables with strong and weak correlations using contour and bubble plots. For gg_ordisurf you may have to try different values of binwidth to get a plot that you like. The default value is calculated as the difference between the minimum and maximum values of the variable divided by 15. Larger values give fewer contours and vice versa.
gg_ordisurf(ord = vare.mds, env.var = varechem$Al, binwidth = 20, pt.size = 1, var.label = "Aluminum")
gg_ordibubble(ord = vare.mds, env.var = varechem$Al, var.label = "Aluminum")
The contours are close to parallel and both plots show that the Al concentration increases from left to right, in the same direction as the vector in the previous plot.
Here we can compare directions because we are using the same ordination data for all of the plots. That is, we are not re-running the ordination each time. If we were, we would have to be more careful in making comparisons with directions. The signs along the axes have no real meaning, and axes can flip unexpectedly. This is especially true of NMDS which is essentially a randomization process. This is why I set the randomization seed, because otherwise the axes could flip when I compiled the document and would not match the text I had written. But axes can flip with other methods, too, if different algorithms are used or samples are added or removed.
Now make the same types of plots for a variable with a weak correlation, % of ground that is bare soil.
gg_ordisurf(ord = vare.mds, env.var = varechem$Baresoil, binwidth = 5, var.label = "Bare Soil")
gg_ordibubble(ord = vare.mds, env.var = varechem$Baresoil, var.label = "Bare Soil")
This relationship is obviously not linear. The value peaks near the center of the first axis (i.e.is monotonic) and there is no clear relationship with the second axis.
I tried making plots for molybdenum (Mo), but the correlation was so bad that the function to fit contours failed. From the bubble plot we can see why.
gg_ordibubble(ord = vare.mds, env.var = varechem$Mo, var.label = "Molybdenum")
gg_ordisurf and gg_ordibubble also silently return lists of the data frames used and the plots made. You may examine, use, and modify these with the methods you have learned above.
