diff --git a/estimation.R b/estimation.R
new file mode 100644
index 0000000000000000000000000000000000000000..5ac35a783d4c51b2f581eb69743c9842eab5d268
--- /dev/null
+++ b/estimation.R
@@ -0,0 +1,407 @@
+library(INLA)
+library(HDInterval)
+library(modeest)
+library(mvtnorm)
+source("load_data.R")
+
+# Required functions
+
+f_G <- function(x){
+ x * dnorm(x) / pnorm(x)
+}
+
+f_z <- function(eta, rho2){
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2))) /
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2) + sqrt(rho2)) + 2 * eta + rho2)
+}
+
+f_sum <- function(x, p = 0.95){
+ x <- x$value
+ tmp <- c(mean(x), sd(x), median(x), meanshift(x), hdi(x, p),
+ quantile(x, probs = c((1 - p)/2, 1 - (1 - p) / 2)))
+ names(tmp)[1:4] <- c("mean", "sd", "median", "mode")
+ return(as_tibble_row(tmp))
+}
+
+# Number of samples from posterior distribution
+# Start off with fewer first!
+sample_size <- 100
+# sample_size <- 10000
+
+# Basic model with all the variance components
+data_ammarnas_inla <- data_ammarnas %>%
+ mutate(species = as.numeric(as.factor(species)),
+ year2 = year)
+
+data_ammarnas_inla %>%
+ ggplot(aes(x = year, y = abundance)) +
+ geom_line(aes(group = species))
+
+m01_am <- inla(abundance ~ f(year,
+ model = "ou",
+ group = species,
+ control.group = list(model = "iid")) +
+ f(species, model = "iid") +
+ f(year2,
+ model = "ou"),
+ family = "poisson",
+ data = data_ammarnas_inla)
+
+# Precision for year and species seem ok!
+summary(m01_am)
+
+tmp_am_par <- inla.hyperpar.sample(sample_size, m01_am) %>%
+ as_tibble() %>%
+ mutate(se2 = 1 / `Precision for year`,
+ gamma = `Phi for year`,
+ gamma_c = `Phi for year2`,
+ sh2 = 1 / `Precision for species`,
+ sc2 = 1 / `Precision for year2`,
+ eta = inla.rmarginal(sample_size,
+ m01_am$marginals.fixed$`(Intercept)`),
+ rho2 = se2 + sh2 + sc2,
+ z = f_z(eta, rho2))
+
+# How many species can there be in total:
+n_years <- n_distinct(data_ammarnas$year)
+n_sp <- 100
+
+
+p0_am = c()
+for (i in 1:sample_size){
+ sim <- tibble(species = 1:n_sp,
+ h_i = rnorm(n_sp, mean = 0,
+ sd = sqrt(tmp_am_par$sh2[i]))) %>%
+ add_column(year = list(tibble(year = 1:n_years,
+ c_t = c(rmvnorm(1,
+ mean = rep(0, n_years),
+ sigma = tmp_am_par$sc2[i] *
+ exp(-tmp_am_par$gamma_c[i] *
+ as.matrix(dist(1:n_years)))))))) %>%
+ unnest(cols = c(year)) %>%
+ arrange(year, species) %>%
+ mutate(s_i_t = c(rmvnorm(n_sp,
+ mean = rep(0, n_years),
+ sigma = tmp_am_par$se2[i] *
+ exp(-tmp_am_par$gamma[i] *
+ as.matrix(dist(1:n_years)))))) %>%
+ add_column(eta = tmp_am_par$eta[i])
+
+ sim <- sim %>%
+ mutate(abundance = rpois(nrow(sim), exp(eta + h_i + c_t + s_i_t)))
+
+ p0_am <- c(p0_am,
+ 1 - n_distinct(select(filter(sim, abundance > 0), species)) / n_sp)
+
+}
+
+tmp_am_par <- tmp_am_par %>%
+ add_column(p0 = p0_am) %>%
+ mutate(ES = n_distinct(data_ammarnas$species) / (1 - p0),
+ I = 1 / (z * ES))
+
+tmp_am_par_sum <- tmp_am_par %>%
+ select(se2:I) %>%
+ pivot_longer(se2:I) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+# write_rds(list(m01_am = m01_am,
+# tmp_am_par = tmp_am_par,
+# tmp_am_par_sum = tmp_am_par_sum),
+# "model_ammarnas_inla.rds")
+
+# Basic model with all the variance components
+data_budal_inla <- data_budal %>%
+ mutate(species = as.numeric(as.factor(species)),
+ year2 = year)
+
+data_budal_inla %>%
+ ggplot(aes(x = year, y = abundance)) +
+ geom_line(aes(group = species))
+
+m01_bu <- inla(abundance ~ f(year,
+ model = "ou",
+ group = species,
+ control.group = list(model = "iid")) +
+ f(species, model = "iid") +
+ f(year2,
+ model = "ou"),
+ family = "poisson",
+ data = data_budal_inla)
+
+# Precision for year and species seem ok!
+summary(m01_bu)
+
+tmp_bu_par <- inla.hyperpar.sample(sample_size, m01_bu) %>%
+ as_tibble() %>%
+ mutate(se2 = 1 / `Precision for year`,
+ gamma = `Phi for year`,
+ gamma_c = `Phi for year2`,
+ sh2 = 1 / `Precision for species`,
+ sc2 = 1 / `Precision for year2`,
+ eta = inla.rmarginal(sample_size,
+ m01_bu$marginals.fixed$`(Intercept)`),
+ rho2 = se2 + sh2 + sc2,
+ z = f_z(eta, rho2))
+
+# How many species can there be in total:
+n_years <- n_distinct(data_budal$year)
+n_sp <- 100
+
+library(mvtnorm)
+
+p0_bu = c()
+for (i in 1:sample_size){
+ sim <- tibble(species = 1:n_sp,
+ h_i = rnorm(n_sp, mean = 0,
+ sd = sqrt(tmp_bu_par$sh2[i]))) %>%
+ add_column(year = list(tibble(year = 1:n_years,
+ c_t = c(rmvnorm(1,
+ mean = rep(0, n_years),
+ sigma = tmp_bu_par$sc2[i] *
+ exp(-tmp_bu_par$gamma_c[i] *
+ as.matrix(dist(1:n_years)))))))) %>%
+ unnest(cols = c(year)) %>%
+ arrange(year, species) %>%
+ mutate(s_i_t = c(rmvnorm(n_sp,
+ mean = rep(0, n_years),
+ sigma = tmp_bu_par$se2[i] *
+ exp(-tmp_bu_par$gamma[i] *
+ as.matrix(dist(1:n_years)))))) %>%
+ add_column(eta = tmp_bu_par$eta[i])
+
+ sim <- sim %>%
+ mutate(abundance = rpois(nrow(sim), exp(eta + h_i + c_t + s_i_t)))
+
+ p0_bu <- c(p0_bu,
+ 1 - n_distinct(select(filter(sim, abundance > 0), species)) / n_sp)
+
+}
+
+tmp_bu_par <- tmp_bu_par %>%
+ add_column(p0 = p0_bu) %>%
+ mutate(ES = n_distinct(data_budal$species) / (1 - p0),
+ I = 1 / (z * ES))
+
+tmp_bu_par_sum <- tmp_bu_par %>%
+ select(se2:I) %>%
+ pivot_longer(se2:I) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+# write_rds(list(m01_bu = m01_bu,
+# tmp_bu_par = tmp_bu_par,
+# tmp_bu_par_sum = tmp_bu_par_sum),
+# "model_budal_inla.rds")
+
+# Basic model with all the variance components
+data_surrey_inla <- data_surrey %>%
+ mutate(species = as.numeric(as.factor(species)),
+ year2 = year)
+
+data_surrey_inla %>%
+ ggplot(aes(x = year, y = abundance)) +
+ geom_line(aes(group = species))
+
+m01_su <- inla(abundance ~ f(year,
+ model = "ou",
+ group = species,
+ control.group = list(model = "iid")) +
+ f(species, model = "iid") +
+ f(year2,
+ model = "ou"),
+ family = "poisson",
+ data = data_surrey_inla)
+
+# Precision for year and species are not ok, so we change the prior of species.
+summary(m01_su)
+
+m01_su <- inla(abundance ~ f(year,
+ model = "ou",
+ group = species,
+ control.group = list(model = "iid")) +
+ f(species, model = "iid",
+ hyper = list(theta = list(prior = "loggamma",
+ param = c(1, 5e-1)))) +
+ f(year2,
+ model = "ou"),
+ family = "poisson",
+ data = data_surrey_inla)
+
+summary(m01_su)
+
+tmp_su_par <- inla.hyperpar.sample(sample_size, m01_su) %>%
+ as_tibble() %>%
+ mutate(se2 = 1 / `Precision for year`,
+ gamma = `Phi for year`,
+ gamma_c = `Phi for year2`,
+ sh2 = 1 / `Precision for species`,
+ sc2 = 1 / `Precision for year2`,
+ eta = inla.rmarginal(sample_size,
+ m01_su$marginals.fixed$`(Intercept)`),
+ rho2 = se2 + sh2 + sc2,
+ z = f_z(eta, rho2))
+
+# How many species can there be in total:
+n_years <- n_distinct(data_surrey$year)
+n_sp <- 100
+
+library(mvtnorm)
+
+p0_su = c()
+for (i in 1:sample_size){
+ sim <- tibble(species = 1:n_sp,
+ h_i = rnorm(n_sp, mean = 0,
+ sd = sqrt(tmp_su_par$sh2[i]))) %>%
+ add_column(year = list(tibble(year = 1:n_years,
+ c_t = c(rmvnorm(1,
+ mean = rep(0, n_years),
+ sigma = tmp_su_par$sc2[i] *
+ exp(-tmp_su_par$gamma_c[i] *
+ as.matrix(dist(1:n_years)))))))) %>%
+ unnest(cols = c(year)) %>%
+ arrange(year, species) %>%
+ mutate(s_i_t = c(rmvnorm(n_sp,
+ mean = rep(0, n_years),
+ sigma = tmp_su_par$se2[i] *
+ exp(-tmp_su_par$gamma[i] *
+ as.matrix(dist(1:n_years)))))) %>%
+ add_column(eta = tmp_su_par$eta[i])
+
+ sim <- sim %>%
+ mutate(abundance = rpois(nrow(sim), exp(eta + h_i + c_t + s_i_t)))
+
+ p0_su <- c(p0_su,
+ 1 - n_distinct(select(filter(sim, abundance > 0), species)) / n_sp)
+
+}
+
+tmp_su_par <- tmp_su_par %>%
+ add_column(p0 = p0_su) %>%
+ mutate(ES = n_distinct(data_surrey$species) / (1 - p0),
+ I = 1 / (z * ES))
+
+tmp_su_par_sum <- tmp_su_par %>%
+ select(se2:I) %>%
+ pivot_longer(se2:I) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+# write_rds(list(m01_su = m01_su,
+# tmp_su_par = tmp_su_par,
+# tmp_su_par_sum = tmp_su_par_sum),
+# "model_surrey_inla.rds")
+
+# Basic model with all the variance components
+data_valley_inla <- data_valley %>%
+ mutate(species = as.numeric(as.factor(species)),
+ year2 = year)
+
+data_valley_inla %>%
+ ggplot(aes(x = year, y = abundance)) +
+ geom_line(aes(group = species))
+
+m01_va <- inla(abundance ~ f(year,
+ model = "ou",
+ group = species,
+ control.group = list(model = "iid")) +
+ f(species, model = "iid") +
+ f(year2,
+ model = "ou"),
+ family = "poisson",
+ data = data_valley_inla)
+
+# Precision for year and species are not ok, so we change the prior of species.
+summary(m01_va)
+
+m01_va <- inla(abundance ~ f(year,
+ model = "ou",
+ group = species,
+ control.group = list(model = "iid")) +
+ f(species, model = "iid",
+ hyper = list(theta = list(prior = "loggamma",
+ param = c(0.5, 5e-1)))) +
+ f(year2, model = "ou"),
+ family = "poisson",
+ data = data_valley_inla)
+
+summary(m01_va)
+
+tmp_va_par <- inla.hyperpar.sample(sample_size, m01_va) %>%
+ as_tibble() %>%
+ mutate(se2 = 1 / `Precision for year`,
+ gamma = `Phi for year`,
+ gamma_c = `Phi for year2`,
+ sh2 = 1 / `Precision for species`,
+ sc2 = 1 / `Precision for year2`,
+ eta = inla.rmarginal(sample_size,
+ m01_va$marginals.fixed$`(Intercept)`),
+ rho2 = se2 + sh2 + sc2,
+ z = f_z(eta, rho2))
+
+# How many species can there be in total:
+n_years <- n_distinct(data_valley$year)
+n_sp <- 100
+
+library(mvtnorm)
+
+p0_va = c()
+for (i in 1:sample_size){
+ sim <- tibble(species = 1:n_sp,
+ h_i = rnorm(n_sp, mean = 0,
+ sd = sqrt(tmp_va_par$sh2[i]))) %>%
+ add_column(year = list(tibble(year = 1:n_years,
+ c_t = c(rmvnorm(1,
+ mean = rep(0, n_years),
+ sigma = tmp_va_par$sc2[i] *
+ exp(-tmp_va_par$gamma_c[i] *
+ as.matrix(dist(1:n_years)))))))) %>%
+ unnest(cols = c(year)) %>%
+ arrange(year, species) %>%
+ mutate(s_i_t = c(rmvnorm(n_sp,
+ mean = rep(0, n_years),
+ sigma = tmp_va_par$se2[i] *
+ exp(-tmp_va_par$gamma[i] *
+ as.matrix(dist(1:n_years)))))) %>%
+ add_column(eta = tmp_va_par$eta[i])
+
+ sim <- sim %>%
+ mutate(abundance = rpois(nrow(sim), exp(eta + h_i + c_t + s_i_t)))
+
+ p0_va <- c(p0_va,
+ 1 - n_distinct(select(filter(sim, abundance > 0), species)) / n_sp)
+
+}
+
+tmp_va_par <- tmp_va_par %>%
+ add_column(p0 = p0_va) %>%
+ mutate(ES = n_distinct(data_valley$species) / (1 - p0),
+ I = 1 / (z * ES))
+
+tmp_va_par_sum <- tmp_va_par %>%
+ select(se2:I) %>%
+ pivot_longer(se2:I) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+# write_rds(list(m01_va = m01_va,
+# tmp_va_par = tmp_va_par,
+# tmp_va_par_sum = tmp_va_par_sum),
+# "model_valley_inla.rds")
diff --git a/fig1.R b/fig1.R
new file mode 100644
index 0000000000000000000000000000000000000000..a859eafc86c760c82402f0f86f02857a80ab6111
--- /dev/null
+++ b/fig1.R
@@ -0,0 +1,174 @@
+# Figure 1 ####
+
+# Required packages
+library(tidyverse)
+library(mvtnorm)
+library(patchwork)
+
+# Parameters for community in (a)
+gamma <- 0.1
+lnK <- 5
+se2 <- 0.1
+sh2 <- 0.3
+sc2 <- 1.1
+n_sp <- 30
+n_years <- 30
+
+# To get the same results
+set.seed(134)
+
+# Correlation matrix
+cor_mat <- exp(-gamma * as.matrix(dist(1:n_years)))
+
+sim01 <- tibble(species = 1:n_sp,
+ h_i = rnorm(n_sp, mean = 0, sd = sqrt(sh2))) %>%
+ add_column(year = list(tibble(year = 1:n_years,
+ c_t = c(rmvnorm(1,
+ mean = rep(0, n_years),
+ sigma = sc2 * cor_mat))))) %>%
+ unnest(cols = c(year)) %>%
+ arrange(year, species) %>%
+ mutate(s_i_t = c(rmvnorm(n_sp,
+ mean = rep(0, n_years),
+ sigma = se2 * cor_mat))) %>%
+ add_column(eta = lnK)
+
+# Parameters in community (b)
+lnK <- 5
+se2 <- 0.6
+sh2 <- 0.3
+sc2 <- 0.6
+
+sim02 <- tibble(species = 1:n_sp,
+ h_i = rnorm(n_sp, mean = 0, sd = sqrt(sh2))) %>%
+ add_column(year = list(tibble(year = 1:n_years,
+ c_t = c(rmvnorm(1,
+ mean = rep(0, n_years),
+ sigma = sc2 * cor_mat))))) %>%
+ unnest(cols = c(year)) %>%
+ arrange(year, species) %>%
+ mutate(s_i_t = c(rmvnorm(n_sp,
+ mean = rep(0, n_years),
+ sigma = se2 * cor_mat))) %>%
+ add_column(eta = lnK)
+
+# Parameters in community (c)
+lnK <- 5
+se2 <- 1.1
+sh2 <- 0.3
+sc2 <- 0.1
+
+sim03 <- tibble(species = 1:n_sp,
+ h_i = rnorm(n_sp, mean = 0, sd = sqrt(sh2))) %>%
+ add_column(year = list(tibble(year = 1:n_years,
+ c_t = c(rmvnorm(1,
+ mean = rep(0, n_years),
+ sigma = sc2 * cor_mat))))) %>%
+ unnest(cols = c(year)) %>%
+ arrange(year, species) %>%
+ mutate(s_i_t = c(rmvnorm(n_sp,
+ mean = rep(0, n_years),
+ sigma = se2 * cor_mat))) %>%
+ add_column(eta = lnK)
+
+fig_data <- bind_rows(add_column(sim01, mod = "common"),
+ add_column(sim02, mod = "mix"),
+ add_column(sim03, mod = "species"))
+
+fig_data <- fig_data %>%
+ mutate(mod = factor(mod, levels = c("species", "common", "mix")))
+
+
+fig_data_a <- fig_data %>%
+ filter(year == 5)
+
+fig_data_b <- fig_data %>%
+ filter(year == 25)
+
+# Illustrate temporal dynamics of communities
+
+fig1a <- fig_data %>%
+ ggplot(aes(x = year, y = eta + h_i + c_t + s_i_t)) +
+ geom_line(data = filter(fig_data,
+ species != 1),
+ aes(group = species), colour = "light grey") +
+ geom_line(data = filter(fig_data,
+ species == 1),
+ aes(group = species), colour = "black") +
+ # geom_point(data = fig_data_a, shape = 21, size = rel(2), fill = "white") +
+ # geom_point(data = fig_data_b, size = rel(2)) +
+ geom_text(data = tibble(mod = factor(c("species", "common", "mix"),
+ levels = c("species", "common", "mix")),
+ x = rep(2, 3),
+ y = rep(7.7, 3),
+ label = c("(a)", "(b)", "(c)")),
+ aes(x = x, y = y, label = label), size = rel(5)) +
+ facet_wrap(~ mod, nrow = 1) +
+ labs(x = "Time", y = "Log abundance") +
+ theme_bw() +
+ theme(strip.background = element_blank(),
+ strip.text = element_blank(),
+ axis.title = element_text(size = rel(1.2)),
+ axis.text = element_text(size = rel(1.2)))
+
+fig_data_c <- fig_data %>%
+ filter(year %in% c(5, 25))
+
+# Illustrate log abundance distribution
+
+fig1b <- fig_data_c %>%
+ ggplot(aes(x = eta + h_i + c_t + s_i_t)) +
+ geom_histogram(aes(fill = as.factor(year)),
+ position = "dodge", colour = "black") +
+ geom_text(data = tibble(mod = factor(c("species", "common", "mix"),
+ levels = c("species", "common", "mix")),
+ x = rep(2.1, 3),
+ y = rep(5.8, 3),
+ label = c("(d)", "(e)", "(f)")),
+ aes(x = x, y = y, label = label), size = rel(5)) +
+ facet_wrap(~ mod, nrow = 1) +
+ scale_fill_manual("Time", values = c("white", "black")) +
+ labs(x = "Log abundance", y = "Count") +
+ theme_bw() +
+ theme(legend.position = c(.23, .99),
+ legend.justification = c("left", "top"),
+ legend.box.just = "right",
+ legend.margin = margin(6, 6, 6, 6),
+ legend.text = element_text(size = rel(1.2)),
+ legend.title = element_text(size = rel(1.2)),
+ strip.background = element_blank(),
+ strip.text = element_blank(),
+ axis.title = element_text(size = rel(1.2)),
+ axis.text = element_text(size = rel(1.2)))
+
+# Illustrate covariation between two time points
+
+fig1c <- fig_data %>%
+ mutate(log_abundance = eta + h_i + c_t + s_i_t) %>%
+ select(-eta, -h_i, -c_t, -s_i_t) %>%
+ pivot_wider(values_from = log_abundance,
+ names_from = year) %>%
+ ggplot(aes(x = `5`,
+ y = `25`)) +
+ geom_abline(slope = 1, intercept = 0, col = "grey") +
+ geom_point(size = rel(2), shape = 21, fill = "grey") +
+ geom_text(data = tibble(mod = factor(c("species", "common", "mix"),
+ levels = c("species", "common", "mix")),
+ x = rep(2.1, 3),
+ y = rep(7.0, 3),
+ label = c("(g)", "(h)", "(i)")),
+ aes(x = x, y = y, label = label), size = rel(5)) +
+ labs(x = "Log abundance at time 5",
+ y = "Log abundance at time 25") +
+ facet_wrap(~ mod, nrow = 1) +
+ theme_bw() +
+ theme(strip.background = element_blank(),
+ strip.text = element_blank(),
+ axis.title = element_text(size = rel(1.2)),
+ axis.text = element_text(size = rel(1.2)))
+
+# Join the figures together
+
+fig1a + fig1b + fig1c + plot_layout(nrow = 3)
+
+# ggsave("Fig1.pdf", height = 8, width = 9)
diff --git a/fig2.R b/fig2.R
new file mode 100644
index 0000000000000000000000000000000000000000..1924940b89496b63b5a1795e370162e2338d7c68
--- /dev/null
+++ b/fig2.R
@@ -0,0 +1,138 @@
+# Figure 2 ####
+
+# Required packages
+library(tidyverse)
+library(patchwork)
+
+# Required functions
+
+f_G <- function(x){
+ x * dnorm(x) / pnorm(x)
+}
+
+f_z <- function(eta, rho2){
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2))) /
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2) + sqrt(rho2)) + 2 * eta + rho2)
+}
+
+
+f_EN <- function(eta, rho2, ES){
+ prop <- pnorm(0, mean = eta,
+ sd = sqrt(rho2), lower.tail = FALSE)
+ EN <- integrate(function(x){exp(x) * ES * dnorm(x,
+ mean = eta,
+ sd = sqrt(rho2)) / prop},
+ lower = 0, upper = 100)$value
+ return(EN)
+}
+
+# Calculate z for each simulation ####
+
+fig2a_data <- tibble(case = c("A", "B", "C"),
+ mu = rep(5, 3),
+ sd = sqrt(c(0.4, 0.9, 1.4))) %>%
+ mutate(x = list(seq(1, 11, length.out = 1000))) %>%
+ unnest(x) %>%
+ mutate(ds = dnorm(x, mean = mu, sd = sd))
+
+fig2a <- fig2a_data %>%
+ ggplot(aes(x = x, y = ds)) +
+ geom_line(aes(colour = case), linewidth = 1) +
+ theme_bw() +
+ scale_x_continuous(expand = c(0, 0)) +
+ scale_y_continuous(expand = c(0, 0), limits = c(0, 0.65)) +
+ scale_color_discrete(name = "Total variance:",
+ breaks = c("A", "B", "C"),
+ labels = c("0.4", "0.9", "1.4")) +
+ labs(x = "Log abundance, x",
+ y = "Density") +
+ theme(axis.title = element_text(size = rel(1.2)),
+ axis.text = element_text(size = rel(1.2)),
+ legend.title = element_text(size = rel(1.2)),
+ legend.text = element_text(size = rel(1.2)),
+ legend.position = c(0.75, 0.75),
+ plot.margin = margin(-5.5, 33, 5.5, 5.5, "pt"))
+
+fig2b_data <- tibble(case = c("A", "B", "C"),
+ mu = rep(5, 3),
+ sd = sqrt(c(0.4, 0.9, 1.4)),
+ z = f_z(mu, sd ^ 2)) %>%
+ mutate(dlnEN = list(seq(0, 4, length.out = 100))) %>%
+ unnest(dlnEN) %>%
+ mutate(dlnES = dlnEN * z)
+
+fig2b <- fig2b_data %>%
+ ggplot(aes(x = dlnEN, y = dlnES)) +
+ geom_line(aes(colour = case), linewidth = 1) +
+ theme_bw() +
+ scale_x_continuous(expand = c(0, 0)) +
+ scale_y_continuous(expand = c(0, 0), limits = c(0, 4)) +
+ labs(x = "Change in log total abundance, dln EN",
+ y = "Change in log number of species, dln ES") +
+ guides(colour = "none") +
+ theme(axis.title = element_text(size = rel(1.2)),
+ axis.text = element_text(size = rel(1.2)),
+ plot.margin = margin(-5.5, 33, 5.5, 5.5, "pt"))
+
+
+fig2c_data <- tibble(case = c("A", "B", "C"),
+ eta = rep(5, 3),
+ ES = c(20, 30, 50)) %>%
+ mutate(rho2 = list(seq(0, 10, length.out = 1000))) %>%
+ unnest(rho2) %>%
+ rowwise() %>%
+ mutate(z = f_z(eta, rho2),
+ EN = f_EN(eta, rho2, ES),
+ I = 1 / (z * ES))
+
+fig2c <- fig2c_data %>%
+ ggplot(aes(x = log(EN), y = z)) +
+ geom_line(aes(linetype = case), linewidth = 1) +
+ geom_point(data = filter(fig2c_data,
+ round(rho2, digits = 2) %in% c(0.61, 1.00, 1.40)),
+ aes(x = log(EN), y = z, colour = as.factor(rho2)),
+ size = rel(3)) +
+ theme_bw() +
+ scale_x_continuous(expand = c(0, 0), limits = c(8, 12)) +
+ scale_y_continuous(expand = c(0, 0), limits = c(0, 1)) +
+ scale_linetype_manual(name = "ES",
+ breaks = c("A", "B", "C"),
+ labels = c("20", "30", "50"),
+ values = c(2, 1, 3)) +
+ labs(x = "Log expected total community abundance, ln EN",
+ y = "Community sensitivity, z") +
+ guides(colour = "none") +
+ theme(axis.title = element_text(size = rel(1.2)),
+ axis.text = element_text(size = rel(1.2)),
+ legend.title = element_text(size = rel(1.2)),
+ legend.text = element_text(size = rel(1.2)),
+ legend.position = c(0.75, 0.75),
+ plot.margin = margin(-5.5, 33, 5.5, 5.5, "pt"))
+
+fig2d <- fig2c_data %>%
+ ggplot(aes(x = log(EN), y = I)) +
+ geom_line(aes(linetype = case), linewidth = 1) +
+ geom_point(data = filter(fig2c_data,
+ round(rho2, digits = 2) %in% c(0.61, 1.00, 1.40)),
+ aes(x = log(EN), y = I, colour = as.factor(rho2)),
+ size = rel(3)) +
+ theme_bw() +
+ scale_x_continuous(expand = c(0, 0), limits = c(8, 12)) +
+ scale_y_continuous(expand = c(0, 0), limits = c(0, 0.4)) +
+ scale_linetype_manual(name = "ES",
+ breaks = c("A", "B", "C"),
+ labels = c("20", "30", "50"),
+ values = c(2, 1, 3)) +
+ labs(x = "Log expected total community abundance, ln EN",
+ y = "Community resistance, I") +
+ guides(linetype = "none",
+ colour = "none") +
+ theme(axis.title = element_text(size = rel(1.2)),
+ axis.text = element_text(size = rel(1.2)),
+ plot.margin = margin(-5.5, 33, 5.5, 5.5, "pt"))
+
+fig2a + fig2b + fig2c + fig2d +
+ plot_layout(ncol = 2) +
+ plot_annotation(tag_levels = "a")
+
+# ggsave("Fig2.pdf", height = 8, width = 9)
diff --git a/fig3.R b/fig3.R
new file mode 100644
index 0000000000000000000000000000000000000000..6fbd94a3c400946bd0542339d91876fdee2b93d1
--- /dev/null
+++ b/fig3.R
@@ -0,0 +1,256 @@
+# Figure 3 ####
+
+# Required packages
+library(tidyverse)
+library(HDInterval)
+library(modeest)
+library(patchwork)
+
+# Required functions
+
+f_G <- function(x){
+ x * dnorm(x) / pnorm(x)
+}
+
+f_z <- function(eta, rho2){
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2))) /
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2) + sqrt(rho2)) + 2 * eta + rho2)
+}
+
+f_sum <- function(x, p = 0.95){
+ x <- x$value
+ tmp <- c(mean(x), sd(x), median(x), meanshift(x), hdi(x, p),
+ quantile(x, probs = c((1 - p)/2, 1 - (1 - p) / 2)))
+ names(tmp)[1:4] <- c("mean", "sd", "median", "mode")
+ return(as_tibble_row(tmp))
+}
+
+# Estimating expected total abundance
+
+# We assume that lambda(x) is approximately normal and use the estimated
+# expected number of species to estimate the total abundance:
+
+f_EN <- function(eta, rho2, ES){
+ prop <- pnorm(0, mean = eta,
+ sd = sqrt(rho2), lower.tail = FALSE)
+ EN <- integrate(function(x){exp(x) * ES * dnorm(x,
+ mean = eta,
+ sd = sqrt(rho2)) / prop},
+ lower = 0, upper = 100)$value
+ return(EN)
+}
+
+model_ammarnas <- read_rds("model_ammarnas_inla.rds")
+tmp_am_par <- model_ammarnas$tmp_am_par
+tmp_am_par_sum <- model_ammarnas$tmp_am_par_sum
+
+tmp_am_par <- tmp_am_par %>%
+ rowwise() %>%
+ mutate(EN = f_EN(eta, rho2, ES),
+ logEN = log(EN),
+ lnK = log(EN / ES),
+ r = lnK * gamma) %>%
+ ungroup()
+
+tmp_am_par_sum <- tmp_am_par %>%
+ select(se2:r) %>%
+ pivot_longer(se2:r) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+model_budal <- read_rds("model_budal_inla.rds")
+tmp_bu_par <- model_budal$tmp_bu_par
+tmp_bu_par_sum <- model_budal$tmp_bu_par_sum
+tmp_bu_par <- tmp_bu_par %>%
+ rowwise() %>%
+ mutate(EN = f_EN(eta, rho2, ES),
+ logEN = log(EN),
+ lnK = log(EN / ES),
+ r = lnK * gamma) %>%
+ ungroup()
+
+tmp_bu_par_sum <- tmp_bu_par %>%
+ select(se2:r) %>%
+ pivot_longer(se2:r) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+model_surrey <- read_rds("model_surrey_inla.rds")
+tmp_su_par <- model_surrey$tmp_su_par
+tmp_su_par_sum <- model_surrey$tmp_su_par_sum
+tmp_su_par <- tmp_su_par %>%
+ rowwise() %>%
+ mutate(EN = f_EN(eta, rho2, ES),
+ logEN = log(EN),
+ lnK = log(EN / ES),
+ r = lnK * gamma) %>%
+ ungroup()
+
+tmp_su_par_sum <- tmp_su_par %>%
+ select(se2:r) %>%
+ pivot_longer(se2:r) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+model_valley <- read_rds("model_valley_inla.rds")
+tmp_va_par <- model_valley$tmp_va_par
+tmp_va_par_sum <- model_valley$tmp_va_par_sum
+tmp_va_par <- tmp_va_par %>%
+ rowwise() %>%
+ mutate(EN = f_EN(eta, rho2, ES),
+ logEN = log(EN),
+ lnK = log(EN / ES),
+ r = lnK * gamma) %>%
+ ungroup()
+
+tmp_va_par_sum <- tmp_va_par %>%
+ select(se2:r) %>%
+ pivot_longer(se2:r) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+locations <- c("Ammarnäs", "Budal", "Eastern Wood", "Birdsong Valley")
+
+fig3_data <- bind_rows(add_column(tmp_am_par,
+ Location = "Ammarnäs"),
+ add_column(tmp_bu_par,
+ Location = "Budal"),
+ add_column(tmp_su_par,
+ Location = "Eastern Wood"),
+ add_column(tmp_va_par,
+ Location = "Birdsong Valley")) %>%
+ # Select posterior samples that are within the upper and lower highest
+ # density interval for community sensitivity (z) and resistance (I)
+ filter((Location == "Ammarnäs" &
+ I > tmp_am_par_sum$lower[11] &
+ I < tmp_am_par_sum$upper[11] &
+ z > tmp_am_par_sum$lower[8] &
+ z < tmp_am_par_sum$upper[8]) |
+ (Location == "Budal" &
+ I > tmp_bu_par_sum$lower[11] &
+ I < tmp_bu_par_sum$upper[11] &
+ z > tmp_bu_par_sum$lower[8] &
+ z < tmp_bu_par_sum$upper[8]) |
+ (Location == "Eastern Wood" &
+ I > tmp_su_par_sum$lower[11] &
+ I < tmp_su_par_sum$upper[11] &
+ z > tmp_su_par_sum$lower[8] &
+ z < tmp_su_par_sum$upper[8]) |
+ (Location == "Birdsong Valley" &
+ I > tmp_va_par_sum$lower[11] &
+ I < tmp_va_par_sum$upper[11] &
+ z > tmp_va_par_sum$lower[8] &
+ z < tmp_va_par_sum$upper[8])) %>%
+ mutate(Location = factor(Location, levels = locations))
+
+asterisk_data <- bind_rows(add_column(tmp_am_par_sum,
+ Location = "Ammarnäs"),
+ add_column(tmp_bu_par_sum,
+ Location = "Budal"),
+ add_column(tmp_su_par_sum,
+ Location = "Eastern Wood"),
+ add_column(tmp_va_par_sum,
+ Location = "Birdsong Valley")) %>%
+ mutate(Location = factor(Location, levels = locations))
+
+asterisk_z_data <- asterisk_data %>%
+ filter(name %in% c("z", "logEN")) %>%
+ select(name, mean, Location) %>%
+ pivot_wider(names_from = name, values_from = mean) %>%
+ pivot_longer(c(z))
+
+asterisk_I_data <- asterisk_data %>%
+ filter(name %in% c("I", "logEN")) %>%
+ select(name, mean, Location) %>%
+ pivot_wider(names_from = name, values_from = mean) %>%
+ pivot_longer(c(I))
+
+fig3a <- fig3_data %>%
+ select(Location, EN, z) %>%
+ pivot_longer(z) %>%
+ ggplot(aes(x = log(EN), y = value)) +
+ geom_point(aes(colour = Location), size = rel(2)) +
+ geom_point(data = asterisk_z_data,
+ aes(x = logEN, y = value), size = rel(4),
+ shape = "asterisk", stroke = 2) +
+ geom_text(data = tibble(Location = factor(locations, levels = locations),
+ logEN = rep(12, 4),
+ z = rep(0.33, 4),
+ label = c("(a)", "(b)", "(c)", "(d)")),
+ aes(x = logEN, y = z, label = label), size = rel(5)) +
+ facet_wrap(~ Location, ncol = 4) +
+ guides(colour = "none") +
+ labs(y = "Community sensitivity, z",
+ x = "ln EN") +
+ theme_bw() +
+ theme(
+ legend.position = c(.99, .99),
+ legend.justification = c("right", "top"),
+ legend.box.just = "right",
+ legend.margin = margin(6, 6, 6, 6),
+ strip.background = element_rect(fill = "white"),
+ strip.text = element_text(size = rel(1.1)),
+ axis.title.y = element_text(size = rel(1.2)),
+ axis.text.y = element_text(size = rel(1.2)),
+ axis.title.x = element_blank(),
+ axis.text.x = element_blank()
+ )
+
+
+fig3b <- fig3_data %>%
+ select(Location, EN, I) %>%
+ pivot_longer(I) %>%
+ ggplot(aes(x = log(EN), y = value)) +
+ geom_point(aes(colour = Location), size = rel(2)) +
+ geom_point(data = asterisk_I_data,
+ aes(x = logEN, y = value), size = rel(4),
+ shape = "asterisk", stroke = 2) +
+ geom_text(data = tibble(Location = factor(locations, levels = locations),
+ logEN = rep(5, 4),
+ I = rep(0.7, 4),
+ label = c("(e)", "(f)", "(g)", "(h)")),
+ aes(x = logEN, y = I, label = label), size = rel(5)) +
+ facet_wrap(~ Location, ncol = 4) +
+ guides(colour = "none") +
+ labs(y = "Community resistance, I",
+ x = "Log expected total community abundance, ln EN") +
+ theme_bw() +
+ theme(
+ legend.position = c(.99, .99),
+ legend.justification = c("right", "top"),
+ legend.box.just = "right",
+ legend.margin = margin(6, 6, 6, 6),
+ strip.background = element_rect(fill = "white"),
+ strip.text = element_blank(),
+ axis.title.y = element_text(size = rel(1.2)),
+ axis.text.y = element_text(size = rel(1.2)),
+ axis.title.x = element_text(size = rel(1.2)),
+ axis.text.x = element_text(size = rel(1.2))
+ )
+
+fig3a + fig3b + plot_layout(nrow = 2)
+
+# ggsave("Fig3.pdf", height = 8, width = 9)
+
+# Tables for SI:
+
+xtable::xtable(tmp_am_par_sum)
+xtable::xtable(tmp_bu_par_sum)
+xtable::xtable(tmp_su_par_sum)
+xtable::xtable(tmp_va_par_sum)
diff --git a/fig4.R b/fig4.R
new file mode 100644
index 0000000000000000000000000000000000000000..c43b5e93335444f151c6586e63dc08df4cc6d1a4
--- /dev/null
+++ b/fig4.R
@@ -0,0 +1,94 @@
+# Figure 4 ####
+
+# Required packages
+library(tidyverse)
+library(patchwork)
+
+source("load_data.R")
+
+f_results <- function(location = "am", type = "single"){
+ tmp <- read_rds(paste(switch(type,
+ "single" = "res_single_",
+ "inc" = "res_inc_year_"), location, ".rds", sep = ""))
+ count <- 0
+ tmp_data <- get(switch(location,
+ "am" = "data_ammarnas",
+ "bu" = "data_budal",
+ "su" = "data_surrey",
+ "va" = "data_valley"))
+ tmp_res <- c()
+
+ range <- ifelse(type == "single", 1, 2):length(unique(tmp_data$year))
+ for (k in unique(tmp_data$year)[range]){
+ count <- count + 1
+ tmp_res <- tmp_res %>%
+ bind_rows(add_column(tmp[[count]][[3]], year = k))
+ }
+
+ tmp_res <- tmp_res %>%
+ add_column(Location = switch(location,
+ "am" = "Ammarnäs",
+ "bu" = "Budal",
+ "su" = "Eastern Wood",
+ "va" = "Birdsong Valley"))
+
+ return(tmp_res)
+}
+
+res_single_am_sum <- f_results(location = "am")
+res_single_bu_sum <- f_results(location = "bu")
+res_single_su_sum <- f_results(location = "su")
+res_single_va_sum <- f_results(location = "va")
+res_inc_year_am_sum <- f_results(location = "am", type = "inc")
+res_inc_year_bu_sum <- f_results(location = "bu", type = "inc")
+res_inc_year_su_sum <- f_results(location = "su", type = "inc")
+res_inc_year_va_sum <- f_results(location = "va", type = "inc")
+
+# How does the increasing number of years affect the estimates?
+
+tmp_data <- bind_rows(res_single_am_sum,
+ res_single_bu_sum,
+ res_single_su_sum,
+ res_single_va_sum) %>%
+ add_column(data = "Single year") %>%
+ bind_rows(add_column(bind_rows(res_inc_year_am_sum,
+ res_inc_year_bu_sum,
+ res_inc_year_su_sum,
+ res_inc_year_va_sum),
+ data = "Increasing years"))
+
+tmp_data <- tmp_data %>%
+ mutate(Location = factor(Location,
+ levels = c("Ammarnäs",
+ "Budal",
+ "Eastern Wood",
+ "Birdsong Valley")))
+
+tbl_name <- tibble(name = c("rho2", "z", "ES", "I"),
+ new_name = factor(c("sigma^2", "z", "ES", "I"),
+ levels = c("sigma^2", "z", "ES", "I")))
+
+tmp_data %>%
+ filter(name %in% c("rho2", "z", "ES", "I")) %>%
+ left_join(tbl_name) %>%
+ filter(mean < 150) %>%
+ ggplot(aes(x = year, y = mean)) +
+ geom_point(aes(colour = data)) +
+ geom_line(aes(colour = data)) +
+ facet_grid(new_name ~ Location,
+ scales = "free",
+ labeller = labeller(new_name = label_parsed)) +
+ scale_color_manual(values = c("orange", "skyblue")) +
+ guides(colour = "none") +
+ theme_bw() +
+ scale_x_continuous(name = "Year", expand = c(0.05,0.05)) +
+ scale_y_continuous(name = "Parameter estimate", expand = c(0.1,0.1)) +
+ theme(legend.position = c(.99, .99),
+ legend.justification = c("right", "top"),
+ legend.box.just = "right",
+ legend.margin = margin(6, 6, 6, 6),
+ strip.background = element_rect(fill = "white"),
+ strip.text = element_text(size = rel(1.2)),
+ axis.title = element_text(size = rel(1.2)),
+ axis.text = element_text(size = rel(1.2))
+ )
diff --git a/fig5.R b/fig5.R
new file mode 100644
index 0000000000000000000000000000000000000000..3413b96dcace8098d3a4c7b62b729c2845e6295a
--- /dev/null
+++ b/fig5.R
@@ -0,0 +1,101 @@
+# Figure 5 ####
+
+# Required packages
+library(tidyverse)
+library(patchwork)
+
+# Required functions
+f_G <- function(x){
+ x * dnorm(x) / pnorm(x)
+}
+
+f_z <- function(eta, rho2){
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2))) /
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2) + sqrt(rho2)) + 2 * eta + rho2)
+}
+
+# Load results
+model_ammarnas <- read_rds("model_ammarnas_inla.rds")
+model_budal <- read_rds("model_budal_inla.rds")
+model_surrey <- read_rds("model_surrey_inla.rds")
+model_valley <- read_rds("model_valley_inla.rds")
+tmp_am_par_sum <- model_ammarnas$tmp_am_par_sum
+tmp_bu_par_sum <- model_budal$tmp_bu_par_sum
+tmp_su_par_sum <- model_surrey$tmp_su_par_sum
+tmp_va_par_sum <- model_valley$tmp_va_par_sum
+
+fig5_data <- add_column(tmp_am_par_sum, Location = "Ammarnäs") %>%
+ bind_rows(add_column(tmp_bu_par_sum, Location = "Budal")) %>%
+ bind_rows(add_column(tmp_su_par_sum, Location = "Eastern Wood")) %>%
+ bind_rows(add_column(tmp_va_par_sum, Location = "Birdsong Valley")) %>%
+ mutate(Location = factor(Location,
+ levels = c("Ammarnäs",
+ "Budal",
+ "Eastern Wood",
+ "Birdsong Valley"))) %>%
+ filter(name %in% c("eta", "rho2", "ES")) %>%
+ select(name, mean, Location) %>%
+ pivot_wider(names_from = name,
+ values_from = mean) %>%
+ mutate(nlog10nu = list(seq(0, 3, length.out = 100))) %>%
+ unnest(cols = c(nlog10nu)) %>%
+ rowwise() %>%
+ mutate(nu = 10 ^ (-nlog10nu),
+ z = f_z(log(exp(eta) * nu), rho2),
+ I = 1 / (ES * z))
+
+fig5a <- fig5_data %>%
+ ggplot(aes(x = nlog10nu, y = z)) +
+ geom_line(aes(colour = Location), linewidth = 1) +
+ labs(y = "Community sensitivity, z",
+ x = "Sampling intensity (proportion sampled)") +
+ theme_bw() +
+ scale_x_continuous(expand = c(0.01, 0.01),
+ breaks = c(0, 1, 2, 3),
+ labels = c("1", "0.1", "0.01", "0.001")) +
+ scale_y_continuous(limits = c(0, 0.6), expand = c(0.01, 0.01)) +
+ theme(
+ legend.position = c(.01, .99),
+ legend.justification = c("left", "top"),
+ legend.box.just = "right",
+ legend.margin = margin(6, 6, 6, 6),
+ legend.title = element_text(size = rel(1.2)),
+ legend.text = element_text(size = rel(1.2)),
+ strip.background = element_rect(fill = "white"),
+ strip.text = element_blank(),
+ axis.title.y = element_text(size = rel(1.2)),
+ axis.text.y = element_text(size = rel(1.2)),
+ axis.title.x = element_text(size = rel(1.2)),
+ axis.text.x = element_text(size = rel(1.2)),
+ plot.margin = margin(5.5, 10, 5.5, 5.5)
+ )
+
+fig5b <- fig5_data %>%
+ ggplot(aes(x = nlog10nu, y = I)) +
+ geom_line(aes(colour = Location), linewidth = 1) +
+ guides(colour = "none") +
+ labs(y = "Community resistance, I",
+ x = "Sampling intensity (proportion sampled)") +
+ theme_bw() +
+ scale_x_continuous(expand = c(0.01, 0.01),
+ breaks = c(0, 1, 2, 3),
+ labels = c("1", "0.1", "0.01", "0.001")) +
+ scale_y_continuous(limits = c(0, 0.6), expand = c(0.01, 0.01)) +
+ theme(
+ legend.position = c(.01, .99),
+ legend.justification = c("left", "top"),
+ legend.box.just = "right",
+ legend.margin = margin(6, 6, 6, 6),
+ legend.title = element_text(size = rel(1.2)),
+ legend.text = element_text(size = rel(1.2)),
+ strip.background = element_rect(fill = "white"),
+ strip.text = element_blank(),
+ axis.title.y = element_text(size = rel(1.2)),
+ axis.text.y = element_text(size = rel(1.2)),
+ axis.title.x = element_text(size = rel(1.2)),
+ axis.text.x = element_text(size = rel(1.2)),
+ plot.margin = margin(5.5, 10, 5.5, 5.5)
+ )
+
+fig5a + fig5b + plot_layout(nrow = 1) + plot_annotation(tag_levels = "a")
+# ggsave("Fig5.pdf", height = 6, width = 8, )
diff --git a/load_data.R b/load_data.R
new file mode 100644
index 0000000000000000000000000000000000000000..9862e331994381423dc28a5cf957ff692c0d7b85
--- /dev/null
+++ b/load_data.R
@@ -0,0 +1,111 @@
+library(tidyverse)
+
+# Data preparation given the initial data files, stored in a folder "datasets"
+
+# Ammarnas ####
+
+bird_ammarnas <- read_delim("datasets/bird_ammarnas.csv",
+ delim = ";",
+ escape_double = FALSE,
+ trim_ws = TRUE,
+ show_col_types = FALSE)
+
+bird_ammarnas_wide <- bird_ammarnas %>%
+ pivot_wider(names_from = Species,
+ values_from = Abundance,
+ values_fill = 0)
+
+bird_ammarnas_long <- bird_ammarnas_wide %>%
+ pivot_longer(cols = c(-Year),
+ names_to = "species",
+ values_to = "abundance") %>%
+ rename(year = Year) %>%
+ add_column(omr_id = "Ammarnas")
+
+data_ammarnas <- bird_ammarnas_long %>%
+ mutate(year = year - min(year))
+
+# Budal ####
+
+bird_budal <- readxl::read_excel("datasets/Budalen_pr_2018_revised.xlsx")
+
+bird_budal_long <- bird_budal %>%
+ dplyr::select(year:buskskv) %>%
+ filter(!(year %in% 1998:1999)) %>%
+ pivot_longer(cols = c(-year),
+ names_to = "species",
+ values_to = "abundance") %>%
+ mutate(abundance = round(as.numeric(ifelse(abundance == "NA",
+ "0",
+ abundance))))
+
+bird_budal_wide <- bird_budal_long %>%
+ pivot_wider(names_from = "species",
+ values_from = "abundance",
+ values_fill = 0)
+
+bird_budal_long <- bird_budal_wide %>%
+ pivot_longer(cols = c(-year),
+ names_to = "species",
+ values_to = "abundance") %>%
+ add_column(omr_id = "Budal")
+
+data_budal <- bird_budal_long %>%
+ mutate(year = year - min(year))
+
+# Eastern Wood (sorry about the wrong names!) ####
+
+bird_surrey <- read_delim("datasets/bird_sussex.txt", delim = ";",
+ escape_double = FALSE,
+ trim_ws = TRUE,
+ show_col_types = FALSE)
+
+bird_surrey_wide <- bird_surrey %>%
+ filter(Year > 1949) %>%
+ pivot_wider(names_from = Species,
+ values_from = Abundance,
+ values_fill = 0)
+
+bird_surrey_long <- bird_surrey_wide %>%
+ pivot_longer(cols = c(-Year),
+ names_to = "species",
+ values_to = "abundance") %>%
+ rename(year = Year) %>%
+ add_column(omr_id = "Surrey") %>%
+ mutate(abundance = round(abundance))
+
+data_surrey <- bird_surrey_long %>%
+ mutate(year = year - min(year))
+
+# Birdsong Valley ####
+
+bird_valley <- readxl::read_excel("datasets/Birdsong Valley.xlsx",
+ col_names = c("species",
+ "year",
+ "abundance",
+ "skip_1",
+ "skip_2",
+ "skip_3"))
+
+bird_valley_wide <- bird_valley %>%
+ dplyr::select(year, species, abundance) %>%
+ mutate(abundance = ifelse(is.na(abundance), 0, abundance)) %>%
+ pivot_wider(names_from = "species",
+ values_from = "abundance",
+ values_fill = 0)
+
+bird_valley_long <- bird_valley_wide %>%
+ pivot_longer(cols = c(-year),
+ names_to = "species",
+ values_to = "abundance") %>%
+ add_column(omr_id = "Birdsong Valley")
+
+data_valley <- bird_valley_long %>%
+ mutate(year = year - min(year))
+
+# Clean-up ####
+
+rm(list = c("bird_budal", "bird_budal_long", "bird_budal_wide",
+ "bird_surrey", "bird_surrey_long", "bird_surrey_wide",
+ "bird_ammarnas", "bird_ammarnas_long", "bird_ammarnas_wide",
+ "bird_valley", "bird_valley_long", "bird_valley_wide"))
\ No newline at end of file
diff --git a/model_ammarnas_inla.rds b/model_ammarnas_inla.rds
new file mode 100644
index 0000000000000000000000000000000000000000..ef481f5f5d1fe1c287e9ad197a92beb2fbb17c43
Binary files /dev/null and b/model_ammarnas_inla.rds differ
diff --git a/model_budal_inla.rds b/model_budal_inla.rds
new file mode 100644
index 0000000000000000000000000000000000000000..eac5b930b601dcb579bbd701b3607be42db57c66
Binary files /dev/null and b/model_budal_inla.rds differ
diff --git a/model_surrey_inla.rds b/model_surrey_inla.rds
new file mode 100644
index 0000000000000000000000000000000000000000..245393b7b591164010dc2583333cb736d2ce422d
Binary files /dev/null and b/model_surrey_inla.rds differ
diff --git a/model_valley_inla.rds b/model_valley_inla.rds
new file mode 100644
index 0000000000000000000000000000000000000000..a7f889ce453f6930749035fa580b3f572b76e364
Binary files /dev/null and b/model_valley_inla.rds differ
diff --git a/tables.R b/tables.R
new file mode 100644
index 0000000000000000000000000000000000000000..cd7173a7e93a17ce2c22593dfbfd6133158d3223
--- /dev/null
+++ b/tables.R
@@ -0,0 +1,162 @@
+# Tables ####
+
+# Required packages
+library(tidyverse)
+library(HDInterval)
+library(modeest)
+
+# Required functions
+
+f_G <- function(x){
+ x * dnorm(x) / pnorm(x)
+}
+
+f_z <- function(eta, rho2){
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2))) /
+ (1 + eta ^ 2 / rho2 + f_G(eta / sqrt(rho2) + sqrt(rho2)) + 2 * eta + rho2)
+}
+
+f_sum <- function(x, p = 0.95){
+ x <- x$value
+ tmp <- c(mean(x), sd(x), median(x), meanshift(x), hdi(x, p),
+ quantile(x, probs = c((1 - p)/2, 1 - (1 - p) / 2)))
+ names(tmp)[1:4] <- c("mean", "sd", "median", "mode")
+ return(as_tibble_row(tmp))
+}
+
+# Estimating expected total abundance
+
+# We assume that lambda(x) is approximately normal and use the estimated
+# expected number of species to estimate the total abundance:
+
+f_EN <- function(eta, rho2, ES){
+ prop <- pnorm(0, mean = eta,
+ sd = sqrt(rho2), lower.tail = FALSE)
+ EN <- integrate(function(x){exp(x) * ES * dnorm(x,
+ mean = eta,
+ sd = sqrt(rho2)) / prop},
+ lower = 0, upper = 100)$value
+ return(EN)
+}
+
+model_ammarnas <- read_rds("model_ammarnas_inla.rds")
+tmp_am_par <- model_ammarnas$tmp_am_par
+tmp_am_par_sum <- model_ammarnas$tmp_am_par_sum
+
+tmp_am_par <- tmp_am_par %>%
+ rowwise() %>%
+ mutate(EN = f_EN(eta, rho2, ES),
+ ss2 = se2 * 2 * gamma,
+ sE2 = sc2 * 2 * gamma_c,
+ sr2 = sh2 * gamma ^ 2,
+ logEN = log(EN),
+ lnK = log(EN / ES),
+ r = lnK * gamma) %>%
+ ungroup()
+
+tmp_am_par_sum <- tmp_am_par %>%
+ select(se2:r) %>%
+ pivot_longer(se2:r) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+model_budal <- read_rds("model_budal_inla.rds")
+tmp_bu_par <- model_budal$tmp_bu_par
+tmp_bu_par_sum <- model_budal$tmp_bu_par_sum
+tmp_bu_par <- tmp_bu_par %>%
+ rowwise() %>%
+ mutate(EN = f_EN(eta, rho2, ES),
+ ss2 = se2 * 2 * gamma,
+ sE2 = sc2 * 2 * gamma_c,
+ sr2 = sh2 * gamma ^ 2,
+ logEN = log(EN),
+ lnK = log(EN / ES),
+ r = lnK * gamma) %>%
+ ungroup()
+
+tmp_bu_par_sum <- tmp_bu_par %>%
+ select(se2:r) %>%
+ pivot_longer(se2:r) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+model_surrey <- read_rds("model_surrey_inla.rds")
+tmp_su_par <- model_surrey$tmp_su_par
+tmp_su_par_sum <- model_surrey$tmp_su_par_sum
+tmp_su_par <- tmp_su_par %>%
+ rowwise() %>%
+ mutate(EN = f_EN(eta, rho2, ES),
+ ss2 = se2 * 2 * gamma,
+ sE2 = sc2 * 2 * gamma_c,
+ sr2 = sh2 * gamma ^ 2,
+ logEN = log(EN),
+ lnK = log(EN / ES),
+ r = lnK * gamma) %>%
+ ungroup()
+
+tmp_su_par_sum <- tmp_su_par %>%
+ select(se2:r) %>%
+ pivot_longer(se2:r) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+model_valley <- read_rds("model_valley_inla.rds")
+tmp_va_par <- model_valley$tmp_va_par
+tmp_va_par_sum <- model_valley$tmp_va_par_sum
+tmp_va_par <- tmp_va_par %>%
+ rowwise() %>%
+ mutate(EN = f_EN(eta, rho2, ES),
+ ss2 = se2 * 2 * gamma,
+ sE2 = sc2 * 2 * gamma_c,
+ sr2 = sh2 * gamma ^ 2,
+ logEN = log(EN),
+ lnK = log(EN / ES),
+ r = lnK * gamma) %>%
+ ungroup()
+
+tmp_va_par_sum <- tmp_va_par %>%
+ select(se2:r) %>%
+ pivot_longer(se2:r) %>%
+ group_by(name) %>%
+ nest() %>%
+ mutate(stats = map(data, f_sum)) %>%
+ unnest(stats) %>%
+ select(-data) %>%
+ ungroup()
+
+# Tables in paper:
+
+res_table <- tibble(Location = c("Ammarnäs", "Budal", "Eastern Wood", "Birdsong Valley")) %>%
+ bind_cols(bind_rows(pivot_wider(select(tmp_am_par_sum, name, mean), values_from = mean, names_from = name),
+ pivot_wider(select(tmp_bu_par_sum, name, mean), values_from = mean, names_from = name),
+ pivot_wider(select(tmp_su_par_sum, name, mean), values_from = mean, names_from = name),
+ pivot_wider(select(tmp_va_par_sum, name, mean), values_from = mean, names_from = name)))
+
+xtable::xtable(res_table[, c("Location", "ss2", "gamma", "sE2", "gamma_c", "r", "sr2")], digits = 3)
+
+res_table <- tibble(Location = c("Ammarnäs", "Budal", "Eastern Wood", "Birdsong Valley")) %>%
+ bind_cols(bind_rows(pivot_wider(select(tmp_am_par_sum, name, mean), values_from = mean, names_from = name),
+ pivot_wider(select(tmp_bu_par_sum, name, mean), values_from = mean, names_from = name),
+ pivot_wider(select(tmp_su_par_sum, name, mean), values_from = mean, names_from = name),
+ pivot_wider(select(tmp_va_par_sum, name, mean), values_from = mean, names_from = name)))
+
+xtable::xtable(res_table[, c("Location", "se2", "sh2", "sc2", "eta", "rho2", "z", "p0", "ES", "I")])
+
+# Tables for SI:
+
+xtable::xtable(tmp_am_par_sum) #, digits = 4) # to get the details on ss2 and sr2)
+xtable::xtable(tmp_bu_par_sum)
+xtable::xtable(tmp_su_par_sum) #, digits = 4) # to get the details on ss2, sE2 and sr2
+xtable::xtable(tmp_va_par_sum) #, digits = 4) # to get the details on ss2, sE2 and sr2