# --- ROBUST LINEAR MIXED MODELS & PAIRWISE CONTRASTS --- #

# Load required packages and dataset
library(robustlmm)
library(emmeans)
library(tidyr)
library(dplyr)
woodpecker.data <- read.csv("/RLMM_and_hedgesg_data.csv")

# PC1 analysis: fit robust model, calculate marginal means, and format pairwise contrasts
LMMPC1 <- robustlmm::rlmer(RC1 ~ spp + (1 | individual2) * (1 | subspecies) * (1|type.2), data = woodpecker.data, REML = TRUE)
EMMPC1 <- emmeans::emmeans(LMMPC1, specs = list(pairwise ~ spp), adjust = "bh")
EMMPC1_contrasts <- as.data.frame(emmeans::contrast(EMMPC1, method = "pairwise", adjust = "bh")) %>%
  separate(col = contrast, into = c("spp_ref", "spp_foc"), sep = " - ", remove = TRUE)

# PC2 analysis: fit robust model, calculate marginal means, and format pairwise contrasts
LMMPC2 <- robustlmm::rlmer(RC2 ~ spp + (1 | individual2) * (1 | subspecies) * (1|type.2), data = woodpecker.data, REML = TRUE)
EMMPC2 <- emmeans::emmeans(LMMPC2, specs = list(pairwise ~ spp), adjust = "bh")
EMMPC2_contrasts <- as.data.frame(emmeans::contrast(EMMPC2, method = "pairwise", adjust = "bh")) %>%
  separate(col = contrast, into = c("spp_ref", "spp_foc"), sep = " - ", remove = TRUE)

# PC3 analysis: fit robust model, calculate marginal means, and format pairwise contrasts
LMMPC3 <- robustlmm::rlmer(RC3 ~ spp + (1 | individual2) * (1 | subspecies) * (1|type.2), data = woodpecker.data, REML = TRUE)
EMMPC3 <- emmeans::emmeans(LMMPC3, specs = list(pairwise ~ spp), adjust = "bh")
EMMPC3_contrasts <- as.data.frame(emmeans::contrast(EMMPC3, method = "pairwise", adjust = "bh")) %>%
  separate(col = contrast, into = c("spp_ref", "spp_foc"), sep = " - ", remove = TRUE)

# Export p-values for heatmap visualization
write.csv(EMMPC1_contrasts, "/pvalues_pc1.csv", row.names = FALSE)
write.csv(EMMPC2_contrasts, "/pvalues_pc2.csv", row.names = FALSE)
write.csv(EMMPC3_contrasts, "/pvalues_pc3.csv", row.names = FALSE)

