sites.long: Helper Functions to Prepare Plotting of Accumulation, Diversity Profile and Ordiplot Results via ggplot2

Description

These functions organize outputs from ordiplot, accumcomp and renyicomp so these can be plotted with ggplot.

Usage

sites.long(x, env.data = NULL)
species.long(x, spec.data = NULL)
centroids.long(y, grouping, FUN = mean, centroids.only = FALSE)
vectorfit.long(z)
ordisurfgrid.long(z)
ordiellipse.long(z, grouping.name = "Grouping")
axis.long(w, choices = c(1, 2), cmdscale.model=FALSE, CAPdiscrim.model=FALSE)
accumcomp.long(x, ci = 2, label.freq = 1)
renyicomp.long(x, label.freq = 1)
renyi.long(x, env.data=NULL, label.freq = 1)

Arguments

Result of ordiplot, accumcomp or renyicomp

env.data

Environmental descriptors for each site.

spec.data

Descriptors for each species.

Result of function sites.long.

grouping

Variable defining the centroids

FUN

A function to compute the summary statistics which can be applied to all data subsets, as in aggregate

centroids.only

Return the coordinates for the centroids

Result of envfit, ordisurf or ordiellipse

grouping.name

Name for the categorical variable, expected as the factor used in the earlier ordiellipse call.

Ordination object from which the ordiplot was obtained, expected to be fitted in vegan.

choices

Ordination axes selected, as in ordiplot

cmdscale.model

Use TRUE is the model was fitted via cmdscale

CAPdiscrim.model

Use TRUE is the model was fitted via CAPdiscrim

Multiplier for confidence intervals as in specaccum. In case 'NA' is provided, then the multiplier is calculated via qt.

label.freq

Frequency of labels along the x-axis (count between labels).

Value

These functions produce data.frames that can subsequentially be plotted via ggplot methods.

Details

Examples for ordiplot results are shown here.

See also: for examples for species accumulation curves and Renyi diversity profiles

References

Kindt, R. & Coe, R. (2005) Tree diversity analysis: A manual and software for common statistical methods for ecological and biodiversity studies.

https://www.worldagroforestry.org/output/tree-diversity-analysis

https://rpubs.com/Roeland-KINDT

Examples

Run this code

# NOT RUN {
# }
# NOT RUN {
# ggplot2 plotting method
library(ggplot2)
library(ggforce)
library(concaveman)
library(ggrepel)
library(ggsci)
library(dplyr)

library(vegan)
data(dune)
data(dune.env)

attach(dune)
attach(dune.env)

Ordination.model1 <- capscale(dune ~ Management, data=dune.env, 
  distance='kulczynski', sqrt.dist=F, add='cailliez')

plot1 <- ordiplot(Ordination.model1, choices=c(1,2), scaling='species')

# obtain 'long' data from the ordiplot object
sites1 <- sites.long(plot1, env.data=dune.env)
species1 <- species.long(plot1)
centroids1 <- centroids.long(sites1, Management, FUN=median)
centroids2 <- centroids.long(sites1, Management, FUN=median, centroids.only=TRUE)

# information on percentage variation from the fitted ordination
axislabs <- axis.long(Ordination.model1, choices=c(1 , 2))

# change the theme
# possibly need for extrafont::loadfonts(device="win") to have Arial
# as alternative, use library(ggThemeAssist)
BioR.theme <- theme(
        panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid = element_blank(),
        axis.line = element_line("gray25"),
        text = element_text(size = 12, family="Arial"),
        axis.text = element_text(size = 10, colour = "gray25"),
        axis.title = element_text(size = 14, colour = "gray25"),
        legend.title = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.key = element_blank())

# no species scores
# centroids calculated directly via the centroids.long function 
plotgg1 <- ggplot() + 
    geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
    geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +  
    xlab(axislabs[1, "label"]) +
    ylab(axislabs[2, "label"]) +  
    scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
    scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +    
    geom_mark_hull(data=sites1, aes(x=axis1, y=axis2, colour = Management), 
        concavity = 0.1, alpha=0.8, size=0.2, show.legend=FALSE) +
    geom_point(data=sites1, aes(x=axis1, y=axis2, colour=Management, shape=Management), 
        size=5) +
#    geom_segment(data=species1, aes(x=0, y=0, xend=axis1*2, yend=axis2*2), 
#        size=1.2, arrow=arrow()) +
#    geom_label_repel(data=species1, aes(x=axis1*2, y=axis2*2, label=labels)) +
    geom_point(data=centroids.long(sites1, grouping=Management, centroids.only=TRUE), 
        aes(x=axis1c, y=axis2c, colour=Centroid, shape=Centroid), size=10, show.legend=FALSE) +
    geom_segment(data=centroids.long(sites1, grouping=Management), 
        aes(x=axis1c, y=axis2c, xend=axis1, yend=axis2, colour=Management), 
        size=1, show.legend=FALSE) +
    BioR.theme +
    ggsci::scale_colour_npg() +
    coord_fixed(ratio=1)

plotgg1

# select species to plot based on goodness of fit
spec.envfit <- envfit(plot1, env=dune)
spec.data1 <- data.frame(r=spec.envfit$vectors$r, p=spec.envfit$vectors$pvals)
species2 <- species.long(plot1, spec.data=spec.data1)
species2 <- species2[species2$r > 0.6, ]

# after_scale introduced in ggplot2 3.3.0
plotgg2 <- ggplot() + 
    geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
    geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +  
    xlab(axislabs[1, "label"]) +
    ylab(axislabs[2, "label"]) +   
    scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
    scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +    
    geom_mark_ellipse(data=sites1, aes(x=axis1, y=axis2, 
        colour=Management, fill=after_scale(alpha(colour, 0.2))), 
        expand=0, size=0.2, show.legend=TRUE) +
    geom_point(data=sites1, aes(x=axis1, y=axis2, colour=Management, shape=Management), 
        size=5) +
    geom_segment(data=centroids.long(sites1, grouping=Management), 
        aes(x=axis1c, y=axis2c, xend=axis1, yend=axis2, colour=Management), 
        size=1, show.legend=FALSE) +
    geom_segment(data=species2, aes(x=0, y=0, xend=axis1*2, yend=axis2*2), 
        size=1.2, arrow=arrow()) +
    geom_label_repel(data=species2, aes(x=axis1*2, y=axis2*2, label=labels)) +
    BioR.theme +
    ggsci::scale_colour_npg() +
    coord_fixed(ratio=1)

plotgg2

# Add contour and vector for a continuous explanatory variable
Ordination.model2 <- capscale(dune ~ Management, data=dune.env, 
  distance='kulczynski', sqrt.dist=F, add='cailliez')

plot2 <- ordiplot(Ordination.model2, choices=c(1,2), scaling='species')

sites2 <- sites.long(plot2, env.data=dune.env)
axislabs <- axis.long(Ordination.model2, choices=c(1 , 2))

dune.envfit <- envfit(plot2, dune.env)
vectors2 <- vectorfit.long(dune.envfit)

A1.surface <- ordisurf(plot2, y=A1)
A1.surface
A1.grid <- ordisurfgrid.long(A1.surface)

plotgg3 <- ggplot() + 
    geom_contour_filled(data=A1.grid, aes(x=x, y=y, z=z)) +
    geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
    geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +  
    xlab(axislabs[1, "label"]) +
    ylab(axislabs[2, "label"]) +  
    scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
    scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +    
    geom_point(data=sites2, aes(x=axis1, y=axis2, size=A1), shape=21, colour="black", fill="red") +
    geom_segment(data=subset(vectors2, vector=A1), aes(x=0, y=0, xend=axis1*2, yend=axis2*2), 
        size=1.2, arrow=arrow()) +
    geom_label_repel(data=subset(vectors2, vector=A1), aes(x=axis1*2, y=axis2*2, 
        label=vector), size=5) +
    BioR.theme +
    scale_fill_viridis_d() +
    scale_size(range=c(1, 20)) +
    labs(fill="A1") +
    coord_fixed(ratio=1)

plotgg3

# after_stat introduced in ggplot2 3.3.0
plotgg4 <- ggplot() + 
    geom_contour(data=A1.grid, aes(x=x, y=y, z=z, colour=factor(after_stat(level))), size=2) +
    geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
    geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +  
    xlab(axislabs[1, "label"]) +
    ylab(axislabs[2, "label"]) +  
    scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
    scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +    
    geom_point(data=sites2, aes(x=axis1, y=axis2, size=A1), shape=21, colour="black", fill="red") +
    geom_label_repel(data=sites2, aes(x=axis1, y=axis2, label=labels), 
        colour='black', size=4) +
    BioR.theme +
    scale_colour_viridis_d() +
    scale_size(range=c(1, 20)) +
    labs(colour="A1") +
    coord_fixed(ratio=1)

plotgg4

# example of Voronoi segmentation
plotgg5 <- ggplot(data=sites2, aes(x=axis1, y=axis2)) +
    geom_voronoi_tile(aes(fill = Management, group=-1L), max.radius=0.2) +
    geom_voronoi_segment(colour="grey50") +
    geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
    geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +  
    xlab(axislabs[1, "label"]) +
    ylab(axislabs[2, "label"]) +  
    scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
    scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +    
    geom_point() +
    BioR.theme +
    ggsci::scale_colour_npg() +
    coord_fixed(ratio=1)

plotgg5

# adding ellipse via ordiellipse

plot3 <- ordiplot(Ordination.model1, choices=c(1,2), scaling='species')
axislabs <- axis.long(Ordination.model1, choices=c(1 , 2))

Management.ellipses <- ordiellipse(plot3, groups=Management, display="sites", kind="sd")
Management.ellipses.long2 <- ordiellipse.long(Management.ellipses, grouping.name="Management")

plotgg6 <- ggplot() + 
    geom_vline(xintercept = c(0), color = "grey70", linetype = 2) +
    geom_hline(yintercept = c(0), color = "grey70", linetype = 2) +  
    xlab(axislabs[1, "label"]) +
    ylab(axislabs[2, "label"]) +  
    scale_x_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
    scale_y_continuous(sec.axis = dup_axis(labels=NULL, name=NULL)) +
    geom_polygon(data=Management.ellipses.long2, 
                   aes(x=axis1, y=axis2, colour=Management, 
                       fill=after_scale(alpha(colour, 0.2))), 
              size=0.2, show.legend=FALSE) +
    geom_point(data=sites1, aes(x=axis1, y=axis2, colour=Management, shape=Management), 
        size=5) +
    geom_segment(data=centroids.long(sites1, grouping=Management), 
        aes(x=axis1c, y=axis2c, xend=axis1, yend=axis2, colour=Management), 
        size=1, show.legend=FALSE) +
    BioR.theme +
    ggsci::scale_colour_npg() +
    coord_fixed(ratio=1)

plotgg6

# }
# NOT RUN {
 # dontrun
# }

Run the code above in your browser using DataLab