### Figure5.R
### Boxplots of root dry mass across experiments and soils at each harvest point


# Load necessary libraries and set up Groundhog for reproducibility
library(groundhog) 
set.groundhog.folder("C:\\Groundhog_R\\R4.4.1_2024-05-14") # Specifies where to store library versions
GroundhogDay <- '2024-05-14' # Set the snapshot date for reproducibility

# Load specific libraries using Groundhog for package management
groundhog.library("tidyverse",GroundhogDay) # Loads tidyverse for data manipulation and visualization
groundhog.library("rstatix",GroundhogDay)   # Loads rstatix for statistical tests
groundhog.library("multcompView",GroundhogDay) # Loads multcompView for multiple comparisons

# Set the working directory to the location of data files
setwd("D:\\Google Drive\\0_SpekboomResearchGroup2.0\\3_Experiments\\NMU22(07)SoilRootingExperiment\\DataForSubmission/")

# --- Load and Prepare Data ---
# Import the dataset from CSV
dat <- read_csv("NMU22(07)SoilRootingExperiment_data.csv") 

# Modify the 'Harvest' column to represent days after planting and create HarvestLabel for plotting
dat <- dat %>%
  mutate(
    Harvest = Harvest * 7 + 7, # Convert harvest from weeks to days after planting
    HarvestLabel = factor(Harvest, 
                          levels = 1:7 * 7 + 7, # Labels for harvest days
                          labels = paste0(1:7 * 7 + 7, " days after planting"))
  ) %>%
  filter(RootDryMass_g < 0.5) %>% # Filter out data with RootDryMass greater than 0.5
  mutate(SoilEnv = factor(paste(SiteCode, Environment, sep = "_"))) # Create a combined variable for Site and Environment

##########################################################################################
# Conduct Levene’s Test for Homogeneity of Variances
# This checks for variance equality across groups
for (i_harvest in unique(dat$HarvestLabel)) {
  dat %>%
    filter(HarvestLabel == i_harvest) %>%
    levene_test(RootDryMass_g ~ SoilEnv) %>% # Test variance across SoilEnv within each harvest group
    print() # Print results of each test
}

# Conduct Shapiro-Wilk Test for Normality
# This checks if data is normally distributed within each group
for (i_harvest in unique(dat$HarvestLabel)) {
  dat %>%
    filter(HarvestLabel == i_harvest) %>%
    group_by(SiteCode, Environment) %>% # Group data by site and environment
    shapiro_test(RootDryMass_g) %>% # Perform Shapiro-Wilk test for normality
    print() # Print results of each test
}

# Since normality assumption is violated, use non-parametric Kruskal-Wallis test

##########################################################################################
# Perform Kruskal-Wallis Test and Post-Hoc Dunn's Test for each harvest day

TK <- NULL # Initialize an empty list to store results
KS <- list() # Initialize list to store Kruskal-Wallis test results
for (i_harvest in unique(dat$HarvestLabel)) {
  # Kruskal-Wallis test to compare RootDryMass across different SoilEnv groups
  dat %>%
    filter(HarvestLabel == i_harvest) %>%
    kruskal_test(RootDryMass_g ~ SoilEnv, data = .) -> ktest # Perform Kruskal-Wallis test
  KS[[i_harvest]] <- ktest # Store Kruskal-Wallis result for each harvest
  
  # Post-hoc Dunn's Test to compare all pairs of groups within each harvest
  dat %>%
    filter(HarvestLabel == i_harvest) %>%
    dunn_test(RootDryMass_g ~ SoilEnv, data = .) -> dtest # Perform Dunn’s test
  p <- dtest$p.adj # Adjusted p-values from Dunn’s test
  names(p) <- paste0(dtest$group1, "-", dtest$group2) # Name adjusted p-values based on group comparisons
  
  cld <- multcompLetters(p) # Conduct letter-based multiple comparison for significance
  print(cld) # Print the letter results
  
  # Summarize maximum root dry mass by SoilEnv and Environment for each harvest
  dat %>%
    filter(HarvestLabel == i_harvest) %>%
    group_by(SoilEnv, Environment) %>%
    summarise(max = max(RootDryMass_g, na.rm = T)) %>% # Calculate max root dry mass
    arrange(desc(max)) -> Tk # Sort by max root dry mass in descending order
  Tk$Letters <- cld$Letters[Tk$SoilEnv] # Add letter-based significance results to data
  Tk$HarvestLabel <- i_harvest # Add harvest label to the data
  TK <- bind_rows(TK, Tk) # Combine results across harvest days
}

# --- Plot Results ---
# Plot boxplots for Root Dry Mass across SoilEnv and Environment, with statistical annotations
dat %>%
  mutate(SoilEnv = factor(SoilEnv,
                          levels = c("Site A_Glasshouse",
                                     "Site B_Glasshouse",
                                     "Site C_Glasshouse",
                                     "Site A_Growth_Chamber",
                                     "Site B_Growth_Chamber",
                                     "Site C_Growth_Chamber"))) %>%
  ggplot(aes(SoilEnv, RootDryMass_g)) +
  geom_boxplot(width = 0.3, aes(fill = Environment), alpha = 100) + # Boxplot of root dry mass
  geom_point(aes(x = as.numeric(SoilEnv) - 0.28, fill = Environment), shape = 22) + # Overlay data points
  scale_fill_manual(values = c("darkgrey", "white"), name = "Experiment") + # Set colors for fill
  geom_text(data = TK, aes(x = SoilEnv, y = max + 0.1, label = Letters)) + # Add letters indicating significance
  theme_bw() + # Set theme to black and white
  facet_wrap(~HarvestLabel, ncol = 2) + # Create separate panels for each harvest label
  geom_vline(xintercept = 3.5, colour = "darkgrey") + # Vertical line to separate groups
  scale_x_discrete("", labels = c("Site A_Glasshouse" = "Site A",
                                  "Site B_Glasshouse" = "Site B",
                                  "Site C_Glasshouse" = "Site C",
                                  "Site A_Growth_Chamber" = "Site A",
                                  "Site B_Growth_Chamber" = "Site B",
                                  "Site C_Growth_Chamber" = "Site C")) + # Custom x-axis labels
  theme(legend.position = c(0.80, 0.05), legend.justification = c(1, 0)) + # Position legend
  ylab("Root dry mass (g)") # Label for y-axis

# Save the plot as a JPEG image
ggsave("Figure5.jpg", width = 8, height = 10.5)