library(tidyverse)
library(vegan)
library(indicspecies)

## 1) Read data and join days_standing ####
df1 <- read.csv("2021a carabids.csv", sep = ";", header = TRUE) %>%
  mutate(DATE = as.Date(DATE)) %>%
  filter(DATE < as.Date("2021-09-01"))

date_mapping_1 <- tibble(
  DATE = as.Date(c("2021-06-07", "2021-06-10", "2021-06-14",
                   "2021-06-17", "2021-07-12", "2021-07-16")),
  days_standing = c(3, 3, 4, 3, 5, 4)
)

df1 <- df1 %>%
  inner_join(date_mapping_1, by = "DATE") %>%
  mutate(TRIAL = "2021.RS")

df2 <- read.csv("2021b carabids.csv", sep = ";", header = TRUE) %>%
  mutate(DATE = as.Date(DATE)) %>%
  filter(DATE < as.Date("2021-09-01"))

date_mapping_2 <- tibble(
  DATE = as.Date(c("2021-06-15","2021-06-18","2021-06-22",
                   "2021-06-25","2021-07-12","2021-07-16")),
  days_standing = c(4, 3, 4, 3, 4, 4)
)

df2 <- df2 %>%
  inner_join(date_mapping_2, by = "DATE") %>%
  mutate(TRIAL = "2021.MS")

df3 <- read.csv("2022 carabids.csv", sep = ";", header = TRUE) %>%
  mutate(DATE = as.Date(DATE))

date_mapping_3 <- tibble(
  DATE = as.Date(c("2022-06-13","2022-06-15","2022-06-30",
                   "2022-07-04","2022-07-07","2022-07-11","2022-07-14")),
  days_standing = c(3, 2, 3, 4, 3, 4, 3)
)

df3 <- df3 %>%
  inner_join(date_mapping_3, by = "DATE") %>%
  mutate(TRIAL = "2022.MS")


## 2) Harmonise columns: add dummy DIRECTION/SURROUNDING for 2022 ####

if (!"DIRECTION"   %in% names(df3)) df3$DIRECTION   <- "std"
if (!"SURROUNDING" %in% names(df3)) df3$SURROUNDING <- "mulch"


## 3) Combine trials and select trap types ####

combined_df <- bind_rows(df1, df2, df3)

# keep only mulch traps for WT and S in 2021; controls keep all traps;
# 2022 has only "mulch"-equivalent traps
combined_df <- combined_df %>%
  filter(
    TREATMENT == "C" |
      (TREATMENT %in% c("WT", "S", "SIL") & SURROUNDING == "mulch")
  )

meta_cols <- c("TRIAL","DATE","BLOCK","TREATMENT","ROW",
               "DIRECTION","SURROUNDING","days_standing")
species_cols <- setdiff(names(combined_df), meta_cols)


## 4) Long format and remove NA trap records ####
combined_long <- combined_df %>%
  pivot_longer(cols = all_of(species_cols),
               names_to = "species",
               values_to = "count") %>%
  filter(!is.na(count))


## 5) Aggregate ####

# sampling effort = sum of trap-days per row
# activity density = total individuals / total trap-days
row_agg <- combined_long %>%
  group_by(TRIAL, TREATMENT, BLOCK, ROW, DATE, DIRECTION, SURROUNDING, species) %>%
  summarise(
    count_day    = sum(count),
    days_standing = first(days_standing),
    .groups = "drop"
  ) %>%
  group_by(TRIAL, TREATMENT, BLOCK, ROW, species) %>%
  summarise(
    total_count    = sum(count_day),
    total_trapdays = sum(days_standing),
    .groups = "drop"
  ) %>%
  filter(total_trapdays > 0) %>%
  mutate(activity_density = total_count / total_trapdays)


## 6) Community matrix for NMDS ####

row_agg <- row_agg %>%
  unite("Group", TRIAL, TREATMENT, BLOCK, ROW, sep = "_", remove = FALSE)

# remove species with zero total activity across all rows
row_agg <- row_agg %>%
  group_by(species) %>%
  filter(sum(activity_density) > 0) %>%
  ungroup()

comm_mat <- row_agg %>%
  dplyr::select(Group, species, activity_density) %>%
  pivot_wider(names_from = species,
              values_from = activity_density,
              values_fill = 0) %>%
  as.data.frame()

lookup_rows <- comm_mat %>%
  dplyr::select(Group) %>%
  separate(Group, into = c("TRIAL","TREATMENT","BLOCK","ROW"),
           sep = "_", remove = FALSE)

rownames(comm_mat) <- comm_mat$Group
comm_mat <- as.matrix(comm_mat[ , -1])


## 7) NMDS + plotting ####

# 2021.RS / 2021a
sel_RS_2021  <- lookup_rows$TRIAL == "2021.RS"
mat_RS_2021  <- comm_mat[sel_RS_2021, ]

dist_RS_2021 <- vegdist(mat_RS_2021, method = "bray")
set.seed(23841)
nmds_RS_2021 <- metaMDS(dist_RS_2021)

scores_RS_2021 <- scores(nmds_RS_2021) %>%
  as_tibble(rownames = "Group") %>%
  inner_join(lookup_rows, by = "Group")

centroid_RS_2021 <- scores_RS_2021 %>%
  group_by(TREATMENT) %>%
  mutate(
    centroid1 = mean(NMDS1),
    centroid2 = mean(NMDS2),
    .groups   = "drop"
  )

scores_RS_2021 %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = TREATMENT, fill = TREATMENT)) +
  theme_light() +
  stat_ellipse(geom = "polygon", level = 0.7, alpha = 0.1, show.legend = FALSE) +
  geom_point() +
  geom_point(data = centroid_RS_2021,
             mapping = aes(x = centroid1, y = centroid2, color = TREATMENT),
             size = 4, shape = 17) +
  labs(x = "NMDS 1", y = "NMDS 2",
       title = "carabids - trial 2021a") +
  theme(
    axis.text.x   = element_text(size = 8),
    axis.text.y   = element_text(size = 8),
    axis.title.x  = element_text(size = 9, colour = "gray30"),
    axis.title.y  = element_text(size = 9, colour = "gray30"),
    plot.title    = element_text(hjust = 0.98, size = 11, vjust = 0),
    panel.border  = element_rect(colour = "lightgray", size = 0.5),
    legend.position = "none"
  ) +
  scale_color_manual(name = "treatment",
                     breaks = c("C", "WT", "S"),
                     labels = c("control", "triticale/vetch", "straw"),
                     values = c("#EE6677", "#CCBB44", "#4477AA")) +
  scale_fill_manual(name = "treatment",
                    breaks = c("C", "WT", "S"),
                    labels = c("control", "triticale/vetch", "straw"),
                    values = c("#EE6677", "#CCBB44", "#4477AA"))
ggsave(filename = "NMDS 2021a - carabids.pdf",
       width = 8, height = 7, units = "cm")


# 2021.MS / 2021b
sel_MS_2021  <- lookup_rows$TRIAL == "2021.MS"
mat_MS_2021  <- comm_mat[sel_MS_2021, ]

dist_MS_2021 <- vegdist(mat_MS_2021, method = "bray")
set.seed(571903)
nmds_MS_2021 <- metaMDS(dist_MS_2021)

scores_MS_2021 <- scores(nmds_MS_2021) %>%
  as_tibble(rownames = "Group") %>%
  inner_join(lookup_rows, by = "Group")

centroid_MS_2021 <- scores_MS_2021 %>%
  group_by(TREATMENT) %>%
  mutate(
    centroid1 = mean(NMDS1),
    centroid2 = mean(NMDS2),
    .groups   = "drop"
  )

scores_MS_2021 %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = TREATMENT, fill = TREATMENT)) +
  theme_light() +
  stat_ellipse(geom = "polygon", level = 0.7, alpha = 0.1, show.legend = FALSE) +
  geom_point() +
  geom_point(data = centroid_MS_2021,
             mapping = aes(x = centroid1, y = centroid2, color = TREATMENT),
             size = 4, shape = 17) +
  labs(x = "NMDS 1", y = "NMDS 2",
       title = "carabids - trial 2021b") +
  theme(
    axis.text.x   = element_text(size = 8),
    axis.text.y   = element_text(size = 8),
    axis.title.x  = element_text(size = 9, colour = "gray30"),
    axis.title.y  = element_text(size = 9, colour = "gray30"),
    plot.title    = element_text(hjust = 0.98, size = 11, vjust = 0),
    panel.border  = element_rect(colour = "lightgray", size = 0.5),
    legend.position = "none"
  ) +
  scale_color_manual(name = "treatment",
                     breaks = c("C", "WT", "S"),
                     labels = c("control", "triticale/vetch", "straw"),
                     values = c("#EE6677", "#CCBB44", "#4477AA")) +
  scale_fill_manual(name = "treatment",
                    breaks = c("C", "WT", "S"),
                    labels = c("control", "triticale/vetch", "straw"),
                    values = c("#EE6677", "#CCBB44", "#4477AA"))
ggsave(filename = "NMDS 2021b - carabids.pdf",
       width = 8, height = 7, units = "cm")

# 2022.MS / 2022
sel_MS_2022  <- lookup_rows$TRIAL == "2022.MS"
mat_MS_2022  <- comm_mat[sel_MS_2022, ]

dist_MS_2022 <- vegdist(mat_MS_2022, method = "bray")
set.seed(98217)
nmds_MS_2022 <- metaMDS(dist_MS_2022)

scores_MS_2022 <- scores(nmds_MS_2022) %>%
  as_tibble(rownames = "Group") %>%
  inner_join(lookup_rows, by = "Group")

centroid_MS_2022 <- scores_MS_2022 %>%
  group_by(TREATMENT) %>%
  mutate(
    centroid1 = mean(NMDS1),
    centroid2 = mean(NMDS2),
    .groups   = "drop"
  )

scores_MS_2022 %>%
  ggplot(aes(x = NMDS1, y = NMDS2, color = TREATMENT, fill = TREATMENT)) +
  theme_light() +
  stat_ellipse(geom = "polygon", level = 0.7, alpha = 0.1, show.legend = FALSE) +
  geom_point() +
  geom_point(data = centroid_MS_2022,
             mapping = aes(x = centroid1, y = centroid2, color = TREATMENT),
             size = 4, shape = 17) +
  labs(x = "NMDS 1", y = "NMDS 2",
       title = "carabids - trial 2022") +
  theme(
    axis.text.x   = element_text(size = 8),
    axis.text.y   = element_text(size = 8),
    axis.title.x  = element_text(size = 9, colour = "gray30"),
    axis.title.y  = element_text(size = 9, colour = "gray30"),
    plot.title    = element_text(hjust = 0.98, size = 11, vjust = 0),
    panel.border  = element_rect(colour = "lightgray", size = 0.5),
    legend.position = "none"
  ) +
  scale_color_manual(name = "treatment",
                     breaks = c("C", "WT", "SIL"),
                     labels = c("control", "triticale/vetch", "silage"),
                     values = c("#EE6677", "#CCBB44", "#228833")) +
  scale_fill_manual(name = "treatment",
                    breaks = c("C", "WT", "SIL"),
                    labels = c("control", "triticale/vetch", "silage"),
                    values = c("#EE6677", "#CCBB44", "#228833"))
ggsave(filename = "NMDS 2022 - carabids.pdf",
       width = 8, height = 7, units = "cm")


# Legend NMDS
dummy_plot <- ggplot(data.frame(x = 1:4, y = 1:4, shape = c("control", "straw", "triticale/vetch", "silage")), 
                     aes(x, y, shape = shape, color = shape)) +
  geom_point(size = 3) +
  scale_shape_manual(name = "treatments:",
                     values = rep(17, 4),
                     labels = c("control", "straw","triticale/vetch", "silage")) +
  scale_color_manual(name = "treatments:",
                     values = c("#EE6677","#4477AA", "#CCBB44", "#228833"),
                     labels = c("control", "straw", "triticale/vetch", "silage")) +
  theme_void() +
  theme(legend.position = "bottom",
        legend.title = element_text(size = 10),
        legend.text = element_text(size = 10))
dummy_plot
leg <- get_legend(dummy_plot)
as_ggplot(leg)
ggsave(filename = "legend_NMDS.pdf",
      width = 4.25, height = 0.2)


## 8) PERMANOVA (adonis2) and dispersion tests ####

pairwise_permanova <- function(dist_mat, meta_df, treat_col = "TREATMENT",
                               block_col = "BLOCK", seed = 446531) {
  treatments <- unique(meta_df[[treat_col]])
  pairs <- combn(treatments, 2, simplify = FALSE)
  
  pairwise_results <- lapply(pairs, function(pair) {
    subset_indices <- meta_df[[treat_col]] %in% pair
    subset_dist <- as.dist(as.matrix(dist_mat)[subset_indices, subset_indices])
    subset_meta <- meta_df[subset_indices, , drop = FALSE]
    
    set.seed(seed)
    result <- adonis2(subset_dist ~ subset_meta[[treat_col]],
                      strata = subset_meta[[block_col]])
    
    data.frame(
      Comparison = paste(pair[1], "vs", pair[2]),
      R2 = result$R2[1],
      F  = result$F[1],
      P  = result$`Pr(>F)`[1]
    )
  })
  
  pairwise_table <- do.call(rbind, pairwise_results)
  pairwise_table$P_adjusted <- p.adjust(pairwise_table$P, method = "bonferroni")
  pairwise_table
}

# 2021.RS / 2021a
dist <- dist_RS_2021
lookup <- scores_RS_2021 %>%
  dplyr::select(Group, TRIAL, TREATMENT, BLOCK, ROW)

adonis2(dist ~ lookup$TREATMENT,
        strata = lookup$BLOCK)

disp <- betadisper(dist, lookup$TREATMENT)
print(disp)
perm_test <- permutest(disp)
print(perm_test)

pairwise_table_RS_2021 <- pairwise_permanova(dist, lookup,
                                             treat_col = "TREATMENT",
                                             block_col = "BLOCK")
print(pairwise_table_RS_2021)


# 2021.MS / 2021b
dist <- dist_MS_2021
lookup <- scores_MS_2021 %>%
  dplyr::select(Group, TRIAL, TREATMENT, BLOCK, ROW)

adonis2(dist ~ lookup$TREATMENT,
        strata = lookup$BLOCK)

disp <- betadisper(dist, lookup$TREATMENT)
print(disp)
perm_test <- permutest(disp)
print(perm_test)

pairwise_table_MS_2021 <- pairwise_permanova(dist, lookup,
                                             treat_col = "TREATMENT",
                                             block_col = "BLOCK")
print(pairwise_table_MS_2021)


# 2022.MS / 2022
dist <- dist_MS_2022
lookup <- scores_MS_2022 %>%
  dplyr::select(Group, TRIAL, TREATMENT, BLOCK, ROW)

adonis2(dist ~ lookup$TREATMENT,
        strata = lookup$BLOCK)

disp <- betadisper(dist, lookup$TREATMENT)
print(disp)
perm_test <- permutest(disp)
print(perm_test)

pairwise_table_MS_2022 <- pairwise_permanova(dist, lookup,
                                             treat_col = "TREATMENT",
                                             block_col = "BLOCK")
print(pairwise_table_MS_2022)



## 9) Diversity calculations ####

comm_df <- as.data.frame(comm_mat) %>%
  rownames_to_column(var = "Group")

total_activity <- comm_df %>%
  mutate(total_activity = rowSums(across(-Group))) %>%
  dplyr::select(Group, total_activity)

diversities <- comm_df %>%
  column_to_rownames("Group") %>%
  as.matrix() %>%
  {
    tibble(
      Group            = rownames(.),
      species_richness = specnumber(.),                  # number of species with activity > 0
      shannon          = diversity(., index = "shannon"),
      simpson          = diversity(., index = "simpson"),
      invsimpson       = 1 / diversity(., index = "simpson"),
      evenness         = shannon / log(species_richness)
    )
  } %>%
  left_join(total_activity, by = "Group")

metadata_diversities <- diversities %>%
  inner_join(lookup_rows, by = "Group")

metadata_diversities$TRIAL <- factor(metadata_diversities$TRIAL, 
                                     levels = c("2021.RS", "2021.MS", "2022.MS"))

trial_labels <- c(
  "2021.RS" = "trial 2021a",
  "2021.MS" = "trial 2021b",
  "2022.MS" = "trial 2022"
)


## 10) Plots: species number and Shannon ####

ggplot(metadata_diversities,
       aes(y = factor(TREATMENT, levels = c("SIL", "S", "WT", "C")),
           color = factor(TREATMENT, levels = c("SIL", "S", "WT", "C")),
           x = species_richness)) +
  geom_jitter(width = 0, height = 0.15, size = 1.5, alpha = 0.5) +
  scale_y_discrete(breaks = c("C", "WT", "SIL", "S"),
                   labels = c("control", "triticale/vetch",
                              "silage", "straw")) +
  scale_x_continuous(breaks = seq(5, 14, by = 2)) + 
  scale_color_manual(values = c("C" = "#EE6677", "WT" = "#CCBB44", 
                                "SIL" = "#228833", "S" = "#4477AA")) +
  stat_summary(fun = mean, geom = "point", 
               shape = 4, size = 2, stroke = 1,  
               color = "black", aes(group = TREATMENT)) +
  theme_light() +
  labs(x = "Number of species", y = NULL) +
  theme(axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 9),
        axis.title.x = element_text(size = 9),
        legend.position = "none") +
  facet_grid(~TRIAL, scales = "free_y", space = "free_y",
             labeller = labeller(TRIAL = trial_labels))

ggsave("Speciesnumber jitter - carabids.pdf",
       width = 17, height = 4.5, units = "cm")

# Shannon jitter
ggplot(metadata_diversities,
       aes(y = factor(TREATMENT, levels = c("SIL", "S", "WT", "C")),
           color = factor(TREATMENT, levels = c("SIL", "S", "WT", "C")),
           x = shannon)) +
  geom_jitter(width = 0, height = 0.15, size = 1.5, alpha = 0.5) +
  scale_y_discrete(breaks = c("C", "WT", "SIL", "S"),
                   labels = c("control", "triticale/vetch",
                              "silage", "straw")) +
  scale_color_manual(values = c("C" = "#EE6677", "WT" = "#CCBB44", 
                                "SIL" = "#228833", "S" = "#4477AA")) +
  stat_summary(fun = mean, geom = "point", 
               shape = 4, size = 2, stroke = 1,  
               color = "black", aes(group = TREATMENT)) +
  theme_light() +
  labs(x = "Shannon index", y = NULL) +
  theme(axis.text.x = element_text(size = 8),
        axis.text.y = element_text(size = 9),
        axis.title.x = element_text(size = 9),
        legend.position = "none") +
  facet_grid(~TRIAL, scales = "free_y", space = "free_y",
             labeller = labeller(TRIAL = trial_labels)) +
  scale_x_continuous(breaks = c(1.0, 1.5, 2.0))

ggsave("Shannon jitter - carabids.pdf",
       width = 17, height = 4.5, units = "cm")


## 3) Indicator species analysis (IndVal) ####

run_indval <- function(trial_id,
                       comm_mat,
                       lookup_rows,
                       p_adj = "BH",
                       seed  = 1) {
  sel       <- lookup_rows$TRIAL == trial_id
  comm_sub  <- comm_mat[sel, , drop = FALSE]
  meta_sub  <- lookup_rows[sel, , drop = FALSE]
  
  df_sub <- as.data.frame(comm_sub)
  df_sub$TREATMENT <- factor(meta_sub$TREATMENT)
  
  set.seed(seed)
  indval <- multipatt(df_sub %>% dplyr::select(-TREATMENT),
                      df_sub$TREATMENT,
                      duleg   = TRUE,
                      control = how(nperm = 999))
  
  p_values   <- indval$sign$p.value
  p_adjusted <- p.adjust(p_values, method = p_adj)
  indval$sign$p.adjusted <- p_adjusted
  
  indval
}

# 2021.RS / 2021a
indval_c_2021_RS <- run_indval("2021.RS",
                               comm_mat    = comm_mat,
                               lookup_rows = lookup_rows,
                               p_adj = "BH",
                               seed  = 235)
summary(indval_c_2021_RS)

# 2021.MS / 2021b
indval_c_2021_MS <- run_indval("2021.MS",
                               comm_mat    = comm_mat,
                               lookup_rows = lookup_rows,
                               p_adj = "BH",
                               seed  = 878556)
summary(indval_c_2021_MS)

# 2022.MS / 2022
indval_c_2022_MS <- run_indval("2022.MS",
                               comm_mat    = comm_mat,
                               lookup_rows = lookup_rows,
                               p_adj = "BH",
                               seed  = 95436)
summary(indval_c_2022_MS)