## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
## ----message=FALSE------------------------------------------------------------
{
library(bayestestR)
library(dplyr)
library(ggplot2)
library(ggdist)
library(ggtext)
library(nichetools)
library(purrr)
library(SIBER)
library(tidyr)
library(viridis)
}
## -----------------------------------------------------------------------------
glimpse(demo.siber.data.2)
## -----------------------------------------------------------------------------
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)
## -----------------------------------------------------------------------------
# ---- 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)
## -----------------------------------------------------------------------------
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", " (‰)")
)
## -----------------------------------------------------------------------------
demo_siber_data <- demo.siber.data.2 %>%
dplyr::select(iso1, iso2, group, community) %>%
arrange(community, group) %>%
as.data.frame()
## -----------------------------------------------------------------------------
siber_example <- createSiberObject(demo_siber_data)
## -----------------------------------------------------------------------------
# options for running jags
parms <- list()
parms$n.iter <- 2 * 10^4 # number of iterations to run the model for
parms$n.burnin <- 1 * 10^3 # discard the first set of values
parms$n.thin <- 10 # thin the posterior by this many
parms$n.chains <- 2 # run this many chains
## -----------------------------------------------------------------------------
# fit the ellipses which uses an Inverse Wishart prior on the
# covariance matrix Sigma, and a vague normal prior on the
# means.
priors <- list()
priors$R <- 1 * diag(2)
priors$k <- 2
priors$tau.mu <- 1.0E-3
## ----message=FALSE, results='hide'--------------------------------------------
ellipses_posterior <- siberMVN(siber_example, parms, priors)
## ----message = FALSE----------------------------------------------------------
df_mu <- extract_mu(ellipses_posterior, pkg = "SIBER",
data_format = "wide",
community_df = cg_names)
## -----------------------------------------------------------------------------
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", " (‰)")
)
## ----message=FALSE------------------------------------------------------------
df_mu_long <- extract_mu(ellipses_posterior, pkg = "SIBER",
community_df = cg_names)
## -----------------------------------------------------------------------------
df_sigma <- extract_sigma(ellipses_posterior, pkg = "SIBER")
## ----message = FALSE----------------------------------------------------------
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)
## -----------------------------------------------------------------------------
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", " (‰)")
)
## -----------------------------------------------------------------------------
sea_b <- siberEllipses(corrected.posteriors = ellipses_posterior)
## -----------------------------------------------------------------------------
seb_convert <- extract_niche_size(data = sea_b,
pkg = "SIBER",
community_df = cg_names)
## -----------------------------------------------------------------------------
group_ml <- groupMetricsML(siber_example)
## -----------------------------------------------------------------------------
group_convert <- extract_group_metrics(data = group_ml,
community_df = cg_names)
## -----------------------------------------------------------------------------
sea_c <- group_convert %>%
filter(metric == "SEAc")
## -----------------------------------------------------------------------------
ggplot() +
geom_violin(data = seb_convert, aes(x = community_name,
y = sea,
fill = group_name)) +
scale_fill_viridis_d(option = "A", begin = 0.25, end = 0.85,
name = "Groups", alpha = 0.75) +
geom_point(data = sea_c, aes(x = community_name,
y = est,
group = group_name,
),
size = 2.5,
fill = "white",
shape = 21,
position = position_dodge(width = 0.9)) +
theme_bw(
base_size = 15
) +
theme(
strip.background = element_blank(),
panel.grid = element_blank(),
legend.position = "inside",
legend.background = element_blank(),
legend.position.inside = c(0.85, 0.8)
) +
labs(
x = "Communities",
y = expression(paste("Niche Size p(", "‰"^2, "| x)"))
)
## -----------------------------------------------------------------------------
cg_names_within_com <- cg_names %>%
create_comparisons(comparison = "within")
## -----------------------------------------------------------------------------
ml_within_overlap <- cg_names_within_com %>%
map(~ maxLikOverlap(.x$cg_1, .x$cg_2, siber_example,
p.interval = 0.95, n = 100))
## -----------------------------------------------------------------------------
ml_95_within_com <- extract_similarities(ml_within_overlap,
type = "ml",
community_df = cg_names)
## -----------------------------------------------------------------------------
bayes95_overlap <- cg_names_within_com %>%
map(~ bayesianOverlap(.x$cg_1, .x$cg_2, ellipses_posterior,
draws = 100, p.interval = 0.95,
n = 100)
)
## -----------------------------------------------------------------------------
bays_95_overlap <- extract_similarities(bayes95_overlap,
type = "bay",
community_df = cg_names)
## -----------------------------------------------------------------------------
viridis_colors_s <- viridis(4, begin = 0.25, end = 0.85,
option = "A",
alpha = 0.75
)
## -----------------------------------------------------------------------------
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)"))
)
## -----------------------------------------------------------------------------
community_ml <- communityMetricsML(siber_example)
## -----------------------------------------------------------------------------
layman_ml <- extract_layman(community_ml,
type = "ml",
community_df = c_names)
## -----------------------------------------------------------------------------
mu_post <- extractPosteriorMeans(siber_example, ellipses_posterior)
## -----------------------------------------------------------------------------
layman_b <- bayesianLayman(mu.post = mu_post)
## -----------------------------------------------------------------------------
layman_be <- extract_layman(layman_b, community_df = c_names)
## -----------------------------------------------------------------------------
viridis_colors <- viridis(2, begin = 0.25, end = 0.85,
option = "G",
alpha = 0.75
)
## -----------------------------------------------------------------------------
ggplot() +
stat_pointinterval(
data = layman_be, aes(x = labels,
y = post_est,
point_fill = community_name),
point_size = 2.5,
interval_colour = "grey60",
position = position_dodge(0.4),
shape = 21
) +
geom_point(data = layman_ml, aes(x = labels,
y = estimate,
group = community_name),
shape = 21,
fill = "white",
alpha = 0.5,
position = position_dodge(0.4)) +
scale_fill_manual(name = "Community",
aesthetics = "point_fill",
values = viridis_colors) +
theme_bw(
base_size = 15
) +
theme(
panel.grid = element_blank(),
axis.text = element_markdown(),
legend.position = "inside",
legend.position.inside = c(0.88, 0.85),
) +
labs(
x = "Community Metrics",
y = expression(paste("p(", "‰", "|X)"))
)