Bring in trophic data
First we will load the necessary packages to preform the analysis and
visualization. We will use {SIBER} and {nichetools} to
preform the analysis. We will use {bayestestR} to
calculate and extract medians and Equal-Tailed Interval (ETI) of
posterior distributions that we want to plot. We will use {dplyr}, {tidyr}, and {purrr} to manipulate data and
iterate processes. Lastly, we will use {ggplot2}, {ggtext}, and {ggdist} to plot and add
labels.
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.
We will first load the packages that were just mentioned
{
library(bayestestR)
library(dplyr)
library(ggplot2)
library(ggdist)
library(ggtext)
library(nichetools)
library(purrr)
library(SIBER)
library(tidyr)
library(viridis)
}
For the purpose of the vignette we will be using the
demo.siber.data.2 data frame that is available within
{SIBER}. We will first look at the structure of this data
frame using the {dplyr} function glimpse().
For your purposes you will need to replace this with your data frame
either by loading a csv, rds, or qs file. You can do this multiple ways,
I prefer using readr::read_csv() for a csv but base R’s
read.csv() works perfectly fine. Note, {SIBER}
functions do not take tibbles or data.tables
so you will have to convert either to data.frame class
prior to running functions in {SIBER}.
glimpse(demo.siber.data.2)
#> Rows: 150
#> Columns: 4
#> $ iso1 <dbl> -11.048902, -8.809360, -9.256212, -8.380952, -10.401561, -7.…
#> $ iso2 <dbl> 1.02044315, 1.97597511, 3.28707260, -1.95750266, -3.11949113…
#> $ group <chr> "city", "city", "city", "city", "city", "city", "city", "cit…
#> $ community <chr> "dublin", "dublin", "dublin", "dublin", "dublin", "dublin", …
You will notice that the community and group column are character
strings that are the actual names of the communities and groups. I
advise changing them into factors, thus allow you to know the order for
each column prior to converting them into a numeric and
then a charcter. The reason why this is important, is that
functions in {SIBER} use for loops based on
the indexing of the SiberObject. If this order does not
match up you can have issues with the names of the communities and
groups you are working with.
Let us change these into a factor, then preserve the
column and create an id column that is the numerical order that will
become the community and group names provided to
createSiberObject().
demo.siber.data.2 <- demo.siber.data.2 %>%
mutate(
group = factor(group),
community = factor(community),
group_id = as.numeric(group) %>%
as.character(),
community_id = as.numeric(community) %>%
as.character()
) %>%
rename(
group_name = group,
community_name = community,
group = group_id,
community = community_id
)
glimpse(demo.siber.data.2)
#> Rows: 150
#> Columns: 6
#> $ iso1 <dbl> -11.048902, -8.809360, -9.256212, -8.380952, -10.401561…
#> $ iso2 <dbl> 1.02044315, 1.97597511, 3.28707260, -1.95750266, -3.119…
#> $ group_name <fct> city, city, city, city, city, city, city, city, city, c…
#> $ community_name <fct> dublin, dublin, dublin, dublin, dublin, dublin, dublin,…
#> $ group <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", …
#> $ community <chr> "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", "1", …
After we have done this, we are going to create two data frames that
are the names of our communities and groups with their associated id
values. We will use these data frames later on to join up the actual
names, allowing us to know what estimates belong to which communities
and groups. We are doing this because it is unlikely you will have
communities and groups named 1, 2, 3 ect. and instead will have actual
names.
# ---- create name with group and community data frame ----
cg_names <- demo.siber.data.2 %>%
distinct(group,
community,
group_name,
community_name) %>%
arrange(community, group)
# ---- create community names data frame ----
c_names <- demo.siber.data.2 %>%
distinct(community,
community_name) %>%
arrange(community)
We will then plot our biplot to confirm we have the correct
structure.
ggplot(data = demo.siber.data.2, aes(x = iso1, y = iso2,
colour = group_name)) +
geom_point() +
facet_wrap(~ community_name) +
scale_colour_viridis_d(option = "A", begin = 0.25, end = 0.85,
name = "Groups", alpha = 0.75) +
theme_bw(
base_size = 15
) +
theme(
strip.background = element_blank(),
panel.grid = element_blank(),
axis.title = element_markdown(),
legend.position = "inside",
legend.position.inside = c(0.65, 0.75)
) +
labs(
x = paste0("\U03B4","<sup>", 13, "</sup>", "C", " (‰)"),
y = paste0("\U03B4","<sup>", 15, "</sup>", "N", " (‰)")
)
Next we will grab the isotopes we need and the community and group
ids that have already been renamed to community and
group. This is important as creatSiberObject()
will 1) only take the following order with the following names
iso1, iso2, group, and
community and 2) we will transform this tibble
into a data.frame as {SIBER} will only work
with a data.frame. In this case we were using
tibbles but we could also be working with
data.table and will need to do the same thing.
demo_siber_data <- demo.siber.data.2 %>%
dplyr::select(iso1, iso2, group, community) %>%
arrange(community, group) %>%
as.data.frame()
Extract posterior distributions for μ and Σ
We will first extract posterior distribution for \(\mu\) using extract_mu() in
{nichetools}. We will need to set the argument
pkg to "SIBER" and the argument
data_format to "wide". This argument takes
"long" or "wide" which dictates whether the
data object returned is in wide or long format. We will also use the
function seperate_wider_delim() from {tidyr}
to separate the community and groups names as they are joined by a
.. We then will use left_join() and the
cg_names data frame we created above to add in the correct
community and group names.
df_mu <- extract_mu(ellipses_posterior, pkg = "SIBER",
data_format = "wide",
community_df = cg_names)
We can confirm that the posterior estimates of \(\mu\) are correct by plotting them with
{ggplot2}.
ggplot() +
geom_point(data = df_mu, aes(x = d13c, y = d15n,
colour = group_name)) +
geom_point(data = demo.siber.data.2, aes(x = iso1, y = iso2,
colour = group_name)) +
facet_wrap( ~ community_name) +
scale_colour_viridis_d(option = "A", begin = 0.25, end = 0.85,
name = "Groups", alpha = 0.75) +
theme_bw(
base_size = 15
) +
theme(
strip.background = element_blank(),
panel.grid = element_blank(),
axis.title = element_markdown(),
legend.position = "inside",
legend.position.inside = c(0.65, 0.75)
) +
labs(
x = paste0("\U03B4","<sup>", 13, "</sup>", "C", " (‰)"),
y = paste0("\U03B4","<sup>", 15, "</sup>", "N", " (‰)")
)
Notice the density of points in the center of the raw data for the
corresponding colours. This density of points are the posterior
estimates from the model for \(\mu\)
and are an indication that the model is iterating over the groups and
communities correctly.
We are also going to extract \(\mu\)
in long format for creating ellipse, as the function that create ellipse
needs the data frame of \(\mu\) to be
in long format.
df_mu_long <- extract_mu(ellipses_posterior, pkg = "SIBER",
community_df = cg_names)
Next we are going to extract posterior estimates of \(\Sigma\) using
extract_sigma().
df_sigma <- extract_sigma(ellipses_posterior, pkg = "SIBER")
Lastly, we are going to feed each \(\mu\) and \(\Sigma\) estimate to
niche_ellipse() to estimate each ellipse. There are few
things to know about this function and they include the following: 1) it
will randomly sample 10 ellipse from the total posterior distribution of
\(\mu\) and \(\Sigma\), this seems quite standard,
however, you can adjust the number of samples by changing the argument
n. 2) to make the function consistently randomly sample the
same set you will need to set set_seed to a numerical
value. If this is not set then it will randomly sample a different set
of 10 ellipses every time. 3) if you would like the function to not
randomly sample set the argument random to
FALSE. 4) by default it will tell you how long it takes to
generate the ellipse and will have progress bars at each step. If you
want to turn this off set message to FALSE. 5)
if you are wanting to change the confidence level of the ellipse you can
do so using the argument p_ell. This value is bound between
0 and 1.
df_el <- niche_ellipse(dat_mu = df_mu_long,
dat_sigma = df_sigma,
set_seed = 4, n = 20) %>%
separate_wider_delim(sample_name, cols_remove = FALSE,
delim = ".", names = c("community",
"group")) %>%
left_join(cg_names)
Now that we have the ellipses created we can plot them.
ggplot(data = df_el,
aes(x = d13c, y = d15n,
group = interaction(sample_number,
sample_name),
colour = group_name)) +
geom_polygon(linewidth = 0.5, fill = NA) +
facet_wrap( ~ community_name) +
scale_colour_viridis_d(option = "A", begin = 0.25, end = 0.85,
name = "Groups", alpha = 0.75) +
theme_bw(
base_size = 15
) +
theme(
strip.background = element_blank(),
panel.grid = element_blank(),
legend.position = "inside",
axis.title = element_markdown(),
legend.background = element_blank(),
legend.position.inside = c(0.65, 0.75)
) +
labs(
x = paste0("\U03B4","<sup>", 13, "</sup>", "C", " (‰)"),
y = paste0("\U03B4","<sup>", 15, "</sup>", "N", " (‰)")
)
Niche Similarties
Now that we have extracted Bayesian estimates of niche size we are
likely wanting to know how much do these niches have in common or are
similar.
We can use maximum-likelihood and Bayesian frameworks to estimate the
percentage of similarity within communities between groups and/or among
communities with groups being consistent. We can use functions in
{SIBER} to create maximum-likelihood and Bayesian estimates
of niche similarity followed by functions in {nichetools}
to extract these similarities.
First we need to create our comparisons that we are wanting to
evaluate. We can use the function create_comparisons() to
generate a list that has two-column data frames of each comparison. You
can change the argument comparison to "among"
to compare communities for the same group, versus "within"
which compares groups within a community.
cg_names_within_com <- cg_names %>%
create_comparisons(comparison = "within")
Next we can feed this listed data frames to either
maxLikOverlap() or bayesianOverlap() using
map() from {purrr}. For this exercise I have
not changed .progress argument of map() but
when working with large data sets I often turn this argument to
TRUE to provide a progress update.
ml_within_overlap <- cg_names_within_com %>%
map(~ maxLikOverlap(.x$cg_1, .x$cg_2, siber_example,
p.interval = 0.95, n = 100))
Next we can extract estimates using
extract_similarities() with the type argument
set to "ml".
ml_95_within_com <- extract_similarities(ml_within_overlap,
type = "ml",
community_df = cg_names)
Now we will repeat the process for bayesianOverlap(),
first supplying the function with our list of data frames, next having
map() iterate over this list.
bayes95_overlap <- cg_names_within_com %>%
map(~ bayesianOverlap(.x$cg_1, .x$cg_2, ellipses_posterior,
draws = 100, p.interval = 0.95,
n = 100)
)
Next we can extract estimates using
extract_similarities() with the type argument
set to "bay".
bays_95_overlap <- extract_similarities(bayes95_overlap,
type = "bay",
community_df = cg_names)
Now that we have extracted maximum-likelihood and Bayesian estimates
we can visualize them.
Prior to creating a point interval plot, we need to create a colour
palette that will be used to identify each community.
viridis_colors_s <- viridis(4, begin = 0.25, end = 0.85,
option = "A",
alpha = 0.75
)
Now we can use point interval plots to visually represent posterior
distributions.
ggplot() +
stat_pointinterval(data = bays_95_overlap,
aes(x = group_1,
y = prop_overlap,
point_fill = group_2),
interval_colour = "grey60",
point_size = 3,
shape = 21,
position = position_dodge(0.4)) +
geom_point(data = ml_95_within_com, aes(x = group_1,
y = prop_overlap,
group = group_2),
shape = 21,
fill = "white",
size = 2,
alpha = 0.5,
position = position_dodge(0.4)) +
scale_fill_manual(name = "Group",
aesthetics = "point_fill",
values = viridis_colors_s) +
theme_bw(
base_size = 15
) +
theme(
panel.grid = element_blank(),
strip.background = element_blank(),
legend.position = "inside",
legend.position.inside = c(0.15, 0.80),
) +
labs(
x = "Group",
y = expression(paste("p(", "‰", "|X)"))
)
