library(readxl)
library(tidyverse)
library(cowplot)
library(RColorBrewer)
library(segmented)
library(MLmetrics)
library(zoo)
library(relaimpo)
library(ggstatsplot)
library(mclust)
library(tcltk)

setwd(getwd()) ##Path to the folder which contains the data files

##################PART I: Constructing model on sN######################
data1 <- read_xlsx("DataS1_N2O_exp.xlsx")
N2O.all <- cbind.data.frame(data1$Lon, data1$Lat, data1$xlon, data1$ylat,
                            data1$`MAP(mm)`, data1$`MAT(K)`, data1$MAP.cv, data1$MAT.cv,
                            data1$`Ndepo(kgN ha-1 yr-1)`, as.numeric(data1$`MAN(kgN ha-1 yr-1)`), data1$MAN.cv,
                            as.factor(data1$Biome), as.factor(data1$Region),
                            data1$`Clay(%)`, data1$`Sand(%)`,
                            data1$`N_Add_Std(kgN ha-1 yr-1)`, as.numeric(data1$FluxYear),
                            data1$`RN2O(kgN ha-1 yr-1)`, data1$Reference)
colnames(N2O.all) <- c("Lon", "Lat", "xlon", "ylat", "MAP", "MAT", "MAP.cv", "MAT.cv", 
                       "Ndepo", "MAN", "MAN.cv", "Biome", "Region", "Clay", "Sand",
                       "N_Add", "Year", "RN2O", "Reference")
N2O.all$Biome <- ordered(N2O.all$Biome, levels = c("Tropical", "Temperate", "Boreal"))
N2O.all$N_Add[is.na(N2O.all$N_Add)] <- 0
N2O.all$N_Input <- N2O.all$N_Add + N2O.all$Ndepo
data2 <- read_xlsx("DataS2_N2O_obs.xlsx")
N2O.valid <- cbind.data.frame(data2$ID, data2$Lon, data2$Lat, data2$xlon, data2$ylat, data2$`MAP(mm)`, data2$`MAT(K)`, 
                              data2$MAP.cv, data2$MAT.cv, data2$`Ndepo(kgN ha-1 yr-1)`, data2$`MAN(kgN ha-1 yr-1)`, data2$MAN.cv, 
                              data2$Biome, data2$`Clay(%)`, data2$`Sand(%)`, data2$FluxYear, data2$`N2O flux (kg N ha-1 y-1)`,
                              data2$Reference)
colnames(N2O.valid) <- c("ID", "Lon", "Lat", "xlon", "ylat", "MAP", "MAT", "MAP.cv", "MAT.cv", 
                         "Ndepo", "MAN", "MAN.cv", "Biome", "Clay", "Sand",
                         "Year", "RN2O", "Reference")
N2O.valid$Biome <- ordered(N2O.valid$Biome, levels = c("Tropical", "Temperate", "Boreal"))
data3 <- read_xlsx("DataS3_Ncycle_exp.xlsx", skip = 1)
Ncycling <- cbind.data.frame(data3$Lon, data3$Lat, data3$Biome, data3$`c1 (kgN kgN-1)`, data3$`c2 (kgN kgN-1)`, data3$Reference)
colnames(Ncycling) <- c("Lon", "Lat", "Biome", "c1", "c2", "Reference")
Ncycling$Biome <- ordered(Ncycling$Biome, levels = c("Tropical", "Temperate", "Boreal"))

old <- theme_set(theme_bw())
theme_update(
  axis.text.x = element_text(angle = 0, hjust = 1, vjust = 1.1, size = 8), # Angle x-axis text
  axis.text.y = element_text(size = 8), # Size y-axis numbers
  axis.ticks = element_blank(), # No ticks
  axis.title.x = element_text(size = 8), # Size x-axis title
  axis.title.y = element_text(size = 8), # Size y-axis title
  panel.grid.minor = element_blank(), # No minor grid lines
  panel.grid.major.x = element_blank(), # No major x-axis grid lines
  panel.grid.major.y = element_blank(),
  text=element_text(size = 8)
)

##Detecting change point in the relationship between N input and RN2O
N2O.all400 <- N2O.all %>%
  filter(N_Input <= 400)
ggplot(data = N2O.all400, aes(x = N_Input, y = RN2O))+
  geom_point(color = "grey90", size = 2) +
  geom_smooth(color = "red", fill = "grey40")
mod0 <- lm(data = N2O.all400 , RN2O ~ N_Input)
seg.mod1 <- segmented(mod0, seg.Z = ~ N_Input, npsi = 1)
summary(seg.mod1)
ggplot(data = N2O.all400, aes(x = N_Input, y = RN2O)) +
  geom_point(color = "grey", size = 2) +
  geom_segment(x = 0, y = seg.mod1$coefficients[1], xend = seg.mod1$psi[1,2], 
               yend = seg.mod1$coefficients[1] + seg.mod1$psi[1,2]*seg.mod1$coefficients[2],
               color = "red")+
  geom_segment(x = seg.mod1$psi[1,2], y = seg.mod1$coefficients[1] + seg.mod1$psi[1,2]*seg.mod1$coefficients[2],
               xend = 400, yend = (seg.mod1$coefficients[2] + seg.mod1$coefficients[3]) * 400 + 
               (seg.mod1$coefficients[1] - seg.mod1$psi[1,2] * seg.mod1$coefficients[3]),
               color = "red")+
  geom_vline(xintercept = seg.mod1$psi[1,2], color = "blue")+
  geom_ribbon(aes(xmin = seg.mod1$psi[1,2] - seg.mod1$psi[1,3], xmax = seg.mod1$psi[1,2] + seg.mod1$psi[1,3]),
              fill = "blue", alpha = 0.2)+
  xlab("N Input Rate (kgN ha-1 yr-1)")+
  ylab("RN2O (kgN ha-1 yr-1)")
#ggsave("plot.segmod2.svg", width = 9, height = 7, units = "cm")
seg.mod2 <- segmented(mod0, seg.Z = ~N_Input, npsi = 4)  ##change npsi to set No. of change points to be estimated
exp((BIC(seg.mod2) - BIC(seg.mod1))/2)  #Calculating bayes factor, 
#BF>150 means we have >99% confidence that seg.mod1 is better

##Calculating sN and R0 for each grid having observations
N2O.LN <- subset(N2O.all, N2O.all$N_Add <= 150)
s <- unique.data.frame(cbind(N2O.LN$xlon, N2O.LN$ylat))
mod.summ <- {}
mod.valid <- {}
for(i in 1:length(s[,1])){
  temp1 <- subset(N2O.LN, N2O.LN$xlon==s[i,1] & N2O.LN$ylat==s[i,2])
  if(any(temp1$N_Add==0)&any(temp1$N_Add>0)){
    mod1 <- lm(temp1$RN2O ~ temp1$N_Input)
    a <- summary(mod1)
    p <- pf(a$fstatistic[1], a$fstatistic[2], a$fstatistic[3], lower.tail = FALSE)
    if(is.na(a$adj.r.squared)| a$adj.r.squared>0){
      mod.summ <- rbind(mod.summ, c(temp1$xlon[1], temp1$ylat[1], temp1$Lon[1], temp1$Lat[1],temp1$Biome[1], 
                                    temp1$MAP[1], temp1$MAP.cv[1], 
                                    temp1$MAT[1], temp1$MAT.cv[1], temp1$Ndepo[1], temp1$MAN[1], temp1$MAN.cv[1], temp1$Clay[1], temp1$Sand[1], #paste(temp1$Reference, collapse = ""), 
                                    a$coefficients[2,1], a$coefficients[2,4], a$coefficients[1,1], a$coefficients[1,4], a$df[2]+2, a$adj.r.squared, p))
    }
  }
}
mod.summ <- as.data.frame(mod.summ)
colnames(mod.summ) <- c("xlon", "ylat", "Lon", "Lat", "Biome", "MAP", "MAP.cv", "MAT", "MAT.cv", "Ndepo", "MAN", "MAN.cv", "Clay", "Sand", #"Reference",
                        "sN", "p_sN", "R0", "p_R0", "n", "R2adjusted", "pvalue")
mod.summ$Biome <- ordered(mod.summ$Biome, levels = c(1,2,3), labels = c("Tropical", "Temperate", "Boreal"))
mod.valid <- N2O.valid %>%
  group_by(xlon, ylat, Year, Biome) %>%
  summarize(across(where(is.numeric), mean))

#Checking normality
dens_sN <- ggplot(data = mod.summ)+
  geom_density(aes(sN))
dens_sqrtsN <- ggplot(data = mod.summ)+
  geom_density(aes(sqrt(sN)))
dens_R0 <- ggplot(data = mod.summ)+
  geom_density(aes(R0))
dens_logR0p <- ggplot(data = mod.summ)+
  geom_density(aes(log(R0+0.5)))
plot_grid(dens_sN, dens_sqrtsN, dens_R0, dens_logR0p, ncol=2)
shapiro.test(mod.summ$sN)
shapiro.test(sqrt(mod.summ$sN))
shapiro.test(mod.summ$R0)
shapiro.test(log(mod.summ$R0+0.5))
#ggsave("normality.svg", height = 12, width = 16, units = "cm")

#Correlation bewteen sN and environmental factors
cor.test(log(mod.summ$sN), log(mod.summ$MAN.cv))
cor.test(log(mod.summ$sN), log(mod.summ$MAP.cv))
cor.test(log(mod.summ$sN), log(mod.summ$MAN))
plot.sN.MAN <- ggplot(data = mod.summ, aes(x = log(sN), y = log(MAN)))+
  geom_point()+
  geom_smooth(method = "lm")
plot.sN.MAN.cv <- ggplot(data = mod.summ, aes(x = log(sN), y = log(MAN.cv)))+
  geom_point()+
  geom_smooth(method = "lm")
plot.sN.MAP.cv <- ggplot(data = mod.summ, aes(x = log(sN), y = log(MAP.cv)))+
  geom_point()+
  geom_smooth(method = "lm")
plot_grid(plot.sN.MAN, plot.sN.MAN.cv, plot.sN.MAP.cv, ncol = 3, nrow = 1, align = "hv")
#ggsave("cor.sN.svg", width = 16, height = 5, units = "cm")

#Generalized linear models on sN and R0
mod1 <- glm(sqrt(sN) ~ log(MAN.cv) * Clay * Sand + MAT * log(MAN) * log(MAP.cv)
            - MAT:log(MAN) - log(MAN.cv) - log(MAP.cv) - log(MAN):log(MAP.cv)
            - MAT:log(MAN):log(MAP.cv) - MAT - MAT:log(MAP.cv)+0
            , data = mod.summ[-c(14),])
summary(mod1)

mod.summ["R0p"] <- mod.summ$R0 + 0.5
mod2 <- glm(R0p ~ MAT * Sand * Clay + MAP * log(MAN.cv) * MAP.cv * log(MAN)
            - MAP:log(MAN.cv):MAP.cv - MAP:log(MAN.cv):MAP.cv:log(MAN)
            - MAT:Sand - Sand:Clay - Sand - log(MAN.cv):MAP.cv
            - log(MAN.cv):MAP.cv:log(MAN) - MAP.cv:log(MAN) - MAP.cv 
            - log(MAN) - MAT - MAP:log(MAN.cv):log(MAN) - MAP:log(MAN)
            - MAP - MAP:MAP.cv - log(MAN.cv):MAP - log(MAN.cv):log(MAN)
            - Clay - MAT:Clay +0 
            ,family = quasipoisson(link = log), data = mod.summ[-c(29,11),])
summary(mod2)

##Model evaluation
mod.test <- data.frame(mod.valid, predict.glm(mod1, newdata = mod.valid, se.fit = T))
mod.test$sN <- (mod.test$fit)^2
mod.test <- data.frame(mod.test, predict.glm(mod2, newdata = mod.valid, se.fit = T))
mod.test$R0 <- exp(mod.test$fit.1) - 0.5
mod.test$RN2O.pred <- mod.test$sN * mod.test$Ndepo + mod.test$R0
RN2O.obs <- mod.test$RN2O
RN2O.pred <- mod.test$RN2O.pred
summary(lm(log(RN2O.pred+1) ~ log(RN2O.obs+1)))
summary(lm(log(RN2O.pred+1) ~ log(RN2O.obs+1)+0))
RMSE(RN2O.pred, RN2O.obs)
cor.test(RN2O.pred, RN2O.obs, method = "s", use = "pairwise.complete.obs")

plot.check <- ggplot(data = mod.test, aes(x = log(RN2O+1), y = log(RN2O.pred+1)))+
  geom_point(size = 1, aes(color = Biome), alpha = 0.8)+
  geom_abline(linetype = "dashed")+
  geom_smooth(method = lm,
              formula = y ~ x + 0, color = "red", se = T, lwd = 0.5) +
  scale_color_brewer(palette = "Set2")+
  xlab("RN2O.obs (kgN ha-1 yr-1)") +
  ylab("RN2O.pred (kgN ha-1 yr-1)") +
  ylim(c(-1, 2))+
  xlim(c(-1, 2))
#ggsave("plotcheck2.svg", height = 6, width = 9, units = "cm")

##################PART II: Calculating sN of global forests######################
data7 <- read.csv("DataS7_GlobalForestEnvFactor.csv")
f.global <- cbind.data.frame(data7$ID, data7$Lon, data7$Lat, data7$MAP, data7$MAT, data7$MAP.cv, 
                             data7$MAT.cv, data7$MAN, data7$MAN.cv, data7$Clay, data7$Sand, data7$Biome, 
                             as.factor(data7$Country), as.factor(data7$CountryName), data7$Region, data7$percentForest)
colnames(f.global) <- c("ID", "Lon", "Lat", "MAP", "MAT", "MAP.cv", "MAT.cv", 
                        "MAN", "MAN.cv", "Clay", "Sand", "Biome", "Country", "CountryName", "Region", "percentForest")
f.global$Biome <- ordered(f.global$Biome, levels = c("Tropical", "Temperate", "Boreal"))
f.global$Area <- 50^2 * cos(abs(f.global$Lat)/180*3.14) *100 * f.global$percentForest/100   #forest cover raster size is 50km (in the equator)

##Calculating sN at the grid level
f.global.cal <- data.frame(f.global, predict.glm(mod1, newdata = f.global, se.fit = T))
f.global.cal$sN <- (f.global.cal$fit)^2
f.global.cal$sN.se <- abs(2 * f.global.cal$fit) * f.global.cal$se.fit
#write.csv(f.global.cal, "GlobalForest_sN.csv")  
###output csv file could be immported to ArcGIS software to produce a rasterized map of sN of global forests

##Latitudinal gradient of sN
sN.lat.mean <- f.global.cal[, c("Lat", "sN", "sN.se")] %>%
  filter(!is.na(sN)) %>%
  group_by(Lat) %>%
  summarize(across(everything(), mean, na.rm = T))
sN.lat.mean.roll <- cbind.data.frame(rollmean(sN.lat.mean$Lat, 2, na.rm = T), rollmean(sN.lat.mean$sN, 2, na.rm = T),
                       rollmean(sN.lat.mean$sN.se, 2, na.rm = T))
colnames(sN.lat.mean.roll) <- c("Lat", "sN", "sN.se")
ggplot(data = sN.lat.mean.roll, aes(Lat, sN)) +
  geom_line() +
  geom_ribbon(aes(ymin = sN - sN.se, ymax = sN + sN.se), color = "grey", alpha = 0.3) +
  scale_x_continuous(breaks = c(-60, -30, 0, 30, 60)) +
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.4, vjust = 1.1, size = 10),
    axis.ticks.y = element_line(),
    panel.grid.major.x = element_line(color = "grey"))
#ggsave("sN.lat.svg", width = 15, height = 6, units = "cm")

##Biome mean c1, c2 and calculation of sN
Ncycling.mean <- Ncycling %>%
  group_by(Lon, Lat, Biome) %>%
  summarise(across(where(is.numeric), mean, na.rm = T)) %>%
  ungroup()
Ncycl.summ <- Ncycling.mean %>%
  group_by(Biome) %>%
  summarize(c1.m = mean(c1, na.rm = T), c1.se = sd(c1, na.rm = T)/sqrt(length(Ncycling.mean$c1[!is.na(Ncycling.mean$c1)])),
            c2.m = mean(c2, na.rm = T), c2.se = sd(c2, na.rm = T)/sqrt(length(Ncycling.mean$c2[!is.na(Ncycling.mean$c2)])))
Ncycl.summ$c3 <- c(0.20, 0.19, 0.19)
Ncycl.summ <- Ncycl.summ %>%
  mutate(sN = c3*(c1.m - c2.m), sN.se = c3*sqrt((c1.se)^2 + (c2.se)^2))
plot.c1 <- ggplot(data = Ncycl.summ)+
  geom_col(aes(x = Biome, y = c1.m, fill = Biome, color = Biome), alpha = 0.5, show.legend = F)+
  geom_errorbar(aes(x = Biome, ymin = c1.m - c1.se, ymax = c1.m + c1.se), width = 0.1) +
  scale_fill_brewer(palette = "Set2")+
  scale_color_brewer(palette = "Set2")+
  ylim(c(0,0.6))+
  ylab("c1") +
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.4, vjust = 1.1, size = 10))
plot.c2 <- ggplot(data = Ncycl.summ)+
  geom_col(aes(x = Biome, y = c2.m, fill = Biome, color = Biome), alpha = 0.5, show.legend = F)+
  geom_errorbar(aes(x = Biome, ymin = c2.m - c2.se, ymax = c2.m + c2.se), width = 0.1) +
  scale_fill_brewer(palette = "Set2")+
  scale_color_brewer(palette = "Set2")+
  ylab("c2") +
  ylim(c(0, 0.5)) +
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.4, vjust = 1.1, size = 10))
sN.modeled <- f.global.cal[, c("Biome", "sN", "sN.se")] %>%
  group_by(Biome) %>%
  summarise(across(everything(), mean, na.rm =T))
sN.calculated <- Ncycl.summ[, c("Biome", "sN", "sN.se")]
sN.all <- merge.data.frame(sN.calculated, sN.modeled, by = "Biome")
colnames(sN.all) <- c("Biome", "Cal.m", "Cal.se", "Mod.m", "Mod.se")
cor.test(sN.all$Cal.m, sN.all$Mod.m)
summary(lm(sN.all$Mod.m ~ sN.all$Cal.m+0))
plot.sN <- ggplot(data = sN.all, aes(x = Cal.m, y= Mod.m))+
  geom_abline(linetype = "dashed")+
#  geom_smooth(method = "lm", formula = "y~x+0", color = "red", se = F, fullrange =T)+
  geom_errorbar(aes(ymin = Mod.m - Mod.se, ymax = Mod.m + Mod.se, color = Biome), width=0.005, show.legend = F)+
  geom_errorbar(aes(xmin = Cal.m - Cal.se, xmax = Cal.m + Cal.se, color = Biome), width=0.003, show.legend = F)+
  geom_point(shape = c("B", "T", "O"), size = 3)+
  scale_color_brewer(palette = "Set2")+
  ylim(c(0, 0.03))
plot.grid1 <- plot_grid(plot.c1, plot.c2, align = "v", ncol = 1, nrow = 2)
plot_grid(plot.grid1, plot.sN, align = "none", rel_widths = c(3.5,6.5))  
#ggsave("meansN.svg", width = 12, height = 8, units = "cm")

##Violin plot of sN by biome
plot.biome <- ggplot(data = f.global.cal, aes(x = Biome, y = sN))+
  geom_violin(stat = "ydensity", aes(fill = Biome, color = Biome), alpha = 0.5, scale = "width", show.legend = F)+
  geom_boxplot(width = 0.2, color = "grey50", outlier.shape = NA)+
  stat_summary(fun = "mean", geom = "point", shape = 23, size = 3)+
  scale_fill_brewer(palette = "Set2")+
  scale_color_brewer(palette = "Set2")+
#  scale_y_continuous(limits = c(-10,5))+
#  ylab("sN(kgN2O-N kgN-1)") +
  coord_flip()+
  theme(
    axis.text.x = element_text(size = 8), 
    axis.text.y = element_text(size = 8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
  )

##Relative importance of environmental factors
s <- levels(f.global.cal$Biome)
output <- data.frame(c(1:8))
for (i in 1:length(s)){
  temp1 <- subset(f.global.cal, f.global.cal$Biome == s[i])
  relimp.mod1 <- calc.relimp(sN ~ MAT + MAP + MAN + Clay + Sand + MAT.cv + MAP.cv + MAN.cv,
                      type = c("lmg"),
                      data = temp1, rank = F, rela = T)
  relim <- as.data.frame(relimp.mod1@lmg)
  colnames(relim) <- s[i]
  output <- cbind(output, relim)
}
output[,1] <- row.names(output)
output2 <- output %>%
  gather("Tropical", "Temperate", "Boreal", key = "Biome", value = "value")
output2 <- as.data.frame(output2)
colnames(output2) <- c("Factor", "Biome", "value")
output2$Factor <- ordered(output2$Factor, levels = c("MAP", "MAP.cv", "MAT", "MAT.cv", "MAN", "MAN.cv", "Sand", "Clay"))
output2$Biome <- ordered(output2$Biome, levels = (c("Tropical", "Temperate", "Boreal")))
plot.factor <- ggplot(data = output2, aes(x = Biome, y = value, fill = Factor))+
  geom_bar(stat = "identity", position = "fill", alpha = 0.8, show.legend = F)+
  scale_fill_manual(values = rev(brewer.pal(11, "Spectral")))+
  coord_flip()+
  theme(
    axis.text.x = element_text(size = 8), 
    axis.text.y = element_text(size = 8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
  )
plot_grid(plot.factor, plot.biome, ncol = 2, nrow = 1)
#ggsave("plot2.svg", width = 16, height = 3.5, units = "cm")

#Comparing means of model-estimated sN of different biomes (bootstrap method)
boot.mean.diff <- function(data, category, value, R){
  s <- unique(data[category])
  s <- t(s)
  pair <- combn(s, 2)
  mean.diff <- matrix(data = NA, nrow = R, ncol = ncol(pair))
  diff.table <- {}
  for(i in 1:ncol(pair)){
    for (j in 1:R){
      data.cat1 <- subset(data[,value], data[,category] == pair[1,i])
      data.cat2 <- subset(data[,value], data[,category] == pair[2,i])
      samp1 <- sample(t(data.cat1), 68, replace = T)
      samp2 <- sample(t(data.cat2), length(t(data.cat2)), replace = T)
      mean.diff[j, i] <- mean(samp1, na.rm = T) - mean(samp2, na.rm = T)
    }
    n.ex <- min(length(subset(mean.diff[,i], mean.diff[,i]<0)),
                length(subset(mean.diff[,i], mean.diff[,i]>0)))
    if(n.ex == 0){
      pvalue <- paste0("<" , 1/(R+1)*2)
    }
    else{pvalue <- (n.ex + 1)/(R+1)*2}
    diff.table <- rbind.data.frame(diff.table, cbind.data.frame(paste0(pair[1,i], " - ", pair[2,i]),
                                                                t(quantile(mean.diff[,i], c(0.0005, 0.005, 0.025, 0.975, 0.995, 0.9995))),
                                                                pvalue))
  }
  colnames(diff.table) <- c("Category", "p0.0005", "p0.005", "p0.025", "p0.975", "p0.995", "p0.9995", "pvalue")
  plot.out <- ggplot(data = diff.table, aes(y = Category)) +
    geom_errorbar(aes(xmin = p0.0005, xmax = p0.9995), width = 0.3) +
    geom_errorbar(aes(xmin = p0.005, xmax = p0.995), size = 10, width=0, color = "yellow", alpha = 0.5) +
    geom_errorbar(aes(xmin = p0.025, xmax = p0.975), size = 5, width=0, color = "red", alpha = 0.5) +
    geom_vline(xintercept = 0, linetype = "dashed")+
    xlim(c(min(diff.table$p0.0005) - (max(diff.table$p0.9995) - min(diff.table$p0.0005))*0.5, 
           max(diff.table$p0.9995)+ (max(diff.table$p0.9995) - min(diff.table$p0.0005))*0.5))+
    xlab(paste0("Differences in mean ", value))+
    theme(
      panel.grid.major.y = element_line(),
      axis.title.y = element_blank())
  list(diff.table, plot.out)
}
boot.mean.diff(f.global.cal, "Biome", "sN", 4999)
#ggsave("meandiff.svg", width = 16, height = 16, units = "cm")

##################PART III: Using sN to classify N-limited and N-saturated forests######################
data4 <- read_xlsx("DataS4_Nleach.xlsx")
Nleach <- cbind.data.frame(data4$ID, data4$Lon, data4$Lat, data4$`N status`, data4$`MAP(mm)`, data4$MAPcv, data4$`MAT(K)`, data4$MATcv,
                           data4$`MAN(kgN ha-1 yr-1)`, data4$MANcv, data4$`Clay(%)`, data4$`Sand(%)`, data4$Biome, data4$Region, 
                           data4$`SOC(%)`, data4$`CEC(cmol/kg)`, data4$pH, as.numeric(data4$CtoN), data4$`Ndepo(kgN ha-1 yr-1)`)
colnames(Nleach) <- c("ID", "Lon", "Lat", "Nstat_obs", "MAP", "MAP.cv", "MAT", "MAT.cv", "MAN", "MAN.cv", "Clay", "Sand",
                      "Biome", "Region", "SOC", "CEC", "pH", "rCN", "Ndepo")
Nleach$sN <- predict.glm(mod1, newdata = Nleach)^2
Nleach$Nstat_obs <- as.factor(Nleach$Nstat_obs)
Nleach$Biome <- ordered(Nleach$Biome, levels = c("Tropical", "Temperate", "Boreal"))
data5 <- read_xlsx("DataS5_NuLi.xlsx", col_types = c(rep("guess", 16), rep("numeric",5), "guess"))
NuLi <- cbind.data.frame(data5$ID, data5$Lat, data5$Lon, data5$`Field tested N limitation`, data5$Vegetype2, 
                         data5$`MAP(mm)`, data5$MAPcv, data5$`MAT(K)`, data5$MATcv, data5$MAN, data5$MANcv, data5$`Clay(%)`, data5$`Sand(%)`, 
                         data5$Biome, data5$Region, data5$`SOC(%)`, data5$pH, data5$`CEC(cmol/kg)`, data5$`Ndepo(kgN ha-1 yr-1)`, data5$CtoN, data5$Reference)
colnames(NuLi) <- c("ID", "Lat", "Lon", "Nstat_obs", "Vegetype", "MAP", "MAP.cv", "MAT", "MAT.cv", "MAN", "MAN.cv", "Clay", "Sand",
                    "Biome", "Region", "SOC", "pH", "CEC", "Ndepo", "rCN", "Reference")
NuLi$sN <- predict.glm(mod1, newdata = NuLi)^2
NuLi$Nstat_obs <- ordered(NuLi$Nstat_obs, levels = c("Limited", "Saturated"), labels = c("limited", "saturated"))
NuLi$Biome <- ordered(NuLi$Biome, levels = c("Tropical", "Temperate", "Boreal"))

##Classifying N-saturated and N-limited forests using Nleach dataset
Nleach$Lon <- round(Nleach$Lon, 2)
Nleach$Lat <- round(Nleach$Lat, 2)
Nleach.stat <- Nleach %>%
  filter(!is.na(sN)) %>%
  group_by(Lon, Lat, Biome, Region) %>%
  summarize(Nstat_obs = majorityVote(Nstat_obs)$majority, sN = mean(sN), 
            SOC = mean(SOC), pH = mean(pH), CEC = mean(CEC), Clay = mean(Clay), rCN = mean(rCN),
            Ndepo = mean(Ndepo))
plot.Nleach <- ggplot(data = Nleach.stat, aes(x = Nstat_obs, y = sN)) +
  geom_boxplot(width=0.4, color="grey30", fill = NA, outlier.alpha=0)+
  geom_violin(width=0.8, fill = NA, scale = "count") +
  geom_jitter(aes(color= Nstat_obs), width = 0.15, size=2.5, alpha=0.5, show.legend = F) +
  stat_summary(fun= mean, geom="point", shape=20, size=6, color="dark red", fill="dark red")+
  stat_summary(fun = mean, geom = "label", col = "black", size = 2.5,
               hjust = -0.6, vjust = 0.5, 
               aes(label = paste("mean =", round(..y.., digits = 3))))+
  scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
  ylim(c(0, 0.035))+
  scale_x_discrete(breaks = c("limited", "saturated"), 
                   labels = c(paste0("limited\n(n =", nrow(Nleach.stat[Nleach.stat$Nstat_obs == "limited",]), ")"), 
                              paste0("saturated\n(n =", nrow(Nleach.stat[Nleach.stat$Nstat_obs == "saturated",]),")")))+
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.4, vjust = 1.1, size = 8), # Angle x-axis text
    axis.text.y = element_text(size = 8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
  )
plot.biome.Nleach <- ggplot(data = Nleach.stat, aes(x = Nstat_obs, y = sN, fill = Nstat_obs))+
  geom_boxplot(alpha = 0.3, show.legend = F)+
  stat_summary(aes(color = Nstat_obs), fun = "mean", geom = "point", position = "identity", 
               shape = 20, size = 5, show.legend = F)+
  facet_wrap(~ Biome)+
  scale_fill_manual(values = c("#4286f4", "#ec2F4B"))+
  scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
  scale_x_discrete(breaks = c("limited", "saturated"), 
                   labels = c("limited", "saturated"))+
  ylim(c(0, 0.035))+
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )
plot_grid(plot.Nleach, plot.biome.Nleach, nrow =2, ncol = 1, rel_heights = c(6,4))
#ggsave("Nleach_plot.svg", width = 12, height = 12, units = "cm")
boot.mean.diff(Nleach.stat, "Nstat_obs", "sN", 19999)
boot.mean.diff(Nleach.stat[Nleach.stat$Biome == "Tropical",], "Nstat_obs", "sN", 19999)
boot.mean.diff(Nleach.stat[Nleach.stat$Biome == "Temperate",], "Nstat_obs", "sN", 19999)
boot.mean.diff(Nleach.stat[Nleach.stat$Biome == "Boreal",], "Nstat_obs", "sN", 19999)

ThreshDetect <- function(dat, VarObs, VarLow, VarHigh, VarClas, precision){
  b1 <- median(t(subset(dat[, VarClas], dat[, VarObs] == VarLow)))
  b2 <- median(t(subset(dat[, VarClas], dat[, VarObs] == VarHigh)))
  tmp <- {}
  for(i in floor(min(b1,b2)/precision) : ceiling(max(b1,b2)/precision)){
    dat.TN <- subset(dat, dat[, VarObs] == VarLow & dat[, VarClas] <= i*precision)
    dat.FP <- subset(dat, dat[, VarObs] == VarLow & dat[, VarClas] > i*precision)
    dat.TP <- subset(dat, dat[, VarObs] == VarHigh & dat[, VarClas] > i*precision)
    dat.FN <- subset(dat, dat[, VarObs] == VarHigh & dat[, VarClas] <= i*precision)
    dat.comp <- rbind.data.frame(dat.TN, dat.FP, dat.TP, dat.FN)
    dat.comp$check <- c(rep("TN", nrow(dat.TN)), rep("FP", nrow(dat.FP)),
                        rep("TP", nrow(dat.TP)), rep("FN", nrow(dat.FN)))
    #constructing binary accuracy table
    n.TN <- length(dat.comp$check[dat.comp$check == "TN"])
    n.FP <- length(dat.comp$check[dat.comp$check == "FP"])
    n.TP <- length(dat.comp$check[dat.comp$check == "TP"])
    n.FN <- length(dat.comp$check[dat.comp$check == "FN"])
    tmp <- rbind.data.frame(tmp, c(i*precision, n.TN/(n.TN + n.FP), n.TP/(n.TP + n.FN),
                                   (n.TN + n.TP)/(n.TP + n.FN + n.TN + n.FP)))
  }
  colnames(tmp) <- c("cutoff", "Acc_Nlimited", "Acc_Nsaturated", "Acc_all")
  Nthreshold <- tmp$cutoff[tmp$Acc_all == max(tmp$Acc_all, na.rm=T) & !is.na(tmp$Acc_all)]
  #  Nthreshold
  list(Nthreshold, tmp)
}
R <- 5000
Nleach.thre <- {}
pb <- tkProgressBar("Progress bar","%", 0, 100) 
for(i in 1:R){
  dat1 <- subset(Nleach.stat, Nleach.stat$Nstat_obs == "limited")
  dat2 <- subset(Nleach.stat, Nleach.stat$Nstat_obs == "saturated")
  dat.cmb <- rbind.data.frame(dat1[sample(1:nrow(dat1), 10, replace = T),], 
                              dat2[sample(1:nrow(dat2), 10, replace = T),])
  thre1 <- ThreshDetect(dat.cmb, "Nstat_obs", "limited", "saturated", "sN", 0.0001)
  Nleach.thre <- c(Nleach.thre, thre1[[1]])
  info<- sprintf("%d%% completed", round(i*100/R))  
  setTkProgressBar(pb, i*100/R, sprintf("Progress (%s)", info),info)  
}

##Classifying N-saturated and N-limited forests using NuLi dataset
NuLi.stat <- NuLi %>%
  filter(Vegetype == "forest" & !is.na(sN))
plot.NuLi <- ggplot(data = NuLi.stat, aes(x = Nstat_obs, y = sN)) +
  geom_boxplot(width=0.4, color="grey30", fill = NA, outlier.alpha=0)+
  geom_violin(width=0.8, fill = NA, scale = "count") +
  geom_jitter(aes(color= Nstat_obs), width = 0.15, size=2.5, alpha=0.5, show.legend = F) +
  stat_summary(fun= mean, geom="point", shape=20, size=6, color="dark red", fill="dark red")+
  stat_summary(fun = mean, geom = "label", col = "black", size = 2.5,
               hjust = -0.6, vjust = 0.5, 
               aes(label = paste("mean =", round(..y.., digits = 3))))+
  scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
  ylim(c(0, 0.03))+
  scale_x_discrete(breaks = c("limited", "saturated"), 
                   labels = c(paste0("limited\n(n =", nrow(NuLi.stat[NuLi.stat$Nstat_obs == "limited",]), ")"), 
                              paste0("saturated\n(n =", nrow(NuLi.stat[NuLi.stat$Nstat_obs == "saturated",]),")")))+
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.4, vjust = 1.1, size = 8), # Angle x-axis text
    axis.text.y = element_text(size = 8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
  )
plot.biome.NuLi <- ggplot(data = NuLi.stat, aes(x = Nstat_obs, y = sN, fill = Nstat_obs))+
  geom_boxplot(alpha = 0.3, show.legend = F)+
  stat_summary(aes(color = Nstat_obs), fun = "mean", geom = "point", position = "identity", 
               shape = 20, size = 5, show.legend = F)+
  facet_wrap(~ Biome)+
  scale_fill_manual(values = c("#4286f4", "#ec2F4B"))+
  scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
  scale_x_discrete(breaks = c("limited", "saturated"), 
                   labels = c("limited", "saturated"))+
  ylim(c(0, 0.03))+
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )
plot_grid(plot.NuLi, plot.biome.NuLi, nrow =2, ncol = 1, rel_heights = c(6,4))
#ggsave("NuLi_plot.svg", width = 12, height = 12, units = "cm")
boot.mean.diff(NuLi.stat, "Nstat_obs", "sN", 19999)
boot.mean.diff(NuLi.stat[NuLi.stat$Biome == "Tropical",], "Nstat_obs", "sN", 19999)
boot.mean.diff(NuLi.stat[NuLi.stat$Biome == "Temperate",], "Nstat_obs", "sN", 19999)
boot.mean.diff(NuLi.stat[NuLi.stat$Biome == "Boreal",], "Nstat_obs", "sN", 19999)
R <- 5000
NuLi.thre <- {}
pb <- tkProgressBar("Progress bar","%", 0, 100) 
for(i in 1:R){
  dat1 <- subset(NuLi.stat, NuLi.stat$Nstat_obs == "limited")
  dat2 <- subset(NuLi.stat, NuLi.stat$Nstat_obs == "saturated")
  dat.cmb <- rbind.data.frame(dat1[sample(1:nrow(dat1), 10, replace = T),], 
                              dat2[sample(1:nrow(dat2), 10, replace = T),])
  Thre2 <- ThreshDetect(dat.cmb, "Nstat_obs", "limited", "saturated", "sN", 0.0001)
  NuLi.thre <- c(NuLi.thre, Thre2[[1]])
  info<- sprintf("%d%% completed", round(i*100/R))  
  setTkProgressBar(pb, i*100/R, sprintf("Progress (%s)", info),info)  
}
sd(Nleach.thre)
sd(NuLi.thre)
Thre.comp <- cbind.data.frame(c(Nleach.thre, NuLi.thre), 
                              c(rep("Nleach", length(Nleach.thre)), rep("NuLi", length(NuLi.thre))))
colnames(Thre.comp) <- c("Thre", "dataset")
Thre.comp.summ <- Thre.comp %>%
  group_by(dataset) %>%
  summarise(means = mean(Thre), sds = sd(Thre))
Thre.comp.summ$ys <- c(250, 280)
plot.thre <- ggplot(data = Thre.comp, aes(x= Thre, group= dataset, fill= dataset, color = dataset)) +
  geom_density(aes(y = ..density..), adjust=1.5, alpha= 0.2, bw = 0.0005) +
  scale_fill_brewer(palette = "Set2") + 
  scale_color_brewer(palette = "Set2") +
  ylim(c(0, 300))+
  geom_point(data = Thre.comp.summ, aes(x = means, y = ys)) +
  geom_segment(aes(y = Thre.comp.summ$ys[1], yend = Thre.comp.summ$ys[1], x = Thre.comp.summ$means[1] - Thre.comp.summ$sds[1],
                   xend = Thre.comp.summ$means[1] + Thre.comp.summ$sds[1]), color = palette("Set2")[1])+
  geom_segment(aes(y = Thre.comp.summ$ys[2], yend = Thre.comp.summ$ys[2], x = Thre.comp.summ$means[2] - Thre.comp.summ$sds[2],
                   xend = Thre.comp.summ$means[2] + Thre.comp.summ$sds[2]), color = palette("Set2")[2])+
  theme(
    axis.title.x = element_blank(),
  )
boot.mean.diff(Thre.comp, "dataset", "Thre", 1999)
#ggsave("plot.svg", width = 10, height = 6, units = "cm")

##Combining datasets Nleach and NuLi, and calculate the threshold value
Nstat.all <- rbind.data.frame(Nleach.stat, NuLi.stat[, colnames(Nleach.stat)])
Nstat.all <- Nstat.all %>%
  group_by(Lon, Lat, Biome, Region) %>%
  summarize(Nstat_obs = majorityVote(Nstat_obs)$majority, sN = mean(sN), 
            SOC = mean(SOC), pH = mean(pH), CEC = mean(CEC), Clay = mean(Clay), rCN = mean(rCN),
            Ndepo = mean(Ndepo))
plot.Nstat <- ggplot(data = Nstat.all, aes(x = Nstat_obs, y = sN)) +
  geom_boxplot(width=0.4, color="grey30", fill = NA, outlier.alpha=0)+
  geom_violin(width=0.8, fill = NA, scale = "count") +
  geom_jitter(aes(color= Nstat_obs), width = 0.15, size=2.5, alpha=0.5, show.legend = F) +
  stat_summary(fun= mean, geom="point", shape=20, size=6, color="dark red", fill="dark red")+
  stat_summary(fun = mean, geom = "label", col = "black", size = 2.5,
               hjust = -0.6, vjust = 0.5, 
               aes(label = paste("mean =", round(..y.., digits = 3))))+
  scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
  ylim(c(0, 0.035))+
  scale_x_discrete(breaks = c("limited", "saturated"), 
                   labels = c(paste0("limited\n(n =", nrow(Nstat.all[Nstat.all$Nstat_obs == "limited",]), ")"), 
                              paste0("saturated\n(n =", nrow(Nstat.all[Nstat.all$Nstat_obs == "saturated",]),")")))+
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.4, vjust = 1.1, size = 8), # Angle x-axis text
    axis.text.y = element_text(size = 8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
  )
plot.biome.Nstat <- ggplot(data = Nstat.all, aes(x = Nstat_obs, y = sN, fill = Nstat_obs))+
  geom_boxplot(alpha = 0.3, show.legend = F)+
  stat_summary(aes(color = Nstat_obs), fun = "mean", geom = "point", position = "identity", 
               shape = 20, size = 5, show.legend = F)+
  facet_wrap(~ Biome)+
  scale_fill_manual(values = c("#4286f4", "#ec2F4B"))+
  scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
  scale_x_discrete(breaks = c("limited", "saturated"), 
                   labels = c("limited", "saturated"))+
  ylim(c(0, 0.035))+
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.8, size = 8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )
plot.region.Nstat <- subset(Nstat.all, Nstat.all$Region == c("Western Europe", "East Asia", "North America")) %>%
  ggplot(aes(x = Nstat_obs, y = sN, fill = Nstat_obs))+
  geom_boxplot(alpha = 0.3, show.legend = F)+
  stat_summary(aes(color = Nstat_obs), fun = "mean", geom = "point", position = "identity", 
               shape = 20, size = 5, show.legend = F)+
  facet_wrap(~ Region)+
  scale_fill_manual(values = c("#4286f4", "#ec2F4B"))+
  scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
  ylim(c(0, 0.035))+
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.8, size = 8),
    axis.title.x = element_blank(),
    axis.title.y = element_blank()
  )
plot_grid(plot.Nstat, plot.biome.Nstat, plot.region.Nstat, nrow =3, ncol = 1, rel_heights = c(6,4,4))
#ggsave("Nstat_plot2.svg", width = 12, height = 16, units = "cm")
boot.mean.diff(Nstat.all, "Nstat_obs", "sN", 19999)
boot.mean.diff(Nstat.all[Nstat.all$Biome == "Tropical",], "Nstat_obs", "sN", 19999)
boot.mean.diff(Nstat.all[Nstat.all$Biome == "Temperate",], "Nstat_obs", "sN", 19999)
boot.mean.diff(Nstat.all[Nstat.all$Biome == "Boreal",], "Nstat_obs", "sN", 19999)
boot.mean.diff(Nstat.all[Nstat.all$Region == "East Asia",], "Nstat_obs", "sN", 19999)
boot.mean.diff(Nstat.all[Nstat.all$Region == "North America",], "Nstat_obs", "sN", 19999)
boot.mean.diff(Nstat.all[Nstat.all$Region == "Western Europe",], "Nstat_obs", "sN", 19999)

Nstat.thre <- {}
R <- 5000
pb <- tkProgressBar("Progress bar","%", 0, 100) 
for(i in 1:R){
  dat1 <- subset(Nstat.all, Nstat.all$Nstat_obs == "limited")
  dat2 <- subset(Nstat.all, Nstat.all$Nstat_obs == "saturated")
  dat.cmb <- rbind.data.frame(dat1[sample(1:nrow(dat1), 10, replace = T),], 
                              dat2[sample(1:nrow(dat2), 10, replace = T),])
  Thre3 <- ThreshDetect(dat.cmb, "Nstat_obs", "limited", "saturated", "sN", 0.0001)
  Nstat.thre <- c(Nstat.thre, Thre3[[1]])
  info<- sprintf("%d%% completed", round(i*100/R))  
  setTkProgressBar(pb, i*100/R, sprintf("Progress (%s)", info),info)  
}
thre <- majorityVote(Nstat.thre)$majority  ##the mode of the detected thresholds was used as the final threshold
Nstat.thre <- as.data.frame(Nstat.thre)
#write.csv(Nstat.thre, "Nstat.thre230206.csv")
ggplot(data = Nstat.thre, aes(x = Nstat.thre))+
  geom_density(aes(y = ..density..),color = "Grey30", fill = "Grey70", alpha= 0.2, bw = 0.0005) +
  geom_vline(aes(xintercept = as.numeric(thre)), color = "red", linetype = "dashed") 
#ggsave("thre3.svg", width = 12, height = 6, units = "cm")

##Contingency table for the accuracy of classification
Nstat.all$Nstat_pred[Nstat.all$sN <= thre] <- "Limited"
Nstat.all$Nstat_pred[Nstat.all$sN > thre] <- "Saturated"
tb.con <- table(Nstat.all[, c("Nstat_pred", "Nstat_obs")])
acc.limited <- tb.con[1,1] / sum(tb.con[,1])
acc.saturated <- tb.con[2,2] / sum(tb.con[,2])
acc.all <- (tb.con[1,1] + tb.con[2,2]) / (sum(tb.con[,1]) + sum(tb.con[,2]))
Nstat.all$Nstat_pred <- NA
Nstat.all$Nstat_pred[Nstat.all$sN <= thre & Nstat.all$Region == "East Asia"] <- "limited"
Nstat.all$Nstat_pred[Nstat.all$sN > thre & Nstat.all$Region == "East Asia"] <- "saturated"
tb.con.ea <- table(Nstat.all[, c("Nstat_pred", "Nstat_obs")])
acc.ea <- (tb.con.ea[1,1] + tb.con.ea[2,2]) / (sum(tb.con.ea[,1]) + sum(tb.con.ea[,2]))
Nstat.all$Nstat_pred <- NA
Nstat.all$Nstat_pred[Nstat.all$sN <= thre & Nstat.all$Region == "North America"] <- "limited"
Nstat.all$Nstat_pred[Nstat.all$sN > thre & Nstat.all$Region == "North America"] <- "saturated"
tb.con.na <- table(Nstat.all[, c("Nstat_pred", "Nstat_obs")])
acc.na <- (tb.con.na[1,1] + tb.con.na[2,2]) / (sum(tb.con.na[,1]) + sum(tb.con.na[,2]))
Nstat.all$Nstat_pred <- NA
Nstat.all$Nstat_pred[Nstat.all$sN <= thre & Nstat.all$Region == "Western Europe"] <- "limited"
Nstat.all$Nstat_pred[Nstat.all$sN > thre & Nstat.all$Region == "Western Europe"] <- "saturated"
tb.con.we <- table(Nstat.all[, c("Nstat_pred", "Nstat_obs")])
acc.we <- (tb.con.we[1,1] + tb.con.we[2,2]) / (sum(tb.con.we[,1]) + sum(tb.con.we[,2]))

##################PART IV: N saturation ratio of global countries/regions######################
f.global$Region[f.global$Region == "NA"] <- NA
f.global$Region[f.global$Country == "NA"] <- NA
f.global$sN <- predict(mod1, newdata = f.global)^2
data6 <- read_xlsx("DataS6_GDPpc.xlsx", skip = 3)
Econ <- data6
Econ$GDP.pc <- apply(Econ[ , as.character(1985:2015)], 1, mean, na.rm= T)

##Summarizing limited and saturated forest area by region
f.global$Nstat[f.global$sN <= thre] <- "limited" 
f.global$Nstat[f.global$sN > thre] <- "saturated"
f.global.Nstat <- f.global %>%
  filter(!is.na(sN)) %>%
  group_by(Nstat) %>%
  summarize(f.area = sum(Area)) %>%
  pivot_wider(names_from = Nstat, values_from = f.area) %>%
  mutate(pSatur = saturated / (saturated + limited))
f.global.Biome <- f.global %>%
  filter(!is.na(sN)) %>%
  group_by(Nstat, Biome) %>%
  summarize(f.area = sum(Area)) %>%
  pivot_wider(names_from = Nstat, values_from = f.area) %>%
  mutate(pSatur = saturated / (saturated + limited))
f.global.region <- f.global %>%
  filter(!is.na(Region), !is.na(sN)) %>%
  group_by(Region, Nstat) %>%
  summarize(f.area = sum(Area))
colnames(f.global.region) <- c("Region", "Nstat", "f.area")
plot.area.region <- ggplot(f.global.region, aes(x = "", y = f.area, fill = Nstat)) +
  geom_bar(stat="identity", position = "fill", width=1, alpha = 0.5, show.legend = F) +
    facet_wrap(~ Region)+
  scale_fill_manual(values = c("#4286f4", "#ec2F4B"))+
  coord_polar("y", start=0) +
  theme(
    axis.text.x = element_blank(),
    axis.title.x = element_blank(),
    axis.title.y = element_blank(),
  )
#ggsave("area_global.svg", height = 10, width = 10, units = "cm")
f.global.region.pivot <- f.global %>%
  filter(!is.na(Region), !is.na(sN)) %>%
  group_by(Region, Nstat) %>%
  summarize(f.area = sum(Area)) %>%
  pivot_wider(names_from = Nstat, values_from = f.area) %>%
  mutate(pSatur = saturated / (saturated + limited))

##Summarizing N-limited and N-saturated forest area by country
f.global.country <- f.global %>%
  filter(!is.na(Country), !is.na(sN)) %>%
  group_by(Country, Nstat) %>%
  summarize(f.area = sum(Area)) %>%
  pivot_wider(names_from = Nstat, values_from = f.area, values_fill = 0) %>%
  mutate(Area = sum(limited, saturated, na.rm = T)/10^6,
         pSatur = saturated / (saturated + limited))
f.global.country.biome <- f.global %>%
  filter(!is.na(Country), !is.na(sN)) %>%
  group_by(Country, Biome) %>%
  summarise(Area = sum(Area, na.rm = T)) %>%
  pivot_wider(names_from = "Biome", values_from = "Area", values_fill = 0)
for (i in 1:nrow(f.global.country.biome)){
  if((f.global.country.biome$Temperate[i] > f.global.country.biome$Tropical[i])
     & (f.global.country.biome$Temperate[i] > f.global.country.biome$Boreal[i])){
    f.global.country.biome$Biome[i] <- "Temperate"
  }
  else if(f.global.country.biome$Tropical[i] > f.global.country.biome$Boreal[i]){
    f.global.country.biome$Biome[i] <- "Tropical"
  }
  else{f.global.country.biome$Biome[i] <- "Boreal"}
}

f.global.country.pivot <- f.global.country %>%
  pivot_longer(cols = c("limited", "saturated"), names_to = "Nstat", values_to = "f.area")
f.global.country2 <- merge.data.frame(f.global.country, f.global.country.biome, by = "Country")
f.global.country3 <- merge.data.frame(f.global.country2, Econ[, c("Country Code", "GDP.pc")], 
                                      by.x = "Country", by.y = "Country Code") %>%
  filter(!is.na(GDP.pc))

cor.test(f.global.country3[f.global.country3$Biome == "Tropical",]$pSatur,
         log(f.global.country3[f.global.country3$Biome == "Tropical",]$GDP.pc))
cor.test(f.global.country3[f.global.country3$Biome == "Temperate",]$pSatur,
         log(f.global.country3[f.global.country3$Biome == "Temperate",]$GDP.pc))
cor.test(f.global.country3[f.global.country3$Biome == "Boreal",]$pSatur,
         log(f.global.country3[f.global.country3$Biome == "Boreal",]$GDP.pc))

summary(lm(pSatur ~ log(GDP.pc), data = f.global.country3[f.global.country3$Biome == "Tropical",]))
summary(lm(pSatur ~ log(GDP.pc), data = f.global.country3[f.global.country3$Biome == "Temperate",]))
summary(lm(pSatur ~ log(GDP.pc), data = f.global.country3[f.global.country3$Biome == "Boreal",]))

plot.GDP.tr <- f.global.country3 %>%
  filter(Biome == "Tropical") %>%
  ggplot(aes(x = log(GDP.pc), y = pSatur*100)) +
  geom_point(color = palette("Set2")[1], alpha = 0.3)+
  geom_smooth(method = "lm", color = palette("Set2")[1], fill = palette("Set2")[1])+
  ylim(c(0,120))+
  xlim(c(6, 11))+
#  geom_text(data = f.global.country3 %>% filter(Biome == "Tropical", pSatur>0, pSatur <0.4, log(GDP.pc) >8),
#            aes(label= Country), size = 2, nudge_x = 0.25, nudge_y = 0.25, check_overlap = F) +
  theme(
    axis.text.x = element_text(hjust = 0.5, vjust = 0.5, size = 8) # Angle x-axis text
  )
plot.GDP.te <- f.global.country3 %>%
  filter(Biome == "Temperate") %>%
  ggplot(aes(x = log(GDP.pc), y = pSatur*100)) +
  geom_point(color = palette("Set2")[2], alpha = 0.3)+
  geom_smooth(method = "lm", color = palette("Set2")[2], fill = palette("Set2")[2])+
  ylim(c(0,120))+
  xlim(c(6, 11))+
#  geom_text(data = f.global.country3 %>% filter(Biome == "Temperate", pSatur >0.8, log(GDP.pc) >10),
#            aes(label= Country), size = 2, nudge_x = 0.25, nudge_y = 0.25, check_overlap = F) +
  theme(
    axis.text.x = element_text(hjust = 0.5, vjust = 0.5, size = 8) # Angle x-axis text
  )
plot.GDP.bo <- f.global.country3 %>%
  filter(Biome == "Boreal") %>%
  ggplot(aes(x = log(GDP.pc), y = pSatur*100)) +
  geom_point(color = palette("Set2")[3], alpha = 0.3)+
  geom_smooth(method = "lm", color = palette("Set2")[3], fill = palette("Set2")[3])+
  ylim(c(0,120))+
  xlim(c(6, 11)) +
#  geom_text(data = f.global.country3 %>% filter(Biome == "Boreal", pSatur >0.7, log(GDP.pc) >10),
#            aes(label= Country), size = 2, nudge_x = 0.25, nudge_y = 0.25, check_overlap = F) +
  theme(
    axis.text.x = element_text(hjust = 0.5, vjust = 0.5, size = 8) # Angle x-axis text
  )
plot_grid(plot.GDP.bo, plot.GDP.te, plot.GDP.tr, ncol = 1, nrow = 3, align = "hv")
#ggsave("econ22.svg", height = 13, width = 5, units = "cm")

##Using other indicators to classify the N-limited and N-saturated forests?
clas.biome <- function(data, indicator, R){
  Data <- data[, c("Nstat_obs", indicator, "Biome")]
  colnames(Data) <- c("Nstat_obs", "indicator", "Biome")
  n.limited <- length(Data$indicator[Data$Nstat_obs == "limited"])
  n.saturated <- length(Data$indicator[Data$Nstat_obs == "saturated"])
  plot.data <- ggplot(data = Data, aes(x = Nstat_obs, y = indicator)) +
    geom_boxplot(width=0.4, color="grey30", fill = NA, outlier.alpha=0)+
    geom_violin(width=0.8, fill = NA, scale = "count") +
    geom_jitter(aes(color= Nstat_obs), width = 0.15, size=2.5, alpha=0.5, show.legend = F) +
    stat_summary(fun= mean, geom="point", shape=20, size=6, color="dark red", fill="dark red")+
    stat_summary(fun = mean, geom = "label", col = "black", size = 2.5,
                 hjust = -0.6, vjust = 0.5, 
                 aes(label = paste("mean =", round(..y.., digits = 3))))+
    scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
    #    ylim(c(0, 0.035))+
    scale_x_discrete(breaks = c("limited", "saturated"), 
                     labels = c(paste0("limited\n(n =", n.limited, ")"), 
                                paste0("saturated\n(n =", n.saturated,")")))+
    theme(
      axis.text.x = element_text(angle = 0, hjust = 0.4, vjust = 1.1, size = 8), # Angle x-axis text
      axis.text.y = element_text(size = 8),
      axis.title.x = element_blank(),
      axis.title.y = element_blank(),
    )
  plot.biome <- ggplot(data = Data, aes(x = Nstat_obs, y = indicator, fill = Nstat_obs))+
    geom_boxplot(alpha = 0.3, show.legend = F)+
    stat_summary(aes(color = Nstat_obs), fun = "mean", geom = "point", position = "identity", 
                 shape = 20, size = 5, show.legend = F)+
    facet_wrap(~ Biome)+
    scale_fill_manual(values = c("#4286f4", "#ec2F4B"))+
    scale_color_manual(values = c("#4286f4", "#ec2F4B"))+
    scale_x_discrete(breaks = c("limited", "saturated"), 
                     labels = c("limited", "saturated"))+
    #ylim(c(0, 0.035))+
    theme(
      axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.8),
      axis.title.x = element_blank(),
      axis.title.y = element_blank()
    )
  plot.out <- plot_grid(plot.data, plot.biome, nrow =2, ncol = 1, rel_heights = c(6,4))
  cmb <- boot.mean.diff(Data, "Nstat_obs", "indicator", R)
  tr <- boot.mean.diff(Data[Data$Biome == "Tropical",], "Nstat_obs", "indicator", R)
  te <- boot.mean.diff(Data[Data$Biome == "Temperate",], "Nstat_obs", "indicator", R)
  bo <- boot.mean.diff(Data[Data$Biome == "Boreal",], "Nstat_obs", "indicator", R)
  pvalue.out <- paste0("pvalues: All ", cmb[[1]]$pvalue, "Tropical ", tr[[1]]$pvalue, 
                       ", Temperate ", te[[1]]$pvalue, ", Boreal ", bo[[1]]$pvalue)
  list(plot.out, pvalue.out)
}
clas.biome(Nstat.all, "sN", 1999)
clas.SOC <- clas.biome(Nstat.all, "SOC", 1999)
clas.Clay <- clas.biome(Nstat.all, "Clay", 1999)
clas.CEC <- clas.biome(Nstat.all, "CEC", 1999)
clas.pH <- clas.biome(Nstat.all, "pH", 1999)
clas.rCN <- clas.biome(Nstat.all, "rCN", 1999)
clas.Ndepo <- clas.biome(Nstat.all, "Ndepo", 1999)
#ggsave("clas_Clay.svg", width = 12, height = 12, units = "cm")

##Testing the sN of different forest types (Deciduous broadleaf vs. Needleleaf)
dat.typ1 <- data7
for(i in 1:nrow(dat.typ1)){
  if(dat.typ1$EDNeedleleaf[i] >= 50){
    dat.typ1$ForestType[i] <- "Needleleaf"
  }else if(dat.typ1$EBroadleaf[i] >= 50){
    dat.typ1$ForestType[i] <- "EBroadleaf"
  }else if(dat.typ1$DBroadleaf[i] >= 50){
    dat.typ1$ForestType[i] <- "DBroadleaf"
  }else{
    dat.typ1$ForestType[i] <- "Mixed"
  }
}
dat.typ2 <- dat.typ1 %>%
  filter(ForestType == "Needleleaf" | ForestType == "DBroadleaf" )
dat.typ2$ForestType <- ordered(dat.typ2$ForestType, levels = c("Needleleaf", "DBroadleaf"))
dat.typ2$Biome <- ordered(dat.typ2$Biome, c("Tropical", "Temperate", "Boreal"))
plot1 <- ggplot(data = dat.typ2, aes(x = ForestType, y = sN)) +
  geom_violin(aes(color = ForestType), width=0.8, show.legend = F) +
  geom_boxplot(aes(fill = ForestType, color = ForestType), width=0.4,  alpha = 0.5, outlier.alpha=0, show.legend = F)+
  stat_summary(fun= mean, geom="point", shape=20, size=6, color="dark red", fill="dark red")+
  stat_summary(fun = mean, geom = "label", col = "black", size = 2.5,
               hjust = -0.6, vjust = 0.5, 
               aes(label = paste("mean =", round(..y.., digits = 3))))+
  ylim(c(0, 0.035))+
  scale_x_discrete(breaks = c("Needleleaf", "DBroadleaf"),
                   labels = c(paste0("Needleleaf\n(n =", nrow(dat.typ2[dat.typ2$ForestType == "Needleleaf",]), ")"), 
                              paste0("DBroadleaf\n(n =", nrow(dat.typ2[dat.typ2$ForestType == "DBroadleaf",]),")")))+
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.4, vjust = 1.1, size = 8), # Angle x-axis text
    axis.text.y = element_text(size = 8),
    axis.title.x = element_blank(),
    #    axis.title.y = element_blank(),
  )
plot2 <- ggplot(data = dat.typ2, aes(x = ForestType, y = sN, fill = ForestType))+
  geom_boxplot(alpha = 0.3, outlier.alpha=0, show.legend = F)+
  stat_summary(aes(color = ForestType), fun = "mean", geom = "point", position = "identity", 
               shape = 20, size = 5, show.legend = F)+
  facet_wrap(~ Biome)+
  ylim(c(0, 0.03))+
  scale_x_discrete(breaks = c("Needleleaf", "DBroadleaf"))+ 
  theme(
    axis.text.x = element_text(angle = 0, hjust = 0.5, vjust = 0.8),
    axis.title.x = element_blank(),
    #    axis.title.y = element_blank()
  )
plot_grid(plot1, plot2, nrow =2, ncol = 1, rel_heights = c(6,4))
boot.mean.diff(dat.typ2, "ForestType", "sN", 1999)

##Comparing temporal trends of N availability indicated by sN and N isotope ratio
ss <- unique.data.frame(cbind(N2O.LN$xlon, N2O.LN$ylat, N2O.LN$Year))
sN.data <- {}
for(i in 1:length(ss[,1])){
  temp1 <- subset(N2O.LN, N2O.LN$xlon==ss[i,1] & N2O.LN$ylat==ss[i,2])
  if(any(temp1$N_Add==0)&any(temp1$N_Add>0)){
    mod1 <- lm(temp1$RN2O ~ temp1$N_Input)
    a <- summary(mod1)
    p <- pf(a$fstatistic[1], a$fstatistic[2], a$fstatistic[3], lower.tail = FALSE)
    if(is.na(a$adj.r.squared)| a$adj.r.squared>0){
      sN.data <- rbind(sN.data, c(temp1$xlon[1], temp1$ylat[1], temp1$Biome[1], temp1$Year[1], temp1$Region[1],
                                  a$coefficients[2,1], a$coefficients[2,4], a$coefficients[1,1], a$coefficients[1,4], a$df[2]+2, a$adj.r.squared, p))
    }
  }
}
sN.data <- as.data.frame(sN.data)
colnames(sN.data) <- c("xlon", "ylat", "Biome", "Year", "Region", "sN", "p_sN", "R0", "p_R0", "n", "R2adjusted", "pvalue")
sN.data$Biome <- ordered(sN.data$Biome, levels = c(1,2,3), labels = c("Tropical", "Temperate", "Boreal"))
sN.data$Region <- factor(sN.data$Region, levels = c(1:6), labels = c("East Asia", "Latin America and Caribbean", "North America",
                                                                     "Southeast Asia and Pacific", "Sub-Saharan Africa", "Western Europe"))
sN.data$Region2 <- sN.data$Region
levels(sN.data$Region2)[levels(sN.data$Region2)%in%c("East Asia", "Southeast Asia and Pacific")] <- "Asia Pacific"
levels(sN.data$Region2)[levels(sN.data$Region2)%in%c("Western Europe")] <- "Europe"
levels(sN.data$Region2)[levels(sN.data$Region2)%in%c("Latin America and Caribbean")] <- "South America"
levels(sN.data$Region2)[levels(sN.data$Region2)%in%c("Sub-Saharan Africa")] <- "Africa"
sN.data$Region2 <- ordered(sN.data$Region2, levels = c("Asia Pacific", "Europe", "North America", "South America", "Africa"))
plot.sN.Year <- sN.data %>%
  ggplot(aes(x = Year, y = sN))+
  geom_point(color = "grey80", size = 0.5)+
  stat_summary(geom = "point", fun = mean, size = 0.7)+
  stat_smooth(method = "lm") +
  facet_wrap(~ Region2, ncol = 5)+
  xlim(c(1985, 2020))+
  theme_update(
    axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0, size = 8))
summary(lm(sN ~ Year, data = subset(sN.data, Region2 == "Asia Pacific")))
summary(lm(sN ~ Year, data = subset(sN.data, Region2 == "Europe")))
summary(lm(sN ~ Year, data = subset(sN.data, Region2 == "North America")))
summary(lm(sN ~ Year, data = subset(sN.data, Region2 == "South America")))

NIso <- read_csv("InputData.csv") #global foliar d15N data were from Craine et al. (2019) (https://doi.org/10.5061/dryad.v2k2607)
NIso$Continent <- ordered(NIso$Continent, levels = c("Asia", "Australia", "Europe", "North America", "South America", "Africa"))
NIso$Region2 <- NIso$Continent
levels(NIso$Region2)[levels(NIso$Region2)%in%c("Asia", "Australia")] <- "Asia Pacific"
NIso2 <- NIso %>%
  filter(!is.na(Region2), Year >1985) %>%
  group_by(Latitude, Longitude, Year, Region2) %>%
  summarize(across(where(is.numeric), mean, na.rm = T))
plot.d15N.Year <-  NIso2 %>%
  ggplot(aes(x = Year, y = Leaf15N)) +
  geom_point(color = "grey80", size = 0.5)+
  stat_summary(geom = "point", fun = mean, size = 0.7)+
  stat_smooth(method = "lm") +
  facet_wrap(~ Region2, ncol = 5)+
  xlim(c(1985, 2020))+
  theme_update(
    axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0, size = 8))
summary(lm(Leaf15N ~ Year, data = subset(NIso2, Region2 == "Asia Pacific" & Year >= 1985)))
summary(lm(Leaf15N ~ Year, data = subset(NIso2, Region2 == "Europe" & Year >= 1985)))
summary(lm(Leaf15N ~ Year, data = subset(NIso2, Region2 == "North America" & Year >= 1985)))
summary(lm(Leaf15N ~ Year, data = subset(NIso2, Region2 == "South America" & Year >= 1985)))
summary(lm(Leaf15N ~ Year, data = subset(NIso2, Region2 == "Africa" & Year >= 1985)))

plot_grid(plot.sN.Year, plot.d15N.Year, ncol = 1, nrow = 2, align = "hv")
#ggsave("plot.Niso.svg", width = 17, height = 9, units = "cm")