### Figure2.R
### PCA and histograms of thicket-wide plot soil parameters showing the three
### sampled soils relative to the available data.

### Script for:
# "Restoring South African subtropical succulent thicket using 
# Portulacaria afra: rooting variation across three soil types"
# by AJ Potts, D Liddell, C Clarke, N Galuszynski

### Libraries
library(groundhog)
set.groundhog.folder("C:\\Groundhog_R\\R4.4.1_2024-05-14") # Sets the folder where Groundhog will store library versions
GroundhogDay <- '2024-05-14' # Specifies the Groundhog snapshot date for reproducibility

groundhog.library("tidyverse", GroundhogDay) # Loads the Tidyverse package for data manipulation and visualization
groundhog.library("vegan", GroundhogDay) # Loads the Vegan package for ecological data analysis

# Set your working directory here
# setwd(")
setwd("D:\\Google Drive\\0_SpekboomResearchGroup2.0\\3_Experiments\\NMU22(07)SoilRootingExperiment\\DataForSubmission/") # Sets the working directory for the project

read_csv("twp_soils.csv") -> dat # Reads the soil data from a CSV file into a dataframe

### FIGURE 2 ###
dat %>%
  select(-ID, -LabNumber, -Soil, -Classification, -Comment, -TValue, -Silt_perc, -Sand_perc, -K_perc, -Na_perc, 
         -Ca_perc, -Mg_perc, -N_perc, -Org_C_perc, -StoneVolume_perc, -Clay_perc) %>% # Removes irrelevant columns for PCA
  na.omit() -> # Removes rows with missing data
  datForPCA

datForPCA %>%
  select(-Plot) %>%
  prcomp(scale = TRUE) -> # Performs PCA, scaling variables for equal weight
  PCA

PCA %>% summary() # Displays PCA results, including variance explained by each component

PCA %>%
  scores() %>% # Extracts PCA scores
  as.matrix() %>%
  as.data.frame() %>%
  bind_cols(datForPCA) %>%
  as_tibble() -> 
  forGG # Combines PCA scores with original data for plotting

PCA$rotation %>%
  as.data.frame() %>%
  rownames_to_column("Variable") %>%
  mutate(Variable = factor(Variable,
                           levels = c("pH_KCl", "C_perc", "PBrayII", "EC_mSm", "K_mgkg", "K_plus", "Ca_2plus", "Mg_2plus", "Na_plus"),
                           labels = c("pH", "Total C", "P Bray II", "EC", "K+", "", "Ca 2+", "Mg 2+", "Na+"))) -> 
  forGGcorrelates # Prepares variable loadings for PCA biplot

# Plot PCA biplot
forGG %>%
  ggplot() +
  geom_vline(xintercept = 0) +
  geom_hline(yintercept = 0) +
  geom_point(aes(PC1, PC2), shape = 1) +
  geom_segment(data = forGGcorrelates,
               aes(x = 0, y = 0, xend = PC1 * 4, yend = PC2 * 4),
               arrow = arrow(length = unit(0.2, 'cm')),
               colour = "darkgrey", linewidth = 1) +
  geom_text(data = forGGcorrelates, aes(PC1 * 5, PC2 * 5, label = Variable), colour = "darkgrey") +
  geom_point(data = filter(forGG, Plot %in% c("SiteA (soil1)", "SiteC (soil2)", "SiteB (soil3)")),
             aes(PC1, PC2), colour = "red", shape = 4, size = 5) +
  geom_label(data = filter(forGG, Plot %in% c("SiteA (soil1)", "SiteC (soil2)", "SiteB (soil3)")),
             aes(PC1, PC2, label = c("A", "C", "B")), colour = "black", size = 5) +
  theme_bw() +
  theme(legend.position = c(.17, .98),
        legend.justification = c("right", "top"),
        legend.box.just = "right",
        legend.margin = margin(6, 6, 6, 6),
        legend.box.background = element_rect(colour = "black")) +
  coord_equal() +
  xlab("PC1 (42.9%)") + # Customizes axis labels
  ylab("PC2 (19.3%)")
ggsave("Figure2.jpg", units = "cm", height = 10, width = 17) # Saves the plot as a JPG