I'm using the SIBER package to analyze some data, and am following this tutorial, but it doesn't seem to be working with my data, and I'm not sure if it's a statistics, coding, or sample size problem (the way I've subdivided my groups, I end up with pretty few datapoints per group). My goal is to produce Bayesian estimates of ellipse dimensions for different combinations of age class, year, and tissue type, but the ellipses I've generated don't seem to align with the data. I've tried increasing the number of iterations and burnin to address any convergence issues, but the ellipses still don't really contain the data. I'm fine with re-grouping my data to answer a different set of questions, but I'd like to confirm that I haven't messed up in a Bayesian or coding way first.
library(tidyr)
library(SIBER)
library(dplyr)
library(ggplot2)
# First we outline the model parameters as a list
parms <- list()
parms$n.iter <- 10 *10^4
parms$n.burnin <- 2 *10^3
parms$n.thin <- 10
parms$n.chains <- 3
# Then we establish priors
# Inverse wishart on covariance matrix of signma
# Vague normal for means
priors <- list()
priors$R <- 1 * diag(2)
priors$k <- 2
priors$tau.mu <- 1.0E-3
# Here is the dataset
siber_df <- structure(list(iso1 = c(-21.86, -22.57, -22.6, -22, -22.68, -22.48,
-23.09, -21.85, -21.68, -22.19, -22.5, -21.59, -21.92, -22.04,
-22.23, -22.4, -20.79, -20.72, -20.76, -21.16, -22.56, -22.76,
-22.29, -21.91, -21.97, -22.61, -22.73, -21.61, -21.91, -23.58,
-23.36, -22.94, -22.34, -22.71, -21.92, -22.09, -24.37, -21.8,
-22.34, -22.17, -22.54, -22.63, -22.43, -20.63, -22.47, -21.55,
-21.03, -22.19, -24.8, -21.9, -21.83, -21, -21.64, -21.77, -21.39,
-20.39, -21.69, -21.18, -21.28, -22.12, -21.95, -21.29, -21.52,
-21.29, -21.61, -21.61, -22.38, -21.8, -21.72, -22.5, -21.26,
-21.45, -22.39, -22.18, -22.19, -22.48, -22.52, -22.15, -22.04,
-20.81, -20.44, -21.29, -22.04, -22.21, -22.18, -21.96, -21.96,
-21.82, -22.12, -22.55, -22.43, -22.38, -22.07, -22.66, -20.76,
-20.56, -20.65, -22.12, -21.97, -22.105, -22.11, -22.37, -22.48,
-21.32, -21.42, -21.2, -20.72, -21.74, -22.1, -21.5, -21.68,
-21.58, -21.63, -21.43, -22.22, -22.33, -20.49, -20.8, -20.69,
-21.19, -21.52, -22.19, -21.8, -21.24, -21.62, -20.96, -20.79,
-21.665, -21.78, -21.76, -21.58, -21.66, -20.93, -20.935, -21.87,
-21.79, -22.33, -21.01, -21.74, -22.08, -20.43, -21.6, -22.1,
-22.07, -21.67, -21.48, -21.67, -21.64, -21.82, -21.89, -21.35,
-21.3, -21.78, -21.57, -21.89, -21.37, -21.97, -22.12, -22.47,
-22.28, -21.14, -20.58, -20.71, -21.17, -20.97), iso2 = c(7.26,
5.74, 5.4, 6.57, 5.99, 6.17, 5.96, 8.81, 6.66, 5.54, 5.55, 6.63,
6.32, 6.71, 6.26, 7.24, 6.08, 6.28, 6.18, 6.02, 7.18, 7.39, 5.89,
5.12, 7.37, 4.08, 7.12, 7.26, 7.49, 5.36, 5.13, 6.75, 6.19, 5.28,
7.15, 7.99, 4.35, 7.23, 6.37, 5.81, 6.95, 5.92, 6.17, 7.15, 7.56,
7.68, 6.97, 6.28, 7.94, 6.58, 5.29, 6.45, 6.47, 6.96, 6.38, 7.85,
6.42, 5.61, 5.57, 5.19, 5.46, 5.85, 5.82, 6.14, 4.5, 4.34, 6.73,
6.81, 6.17, 6.87, 7.3, 7.72, 7.7, 2.67, 3.4, 5.29, 4.92, 3.43,
3.01, 4.95, 5.24, 4.86, 3.21, 2.76, 3.26, 4.61, 4.45, 4.41, 5.44,
2.72, 4.07, 3.85, 5.37, 2.58, 4.45, 4.96, 5.06, 5.04, 5.03, 3.42,
3.46, 5.08, 4.98, 6.62, 6.68, 5.38, 5.95, 4.68, 4.37, 3.6, 3.65,
5.14, 5.03, 4.98, 3.4, 3.59, 5.05, 4.89, 5.07, 6.43, 6.39, 4.67,
4.84, 6.73, 6.3, 6.41, 6.42, 4.48, 4.13, 5.89, 6.43, 4.8, 6.99,
6.42, 6.14, 6.19, 5.75, 4.97, 5.91, 5.59, 6.53, 5.52, 6.15, 6.54,
5.88, 5.38, 5.73, 5.57, 6.09, 5.21, 6.04, 5.88, 5.95, 6.1, 5.87,
6.01, 5.52, 5.4, 5, 5.3, 6.04, 6.27, 6.33, 6.39, 6.22), group = c("ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "NL", "NL", "NL", "NL", "NL", "NL", "ADULT", "FL", "FL",
"FL", "FL", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "HY", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "HY", "HY", "HY", "HY", "HY", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "HY", "ADULT", "FL",
"FL", "ADULT", "ADULT", "ADULT", "FL", "FL", "NL", "NL", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "FL", "FL",
"FL", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "NL", "NL", "NL", "NL", "NL", "NL",
"NL", "ADULT", "ADULT", "FL", "FL", "FL", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "HY", "HY", "ADULT", "ADULT",
"ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "ADULT", "HY", "HY",
"HY", "HY", "HY", "HY", "HY", "HY", "ADULT", "ADULT", "ADULT",
"ADULT", "ADULT", "HY", "HY"), community = c("2019_feathers",
"2021_feathers", "2021_feathers", "2019_feathers", "2021_feathers",
"2021_feathers", "2021_feathers", "2019_feathers", "2019_feathers",
"2019_feathers", "2019_feathers", "2019_feathers", "2019_feathers",
"2019_feathers", "2019_feathers", "2019_feathers", "2019_feathers",
"2019_feathers", "2019_feathers", "2019_feathers", "2019_feathers",
"2019_feathers", "2021_feathers", "2021_feathers", "2019_feathers",
"2019_feathers", "2019_feathers", "2019_feathers", "2019_feathers",
"2021_feathers", "2021_feathers", "2021_feathers", "2021_feathers",
"2021_feathers", "2021_feathers", "2021_feathers", "2021_feathers",
"2021_feathers", "2021_feathers", "2021_feathers", "2021_feathers",
"2021_feathers", "2021_feathers", "2021_feathers", "2021_feathers",
"2021_feathers", "2021_feathers", "2021_feathers", "2021_feathers",
"2021_feathers", "2021_feathers", "2021_feathers", "2021_feathers",
"2021_feathers", "2021_feathers", "2021_feathers", "2019_feathers",
"2019_feathers", "2019_feathers", "2019_feathers", "2019_feathers",
"2019_feathers", "2019_feathers", "2019_feathers", "2019_feathers",
"2019_feathers", "2019_feathers", "2019_feathers", "2019_feathers",
"2019_feathers", "2019_feathers", "2019_feathers", "2019_feathers",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2021_claws",
"2021_claws", "2019_claws", "2019_claws", "2021_claws", "2021_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2019_claws", "2021_claws", "2021_claws", "2019_claws",
"2019_claws", "2019_claws", "2019_claws", "2019_claws", "2019_claws",
"2019_claws", "2021_claws", "2021_claws", "2021_claws", "2021_claws",
"2021_claws", "2021_claws", "2021_claws", "2021_claws", "2021_claws",
"2021_claws", "2021_claws", "2021_claws", "2021_claws", "2021_claws",
"2021_claws", "2021_claws", "2021_claws", "2021_claws", "2021_claws",
"2021_claws", "2021_claws", "2021_claws", "2021_claws", "2021_claws",
"2021_claws", "2021_claws", "2021_claws", "2021_claws", "2021_claws",
"2021_claws", "2021_claws")), class = "data.frame", row.names = c(NA,
-165L))
# Converting the whole dataframe into a siber object
siber_object <- createSiberObject(siber_df)
# Calculating group metrics (basically only so I can use the labels for plotting later)
group.ML <- groupMetricsML(siber_object)
# Labeling the ellipses for the model to identify
# These must have periods or they will not be recognized by siber
claws_19_ad <- "2019_claws.ADULT"
claws_19_fl <- "2019_claws.FL"
claws_19_nl <- "2019_claws.NL"
feath_19_ad <- "2019_feathers.ADULT"
feath_19_nl <- "2019_feathers.NL"
feath_19_fl <- "2019_feathers.FL"
claws_21_ad <- "2021_claws.ADULT"
claws_21_hy <- "2021_claws.HY"
feath_21_ad <- "2021_feathers.ADULT"
feath_21_hy <- "2021_feathers.HY"
# Getting a posterior distribution
community_posterior <- siberMVN(siber_object, parms, priors)
SEA.B <- siberEllipses(community_posterior)
siberDensityPlot(SEA.B, xticklabels = colnames(group.ML),
xlab = c("Community | Group"),
ylab = expression("Standard Ellipse Area " ('permille' ^2) ),
bty = "L",
las = 1,
main = "SIBER ellipses on each group"
)
# how many of the posterior draws do you want?
n.posts <- 10
# decide how big an ellipse you want to draw
p.ell <- 0.95
# for a standard ellipse use
# p.ell <- pchisq(1,2)
# a list to store the results
all_ellipses <- list()
# loop over groups
for (i in 1:length(community_posterior)){
# a dummy variable to build in the loop
ell <- NULL
post.id <- NULL
for ( j in 1:n.posts){
# covariance matrix
Sigma <- matrix(community_posterior[[i]][j,1:4], 2, 2)
# mean
mu <- community_posterior[[i]][j,5:6]
# ellipse points
out <- ellipse::ellipse(Sigma, centre = mu , level = p.ell)
ell <- rbind(ell, out)
post.id <- c(post.id, rep(j, nrow(out)))
}
ell <- as.data.frame(ell)
ell$rep <- post.id
all_ellipses[[i]] <- ell
}
ellipse_df <- bind_rows(all_ellipses, .id = "id")
# now we need the group and community names
# extract them from the community_posterior list
group_comm_names <- names(community_posterior)[as.numeric(ellipse_df$id)]
# split them and conver to a matrix, NB byrow = T
split_group_comm <- matrix(unlist(strsplit(group_comm_names, "[.]")),
nrow(ellipse_df), 2, byrow = TRUE)
ellipse_df$community <- split_group_comm[,1]
ellipse_df$group <- split_group_comm[,2]
ellipse_df <- dplyr::rename(ellipse_df, iso1 = x, iso2 = y)
first.plot <- ggplot(data = siber_df, aes(iso1, iso2)) +
geom_point(aes(color = factor(group):factor(community)), size = 2)+
ylab(expression(paste(delta^{15}, "N (permille)")))+
xlab(expression(paste(delta^{13}, "C (permille)"))) +
theme(text = element_text(size=15))
print(first.plot)
second.plot <- first.plot + facet_wrap(~factor(group):factor(community))
print(second.plot)
third.plot <- second.plot +
geom_polygon(data = ellipse_df,
mapping = aes(iso1, iso2,
group = rep,
color = factor(group):factor(community),
fill = NULL),
fill = NA,
alpha = 0.2)
print(third.plot)