######################################################################################################
# Code for manuscript: Chia et al. Enclosed bird nests driven by predation and thermoregulation
# 2024 November
# This script performs PCA on climatic variables and pairwise correlations among traits

## INPUTS
# Species trait/environmental data "data/all_traits.csv"

## OUTPUTS
# Species trait/environmental data with climatic PCA "data/all_traits_pca.csv"
# PCA and pairwise correlation plots and tables: Table S1 & S2, Figure S1 & S2
######################################################################################################

library(factoextra) # PCA visualization
library(PerformanceAnalytics) # pairwise correlations visualization

# import dataset
dt <- read.csv("data/all_traits.csv", row.names = 1)

## PCA for climatic variables
# (variable explanation: https://chelsa-climate.org/exchelsa-extended-bioclim/)
pca <- prcomp(as.matrix(dt[,c("tmean","dtr","hurs","rsds","vpd","precp","sfcWind")]), scale = T)

# Make necessary change of signs of PC1 & PC2 loadings
# (-hurs for PC1 (aridity), +tmean for PC2 (temperature))
if (pca$rotation["hurs", 1] > 0) {
  pca$x[, "PC1"] <- -pca$x[,"PC1"]
  pca$rotation[, 1] <- -pca$rotation[,1]
}
if (pca$rotation["tmean", 2] < 0) {
  pca$x[, "PC2"] <- -pca$x[, "PC2"]
  pca$rotation[, 2] <- -pca$rotation[, 2]
}

# PCA summary table (loadings and variance explained) (Table S1)
pca.summary <- rbind(pca$rotation, (pca$sdev)^2 / sum((pca$sdev)^2) * 100)
pca.summary

# Add PC1 & PC2 into the trait dataset and save
dt$PC1 <- pca$x[, "PC1"] # aridity
dt$PC2 <- pca$x[, "PC2"] # temperature

write.csv(dt, "data/all_traits_pca.csv")

# Pairwise correlation between life-history traits (Table S2)
cor(dt[, c("Enclosed","Ground","Cooperative","Egg","Clutch")])

# Plot PCA biplot (Figure S1)
rownames(pca$rotation) <- c("Mean tmperature","DTR","Humidity","Solar radiation","Vapor pressure deficit","Precipitation","Wind")
fviz_pca_biplot(pca, repel = TRUE, label = c("var","quali"), axes = c(1, 2),
                col.var = "black", col.ind = "cornflowerblue", alpha.ind = 0.4, labelsize = 4, pointsize = 1, stroke = 0,
                xlab = paste0("PC1 (", round(pca.summary[nrow(pca.summary), 1], 1),"%)"), 
                ylab = paste0("PC2 (", round(pca.summary[nrow(pca.summary), 2], 1),"%)")) +
  theme(panel.grid = element_blank())

# Pairwise correlation between environmental variables (Figure S2)
chart.Correlation(dt[, c("PC1","PC2","npp","tmean")], histogram = T, method = "pearson", pch = 19)
