# 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