# loading Data,sheet=1 is rawdata, sheet=1 is NIE data
rawdata<-readxl::read_excel("./Raw data.xlsx",sheet=1,na=c("NA","#DIV/0!"))


str(rawdata)
rawdata$Year<-as.factor(rawdata$Year)
rawdata$Block<-as.factor(rawdata$Block)
rawdata$Plot<-as.factor(rawdata$Plot)
rawdata$Treatment<-as.factor(rawdata$Treatment)
rawdata$N<-as.factor(rawdata$N)
rawdata$P<-as.factor(rawdata$P)

str(rawdata)


# repeated two-way ANOVA 
library(ez)
library(car)

# repeated two-way ANOVA for ANPP
rm_anova <- ezANOVA(
  data = rawdata,
  dv = .(ANPP),
  wid = .(Plot),
  within = .(Year),
  between = .(N, P),
  type = 3,  
  detailed = TRUE,
  return_aov = TRUE
)

# print results
print(rm_anova)
# Extract the data of the ANOVA table
anova_data <- rm_anova$ANOVA

# Define the significance marking function
get_significance_star <- function(p_value) {
  if (p_value < 0.001) {
    return("***")
  } else if (p_value < 0.01) {
    return("**")
  } else if (p_value < 0.05) {
    return("*")
  } else if (p_value < 0.1) {
    return("^")
  } else {
    return("ns")
  }
}

# Apply thesignificance marking function to the ANOVA table
anova_data$significance <- sapply(anova_data$p, get_significance_star)
anova_data

# Ensure that the F value is numerical and formatted as a common value (not a scientific count).
library(knitr)
anova_data$F <- format(anova_data$F, scientific = FALSE, digits = 4)  
anova_data

# repeated two-way ANOVA for ANPP of carbon uptake rate (CUR)
#the CUE data in 2022-2024 is NA, so susbset it
str(rawdata)

# Select specified columns from rawdata and remove rows with NA values
library(dplyr)
Carbondata <- rawdata %>% 
  select(Year, Block, Plot, Treatment, N, P, CUR, LUE, WUE) %>% 
  na.omit()  # Remove rows containing any NA values

# Verify the structure of the processed data
str(Carbondata)


rm_anova <- ezANOVA(
  data = Carbondata,
  dv = .(CUR),
  wid = .(Plot),
  within = .(Year),
  between = .(N, P),
  type = 3,  
  detailed = TRUE,
  return_aov = TRUE
)

# print results
print(rm_anova)
# Extract the data of the ANOVA table
anova_data <- rm_anova$ANOVA

# Define the significance marking function
get_significance_star <- function(p_value) {
  if (p_value < 0.001) {
    return("***")
  } else if (p_value < 0.01) {
    return("**")
  } else if (p_value < 0.05) {
    return("*")
  } else if (p_value < 0.1) {
    return("^")
  } else {
    return("ns")
  }
}

# Apply thesignificance marking function to the ANOVA table
anova_data$significance <- sapply(anova_data$p, get_significance_star)
anova_data

# Ensure that the F value is numerical and formatted as a common value (not a scientific count).
library(knitr)
anova_data$F <- format(anova_data$F, scientific = FALSE, digits = 4)  
anova_data


# repeated two-way ANOVA for ANPP of light use efficiency (LUE)
str(Carbondata)
rm_anova <- ezANOVA(
  data = Carbondata,
  dv = .(LUE),
  wid = .(Plot),
  within = .(Year),
  between = .(N, P),
  type = 3,  
  detailed = TRUE,
  return_aov = TRUE
)

# print results
print(rm_anova)
# Extract the data of the ANOVA table
anova_data <- rm_anova$ANOVA

# Define the significance marking function
get_significance_star <- function(p_value) {
  if (p_value < 0.001) {
    return("***")
  } else if (p_value < 0.01) {
    return("**")
  } else if (p_value < 0.05) {
    return("*")
  } else if (p_value < 0.1) {
    return("^")
  } else {
    return("ns")
  }
}

# Apply thesignificance marking function to the ANOVA table
anova_data$significance <- sapply(anova_data$p, get_significance_star)
anova_data

# Ensure that the F value is numerical and formatted as a common value (not a scientific count).
library(knitr)
anova_data$F <- format(anova_data$F, scientific = FALSE, digits = 4)  
anova_data



# repeated two-way ANOVA for ANPP of water use efficiency (WUE)
str(Carbondata)
rm_anova <- ezANOVA(
  data = Carbondata,
  dv = .(WUE),
  wid = .(Plot),
  within = .(Year),
  between = .(N, P),
  type = 3,  
  detailed = TRUE,
  return_aov = TRUE
)

# print results
print(rm_anova)
# Extract the data of the ANOVA table
anova_data <- rm_anova$ANOVA

# Define the significance marking function
get_significance_star <- function(p_value) {
  if (p_value < 0.001) {
    return("***")
  } else if (p_value < 0.01) {
    return("**")
  } else if (p_value < 0.05) {
    return("*")
  } else if (p_value < 0.1) {
    return("^")
  } else {
    return("ns")
  }
}

# Apply thesignificance marking function to the ANOVA table
anova_data$significance <- sapply(anova_data$p, get_significance_star)
anova_data

# Ensure that the F value is numerical and formatted as a common value (not a scientific count).
library(knitr)
anova_data$F <- format(anova_data$F, scientific = FALSE, digits = 4)  
anova_data





# repeated two-way ANOVA for ANPP of rain use efficiency (RUE)
str(rawdata)
rm_anova <- ezANOVA(
  data = rawdata,
  dv = .(RUE),
  wid = .(Plot),
  within = .(Year),
  between = .(N, P),
  type = 3,  
  detailed = TRUE,
  return_aov = TRUE
)

# print results
print(rm_anova)
# Extract the data of the ANOVA table
anova_data <- rm_anova$ANOVA

# Define the significance marking function
get_significance_star <- function(p_value) {
  if (p_value < 0.001) {
    return("***")
  } else if (p_value < 0.01) {
    return("**")
  } else if (p_value < 0.05) {
    return("*")
  } else if (p_value < 0.1) {
    return("^")
  } else {
    return("ns")
  }
}

# Apply thesignificance marking function to the ANOVA table
anova_data$significance <- sapply(anova_data$p, get_significance_star)
anova_data

# Ensure that the F value is numerical and formatted as a common value (not a scientific count).
library(knitr)
anova_data$F <- format(anova_data$F, scientific = FALSE, digits = 4)  
anova_data




# repeated two-way ANOVA for ANPP of N use efficiency (NUE)

# Select specified columns from rawdata and remove rows with NA values
library(dplyr)
Nut_data <- rawdata %>% 
  select(Year, Block, Plot, Treatment, N, P, NUE_ANPP, PUE_ANPP) %>% 
  na.omit()  # Remove rows containing any NA values

# Verify the structure of the processed data
str(Nut_data)

rm_anova <- ezANOVA(
  data = Nut_data,
  dv = .(NUE_ANPP),
  wid = .(Plot),
  within = .(Year),
  between = .(N, P),
  type = 3,  
  detailed = TRUE,
  return_aov = TRUE
)

# run error: Error in ezANOVA_main(data = data, dv = dv, wid = wid, within = within,  : 
#One or more cells is missing data. Try using ezDesign() to check your data.

# because 4 data in Block 5 in 2020 is NA, so choose the Linear mixed - effects model


#Linear mixed - effects model for NUE
library(lme4)
library(lmerTest)

model1 <- lmer(NUE_ANPP ~ Year*N*P+(1|Plot), data = rawdata)
anova(model1)


model2 <- lmer(PUE_ANPP ~ Year*N*P+(1|Plot), data = rawdata)
anova(model2)

############################# code for Figure 1
#set the plot theme
library(ggplot2)
theme_custom<- function() {
  theme(
    axis.title.x=element_text(vjust=1,size=40) , 
    axis.title.y = element_text(size = 40, color = "black", margin = margin(r = 0.25, unit = "cm")),  
    axis.text.x= element_text(size=35,margin = unit(c(0.5, 0.5, 0.5, 0.5), 'cm'),angle=0),
    axis.text.y= element_text(size=35,margin = unit(c(0.5, 0.5, 0.5, 0.5), 'cm'),angle=0),
    text = element_text(family = "Times New Roman"),
    axis.line = element_line(color = "black", size = 1 ), 
    axis.title = element_text( color = "black"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black", size =1),
    axis.ticks.length.y = unit(-0.5, 'cm'), 
    axis.ticks.length.x = unit(0, 'cm'), 
    plot.title = element_text(color = "black"),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    line = element_line(size = 2) 
  )
}


str(rawdata)
rawdata$Year<-as.factor(rawdata$Year)
rawdata$Block<-as.factor(rawdata$Block)
rawdata$Plot<-as.factor(rawdata$Plot)
rawdata$Treatment<-as.factor(rawdata$Treatment)
rawdata$N<-as.factor(rawdata$N)
rawdata$P<-as.factor(rawdata$P)
str(rawdata)

#subset data
data20<-rawdata[which((rawdata$Year== '2020')), ]
data21<-rawdata[which((rawdata$Year== '2021')), ]
data22<-rawdata[which((rawdata$Year== '2022')), ]
data23<-rawdata[which((rawdata$Year== '2023')), ]
data24<-rawdata[which((rawdata$Year== '2024')), ]



#Multiple comparisons for ANPP in 2020
str(data20)
#1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
model <- lmer(ANPP~ Treatment+(1|Block), data = data20)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#2-Multiple comparisons
#2.1-means of  Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
#2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha =0.05, reverse = TRUE)
print(cld_result5)



# Merge the mean sum and letter tags
library(Rmisc) 
Mean<- summarySE(data20, measurevar="ANPP", groupvars=c("Year","Treatment"))
Mean

cld_result5
#Merge
Mean_20<- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_20



#Multiple comparisons for ANPP in 2021
str(data21)
#1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
model <- lmer(ANPP~ Treatment+(1|Block), data = data21)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#2-Multiple comparisons
#2.1-means of  Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
#2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha =0.05, reverse = TRUE)
print(cld_result5)



# Merge the mean sum and letter tags
library(Rmisc) 
Mean<- summarySE(data21, measurevar="ANPP", groupvars=c("Year","Treatment"))
Mean

cld_result5
#Merge
Mean_21<- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_21





#Multiple comparisons for ANPP in 2023
str(data22)
#1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
model <- lmer(ANPP~ Treatment+(1|Block), data = data22)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#2-Multiple comparisons
#2.1-means of  Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
#2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha =0.05, reverse = TRUE)
print(cld_result5)



# Merge the mean sum and letter tags
library(Rmisc) 
Mean<- summarySE(data22, measurevar="ANPP", groupvars=c("Year","Treatment"))
Mean

cld_result5
#Merge
Mean_22<- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_22




#Multiple comparisons for ANPP in 2023
str(data23)
#1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
model <- lmer(ANPP~ Treatment+(1|Block), data = data23)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#2-Multiple comparisons
#2.1-means of  Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
#2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha =0.05, reverse = TRUE)
print(cld_result5)



# Merge the mean sum and letter tags
library(Rmisc) 
Mean<- summarySE(data23, measurevar="ANPP", groupvars=c("Year","Treatment"))
Mean

cld_result5
#Merge
Mean_23<- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_23



#Multiple comparisons for ANPP in 2024
str(data24)
#1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
model <- lmer(ANPP~ Treatment+(1|Block), data = data24)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#2-Multiple comparisons
#2.1-means of  Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
#2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha =0.05, reverse = TRUE)
print(cld_result5)



# Merge the mean sum and letter tags
library(Rmisc) 
Mean<- summarySE(data24, measurevar="ANPP", groupvars=c("Year","Treatment"))
Mean

cld_result5
#Merge
Mean_24<- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_24


#Multiple comparisons for ANPP in 2020-2024
str(rawdata)
#1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
model <- lmer(ANPP~ Treatment+(1|Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#Test for normality of model residuals(p-value = 0.000531), so log the data
model_log <- lmer(log(ANPP)~ Treatment+(1|Block), data = rawdata)
anova(model_log )
# Extract the model residuals
residuals <- residuals(model_log , type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#2-Multiple comparisons
#2.1-means of  Treatment
library(emmeans)
emm <- emmeans(model_log , ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
#2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha =0.05, reverse = TRUE)
print(cld_result5)



# Merge the mean sum and letter tags
library(Rmisc) 
Mean<- summarySE(rawdata, measurevar="ANPP", groupvars=c("Treatment"))
Mean
#Add a new column called "Overall"
Mean$Year<-rep("Overall",length(Mean$Treatment))
Mean

cld_result5
#Merge
Mean_Total<- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_Total


#Merge all the data
Mean_Total
Mean_20
Mean_21
Mean_22
Mean_23
Mean_24


A_ANPPSE1<-rbind(Mean_Total,Mean_20, Mean_21,Mean_22,Mean_23,Mean_24)
# Remove the Spaces before and after the letters in the.group column
A_ANPPSE1$.group <- trimws(A_ANPPSE1$.group)
A_ANPPSE1


#Determine the interaction type between N and P by NIE

## loading Data,sheet=1 is rawdata, sheet=1 is NIE data
NIEdata<-readxl::read_excel("./Raw data.xlsx",sheet=2,na=c("NA","#DIV/0!"))


str(NIEdata)
NIEdata$Year<-as.factor(NIEdata$Year)
NIEdata$Block<-as.factor(NIEdata$Block)



#subset data
data20<-NIEdata[which((NIEdata$Year== '2020')), ]
data21<-NIEdata[which((NIEdata$Year== '2021')), ]
data22<-NIEdata[which((NIEdata$Year== '2022')), ]
data23<-NIEdata[which((NIEdata$Year== '2023')), ]
data24<-NIEdata[which((NIEdata$Year== '2024')), ]

# Create a function to convert the t-test results into a data frame when the data is normal
tidy_t_test <- function(t_test_result) {
  data.frame(
    test_type = "One Sample t-test",
    t_value = t_test_result$statistic,
    df = t_test_result$parameter,
    p_value = t_test_result$p.value,
    alternative = t_test_result$alternative,
    conf_low = t_test_result$conf.int[1],
    conf_high = t_test_result$conf.int[2],
    conf_level = attr(t_test_result$conf.int, "conf.level"),
    mean = t_test_result$estimate,
    stringsAsFactors = FALSE
  )
}


# Create a function to convert the wilcox_test results into a data frame when the data is  not normal
tidy_wilcox_test <- function(wilcox_result) {
  data.frame(
    test_type = ifelse(is.null(wilcox_result$method), 
                       "Wilcoxon Test", 
                       wilcox_result$method),
    statistic = wilcox_result$statistic,
    p_value = wilcox_result$p.value,
    alternative = wilcox_result$alternative,
    stringsAsFactors = FALSE
  )
}


# NIE of ANPP in 2020
#Normality test
shapiro.test(data20$ANPP)
#T-test
t_result20<-t.test(data20$ANPP,mu=0)
#convert the t-test results into a data frame using the function
result20 <- tidy_t_test(t_result20)
result20
#Add a Year column 
result20$Year<-rep("2020",length(result20$mean))
result20


# NIE of ANPP in 2021
#Normality test
shapiro.test(data21$ANPP)
#T-test
t_result21<-t.test(data21$ANPP,mu=0)
#convert the t-test results into a data frame using the function
result21 <- tidy_t_test(t_result21)
result21
#Add a Year column 
result21$Year<-rep("2021",length(result21$mean))
result21





# NIE of ANPP in 2022
#Normality test
shapiro.test(data22$ANPP)
#T-test
t_result22<-t.test(data22$ANPP,mu=0)
#convert the t-test results into a data frame using the function
result22 <- tidy_t_test(t_result22)
result22
#Add a Year column 
result22$Year<-rep("2022",length(result22$mean))
result22

# NIE of ANPP in 2023
#Normality test
shapiro.test(data23$ANPP)
#T-test
t_result23<-t.test(data23$ANPP,mu=0)
#convert the t-test results into a data frame using the function
result23 <- tidy_t_test(t_result23)
result23
#Add a Year column 
result23$Year<-rep("2023",length(result23$mean))
result23


# NIE of ANPP in 2024
#Normality test
shapiro.test(data24$ANPP)
#T-test
t_result24<-t.test(data24$ANPP,mu=0)
#convert the t-test results into a data frame using the function
result24 <- tidy_t_test(t_result24)
result24
#Add a Year column 
result24$Year<-rep("2024",length(result24$mean))
result24


# NIE of ANPP in 2020-2024
#Normality test
shapiro.test(NIEdata$ANPP)
#wilcox_test 
wilcox_result<-wilcox.test(NIEdata$ANPP, mu = 0)
#convert the wilcox_test  results into a data frame using the function
wilcox_result <- tidy_wilcox_test(wilcox_result)
wilcox_result
#Add a Year column 
wilcox_result$Year<-rep("Overall",length(wilcox_result$p_value))
wilcox_result
#calculate mean and SD in 2020-2024
library(Rmisc)  
mean_whole_year<- summarySE(NIEdata, measurevar="ANPP")
mean_whole_year

wilcox_result$mean<-mean_whole_year$ANPP
wilcox_result
#merge the data
wilcox_result
result24
result23
result22
result21
result20


# Directly combine all result data frames including p-values
A_ANPP_NIE <- rbind(
  # Extract four columns from the Wilcoxon test results
  data.frame(
    test_type = wilcox_result$test_type,
    mean = wilcox_result$mean,
    Year = wilcox_result$Year,
    p_value = wilcox_result$p_value
  ),
  # Extract four columns from the 2024 t-test results
  data.frame(
    test_type = result24$test_type,
    mean = result24$mean,
    Year = result24$Year,
    p_value = result24$p_value
  ),
  # Extract four columns from the 2023 t-test results
  data.frame(
    test_type = result23$test_type,
    mean = result23$mean,
    Year = result23$Year,
    p_value = result23$p_value
  ),
  # Extract four columns from the 2022 t-test results
  data.frame(
    test_type = result22$test_type,
    mean = result22$mean,
    Year = result22$Year,
    p_value = result22$p_value
  ),
  # Extract four columns from the 2021 t-test results
  data.frame(
    test_type = result21$test_type,
    mean = result21$mean,
    Year = result21$Year,
    p_value = result21$p_value
  ),
  # Extract four columns from the 2020 t-test results
  data.frame(
    test_type = result20$test_type,
    mean = result20$mean,
    Year = result20$Year,
    p_value = result20$p_value
  )
)

# Reset row indices
rownames(A_ANPP_NIE) <- NULL

# View the combined results
print(A_ANPP_NIE)



# 1.定义提取显著性标记
get_significance_star <- function(p_value) {
  if (p_value < 0.001) {
    return("***")
  } else if (p_value < 0.01) {
    return("**")
  } else if (p_value < 0.05) {
    return("*")
  } else if (p_value < 0.1) {
    return("^")
  } else {
    return("ns")
  }
}


# 3.获取显著性星号
A_ANPP_NIE$p_value_mark<- sapply(A_ANPP_NIE$p_value, get_significance_star)
A_ANPP_NIE




str(A_ANPP_NIE)#resluts of NIE (interactive type)
str(A_ANPPSE1)#means of ANPP 

library(dplyr) 

ALl_merged_data <- left_join(A_ANPPSE1, A_ANPP_NIE %>% 
                               dplyr::select(Year, p_value_mark), by = c("Year" = "Year"))


# make the Fig 1
str(ALl_merged_data )
ALl_merged_data$Treatment
ALl_merged_data$Treatment= factor(ALl_merged_data$Treatment, levels=c('Ctrl','N','P','NP'))

#paste
ALl_merged_data$Two_way<-paste(ALl_merged_data$Treatment,ALl_merged_data$Significance,sep = "")
unique(ALl_merged_data$Two_way)
ALl_merged_data$Two_way <- gsub("NP", "N × P",ALl_merged_data$Two_way)


str(ALl_merged_data)

ALl_merged_data <- ALl_merged_data %>%
  mutate(
    Total_mark = case_when(
      Treatment == "Ctrl" ~ paste("Synergistic (NIE > 0", p_value_mark, ")", sep = ""),
      TRUE ~ Two_way  # 非Ctrl时使用Two_way列的值
    )
  )




# Sort the data frames according to the order of the treatment factor levels
library(dplyr)
ALl_merged_data<- ALl_merged_data %>% arrange(Treatment)




mycolor<-c("#A7C2B1","#6C91CB","#B0E0E6","#E69091")


library(ggplot2)
str(ALl_merged_data)#data for Fig 1

unique(ALl_merged_data$Year)

#make plot of ANPP in 2020

# Create a text label column that only contains the Total_mark content of Ctrl
Fig_data_20 <- ALl_merged_data %>% 
  filter(Year == "2020") %>%
  mutate(
    # Use Ctrl's Total_mark directly as the tag
    label_text = ifelse(Treatment == "Ctrl", Total_mark, NA)
  )

# check data
Fig_data_20

# make the plot
plot_ANPP20 <- ggplot(Fig_data_20, aes(x = Treatment, y = ANPP, fill = Treatment)) +
  geom_bar(stat = "identity", position = position_dodge(0.8), width = 0.8, alpha = 0.9) +
  geom_errorbar(aes(ymax = ANPP + se, ymin = ANPP - se), position = position_dodge(0.8), width = 0.15, size = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(aes(y = ANPP + 3*se, label = .group), position = position_dodge(0.8), hjust = 0.5, vjust = 0, size = 8) +
  geom_text(
    data = Fig_data_20 %>% filter(Treatment == "Ctrl"),
    aes(x = 0.8, y = 600, label = label_text),
    hjust = 0, vjust = 1.5, size = 8, family = "Times New Roman", lineheight = 1.2, inherit.aes = FALSE
  ) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 600, 150), limits = c(0, 600)) +
  xlab("Treatment") +
  ylab("ANPP in 2020 (g m⁻²)") +
  theme(legend.position = "none") +
  theme_custom()


plot_ANPP20


## Make plot of ANPP in 2021
# 1. Create a text label column, only including Ctrl's Total_mark content
Fig_data_21 <- ALl_merged_data %>% filter(Year == "2021") %>%
  mutate(
    label_text = ifelse(Treatment == "Ctrl", Total_mark, NA)
  )

plot_ANPP21 <- ggplot(Fig_data_21, aes(x = Treatment, y = ANPP, fill = Treatment)) +
  geom_bar(stat = "identity", position = position_dodge(0.8), width = 0.8, alpha = 0.9) +
  geom_errorbar(aes(ymax = ANPP + se, ymin = ANPP - se), position = position_dodge(0.8), width = 0.15, size = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(aes(y = ANPP + 3*se, label = .group), position = position_dodge(0.8), hjust = 0.5, vjust = 0, size = 8) +
  geom_text(
    data = Fig_data_21 %>% filter(Treatment == "Ctrl"),
    aes(x = 0.8, y = 600, label = label_text),
    hjust = 0, vjust = 1.5, size = 8, family = "Times New Roman", lineheight = 1.2, inherit.aes = FALSE
  ) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 600, 150), limits = c(0, 600)) +
  xlab("Treatment") +
  ylab("ANPP in 2021 (g m⁻²)") +
  theme(legend.position = "none") +
  theme_custom()

## Make plot of ANPP in 2022
# 1. Create a text label column, only including Ctrl's Total_mark content
Fig_data_22 <- ALl_merged_data %>% filter(Year == "2022") %>%
  mutate(
    label_text = ifelse(Treatment == "Ctrl", Total_mark, NA)
  )

plot_ANPP22 <- ggplot(Fig_data_22, aes(x = Treatment, y = ANPP, fill = Treatment)) +
  geom_bar(stat = "identity", position = position_dodge(0.8), width = 0.8, alpha = 0.9) +
  geom_errorbar(aes(ymax = ANPP + se, ymin = ANPP - se), position = position_dodge(0.8), width = 0.15, size = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(aes(y = ANPP + 3*se, label = .group), position = position_dodge(0.8), hjust = 0.5, vjust = 0, size = 8) +
  geom_text(
    data = Fig_data_22 %>% filter(Treatment == "Ctrl"),
    aes(x = 0.8, y = 600, label = label_text),
    hjust = 0, vjust = 1.5, size = 8, family = "Times New Roman", lineheight = 1.2, inherit.aes = FALSE
  ) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 600, 150), limits = c(0, 600)) +
  xlab("Treatment") +
  ylab("ANPP in 2022 (g m⁻²)") +
  theme(legend.position = "none") +
  theme_custom()

## Make plot of ANPP in 2023
# 1. Create a text label column, only including Ctrl's Total_mark content
Fig_data_23 <- ALl_merged_data %>% filter(Year == "2023") %>%
  mutate(
    label_text = ifelse(Treatment == "Ctrl", Total_mark, NA)
  )

plot_ANPP23 <- ggplot(Fig_data_23, aes(x = Treatment, y = ANPP, fill = Treatment)) +
  geom_bar(stat = "identity", position = position_dodge(0.8), width = 0.8, alpha = 0.9) +
  geom_errorbar(aes(ymax = ANPP + se, ymin = ANPP - se), position = position_dodge(0.8), width = 0.15, size = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(aes(y = ANPP + 3*se, label = .group), position = position_dodge(0.8), hjust = 0.5, vjust = 0, size = 8) +
  geom_text(
    data = Fig_data_23 %>% filter(Treatment == "Ctrl"),
    aes(x = 0.8, y = 600, label = label_text),
    hjust = 0, vjust = 1.5, size = 8, family = "Times New Roman", lineheight = 1.2, inherit.aes = FALSE
  ) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 600, 150), limits = c(0, 600)) +
  xlab("Treatment") +
  ylab("ANPP in 2023 (g m⁻²)") +
  theme(legend.position = "none") +
  theme_custom()

## Make plot of ANPP in 2024
# 1. Create a text label column, only including Ctrl's Total_mark content
Fig_data_24 <- ALl_merged_data %>% filter(Year == "2024") %>%
  mutate(
    label_text = ifelse(Treatment == "Ctrl", Total_mark, NA)
  )

plot_ANPP24 <- ggplot(Fig_data_24, aes(x = Treatment, y = ANPP, fill = Treatment)) +
  geom_bar(stat = "identity", position = position_dodge(0.8), width = 0.8, alpha = 0.9) +
  geom_errorbar(aes(ymax = ANPP + se, ymin = ANPP - se), position = position_dodge(0.8), width = 0.15, size = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(aes(y = ANPP + 3*se, label = .group), position = position_dodge(0.8), hjust = 0.5, vjust = 0, size = 8) +
  geom_text(
    data = Fig_data_24 %>% filter(Treatment == "Ctrl"),
    aes(x = 0.8, y = 600, label = label_text),
    hjust = 0, vjust = 1.5, size = 8, family = "Times New Roman", lineheight = 1.2, inherit.aes = FALSE
  ) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 600, 150), limits = c(0, 600)) +
  xlab("Treatment") +
  ylab("ANPP in 2024 (g m⁻²)") +
  theme(legend.position = "none") +
  theme_custom()

## Make plot of ANPP in 2020-2024 (Overall)
# 1. Create a text label column, only including Ctrl's Total_mark content
Fig_data_Overall <- ALl_merged_data %>% filter(Year == "Overall") %>%
  mutate(
    label_text = ifelse(Treatment == "Ctrl", Total_mark, NA)
  )

plot_ANPP_Overall <- ggplot(Fig_data_Overall, aes(x = Treatment, y = ANPP, fill = Treatment)) +
  geom_bar(stat = "identity", position = position_dodge(0.8), width = 0.8, alpha = 0.9) +
  geom_errorbar(aes(ymax = ANPP + se, ymin = ANPP - se), position = position_dodge(0.8), width = 0.15, size = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(aes(y = ANPP + 3*se, label = .group), position = position_dodge(0.8), hjust = 0.5, vjust = 0, size = 8) +
  geom_text(
    data = Fig_data_Overall %>% filter(Treatment == "Ctrl"),
    aes(x = 0.8, y = 600, label = label_text),
    hjust = 0, vjust = 1.5, size = 8, family = "Times New Roman", lineheight = 1.2, inherit.aes = FALSE
  ) +
  scale_y_continuous(expand = c(0, 0), breaks = seq(0, 600, 150), limits = c(0, 600)) +
  xlab("Treatment") +
  ylab("ANPP (Overall) (g m⁻²)") +
  theme(legend.position = "none") +
  theme_custom()






#show the plot
plot_ANPP20
plot_ANPP21
plot_ANPP22
plot_ANPP23
plot_ANPP24
plot_ANPP_Overall

library(cowplot)
library(eoffice)
library(ggplot2)

# Set global font and unified theme
custom_theme <- theme(
  # Set font to Times New Roman, non-bold
  text = element_text(family = "Times New Roman", face = "plain"),
  # Uniform margin between axis titles and the plot
  axis.title.y = element_text(margin = margin(r = 45)),  # 45pt right margin for Y-axis title
  axis.title.x = element_text(margin = margin(t = 55)),  # 55pt top margin for X-axis title
  # Uniform plot margins (critical! Ensure consistent height and width)
  plot.margin = margin(15, 15, 15, 15, "pt"),
  # Remove default theme distractions
  panel.grid = element_blank(),
  plot.background = element_rect(fill = "white"),
  panel.background = element_rect(fill = "white")
)

# Apply the unified theme to all plots
apply_theme <- function(plot) {
  plot + custom_theme
}

# Re-apply the theme to all subplots
plot_ANPP20<- apply_theme(plot_ANPP20)
plot_ANPP21<- apply_theme(plot_ANPP21)
plot_ANPP22<- apply_theme(plot_ANPP22)
plot_ANPP23<- apply_theme(plot_ANPP23)
plot_ANPP24<- apply_theme(plot_ANPP24)
plot_ANPP_Overall <- apply_theme(plot_ANPP_Overall)


# Force alignment of all subplot axes (critical step)
aligned_plots <- align_plots(
  plot_ANPP20, plot_ANPP22,  plot_ANPP24,  
  plot_ANPP21, plot_ANPP23, plot_ANPP_Overall,
  align = "hv",  # Align both horizontally and vertically
  axis = "lr"    # Align left and right axes
)

# Arrange the aligned plots in a 3-row, 2-column layout
plot <- plot_grid(
  aligned_plots[[1]], aligned_plots[[4]],  # First row: A/B
  aligned_plots[[2]], aligned_plots[[5]],  # Second row: C/D
  aligned_plots[[3]], aligned_plots[[6]],  # Third row: E/F
  nrow = 3, ncol = 2,
  rel_heights = rep(1, 3),  # All rows have equal height
  rel_widths = c(1, 1),    # All columns have equal width
  # Force alignment and axis alignment
  align = "v",
  axis = "lr",
  hjust = 0, vjust = 0,    # Align to top-left
  # Add labels
  labels = c("A", "B", "C", "D", "E", "F"),
  label_size = 42,
  label_fontface = "plain",
  label_x = 0.01,  # Left-align labels
  label_y = 1.05   # Slightly shift up to avoid overlap
)

# Create a background plot with white background and margins
background_plot <- ggplot() +
  theme_void() +  # Remove default axes and grids
  theme(
    plot.background = element_rect(fill = "white", colour = NA),  # Set white background
    plot.margin = margin(5, 5, 5, 5, unit = "pt")  # 5pt margins on all sides
  )

# Add the combined plot as an annotation to the background layer
final_plot <- background_plot + 
  annotation_custom(
    ggplotGrob(plot),
    xmin = 0.02,  # 2% inner margin on the left
    xmax = 0.98,  # 2% inner margin on the right
    ymin = 0.02,  # 2% inner margin at the bottom
    ymax = 0.98   # 2% inner margin at the top
  )

# Output the plot
print(final_plot)
topptx(final_plot, filename = "Figure1.pptx", width = 20, height = 20*1.4)

#############################code for Figure 2


# loading Data,sheet=1 is rawdata, sheet=1 is NIE data
rawdata<-readxl::read_excel("./Raw data.xlsx",sheet=1,na=c("NA","#DIV/0!"))
NIEdata<-readxl::read_excel("./Raw data.xlsx",sheet=2,na=c("NA","#DIV/0!"))


str(rawdata)
rawdata$Year<-as.factor(rawdata$Year)
rawdata$Block<-as.factor(rawdata$Block)
rawdata$Plot<-as.factor(rawdata$Plot)
rawdata$Treatment<-as.factor(rawdata$Treatment)
rawdata$N<-as.factor(rawdata$N)
rawdata$P<-as.factor(rawdata$P)
str(rawdata)

str(NIEdata)
NIEdata$Year<-as.factor(NIEdata$Year)
NIEdata$Block<-as.factor(NIEdata$Block)
str(NIEdata)


#set the plot theme
library(ggplot2)
theme_custom_Bar <- function() {
  theme(
    axis.title.x=element_text(vjust=1,size=40) , 
    axis.title.y = element_text(size = 40, color = "black", margin = margin(r = 0.25, unit = "cm")),  
    axis.text.x= element_text(size=35,margin = unit(c(0.5, 0.5, 0.5, 0.5), 'cm'),angle=0),
    axis.text.y= element_text(size=35,margin = unit(c(0.5, 0.5, 0.5, 0.5), 'cm'),angle=0),
    text = element_text(family = "Times New Roman"),
    axis.line = element_line(color = "black", size = 1 ), 
    axis.title = element_text( color = "black"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black", size =1),
    axis.ticks.length.y = unit(-0.5, 'cm'), 
    axis.ticks.length.x = unit(0, 'cm'), 
    plot.title = element_text(color = "black"),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    line = element_line(size = 2) 
  )
}


#Multiple comparisons for CUR
str(rawdata)
#1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
library(lmerTest)
model <- lmer(CUR~ Year*Treatment+(1|Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#2-Multiple comparisons
#2.1-means of  Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
#2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha =0.05, reverse = TRUE)
print(cld_result5)



# Merge the mean sum and letter tags
library(Rmisc) 
Mean<- summarySE(rawdata, measurevar="CUR", groupvars=c("Treatment"),na=T)
Mean

cld_result5
#Merge
Mean_CUR<- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_CUR

# Remove the Spaces before and after the letters in the.group column
Mean_CUR$.group <- trimws(Mean_CUR$.group)
Mean_CUR

## NIE of CUR in 2020-2021
str(NIEdata)
#Normality test
shapiro.test(NIEdata$CUR)
#T-test
t.test(NIEdata$CUR,mu=0)
# Wilcoxon symbolic rank test
wilcox.test(NIEdata$CUR, mu = 0)

#p-value = 0.0007351, NIE>0***

#Draw a bar chart
str(Mean_CUR)
unique(Mean_CUR$Treatment)

Mean_CUR$Treatment= factor(Mean_CUR$Treatment, levels=c('Ctrl','N','P','NP'))

mycolor<-c("#A7C2B1","#6C91CB","#B0E0E6","#E69091")


# Bar plot of CUR
plot_CUR<-ggplot(Mean_CUR, aes(x = Treatment, y = CUR, fill = Treatment)) +
  #Bar 
  geom_bar(stat="identity",size=0.2,position="dodge",width=0.8,alpha = 0.9)+
  #Error bar
  geom_errorbar(aes(ymax = CUR+se,ymin =CUR-se),
                position = position_dodge(0.8), width = 0.15,linewidth = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(
    aes(x = Treatment, y = CUR + 2 * se, label = .group),
    size = 8,
    position = position_dodge(0.8),
    hjust = 0.5,
    vjust = 0)+
  annotate("text",   x = 0.8, y = 16 - 4/2, 
           hjust = 0,  vjust = 0.2,   label = "Syneigistic (NIE > 0***)",  size = 8, 
           color = "black",   family = "Times New Roman") +
  scale_y_continuous( expand = c(0, 0),breaks = seq(0,16,4),limits = c(0,16))+
  xlab("Treatment")+
  ylab("CUR (μmol CO2 m–2 s–1)")+
  theme(legend.position="none")+
  theme_custom_Bar ()


plot_CUR

# Multiple comparisons for LUE
str(rawdata)
# 1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
library(lmerTest)
model <- lmer(LUE ~ Year * Treatment + (1 | Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

# 2-Multiple comparisons
# 2.1-means of Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
# 2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha = 0.05, reverse = TRUE)
print(cld_result5)

# Merge the mean sum and letter tags
library(Rmisc) 
Mean <- summarySE(rawdata, measurevar = "LUE", groupvars = c("Treatment"), na = T)
Mean

cld_result5
# Merge
Mean_LUE <- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_LUE

# Remove the Spaces before and after the letters in the.group column
Mean_LUE$.group <- trimws(Mean_LUE$.group)
Mean_LUE

## NIE of LUE in 2020 - 2021
str(NIEdata)
# Normality test
shapiro.test(NIEdata$LUE)
# T-test
t.test(NIEdata$LUE, mu = 0)
# Wilcoxon symbolic rank test
wilcox.test(NIEdata$LUE, mu = 0)


#p-value = 0.0008033，NIE>0***

# Draw a bar chart
str(Mean_LUE)
unique(Mean_LUE$Treatment)

Mean_LUE$Treatment <- factor(Mean_LUE$Treatment, levels = c('Ctrl', 'N', 'P', 'NP'))

mycolor <- c("#A7C2B1", "#6C91CB", "#B0E0E6", "#E69091")

# Bar plot of LUE
plot_LUE <- ggplot(Mean_LUE, aes(x = Treatment, y = LUE, fill = Treatment)) +
  # Bar 
  geom_bar(stat = "identity", size = 0.2, position = "dodge", width = 0.8, alpha = 0.9) +
  # Error bar
  geom_errorbar(aes(ymax = LUE + se, ymin = LUE - se),
                position = position_dodge(0.8), width = 0.15, linewidth = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(
    aes(x = Treatment, y = LUE + 2 * se, label =.group),
    size = 8,
    position = position_dodge(0.8),
    hjust = 0.5,
    vjust = 0) +
  annotate("text", x = 0.8, y = 10 - 2.5 / 2, 
           hjust = 0, vjust = 0.2, label = "Syneigistic (NIE > 0***)", size = 8, 
           color = "black", family = "Times New Roman") +
   scale_y_continuous( expand = c(0, 0),breaks = seq(0,10, 2.5),limits = c(0,10))+
  xlab("Treatment")+
  ylab("LUE (μmol CO2 mmol-1 photon)") +
  theme(legend.position="none")+
  theme_custom_Bar ()

plot_LUE 


# Multiple comparisons for WUE
str(rawdata)
# 1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
library(lmerTest)
model <- lmer(WUE ~ Year * Treatment + (1 | Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

# 2-Multiple comparisons
# 2.1-means of Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
# 2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha = 0.05, reverse = TRUE)
print(cld_result5)

# Merge the mean sum and letter tags
library(Rmisc) 
Mean <- summarySE(rawdata, measurevar = "WUE", groupvars = c("Treatment"), na = T)
Mean

cld_result5
# Merge
Mean_WUE <- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_WUE

# Remove the Spaces before and after the letters in the.group column
Mean_WUE$.group <- trimws(Mean_WUE$.group)
Mean_WUE

## NIE of WUE in 2020 - 2021
str(NIEdata)
# Normality test
shapiro.test(NIEdata$WUE)
# T-test
t.test(NIEdata$WUE, mu = 0)
# Wilcoxon symbolic rank test
wilcox.test(NIEdata$WUE, mu = 0)


#p-value = 0.001256,NIE>0**

# Draw a bar chart
str(Mean_WUE)
unique(Mean_WUE$Treatment)

Mean_WUE$Treatment <- factor(Mean_WUE$Treatment, levels = c('Ctrl', 'N', 'P', 'NP'))

mycolor <- c("#A7C2B1", "#6C91CB", "#B0E0E6", "#E69091")

# Bar plot of WUE
plot_WUE <- ggplot(Mean_WUE, aes(x = Treatment, y = WUE, fill = Treatment)) +
  # Bar 
  geom_bar(stat = "identity", size = 0.2, position = "dodge", width = 0.8, alpha = 0.9) +
  # Error bar
  geom_errorbar(aes(ymax = WUE + se, ymin = WUE - se),
                position = position_dodge(0.8), width = 0.15, linewidth = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(
    aes(x = Treatment, y = WUE + 2 * se, label =.group),
    size = 8,
    position = position_dodge(0.8),
    hjust = 0.5,
    vjust = 0) +
  annotate("text", x = 0.8, y = 16 - 4 / 2, 
           hjust = 0, vjust = 0.2, label = "Syneigistic (NIE > 0**)", size = 8, 
           color = "black", family = "Times New Roman") +
  scale_y_continuous( expand = c(0, 0),breaks = seq(0,16,4),limits = c(0,16))+
  xlab("Treatment") +
  ylab("WUE (μmol CO2 mmol-1 H2O)") +
  theme(legend.position="none")+
  theme_custom_Bar()
plot_WUE


# Multiple comparisons for RUE
str(rawdata)
# 1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
library(lmerTest)
model <- lmer(RUE ~ Year * Treatment + (1 | Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

# p-value = 0.01809 means the residuals is not normal , so log the data
model <- lmer(log(RUE) ~ Year * Treatment + (1 | Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

# 2-Multiple comparisons
# 2.1-means of Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
# 2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha = 0.05, reverse = TRUE)
print(cld_result5)

# Merge the mean sum and letter tags
library(Rmisc) 
Mean <- summarySE(rawdata, measurevar = "RUE", groupvars = c("Treatment"), na = T)
Mean

cld_result5
# Merge
Mean_RUE <- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_RUE

# Remove the Spaces before and after the letters in the.group column
Mean_RUE$.group <- trimws(Mean_RUE$.group)
Mean_RUE

## NIE of RUE in 2020 - 2021
str(NIEdata)
# Normality test
shapiro.test(NIEdata$RUE)
# T-test
t.test(NIEdata$RUE, mu = 0)
# Wilcoxon symbolic rank test
wilcox.test(NIEdata$RUE, mu = 0)

#p-value = 3.782e-05,NIE>0***

# Draw a bar chart
str(Mean_RUE)
unique(Mean_RUE$Treatment)

Mean_RUE$Treatment <- factor(Mean_RUE$Treatment, levels = c('Ctrl', 'N', 'P', 'NP'))

mycolor <- c("#A7C2B1", "#6C91CB", "#B0E0E6", "#E69091")

# Bar plot of RUE
plot_RUE <- ggplot(Mean_RUE, aes(x = Treatment, y = RUE, fill = Treatment)) +
  # Bar 
  geom_bar(stat = "identity", size = 0.2, position = "dodge", width = 0.8, alpha = 0.9) +
  # Error bar
  geom_errorbar(aes(ymax = RUE + se, ymin = RUE - se),
                position = position_dodge(0.8), width = 0.15, linewidth = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(
    aes(x = Treatment, y = RUE + 2 * se, label =.group),
    size = 8,
    position = position_dodge(0.8),
    hjust = 0.5,
    vjust = 0) +
  annotate("text", x = 0.8, y = 2.4 - 0.6 / 2, 
           hjust = 0, vjust = 0.2, label = "Syneigistic (NIE > 0***)", size = 8, 
           color = "black", family = "Times New Roman") +
  scale_y_continuous( expand = c(0, 0),breaks = seq(0,2.4, 0.6),limits = c(0,2.4))+
  xlab("Treatment") +
  ylab("RUE (g NPP mm-1 rainfall)")+
  theme(legend.position="none")+
  theme_custom_Bar()
plot_RUE



# Multiple comparisons for NUE_ANPP
str(rawdata)
# 1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
library(lmerTest)
model <- lmer(NUE_ANPP ~ Year * Treatment + (1 | Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

# 2-Multiple comparisons
# 2.1-means of Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
# 2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha = 0.05, reverse = TRUE)
print(cld_result5)

# Merge the mean sum and letter tags
library(Rmisc) 
Mean <- summarySE(rawdata, measurevar = "NUE_ANPP", groupvars = c("Treatment"), na = T)
Mean

cld_result5
# Merge
Mean_NUE_ANPP <- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_NUE_ANPP

# Remove the Spaces before and after the letters in the.group column
Mean_NUE_ANPP$.group <- trimws(Mean_NUE_ANPP$.group)
Mean_NUE_ANPP

## NIE of NUE_ANPP in 2020 - 2021
str(NIEdata)
# Normality test
shapiro.test(NIEdata$NUE_ANPP)
# T-test
t.test(NIEdata$NUE_ANPP, mu = 0)
# Wilcoxon symbolic rank test
wilcox.test(NIEdata$NUE_ANPP, mu = 0)

#p-value = 0.2258,NIE=0

# Draw a bar chart
str(Mean_NUE_ANPP)
unique(Mean_NUE_ANPP$Treatment)

Mean_NUE_ANPP$Treatment <- factor(Mean_NUE_ANPP$Treatment, levels = c('Ctrl', 'N', 'P', 'NP'))

mycolor <- c("#A7C2B1", "#6C91CB", "#B0E0E6", "#E69091")

# Bar plot of NUE_ANPP
plot_NUE_ANPP <- ggplot(Mean_NUE_ANPP, aes(x = Treatment, y = NUE_ANPP, fill = Treatment)) +
  # Bar 
  geom_bar(stat = "identity", size = 0.2, position = "dodge", width = 0.8, alpha = 0.9) +
  # Error bar
  geom_errorbar(aes(ymax = NUE_ANPP + se, ymin = NUE_ANPP - se),
                position = position_dodge(0.8), width = 0.15, linewidth = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(
    aes(x = Treatment, y = NUE_ANPP + 2 * se, label =.group),
    size = 8,
    position = position_dodge(0.8),
    hjust = 0.5,
    vjust = 0) +
  annotate("text", x = 0.8, y = 160 - 40 / 2, 
           hjust = 0, vjust = 0.2, label = "Additive (NIE = 0)", size = 8, 
           color = "black", family = "Times New Roman") +
  scale_y_continuous( expand = c(0, 0),breaks = seq(0,160, 40),limits = c(0,160))+
  xlab("Treatment")+
  ylab("NUE (g ANPP g-1 N)")+
  theme(legend.position="none")+
  theme_custom_Bar()

plot_NUE_ANPP

# Multiple comparisons for PUE_ANPP
str(rawdata)
# 1- Effect of Treatment (Ctrl, N, P and NP)
library(lme4)
library(lmerTest)
model <- lmer(PUE_ANPP ~ Year * Treatment + (1 | Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)

#p-value = 2.517e-09 means residuals is not normal, so log the data
model <- lmer(log(PUE_ANPP) ~ Year * Treatment + (1 | Block), data = rawdata)
anova(model)
# Extract the model residuals
residuals <- residuals(model, type = "response")  
# Test for normality of model residuals
shapiro.test(residuals)


# 2-Multiple comparisons
# 2.1-means of Treatment
library(emmeans)
emm <- emmeans(model, ~ Treatment)
print(emm)
# 2.2-Tukey multiple comparisons
tukey_result <- pairs(emm, adjust = "tukey")
print(tukey_result)
# 2.3-Generate letter tags
library(multcomp)
cld_result5 <- cld(emm, adjust = "tukey", Letters = letters, alpha = 0.05, reverse = TRUE)
print(cld_result5)

# Merge the mean sum and letter tags
library(Rmisc) 
Mean <- summarySE(rawdata, measurevar = "PUE_ANPP", groupvars = c("Treatment"), na = T)
Mean

cld_result5
# Merge
Mean_PUE_ANPP <- merge(Mean, cld_result5[, c("Treatment", "emmean", "SE", ".group")], by = "Treatment")
Mean_PUE_ANPP

# Remove the Spaces before and after the letters in the.group column
Mean_PUE_ANPP$.group <- trimws(Mean_PUE_ANPP$.group)
Mean_PUE_ANPP

## NIE of PUE_ANPP in 2020 - 2021
str(NIEdata)
# Normality test
shapiro.test(NIEdata$PUE_ANPP)
# T-test
t.test(NIEdata$PUE_ANPP, mu = 0)
# Wilcoxon symbolic rank test
wilcox.test(NIEdata$PUE_ANPP, mu = 0)

#p-value = 0.4497,NIE=0

# Draw a bar chart
str(Mean_PUE_ANPP)
unique(Mean_PUE_ANPP$Treatment)

Mean_PUE_ANPP$Treatment <- factor(Mean_PUE_ANPP$Treatment, levels = c('Ctrl', 'N', 'P', 'NP'))

mycolor <- c("#A7C2B1", "#6C91CB", "#B0E0E6", "#E69091")

# Bar plot of NUE_ANPP
plot_PUE_ANPP <- ggplot(Mean_PUE_ANPP, aes(x = Treatment, y = PUE_ANPP, fill = Treatment)) +
  # Bar 
  geom_bar(stat = "identity", size = 0.2, position = "dodge", width = 0.8, alpha = 0.9) +
  # Error bar
  geom_errorbar(aes(ymax = PUE_ANPP + se, ymin = PUE_ANPP - se),
                position = position_dodge(0.8), width = 0.15, linewidth = 0.5) +
  scale_fill_manual(values = mycolor) +
  geom_text(
    aes(x = Treatment, y = PUE_ANPP + 2 * se, label =.group),
    size = 8,
    position = position_dodge(0.8),
    hjust = 0.5,
    vjust = 0) +
  annotate("text", x = 0.8, y = 1600 - 400 / 2, 
           hjust = 0, vjust = 0.2, label = "Additive (NIE = 0)", size = 8, 
           color = "black", family = "Times New Roman") +
  scale_y_continuous( expand = c(0, 0),breaks = seq(0,1600, 400),limits = c(0,1600))+
  xlab("Treatment")+
  ylab("PUE (g ANPP g-1 P)")+
  theme(legend.position="none")+
  theme_custom_Bar()

plot_PUE_ANPP


plot_CUR
plot_LUE 
plot_WUE
plot_RUE
plot_NUE_ANPP
plot_PUE_ANPP



library(cowplot)
library(eoffice)
library(ggplot2)

# Set global font and unified theme
custom_theme <- theme(
  # Set font to Times New Roman, non-bold
  text = element_text(family = "Times New Roman", face = "plain"),
  # Uniform margin between axis titles and the plot
  axis.title.y = element_text(margin = margin(r = 45)),  # 45pt right margin for Y-axis title
  axis.title.x = element_text(margin = margin(t = 55)),  # 55pt top margin for X-axis title
  # Uniform plot margins (critical! Ensure consistent height and width)
  plot.margin = margin(15, 15, 15, 15, "pt"),
  # Remove default theme distractions
  panel.grid = element_blank(),
  plot.background = element_rect(fill = "white"),
  panel.background = element_rect(fill = "white")
)

# Apply the unified theme to all plots
apply_theme <- function(plot) {
  plot + custom_theme
}


# Re-apply the theme to all subplots
plot_CUR<- apply_theme(plot_CUR)
plot_LUE <- apply_theme(plot_LUE)
plot_WUE <- apply_theme(plot_WUE)
plot_RUE <- apply_theme(plot_RUE)
plot_NUE<- apply_theme(plot_NUE_ANPP)
plot_PUE<- apply_theme(plot_PUE_ANPP)

# Force alignment of all subplot axes (critical step)
aligned_plots <- align_plots(
  plot_CUR, plot_WUE,plot_NUE, 
  plot_LUE , plot_RUE, plot_PUE,
  align = "hv",  # Align both horizontally and vertically
  axis = "lr"    # Align left and right axes
)

# Arrange the aligned plots in a 3-row, 2-column layout
plot <- plot_grid(
  aligned_plots[[1]], aligned_plots[[4]], 
  aligned_plots[[2]], aligned_plots[[5]],  
  aligned_plots[[3]], aligned_plots[[6]],  
  nrow = 3, ncol = 2,
  rel_heights = rep(1, 3),  # All rows have equal height
  rel_widths = c(1, 1),    # All columns have equal width
  # Force alignment and axis alignment
  align = "v",
  axis = "lr",
  hjust = 0, vjust = 0,    # Align to top-left
  # Add labels
  labels = c("A", "B", "C", "D", "E", "F"),
  label_size = 42,
  label_fontface = "plain",
  label_x = 0.01,  # Left-align labels
  label_y = 1.05   # Slightly shift up to avoid overlap
)

# Create a background plot with white background and margins
background_plot <- ggplot() +
  theme_void() +  # Remove default axes and grids
  theme(
    plot.background = element_rect(fill = "white", colour = NA),  # Set white background
    plot.margin = margin(5, 5, 5, 5, unit = "pt")  # 5pt margins on all sides
  )

# Add the combined plot as an annotation to the background layer
final_plot <- background_plot + 
  annotation_custom(
    ggplotGrob(plot),
    xmin = 0.02,  # 2% inner margin on the left
    xmax = 0.98,  # 2% inner margin on the right
    ymin = 0.02,  # 2% inner margin at the bottom
    ymax = 0.98   # 2% inner margin at the top
  )

# Output the plot
print(final_plot)
topptx(final_plot, filename = "Figure2.pptx", width =20, height = 20*1.3)

########################## code for Figure 3

# loading Data,sheet=1 is rawdata, sheet=1 is NIE data
NIEdata<-readxl::read_excel("./Raw data.xlsx",sheet=2,na=c("NA","#DIV/0!"))


str(NIEdata)
NIEdata$Year<-as.factor(NIEdata$Year)
NIEdata$Block<-as.factor(NIEdata$Block)
str(NIEdata)




#Darw the line plot

#set theme
library(ggplot2)
theme_custom_line<- function() {
  theme(
    axis.title.x=element_text(vjust=1,size=40) , 
    axis.title.y = element_text(size = 40, color = "black", margin = margin(r = 0.25, unit = "cm")),  
    axis.text.x= element_text(size=35,margin = unit(c(0.5, 0.5, 0.5, 0.5), 'cm'),angle=0),
    axis.text.y= element_text(size=35,margin = unit(c(0.5, 0.5, 0.5, 0.5), 'cm'),angle=0),
    text = element_text(family = "Times New Roman"),
    axis.line = element_line(color = "black", size = 1 ), 
    axis.title = element_text( color = "black"),
    axis.text = element_text(color = "black"),
    axis.ticks = element_line(color = "black", size =1),
    axis.ticks.length.y = unit(-0.5, 'cm'), 
    axis.ticks.length.x = unit(-0.5, 'cm'), 
    plot.title = element_text(color = "black"),
    panel.border = element_blank(),
    panel.background = element_blank(),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    line = element_line(size = 2) 
  )
}

library(nlme)
library(cowplot)
library(MuMIn)

# add "NIE" 
x_var <- "NIE"  # add "NIE" for x_var
y_var <- "NIE"  # add "NIE" for y_var
title <- ""     #no add anything for title 

# Define a function for analysis and plotting
analyze_and_plot <- function(y_var, x_var, title) {
  
  formula <- as.formula(paste(y_var, "~", x_var))
  total_model <- lme(formula, random = ~1 | Year, data = subdata)
  total_p_value <- summary(total_model)$tTable[2, 5]
  
  # R2
  r_squared <- r.squaredGLMM(total_model)
  total_r_squared_marginal <- r_squared[1] # marginal R2
  total_r_squared_conditional <- r_squared[2]# conditionalR2
  
  # The line type of the total fitting line is determined based on the total P-value
  total_linetype <- ifelse(total_p_value < 0.05, "solid", "dashed")
  
  # Create the basic drawing object
  p <- ggplot(subdata, aes_string(x = x_var, y = y_var)) +
    geom_point(shape = 1, size = 4.5, color = "#000000", stroke = 1) +  #  stroke to thicken the lines
    labs(x = paste("NIE of", x_var), y = paste("NIE of", y_var), title = title) +
    theme_custom_line() 
  
  # Add the total fitting line
  p <- p + geom_smooth(method = "lm", se = T,size =2, linetype = total_linetype,color = "#000000", fill = "#ADD8E6")
  r_squared_formatted <- sprintf("%.2f", total_r_squared_marginal)
  
  # add R² and P
  text_grob <- grid::textGrob(
    label = paste0("R² = ", r_squared_formatted, 
                   "\nP = ", format(round(total_p_value, 3), nsmall = 3)),#Format the R² value to ensure there are three decimal places
    x = grid::unit(0.25, "npc"),
    y = grid::unit(0.75, "npc"),
    hjust = 0,
    vjust = 1,
    gp = grid::gpar(fontsize = 12 *.pt, family = "Times New Roman")
  )
  
  # Add annotations
  p <- p + annotation_custom(grob = text_grob)
  
  return(p)
}


# Select specified columns from NIEdata and remove rows with NA values
library(dplyr)
subdata<- NIEdata %>% 
  dplyr::select(Year, Block,  ANPP,CUR, LUE, WUE) %>% 
  na.omit()  # Remove rows containing any NA values

str(subdata)
# Draw line plot for CUR
# 分别绘制三个图
line_CUR <- analyze_and_plot("ANPP", "CUR", "") +
  scale_x_continuous(expand = c(0, 0),limits = c(0, 3.0), breaks = seq(0, 3.0, 0.75)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-10, 110), breaks = seq(-10, 110, 30))
line_CUR

line_LUE <- analyze_and_plot("ANPP", "LUE", "") +
  scale_x_continuous(expand = c(0, 0),limits = c(0, 1.6), breaks = seq(0, 1.6, 0.4)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-10, 110), breaks = seq(-10, 110, 30))

line_LUE

line_WUE<- analyze_and_plot("ANPP", "WUE", "") +
  scale_x_continuous(expand = c(0, 0),limits = c(0, 10), breaks = seq(0, 10, 2.5)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-10, 110), breaks = seq(-10, 110, 30))

line_WUE

# Select specified columns from NIEdata and remove rows with NA values
library(dplyr)
subdata<- NIEdata %>% 
  dplyr::select(Year, Block,  ANPP,RUE) %>% 
  na.omit()  # Remove rows containing any NA values

str(subdata)

line_RUE <- analyze_and_plot("ANPP", "RUE", "")   +
  scale_x_continuous(expand = c(0, 0),limits = c(-0.6, 1.0), breaks = seq(-0.6,1.0, 0.4)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-50, 190), breaks = seq(-50, 190, 60))

line_RUE


# Select specified columns from NIEdata and remove rows with NA values
library(dplyr)
subdata<- NIEdata %>% 
  dplyr::select(Year, Block,  ANPP, NUE_ANPP,PUE_ANPP) %>% 
  na.omit()  # Remove rows containing any NA values

str(subdata)

line_NUE <- analyze_and_plot("ANPP", "NUE_ANPP", "")   +
  scale_x_continuous(expand = c(0, 0),limits = c(-50, 30), breaks = seq(-50, 30, 20)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-20, 180), breaks = seq(-20, 180, 50))

line_NUE


line_PUE <- analyze_and_plot("ANPP", "PUE_ANPP", "")   +
  scale_x_continuous(expand = c(0, 0),limits = c(-600, 400), breaks = seq(-600, 400, 250)) +
  scale_y_continuous(expand = c(0, 0),limits = c(-20, 180), breaks = seq(-20, 180, 50))

line_PUE



library(cowplot)
library(eoffice)
library(ggplot2)

# Set global font and unified theme
custom_theme <- theme(
  # Set font to Times New Roman, non-bold
  text = element_text(family = "Times New Roman", face = "plain"),
  # Uniform margin between axis titles and the plot
  axis.title.y = element_text(margin = margin(r = 45)),  # 45pt right margin for Y-axis title
  axis.title.x = element_text(margin = margin(t = 55)),  # 55pt top margin for X-axis title
  # Uniform plot margins (critical! Ensure consistent height and width)
  plot.margin = margin(15, 15, 15, 15, "pt"),
  # Remove default theme distractions
  panel.grid = element_blank(),
  plot.background = element_rect(fill = "white"),
  panel.background = element_rect(fill = "white")
)

# Apply the unified theme to all plots
apply_theme <- function(plot) {
  plot + custom_theme
}

# Re-apply the theme to all subplots
plot_CUR <- apply_theme(line_CUR)
plot_LUE <- apply_theme(line_LUE)
plot_WUE <- apply_theme(line_WUE)
plot_RUE <- apply_theme(line_RUE)
plot_NUE <- apply_theme(line_NUE)
plot_PUE <- apply_theme(line_PUE)


# Force alignment of all subplot axes (critical step)
aligned_plots <- align_plots(
  plot_CUR, plot_WUE, plot_NUE, 
  plot_LUE, plot_RUE, plot_PUE,
  align = "hv",  # Align both horizontally and vertically
  axis = "lr"    # Align left and right axes
)

# Arrange the aligned plots in a 3-row, 2-column layout
plot <- plot_grid(
  aligned_plots[[1]], aligned_plots[[4]],  # First row: A/B
  aligned_plots[[2]], aligned_plots[[5]],  # Second row: C/D
  aligned_plots[[3]], aligned_plots[[6]],  # Third row: E/F
  nrow = 3, ncol = 2,
  rel_heights = rep(1, 3),  # All rows have equal height
  rel_widths = c(1, 1),    # All columns have equal width
  # Force alignment and axis alignment
  align = "v",
  axis = "lr",
  hjust = 0, vjust = 0,    # Align to top-left
  # Add labels
  labels = c("A", "B", "C", "D", "E", "F"),
  label_size = 42,
  label_fontface = "plain",
  label_x = 0.01,  # Left-align labels
  label_y = 1.05   # Slightly shift up to avoid overlap
)

# Create a background plot with white background and margins
background_plot <- ggplot() +
  theme_void() +  # Remove default axes and grids
  theme(
    plot.background = element_rect(fill = "white", colour = NA),  # Set white background
    plot.margin = margin(5, 5, 5, 5, unit = "pt")  # 5pt margins on all sides
  )

# Add the combined plot as an annotation to the background layer
final_plot <- background_plot + 
  annotation_custom(
    ggplotGrob(plot),
    xmin = 0.02,  # 2% inner margin on the left
    xmax = 0.98,  # 2% inner margin on the right
    ymin = 0.02,  # 2% inner margin at the bottom
    ymax = 0.98   # 2% inner margin at the top
  )

# Output the plot
print(final_plot)
topptx(final_plot, filename = "Figure3.pptx", width = 20, height = 20*1.4)