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.08414 0.99645 0.2245 0.0627 .
## P -0.79777 -0.60296 0.1891 0.1130
## K -0.91838 -0.39571 0.1549 0.1676
## Ca -0.87356 -0.48671 0.2836 0.0287 *
## Mg -0.95775 -0.28760 0.2134 0.0810 .
## S -0.46058 -0.88762 0.0502 0.5873
## Al 1.00000 0.00047 0.4819 0.0012 **
## Fe 0.99078 0.13548 0.4073 0.0032 **
## Mn -1.00000 0.00256 0.4182 0.0026 **
## Zn -0.99765 -0.06848 0.1390 0.2079
## Mo 0.77138 0.63638 0.0482 0.5980
## Baresoil -0.99753 0.07023 0.1641 0.1562
## Humdepth -0.91408 -0.40553 0.4759 0.0019 **
## pH 0.99417 -0.10786 0.2063 0.0925 .
## ---
## 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")
## Warning: The dot-dot notation (`..level..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(level)` instead.
## ℹ The deprecated feature was likely used in the ggordiplots package.
## Please report the issue at <https://github.com/jfq3/ggordiplots/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
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.
