dynamicSAD:

# Simulation examples ####

# The current file assumes the estimations have already been done and stored
# as est_res_ex0X.rds, where
# X = 1: assuming temporal correlation and species heterogeneity
# X = 2: assuming no species heterogeneity
# X = 3: assuming no temporal correlation

# First we simulate 100 data sets and analyse those data sets based on
# different assumptions

# Parameters ####

input_ex01 <- list(n_sp = 100, # number of species
                   n_year = 20, # number of time points
                   n_loc = 1, # number of locations
                   n_obs = 2, # number of observers, per space-time point
                   mu = 1, # mean log abundance
                   sigma_r_sq = 2, # species heterogeneity
                   sigma_s_sq = 1, # environmental variance
                   sigma_d_sq = 0.5, # sampling variance
                   delta = 0.2, # "strength of density regulation"
                   eta = 0.2, # 1 / "spatial scaling of noise"
                   pos_x = NULL, # x-coordinates of locations
                   pos_y = NULL) # y-coordinates of locations

sim_ex01 <- list()

for (i in 1:1){

  sim_ex01[[i]] <- f_sim_sad_fast(input_ex01)

  cat(i, "%", " \r")

}
# rm(i)

sim_ex01[[1]] %>%
  ggplot(aes(x = interaction(year_m, pos_m, observer), y = species)) +
  geom_tile(aes(fill = log(abundance)), colour = "black") +
  geom_tile(data = filter(sim_ex01[[1]], abundance == 0),
              aes(fill = abundance), fill = "white", colour = "black") +
  scale_fill_viridis_c() +
  theme(axis.text.x = element_text(hjust = 1, angle = 45))


# Assumption 1: temporal correlation, heterogeneity, sampling variation ####

# est_ex01_mt <- list()
#
# for (i in 1:100){
#
#   est_ex01_mt[[i]] <- try(f_est_sad_mt(data = filter(sim_ex01[[i]],
#                                                      abundance > 0),
#                                        n_boot = 10,
#                                        n_rep = 1,
#                                        n_sp_mult = 1))
#
#   cat(i, "%", "c", " \r")
#
# }

# est_res_ex01 <- c()
#
# for (i in 1:100){
#
#   est_res_ex01 <- rbind(est_res_ex01,
#                         add_column(est_ex01_mt[[i]],
#                                    replicate = i,
#                                    iteration = 1:nrow(est_ex01_mt[[i]])))
#
#   cat(i, "%", " \r")
#
# }

# Assumption 2: temporal correlation, no heterogeneity, sampling variation ####

# est_ex02_mt <- list()
#
# for (i in 1:100){
#
#   est_ex02_mt[[i]] <- try(f_est_sad_mt_nh(data = filter(sim_ex01[[i]],
#                                                      abundance > 0),
#                                        n_boot = 10,
#                                        n_rep = 1,
#                                        n_sp_mult = 1))
#
#   cat(i, "%", "c", " \r")
#
# }


# est_res_ex02 <- c()
#
# for (i in 1:100){
#
#   est_res_ex02 <- rbind(est_res_ex02,
#                         add_column(est_ex02_mt[[i]],
#                                    replicate = i,
#                                    iteration = 1:nrow(est_ex02_mt[[i]])))
#
#   cat(i, "%", " \r")
#
# }

# Assumption 3: no temporal correlation, heterogeneity, sampling variation ####

# est_ex03_mt <- list()
#
# for (i in 1:100){
#
#   est_ex03_mt[[i]] <- try(f_est_sad_mt_nc(data = filter(sim_ex01[[i]],
#                                                         abundance > 0),
#                                           n_boot = 10,
#                                           n_rep = 1,
#                                           n_sp_mult = 1))
#   cat(i, "%", " \r")
#
# }

# est_res_ex03 <- c()
#
# for (i in 1:100){
#
#   est_res_ex03 <- rbind(est_res_ex03,
#                         add_column(est_ex03_mt[[i]],
#                                    replicate = i,
#                                    iteration = 1:nrow(est_ex03_mt[[i]])))
#
#   cat(i, "%", " \r")
#
# }

# Store results ####

# write_rds(est_res_ex01, "est_res_ex01.rds")
# write_rds(est_res_ex02, "est_res_ex02.rds")
# write_rds(est_res_ex03, "est_res_ex03.rds")


# Load results (if available)
est_res_ex01 <- read_rds("est_res_ex01.rds")
est_res_ex02 <- read_rds("est_res_ex02.rds")
est_res_ex03 <- read_rds("est_res_ex03.rds")

# Final figure ####

est_res_sum <- bind_rows(est_res_ex01,
                         est_res_ex02,
                         est_res_ex03) %>%
  add_column(method = rep(c("Temporal correlation and heterogeneity",
                            "No heterogeneity",
                            "No temporal correlation"),
                          c(nrow(est_res_ex01),
                            nrow(est_res_ex01),
                            nrow(est_res_ex01))))

est_res_sum_long <- est_res_sum %>%
  mutate("sigma[s]^2" = exp(lsigma_e)^2,
         delta = ifelse(is.na(ldelta), NA, exp(ldelta)),
         "sigma[r]^2" = ifelse(is.na(lsigma_r), NA, exp(lsigma_r)^2),
         "theta^2" = exp(lsigma_d)^2) %>%
  pivot_longer(cols = c(mu, "sigma[s]^2":"theta^2"))

true_values <- tibble(name = rep(c("mu",
                                   "delta",
                                   "sigma[r]^2",
                                   "theta^2",
                                   "sigma[s]^2"), 3),
                      value = rep(c(1, 0.2, 2, 0.5, 1), 3),
                      method = rep(c("No heterogeneity",
                                     "Temporal correlation and heterogeneity",
                                     "No temporal correlation"), each = 5))

est_res_sum_long %>%
  filter(iteration == 11) %>%
  ggplot(aes(x = value)) +
  geom_histogram(bins = 20, colour = "black", fill = "light grey") +
  geom_vline(data = true_values, aes(xintercept = value)) +
  theme_bw() +
  facet_grid(method ~ name,
             labeller = labeller(method = label_wrap_gen(30),
                                 name = label_parsed),
             scales = "free_x") +
  labs(y = "Count", x = "Estimated value")

# ggsave("example_assumptions.pdf", width = 9, height = 6)