Bring in trophic niche data
First we will load the necessary packages to preform the analysis and
visualization. We will use {nicheROVER}
and {nichetools} to
preform the analysis. We will use {dplyr}, {tidyr}, and {purrr} to manipulate data and
iterate processes. Lastly, we will use {ggplot2}, {ggtext}, and {patchwork} to plot,
add labels, and arrange plots.
I will add that many of the {dplyr} and
{tidyr} functions and processes can be replaced using {data.table}
which is great when working with large data sets.
{
library(dplyr)
library(ggplot2)
library(ggtext)
library(ggh4x)
library(nicheROVER)
library(nichetools)
library(patchwork)
library(purrr)
library(stringr)
library(tidyr)
}
For the purpose of the vignette we will be using the
fish data frame that is available within
{nicheROVER}. We will remove \(\delta\)34S for simplicity of
the vignette. If more than two isotopes or metrics are being used to
compare niche sizes and similarities, please use the functions for each
pairing. Right now some functions (i.e., niche_ellipse())
in {nichetools} doesn’t have the ability to work with more
than two isotopes. This will become a feature at some point but for now.
Please be patient and use the functions for each pairing you have.
We will first use the function janitor::clean_names() to
clean up column names. For your purposes you will need to replace fish
with your data frame either by loading a csv, rds, or qs, with your
data. You can do this multiple ways, I prefer using
readr::read_csv() but base R’s read.csv()
works perfectly fine.
df <- fish %>%
janitor::clean_names()
If there are any isotopic values that did not run and are
NA, they will need to be removed because
{nicheROVER}’s functions will not accommodate values of
NA.
Estimate posterior distribution with Normal-Inverse-Wishart (NIW)
priors.
We will take 1,000 posterior samples for each group. You can change
this but suggest nothing less than 1,000.
We will then split the data frame into a list with each species as a
data frame object within the list, We will then iterate over the list,
using map() from {purrr}, to estimate posterior
distribution using Normal-Inverse-Wishart (NIW) priors.
fish_par <- df %>%
split(.$species) %>%
map(~ select(., d13c, d15n)) %>%
map(~ niw.post(nsample = nsample, X = .))
Plot posterior distrubtion of μ and Σ
For most plotting within this vignette, I will split()
the data frame by isotope, creating a list that I will then use
imap() to iterate over the list to create plots. We will
use geom_density() to represent densities for both \(\mu\) and \(\Sigma\). Plot objects will then be stored
in a list.
First we will plot \(\mu\) for each
isotope. We will use {patchwork} to
configure plots for multi-panel figures. This package is phenomenal and
uses math operators to configure and manipulate the plots to create
multi-panel figures.
For labeling we are also going to use element_markdown()
from {ggtext} to work with
the labels that are needed to correctly display the isotopic signature.
If you are working other data please replace.
posterior_plots <- df_mu %>%
split(.$isotope) %>%
imap(
~ ggplot(data = ., aes(x = mu_est)) +
geom_density(aes(fill = sample_name), alpha = 0.5) +
scale_fill_viridis_d(begin = 0.25, end = 0.75,
option = "D", name = "Species") +
theme_bw() +
theme(panel.grid = element_blank(),
axis.title.x = element_markdown(),
axis.title.y = element_markdown(),
legend.position = "none",
legend.background = element_blank()
) +
labs(
x = paste("\u00b5<sub>\U03B4</sub>", "<sub><sup>",
unique(.$neutron), "</sup></sub>",
"<sub>",unique(.$element), "</sub>", sep = ""),
y = paste0("p(\u00b5 <sub>\U03B4</sub>","<sub><sup>",
unique(.$neutron), "</sub></sup>",
"<sub>",unique(.$element),"</sub>",
" | X)"), sep = "")
)
posterior_plots$d15n +
theme(legend.position = c(0.18, 0.82)) +
posterior_plots$d13c
For labeling purposes we need to add columns that are the element
abbreviation and neutron number. I do this by using
case_when() which are vectorized if else statements.
df_sigma_cn <- df_sigma_cn %>%
mutate(
element_id = case_when(
id == "d15n" ~ "N",
id == "d13c" ~ "C",
),
neutron_id = case_when(
id == "d15n" ~ 15,
id == "d13c" ~ 13,
),
element_iso = case_when(
isotope == "d15n" ~ "N",
isotope == "d13c" ~ "C",
),
neutron_iso = case_when(
isotope == "d15n" ~ 15,
isotope == "d13c" ~ 13,
)
)
Next we will plot the posteriors for \(\Sigma\).
sigma_plots <- df_sigma_cn %>%
group_split(id, isotope) %>%
imap(
~ ggplot(data = ., aes(x = post_sample)) +
geom_density(aes(fill = sample_name), alpha = 0.5) +
scale_fill_viridis_d(begin = 0.25, end = 0.75,
option = "D", name = "Species") +
theme_bw() +
theme(panel.grid = element_blank(),
axis.title.x = element_markdown(),
axis.title.y = element_markdown(),
legend.position = "none"
) +
labs(
x = paste("\U03A3","<sub>\U03B4</sub>",
"<sub><sup>", unique(.$neutron_id), "</sub></sup>",
"<sub>",unique(.$element_id),"</sub>"," ",
"<sub>\U03B4</sub>",
"<sub><sup>", unique(.$neutron_iso), "</sub></sup>",
"<sub>",unique(.$element_iso),"</sub>", sep = ""),
y = paste("p(", "\U03A3","<sub>\U03B4</sub>",
"<sub><sup>", unique(.$neutron_id), "</sub></sup>",
"<sub>",unique(.$element_id),"</sub>"," ",
"<sub>\U03B4</sub>",
"<sub><sup>", unique(.$neutron_iso), "</sub></sup>",
"<sub>",unique(.$element_iso),"</sub>", " | X)", sep = ""),
)
)
sigma_plots[[1]] +
theme(legend.position = c(0.1, 0.82))
Estimate niche ellipse
We then will use niche_ellipse() to easily extract
ellipse for each \(\Sigma\) estimate
(i.e., 1000). If you are to have additional isotopes or metrics, you
will need to create mu and sigma objects for each pairing, as currently
this function only handles two isotopes. In the future, there likely
will be the ability to specify the number of isotopes you have with the
default being two. The reason for the lack of functionality is
ellipse::ellipse() can only work within two-dimensions, not
three, so you will have to create multiple ellipse() calls
for each combination of isotopes or metrics and I haven’t had the time
to implement this. The function will also tell you how long it took to
process as with large sets of isotope data it is nice to know the time
it takes for the function to work. The function also has the argument
random which by default is set to TRUE. This
argument randomly subsamples and returns 10 ellipse estimates out of the
total number of samples taken in this case it is 1,000. The
set_seed argument allows you to change the
set.seed value by giving a numerical value to make the
results of the function reproducable, default is a random value. I
highly suggest using set_seed otherwise it will subsample different
values, there is no default value because CRAN will not allow for a
default value. You can change the number of subsamples but 10 seems
pretty standard. If you’d like all 1,000 ellipses you can set
random to FALSE.
ellipse_df <- niche_ellipse(dat_mu = df_mu, dat_sigma = df_sigma,
set_seed = 4)
#> → Total time processing was 0.03 secs
Plot ellipses, densities of each istope, and isotope biplot
We will first plot the ellipse for each sample_name
ellipse_plots <- ggplot() +
geom_polygon(data = ellipse_df,
mapping = aes(x = d13c, y = d15n,
group = interaction(sample_number, sample_name),
color = sample_name),
fill = NA,
linewidth = 0.5) +
scale_colour_viridis_d(begin = 0.25, end = 0.75,
option = "D", name = "species",
) +
scale_x_continuous(breaks = rev(seq(-20, -40, -2))) +
scale_y_continuous(breaks = seq(6, 16, 2)) +
theme_bw(base_size = 10) +
theme(axis.text = element_text(colour = "black"),
panel.grid = element_blank(),
legend.position = "none",
legend.title = element_text(hjust = 0.5),
legend.background = element_blank()) +
labs(x = expression(paste(delta ^ 13, "C")),
y = expression(paste(delta ^ 15, "N")))
We need to turn df into long format to iterate over
using imap() to easily create density plots. You will
notice that I again use case_when() to make columns of
element abbreviations and neutron numbers that will be used in plot
labeling.
iso_long <- df %>%
pivot_longer(cols = -species,
names_to = "isotope",
values_to = "value") %>%
mutate(
element = case_when(
isotope == "d15n" ~ "N",
isotope == "d13c" ~ "C",
),
neutron = case_when(
isotope == "d15n" ~ 15,
isotope == "d13c" ~ 13,
)
)
We will then make density plots for each isotope using
geom_density()
iso_density <- iso_long %>%
group_split(isotope) %>%
imap(
~ ggplot(data = .) +
geom_density(aes(x = value,
fill = species),
alpha = 0.35,
linewidth = 0.8) +
scale_fill_viridis_d(begin = 0.25, end = 0.75,
option = "D", name = "Species") +
theme_bw(base_size = 10) +
theme(axis.text = element_text(colour = "black"),
panel.grid = element_blank(),
legend.position = c(0.15, 0.55),
legend.background = element_blank(),
axis.title.x = element_markdown(family = "sans")) +
labs(x = paste("\U03B4",
"<sup>", unique(.$neutron), "</sup>",unique(.$element),
sep = ""),
y = "Density")
)
d13c_density <- iso_density[[1]] +
scale_x_continuous(breaks = rev(seq(-20, -34, -2)),
limits = rev(c(-20, -34)))
d15n_density <- iso_density[[2]] +
scale_x_continuous(breaks = seq(5, 15, 2.5),
limits = c(5, 15)) +
theme(
legend.position = "none"
)
Lastly we will use geom_point() to make isotopic
biplot.
iso_biplot <- ggplot() +
geom_point(data = df, aes(x = d13c, y = d15n,
fill = species),
shape = 21, colour = "black",
stroke = 0.8,
size = 3, alpha = 0.70) +
scale_fill_viridis_d(begin = 0.25, end = 0.75,
option = "D", name = "species") +
scale_x_continuous(breaks = rev(seq(-20, -39, -1))) +
scale_y_continuous(breaks = seq(5, 17, 1)) +
theme_bw(base_size = 10) +
theme(axis.text = element_text(colour = "black"),
panel.grid = element_blank(),
legend.position = "none",
legend.background = element_blank()) +
labs(x = expression(paste(delta ^ 13, "C")),
y = expression(paste(delta ^ 15, "N")))
Determine the 95% niche similarties for each species
We will use the overlap() function from {nicheROVER})
to estimate the percentage of similarity among species. We will set
overlap to assess based on 95% similarities.
over_stat <- overlap(fish_par, nreps = nsample, nprob = 1000,
alpha = 0.95)
We then are going transform this output to a data frame using
extract_overlap() plotting so we can assess overall
similarities among species.
over_stat_df <- extract_overlap(data = over_stat)
We then are going to take our newly made data frame and extract out
the median percentage of similarities and the 2.5% and 97.5%
quantiles.
over_sum <- over_stat_df %>%
group_by(sample_name_a, sample_name_b) %>%
summarise(
median_niche_overlap = round(median(niche_overlap_perc), digits = 2),
qual_2.5 = round(quantile(niche_overlap_perc,
probs = 0.025, na.rm = TRUE), digits = 2),
qual_97.5 = round(quantile(niche_overlap_perc,
probs = 0.975, na.rm = TRUE), digits = 2)
) %>%
ungroup() %>%
pivot_longer(cols = -c(sample_name_a, sample_name_b, median_niche_overlap),
names_to = "percentage",
values_to = "niche_overlap_qual") %>%
mutate(
percentage = as.numeric(str_remove(percentage, "qual_"))
)
We are now going to use ggplot(),
geom_violin(), and stat_summary() to represent
the posterior distributions and the median of the posterior
distributions. Representations of posterior distributions should either
use the mode or median, not the mean. See the following publication as a
guide to representation of posterior distributions in vizialations and
in text. You can extract credible intervals or the equal-tailed
intervals using {bayestestR}.
ggplot(data = over_stat_df, aes(x = sample_name_a,
y = niche_overlap_perc,
fill = sample_name_b)) +
geom_violin() +
stat_summary(fun.y = median, geom = "point",
size = 3,
position = position_dodge(width = 0.9)) +
geom_vline(xintercept = seq(1.5, 3.5, 1),
linetype = 2) +
scale_fill_viridis_d(begin = 0.25, end = 0.75,
option = "D", name = "Species",
alpha = 0.35) +
theme_bw() +
theme(
panel.grid = element_blank(),
axis.text = element_text(colour = "black"),
legend.background = element_blank(),
strip.background = element_blank()
) +
labs(x = paste("Overlap Probability (%)", "\u2013",
"Niche Region Size: 95%"),
y = "p(Percent Overlap | X)")
Plot niche size
We will now use geom_violin(),
geom_point(), and geom_errorbar() to plot the
distribution for niche size for each species.
ggplot(data = niche_size,
aes(x = sample_name, y = niche_size)) +
geom_violin(
width = 0.2) +
stat_summary(fun.y = median, geom = "point",
size = 3,
position = position_dodge(width = 0.9)) +
theme_bw(base_size = 15) +
theme(panel.grid = element_blank(),
axis.text = element_text(colour = "black")) +
labs(x = "Species",
y = "Niche Size")
Now that we have our niche sizes and similarities determined we can
make inferences about the species, trophic similarities, and the
ecosystem.
