# load packages
library(openxlsx)
library(tidyverse)
library(nlme)
library(minpack.lm)
library(randomForest)
library(rfUtilities)
library(rfPermute)
library(glmm.hp)
library(ggpubr)
library(lme4)
library(lmerTest)

# read data
data <- read.xlsx("data.xlsx")

grain_biomass_c <- data[data$`Crop Component` == "Grain", ] 
straw_biomass_c <- data[data$`Crop Component` == "Straw", ]

# normal distribution test
shapiro.test(grain_biomass_c$`Crop biomass carbon`)
shapiro.test(straw_biomass_c$`Crop biomass carbon`)

# random forest
Straw_rf <- randomForest(`Crop biomass carbon` ~ ., data = straw_biomass_c, importance = TRUE, proximity = TRUE)
Straw_rf_im <- rf.significance(Straw_rf, straw_biomass_c, nperm = 100, ntree = 500)
Straw_rf_im_var <- rfPermute(`Crop biomass carbon` ~ ., data = straw_biomass_c, ntree = 500, nrep = 100)

Grain_rf <- randomForest(`Crop biomass carbon` ~ ., data = grain_biomass_c, importance = TRUE, proximity = TRUE)
Grain_rf_im <- rf.significance(Grain_rf, grain_biomass_c, nperm = 100, ntree = 500)
Grain_rf_im_var <- rfPermute(`Crop biomass carbon` ~ ., data = grain_biomass_c, ntree = 500, nrep = 100)

# hierarchical proportion
Straw_lm <- lm(`Crop biomass carbon` ~ MAT + MAST + MGT + MGST + MWT + MWST, data = straw_biomass_c)
glmm.hp(Straw_lm)

Grain_lm <- lm(`Crop biomass carbon` ~ MAT + MAST + MGT + MGST + MWT + MWST, data = grain_biomass_c)
glmm.hp(Grain_lm)

# difference of crop biomass carbon between years in each climate and crop type
diff_between_ldgroups <- function (data, treat_reserve = NULL, multi_vars = NULL, 
                                   y = NULL, x = NULL, class_name = NULL, 
                                   var_name = NULL, FUN = "-", sep = "--") 
{
  
  if (is.null(class_name)) {
    stop("`class_name` must be specified.")
  }
  
  if (!is.null(multi_vars) && !is.null(x) && !is.null(y)) {
    stop("Both `multi_vars` parameter and `x` / `y` cannot be specified simultaneously.")
  }
  
  if ((!is.null(x) && is.null(y) | is.null(x) && !is.null(y))) {
    stop("Both `x` and `y` must be specified simultaneously.")
  }
  
  if (!is.null(multi_vars) && !is.null(var_name)) {
    stop("Both `multi_vars` and `var_name` cannot be specified simultaneously.")
  }
  
  require(tidyr)
  require(tidyverse)
  require(dplyr)
  
  results = NULL
  treat_names = NULL
  class_names = NULL
  var_names = NULL
  
  data = data[order(data[, class_name]), ]
  df = data %>% tidyr::gather(key = "vars", value = "values", 
                              multi_vars, na.rm = TRUE)
  x1 = sapply(unique(df[[class_name]]), 
              function(x) 1:nrow(df[df[[class_name]] == x,]))
  x2 = unlist(sapply(x1, function(x) x))
  
  if (is.null(df$vars)) {
    df = df %>% 
      rename(values = var_name) %>% 
      cbind(vars = var_name) %>% 
      cbind(x2) %>% 
      rename(x = x2) %>% 
      dplyr::select(x, all_of(treat_reserve), all_of(class_name), vars, values)
  } else {
    df = df[order(df[, class_name]), ]
    df = df %>% 
      cbind(x2) %>% 
      rename(x = x2) %>% 
      dplyr::select(x, all_of(treat_reserve), all_of(class_name), vars, values)
  }
  
  if (is.null(treat_reserve)) {
    rownames(df) <- 1:nrow(df)
    df1 = tidyr::spread(df, class_name, values)
    df2 = as.data.frame(df1[, -1])
  } else {
    total = tidyr::unite(df, treat_names, unique(treat_reserve), sep = sep)
    rownames(total) <- 1:nrow(total)
    df1 = tidyr::pivot_wider(total, names_from = class_name, values_from = values) %>% 
      arrange(vars)
    df2 = as.data.frame(df1[, -1])
  }
  
  if(is.null(treat_reserve) & !is.null(var_name)) {
    if( is.null(x) & is.null(y)){
      for (i in 1:(ncol(df2) - 1)) {
        for (j in 1:(ncol(df2) - 1)) {
          if (i <= j) next
          values <- c(outer((df2 %>% dplyr::select(-vars))[, i], 
                            (df2 %>% dplyr::select(-vars))[, j], 
                            FUN))
          results <- c(results, values)
          class_labels <- rep(paste(colnames((df2 %>% 
                                                dplyr::select(-vars)))[i], 
                                    colnames((df2 %>% 
                                                dplyr::select(-vars)))[j], sep = sep), 
                              times = length(values))
          class_names <- c(class_names, class_labels)
          var_labels <- rep(var_name, times = length(values))
          var_names <- c(var_names, var_labels)
          final <- data.frame(class_names, var_names, results)
          final <- na.omit(final)
        }
      }
      return(final)
    }  
    if (!is.null(x) & !is.null(y)) {
      values <- c(outer(df2[, y], df2[, x], FUN))
      class_names <- rep(paste(y, x, sep = sep), times = length(values))
      var_labels <- rep(var_name, times = length(values))
      var_names <- c(var_names, var_labels)
      final <- data.frame(class_names, var_names, values)
      final <- na.omit(final)
    }
    return(final)
  }  
  
  if (is.null(treat_reserve) & !is.null(multi_vars)){
    if (is.null(x) & is.null(y)) {
      for (i in 1:(ncol(df2) - 1)) {
        for (j in 1:(ncol(df2) - 1)) {
          for (k in unique(df2$vars)) {
            if (i <= j)  next
            values <- c(outer((df2 %>% dplyr::filter(vars == k) %>% 
                                 dplyr::select(-vars))[, i], 
                              (df2 %>% dplyr::filter(vars == k) %>% 
                                 dplyr::select(-vars))[, j], FUN))
            results <- c(results, values)
            class_labels <- rep(paste(colnames(df2 %>% 
                                                 dplyr::filter(vars == k) %>% 
                                                 dplyr::select(-vars))[i], 
                                      colnames(df2 %>% 
                                                 dplyr::filter(vars == k) %>% 
                                                 dplyr::select(-vars))[j], sep = sep), 
                                times = length(values))
            class_names <- c(class_names, class_labels)
            var_labels <- rep(k, times = length(values))
            var_names <- c(var_names, var_labels)
            fianl <- data.frame(class_names, var_names, results) %>% arrange(var_names)
            final <- na.omit(final)
          }
        }
      }
      return(fianl)
    }
    
    if (!is.null(x) & !is.null(y)) {
      for (k in 1:unique(df2$vars)) {
        values <- c(outer((df2 %>% 
                             dplyr::filter(vars == k) %>% 
                             dplyr::select(-vars))[, y], 
                          (df2 %>% dplyr::filter(vars == k) %>% 
                             dplyr::select(-vars))[, x], FUN))
        results <- c(results, values)
        class_labels <- rep(paste(y, x, sep = sep), 
                            times = length(values))
        class_names <- c(class_names, class_labels)
        var_labels <- rep(k, times = length(values))
        var_names <- c(var_names, var_labels)
        final <- data.frame(class_names, var_names, results) %>% arrange(var_names)
        final <- na.omit(final)
      }
      return(final)
    }
  }
  
  if (is.null(data$vars) & (!is.null(multi_vars) | !is.null(var_name)) & 
      !is.null(treat_reserve) & !is.null(class_name)) {
    if (is.null(x) & is.null(y)) {
      for (i in 1:(ncol(df2) - 2)) {
        for(j in 1:(ncol(df2) - 2)){
          for (k in unique(df2$treat_names)) {
            for (m in unique(df2$vars)) {
              if (i <= j) next
              values <- c(outer((df2 %>%
                                   dplyr::filter(treat_names == k & vars == m) %>% 
                                   dplyr::select(-vars, -treat_names))[, i], 
                                (df2 %>% dplyr::filter(treat_names == k & vars == m) %>%
                                   dplyr::select(-vars,-treat_names))[, j], FUN))
              results = c(results, values)
              class_labels = rep(paste(colnames(df2 %>% 
                                                  dplyr::filter(treat_names == k & vars == m) %>%
                                                  dplyr::select(-vars, -treat_names))[i], 
                                       colnames(df2 %>% dplyr::filter(treat_names == k & vars == m) %>% 
                                                  dplyr::select(-vars, -treat_names))[j], sep = sep), 
                                 times = length(values))
              class_names = c(class_names, class_labels)
              treat_labels = rep(k, times = length(values))
              treat_names = c(treat_names, treat_labels)
              var_labels = rep(m, times = length(values))
              var_names <- c(var_names, var_labels)
              final = data.frame(class_names, var_names, treat_names, results) %>% 
                arrange(treat_names, var_names)
              final <- na.omit(final)
            }
          }
        }
      }
      final <- tidyr::separate(final, col = treat_names, into = unique(treat_reserve), sep = sep)
      return(final)
    }
    
    if (!is.null(x) & !is.null(y)) {
      for (k in unique(df2$treat_names)) {
        for (m in unique(df2$vars)) {
          values <- c(outer((df2 %>%
                               dplyr::filter(treat_names ==  k & 
                                               vars == m) %>% 
                               dplyr::select(-vars, -treat_names))[, y], 
                            (df2 %>% dplyr::filter(treat_names == k & 
                                                     vars == m) %>% 
                               dplyr::select(-vars, -treat_names))[, x], FUN))
          results <- c(results, values)
          class_labels <- rep(paste(y, x, sep = sep), 
                              times = length(values))
          class_names <- c(class_names, class_labels)
          treat_labels <- rep(k, 
                              times = length(values))
          treat_names <- c(treat_names, treat_labels)
          var_labels <- rep(m, times = length(values))
          var_names <- c(var_names, var_labels)
          final <- data.frame(class_names, var_names, treat_names, results) %>% 
            arrange(treat_names, var_names)
          final <- na.omit(final)
        }
      }
      final <- tidyr::separate(final, col = treat_names, into = unique(treat_reserve), sep = sep)
      return(final)
    }
  }
}

Diff_crop_biomass_c <- diff_between_ldgroups(grain_biomass_c, treat_reserve = c("Article ID", "Crop type", "Climate type"), 
                                             multi_vars = c(11:16, 33),
                                             y = NULL, x = NULL, class_name = "Time",
                                             var_name = NULL, FUN = "-", sep = "--")


# model determination
crop_mwst_glm <- glm(`Crop biomass carbon` ~ MWST,
                      data = Diff_crop_biomass_c,
                      family = gaussian)

crop_mwst_nls1 <- nlsLM(`Crop biomass carbon` ~ a * MWST^2 + b * MWST + c,
                         data = Diff_crop_biomass_c,
                         start = c(a = 1, b = 1, c = 1))

crop_mwst_nls2 <- nlsLM(`Crop biomass carbon` ~ a*exp(b*MWST),
                         data = Diff_crop_biomass_c,
                         start = c(a = 1, b = 1))

crop_mwst_nls3 <- nlsLM(`Crop biomass carbon` ~ a*log(MWST)+b,
                         data = Diff_crop_biomass_c,
                         start = c(a = 1, b = 1))

crop_mwst_nls4 <- nlsLM(`Crop biomass carbon` ~ a*MWST^b,
                         data = Diff_crop_biomass_c,
                         start = c(a = 1, b = 1))

crop_mwst_lme <- lme(`Crop biomass carbon` ~ MWST,
                      data = Diff_crop_biomass_c,
                      random = list(Climate = ~ 1, 
                                    Crop = ~ 1),
                      method = "REML",
                      control = lmeControl(maxIter = 100))

crop_mwst_nlme1 <- nlme(`Crop biomass carbon` ~ a * MWST^2 + b * MWST + c,
                         fixed = a + b + c ~ 1,
                         data = Diff_crop_biomass_c,
                         random = list(Climate = a ~ 1, 
                                       Crop = a ~ 1),
                         method = "REML", 
                         start = c(a = 1, b = 1, c = 1),
                         control = lmeControl(maxIter = 100))

crop_mwst_nlme2 <- nlme(`Crop biomass carbon` ~ a * MWST^b,
                         fixed = a + b ~ 1,
                         data = Diff_crop_biomass_c,
                         random = list(Climate = a ~ 1, 
                                       Crop = a ~ 1),
                         method = "REML", 
                         start = c(a = 1, b = 2),
                         control = lmeControl(maxIter = 100))

crop_mwst_nlme3 <- nlme(`Crop biomass carbon` ~  a * exp(b * MWST),
                         fixed = a + b ~ 1,
                         data = Diff_crop_biomass_c,
                         random = list(Climate = a ~ 1, 
                                       Crop = a ~ 1),
                         method = "REML", 
                         start = c(a = 1, b = 1),
                         control = lmeControl(maxIter = 100))

crop_mwst_nlme4 <- nlme(`Crop biomass carbon` ~  a*log(MWST)+b,
                         fixed = a + b ~ 1,
                         data = Diff_crop_biomass_c,
                         random = list(Climate = a ~ 1, 
                                       Crop = a ~ 1),
                         method = "REML", 
                         start = c(a = 1, b = 1),
                         control = lmeControl(maxIter = 100))

AIC(logLik(crop_mwst_glm))
AIC(logLik(crop_mwst_nls1))
AIC(logLik(crop_mwst_nls2))
AIC(logLik(crop_mwst_nls3))
AIC(logLik(crop_mwst_nls4))
AIC(logLik(crop_mwst_lme))
AIC(logLik(crop_mwst_nlme1))
AIC(logLik(crop_mwst_nlme2))
AIC(logLik(crop_mwst_nlme3))
AIC(logLik(crop_mwst_nlme4))

BIC(logLik(crop_mwst_glm))
BIC(logLik(crop_mwst_nls1))
BIC(logLik(crop_mwst_nls2))
BIC(logLik(crop_mwst_nls3))
BIC(logLik(crop_mwst_nls4))
BIC(logLik(crop_mwst_lme))
BIC(logLik(crop_mwst_nlme1))
BIC(logLik(crop_mwst_nlme2))
BIC(logLik(crop_mwst_nlme3))
BIC(logLik(crop_mwst_nlme4))

# AIC calculations
null1 <- lmer(`Crop biomass carbon` ~ (1 | Climate ) , data = Diff_crop_biomass_c)
AIC(null1)

null11 <- lmer(`Crop biomass carbon` ~ (1 | Crop ) , data = Diff_crop_biomass_c)
AIC(null11)

null2 <- lmer(`Crop biomass carbon` ~ (0 + MWST | Climate) , data = Diff_crop_biomass_c)
AIC(null2)

null21 <- lmer(`Crop biomass carbon` ~ (0 + MWST | Crop) , data = Diff_crop_biomass_c)
AIC(null21)

null3 <- lmer(`Crop biomass carbon` ~ (MWST | Climate), data = Diff_crop_biomass_c)
AIC(null3)

null31 <- lmer(`Crop biomass carbon` ~ (MWST | Crop), data = Diff_crop_biomass_c)
AIC(null31)

null41 <- lmer(`Crop biomass carbon` ~ (1 | Climate ) + (1 | Crop) , data = Diff_crop_biomass_c)
AIC(null41)

null42 <- lmer(`Crop biomass carbon` ~ (0 + MWST | Climate) + (0 + MWST | Crop) , 
               data = Diff_crop_biomass_c)
AIC(null42)

null43 <- lmer(`Crop biomass carbon` ~ (0 + MWST | Climate) + (MWST | Crop), data = Diff_crop_biomass_c)
AIC(null43)

null44 <- lmer(`Crop biomass carbon` ~ ( MWST | Climate) + (0 + MWST | Crop), data = Diff_crop_biomass_c)
AIC(null44)

null45 <- lmer(`Crop biomass carbon` ~ (0 + MWST | Climate) + (1 | Crop), data = Diff_crop_biomass_c)
AIC(null45)

null46 <- lmer(`Crop biomass carbon` ~ ( 1 | Climate) + (0 + MWST | Crop), data = Diff_crop_biomass_c)
AIC(null46)

null47 <- lmer(`Crop biomass carbon` ~ (MWST | Climate) + (1 | Crop), data = Diff_crop_biomass_c)
AIC(null47)

null48 <- lmer(`Crop biomass carbon` ~ (1 | Climate) + (MWST | Crop), data = Diff_crop_biomass_c)
AIC(null48)

null49 <- lmer(`Crop biomass carbon` ~ (MWST | Climate) + (MWST | Crop), data = Diff_crop_biomass_c)
AIC(null49)

aov_straw <- anova(null1, null11, null2, null21, null3, null31, null41, null42, null43,
                   null44, null45, null46, null47, null48, null49)

# build up the linear mix-effects model
# lme4
# total
straw.total <- lmer(`Crop biomass carbon` ~ MWST + (MWST | Climate) + (MWST | Crop), data = Diff_crop_biomass_c %>% 
                      filter(CT == 'Straw'), REML = T, 
                    lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
anova(straw.total, type = "I")
summary(straw.total)


grain.total <- lmer(`Crop biomass carbon` ~ MWST + (MWST | Climate) + (MWST | Crop), 
                    data = Diff_crop_biomass_c %>% 
                      filter(CT == 'Grain'), REML = T, 
                    lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
anova(grain.total, type = "I")
summary(grain.total)

# predict confidence interval
# straw
straw.mwst <- data.frame(MWST = seq(-20, 20, length = 100))
straw.mwst.pred <- predict(straw.total, newdata = straw.mwst, re.form = ~ 0)
straw.mwst.ci <- bootMer(straw.total, FUN = function(.)
  predict(.,newdata = straw.mwst,
          re.form = ~0), nsim = 999, ncpus = 12)
straw.mwst.ci.pred <- apply(straw.mwst.ci$t, 2, function(x) x[order(x)][c(25, 975)])
straw.mwst.lb <- straw.mwst.ci.pred[1,]
straw.mwst.ub <- straw.mwst.ci.pred[2,]

straw.mwst$pred <- straw.mwst.pred
straw.mwst$lci <- straw.mwst.lb
straw.mwst$uci <- straw.mwst.ub

straw.mwst$type <- c("Straw C")
straw.mwst$TT <- "MWST"
straw.mwst$LT <- "Mid-high Latitude"
colnames(straw.mwst)[1] <- "Tem value"


# grain-----
grain.mwst <- data.frame(MWST = seq(-20, 20, length = 100))
grain.mwst.pred <- predict(grain.total, newdata = grain.mwst, re.form = ~ 0)
grain.mwst.ci <- bootMer(grain.total, FUN = function(.)
  predict(.,newdata = grain.mwst,
          re.form = ~0), nsim = 999)
grain.mwst.ci.pred <- apply(grain.mwst.ci$t, 2, function(x) x[order(x)][c(25, 975)])
grain.mwst.lb <- grain.mwst.ci.pred[1,]
grain.mwst.ub <- grain.mwst.ci.pred[2,]

grain.mwst$pred <- grain.mwst.pred
grain.mwst$lci <- grain.mwst.lb
grain.mwst$uci <- grain.mwst.ub

grain.mwst$type <- c("grain C")
grain.mwst$TT <- "MWST"
grain.mwst$LT <- "Mid-high Latitude"
colnames(grain.mwst)[1] <- "Tem value"


# mid-high latitude and low latitude
straw.MWST.h.lme <- lmer(`Crop biomass carbon` ~ MWST + (MWST | Climate) + (MWST | Crop), data = Diff_crop_biomass_c %>% 
                           filter(CT == 'Straw' & LT == "Mid-high latitude"), REML = F, 
                         lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
anova(straw.MWST.h.lme, type = "I")
summary(straw.MWST.h.lme)


straw.MWST.l.lme <- lmer(`Crop biomass carbon` ~ MWST + (MWST | Climate) + (MWST | Crop), data = Diff_crop_biomass_c %>% 
                           filter(CT == 'Straw' & LT == "Low latitude" ),
                         REML = F, lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
anova(straw.MWST.l.lme, type = "I")
summary(straw.MWST.l.lme)


grain.MWST.h.lme <- lmer(`Crop biomass carbon` ~ MWST + (MWST | Climate) + (MWST | Crop) , data = Diff_crop_biomass_c %>% 
                           filter(CT == 'Grain' & LT == "Mid-high latitude"),
                         REML = F, lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb'))) 
anova(grain.MWST.h.lme, type="I")
summary(grain.MWST.h.lme)


grain.MWST.l.lme <- lmer(`Crop biomass carbon` ~ MWST + (MWST | Climate) + (MWST | Crop), data = Diff_crop_biomass_c %>% 
                           filter(CT == 'Grain' & LT == "Low latitude"),
                         REML = F, lmerControl(optimizer ='optimx', optCtrl=list(method='nlminb')))
anova(grain.MWST.l.lme, type="I")
summary(grain.MWST.l.lme)
confint(grain.MWST.l.lme)


# straw-mwst-high & low latitude-----
straw.mwst.h <- data.frame(MWST = seq(-20, 20, length = 100))
straw.mwst.pred.h <- predict(straw.MWST.h.lme, newdata = straw.mwst.h, re.form = ~ 0)
straw.mwst.h.ci <- bootMer(straw.MWST.h.lme, FUN = function(.)
  predict(.,newdata = straw.mwst.h,
          re.form = ~0), nsim = 999, ncpus = getOption("boot.ncpus", 12L))
straw.mwst.h.ci.pred <- apply(straw.mwst.h.ci$t, 2, function(x) x[order(x)][c(25, 975)])
straw.mwst.h.lb <- straw.mwst.h.ci.pred[1,]
straw.mwst.h.ub <- straw.mwst.h.ci.pred[2,]

straw.mwst.h$pred <- straw.mwst.pred.h
straw.mwst.h$lci <- straw.mwst.h.lb
straw.mwst.h$uci <- straw.mwst.h.ub

straw.mwst.h$type <- c("Straw C")
straw.mwst.h$TT <- "MWST"
straw.mwst.h$LT <- "Mid-high Latitude"
colnames(straw.mwst.h)[1] <- "Tem value"


straw.mwst.l <- data.frame(MWST = seq(-20, 20, length = 100))
straw.mwst.pred.l <- predict(straw.MWST.l.lme, newdata = straw.mwst.l, re.form = ~ 0)
straw.mwst.l.ci <- bootMer(straw.MWST.l.lme, FUN = function(.)
  predict(.,newdata = straw.mwst.l,
          re.form = ~0), nsim = 999, ncpus = getOption("boot.ncpus", 12L))
straw.mwst.l.ci.pred <- apply(straw.mwst.l.ci$t, 2, function(x) x[order(x)][c(25, 975)])
straw.mwst.l.lb <- straw.mwst.l.ci.pred[1,]
straw.mwst.l.ub <- straw.mwst.l.ci.pred[2,]

straw.mwst.l$pred <- straw.mwst.pred.l
straw.mwst.l$lci <- straw.mwst.l.lb
straw.mwst.l$uci <- straw.mwst.l.ub

straw.mwst.l$type <- c("Straw C")
straw.mwst.l$TT <- "MWST"
straw.mwst.l$LT <- "Low Latitude"
colnames(straw.mwst.l)[1] <- "Tem value"


# grain-mwst-high & low latitude-----
grain.mwst.h <- data.frame(MWST = seq(-20, 20, length = 100))
grain.mwst.pred.h <- predict(grain.MWST.h.lme, newdata = grain.mwst.h, re.form = ~ 0)
grain.mwst.h.ci <- bootMer(grain.MWST.h.lme, FUN = function(.)
  predict(.,newdata = grain.mwst.h,
          re.form = ~0),nsim = 999)
grain.mwst.h.ci.pred <- apply(grain.mwst.h.ci$t, 2, function(x) x[order(x)][c(25, 975)])
grain.mwst.h.lb <- grain.mwst.h.ci.pred[1,]
grain.mwst.h.ub <- grain.mwst.h.ci.pred[2,]

grain.mwst.h$pred <- grain.mwst.pred.h
grain.mwst.h$lci <- grain.mwst.h.lb
grain.mwst.h$uci <- grain.mwst.h.ub

grain.mwst.h$type <- c("Grain C")
grain.mwst.h$TT <- "MWST"
grain.mwst.h$LT <- "Mid-high Latitude"
colnames(grain.mwst.h)[1] <- "Tem value"



grain.mwst.l <- data.frame(MWST = seq(-20, 20, length = 100))
grain.mwst.pred.l <- predict(grain.MWST.l.lme, newdata = grain.mwst.l, re.form = ~ 0)
grain.mwst.l.ci <- bootMer(grain.MWST.l.lme, FUN = function(.)
  predict(.,newdata = grain.mwst.l,
          re.form = ~0),nsim = 999, parallel = "multicore", ncpus = getOption("boot.ncpus", 10L))
grain.mwst.l.ci.pred <- apply(grain.mwst.l.ci$t, 2, function(x) x[order(x)][c(25, 975)])
grain.mwst.l.lb <- grain.mwst.l.ci.pred[1,]
grain.mwst.l.ub <- grain.mwst.l.ci.pred[2,]

grain.mwst.l$pred <- grain.mwst.pred.l
grain.mwst.l$lci <- grain.mwst.l.lb
grain.mwst.l$uci <- grain.mwst.l.ub

grain.mwst.l$type <- c("Grain C")
grain.mwst.l$TT <- "MWST"
grain.mwst.l$LT <- "Low Latitude"
colnames(grain.mwst.l)[1] <- "Tem value"

# nlme
# straw
total_straw_fd <- Diff_crop_biomass_c[Diff_crop_biomass_c$CT == 'Straw', ]
straw_mwst <- lme(fixed = `Crop biomass carbon` ~ MWST, 
                  data = total_straw_fd,
                  random = list(Climate = ~ MWST, Crop = ~ MWST),
                  na.action = "na.omit",
                  method = "ML", 
                  weights = varIdent(form =  ~ MWST | Climate),
                  control = lmeControl(opt = "optim", 
                                       msMaxIter = 500, 
                                       msTol = 1e-4))

summary(straw_mwst)
anova(straw_mwst)
coef(straw_mwst)
AIC(logLik( straw_mwst))
intervals(straw_mwst, which = "fixed")  

# residuals test
# 1. Residuals vs Fitted
residual <- ggplot(data.frame(Fitted = fitted(straw_mwst), Residuals = residuals(straw_mwst)), aes(x = Fitted, y = Residuals)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "#1874CD", fill = "#1874CD", alpha = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_y_continuous(limits = c(-100, 100))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 14, vjust = 1),
    axis.text.y = element_text(size = 14, hjust = 1),
    axis.title.y = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    panel.border = element_rect(color = "black", size = 1),
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.margin = ggplot2::margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = "Fitted Values", y = "Residuals", title = "Residual vs Fitted Plot");residual

# 2. Normality Q-Q plot
qq <- ggplot(data.frame(Standardized_Residuals = residuals(straw_mwst, type = "normalized")), aes(sample = Standardized_Residuals)) +
  stat_qq(color = "#1874CD") +
  stat_qq_line() +
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 14, vjust = 1),
    axis.text.y = element_text(size = 14, hjust = 1),
    axis.title.y = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    panel.border = element_rect(color = "black", size = 1),
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.margin = ggplot2::margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = "Theoretical Quantiles", y = "Observed Quantiles", title = "Normal Q-Q Plot");qq

# 3. Normalized residual histogram
hist <- ggplot(data.frame(Standardized_Residuals = residuals(straw_mwst, type = "normalized")), aes(x = Standardized_Residuals)) +
  geom_histogram(binwidth = 0.5, fill = "#1874CD", color = "black") +
  scale_y_continuous(limits = c(0, 600))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 14, vjust = 1),
    axis.text.y = element_text(size = 14, hjust = 1),
    axis.title.y = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    panel.border = element_rect(color = "black", size = 1),
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.margin = ggplot2::margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = "Standardized Residuals", y = "Frequency", title = "Histogram of Standardized Residuals");hist

# predict by bootstrap
s_mwst_total_pfd <- data.frame(MWST = seq(min(total_straw_fd$MWST), 
                                          max(total_straw_fd$MWST), 
                                          length.out = 100))


constructCIRibbon <- function(newdata, model) {
  df <- newdata %>%
    mutate(Predict = predict(model, newdata = ., level = 0))
  
  mm <- model.matrix(eval(eval(model$call$fixed)[-2]), data = df)
  
  vars <- mm %*% vcov(model) %*% t(mm)
  sds <- sqrt(diag(vars))
  
  df %>% mutate(
    lowCI = Predict - 1.96 * sds,
    HiCI = Predict + 1.96 * sds
  )
}

s_mwst_total_pfd1 <- constructCIRibbon(s_mwst_total_pfd, straw_mwst)
s_mwst_total_pfd1$type <- c("Straw C")
s_mwst_total_pfd1$TT <- "MWST"


# straw-h-mwst
straw_mwst_h_fd <- Diff_crop_biomass_c %>% 
  filter(CT == 'Straw' & LT == "Mid-high latitude")

straw_mwst_h_model <- lme(fixed = `Crop biomass carbon` ~ MWST, 
                          data = straw_mwst_h_fd,
                          random = list(Climate = ~ MWST, Crop = ~ MWST),
                          na.action = "na.omit",
                          method = "ML", 
                          weights = varIdent(form =  ~ MWST | Climate),
                          control = lmeControl(opt = "optim", 
                                               msMaxIter = 500, 
                                               msTol = 1e-4))
summary(straw_mwst_h_model)
anova(straw_mwst_h_model)
coef(straw_mwst_h_model)
AIC(logLik( straw_mwst_h_model))
intervals(straw_mwst_h_model, which = "fixed")  


s_mwst_h_pfd <- data.frame(MWST = seq(min(straw_mwst_h_fd$MWST), max(straw_mwst_h_fd$MWST), 
                                      length.out = 100))


s_mwst_h_pfd1 <- constructCIRibbon(s_mwst_h_pfd, straw_mwst_h_model)
s_mwst_h_pfd1$type <- c("Straw C")
s_mwst_h_pfd1$TT <- "MWST"
s_mwst_h_pfd1$LT <- "Mid-high Latitude"


# straw_mwst_l
straw_mwst_l_fd <- Diff_crop_biomass_c %>% 
  filter(CT == 'Straw' & LT == "Low latitude")

straw_mwst_l_model <- lme(fixed = `Crop biomass carbon` ~ MWST, 
                          data = straw_mwst_l_fd,
                          random = list(Climate = ~ MWST, Crop = ~ MWST),
                          na.action = "na.omit",
                          method = "ML", 
                          weights = varIdent(form =  ~ MWST | Climate),
                          control = lmeControl(opt = "optim", 
                                               msMaxIter = 500, 
                                               msTol = 1e-4))

summary(straw_mwst_l_model)
anova(straw_mwst_l_model)
AIC(logLik( straw_mwst_l_model))
intervals(straw_mwst_l_model, which = "fixed")


s_mwst_l_pfd <- data.frame(MWST = seq(min(straw_mwst_l_fd$MWST), max(straw_mwst_l_fd$MWST), 
                                      length.out = 100))
s_mwst_l_pfd1 <- constructCIRibbon(s_mwst_l_pfd, straw_mwst_l_model)
s_mwst_l_pfd1$type <- c("Straw C")
s_mwst_l_pfd1$TT <- "MWST"
s_mwst_l_pfd1$LT <- "Low Latitude"


# grain
total_grain_fd <- Diff_crop_biomass_c[Diff_crop_biomass_c$CT == 'Grain', ]
grain_mwst <- lme(fixed = `Crop biomass carbon` ~ MWST, 
                  data = total_grain_fd,
                  random = list(Climate = ~ 1, Crop = ~ 1),
                  na.action = "na.omit",
                  method = "ML", 
                  weights = varIdent(form =  ~ MWST | Climate),
                  control = lmeControl(opt = "optim", 
                                       msMaxIter = 500, 
                                       msTol = 1e-4))
summary(grain_mwst)
anova(grain_mwst)
coef(grain_mwst)
AIC(logLik(grain_mwst))
intervals(grain_mwst, which = "fixed")  


# 1. Residuals vs Fitted
residual_g <- ggplot(data.frame(Fitted = fitted(grain_mwst), Residuals = residuals(grain_mwst)), aes(x = Fitted, y = Residuals)) +
  geom_point(alpha = 0.3) +
  geom_smooth(method = "lm", color = "#1874CD", fill = "#1874CD", alpha = 0.2) +
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_y_continuous(limits = c(-200, 200))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 14, vjust = 1),
    axis.text.y = element_text(size = 14, hjust = 1),
    axis.title.y = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    panel.border = element_rect(color = "black", size = 1),
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.margin = ggplot2::margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = "Fitted Values", y = "Residuals", title = "Residual vs Fitted Plot");residual_g

# 2. Normality Q-Q plot
qq_g <- ggplot(data.frame(Standardized_Residuals = residuals(grain_mwst, type = "normalized")), 
               aes(sample = Standardized_Residuals)) +
  stat_qq(color = "#1874CD") +
  stat_qq_line() +
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 14, vjust = 1),
    axis.text.y = element_text(size = 14, hjust = 1),
    axis.title.y = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    panel.border = element_rect(color = "black", size = 1),
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.margin = ggplot2::margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = "Theoretical Quantiles", y = "Observed Quantiles", title = "Normal Q-Q Plot");qq_g

# 3. Normalized residual histogram
hist_g <- ggplot(data.frame(Standardized_Residuals = residuals(grain_mwst, type = "normalized")), 
                 aes(x = Standardized_Residuals)) +
  geom_histogram(binwidth = 0.5, fill = "#1874CD", color = "black") +
  scale_y_continuous(limits = c(0, 15000))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 14, vjust = 1),
    axis.text.y = element_text(size = 14, hjust = 1),
    axis.title.y = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    panel.border = element_rect(color = "black", size = 1),
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.margin = ggplot2::margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = "Standardized Residuals", y = "Frequency", 
       title = "Histogram of Standardized Residuals");hist_g


g_mwst_total_pfd <- data.frame(MWST = seq(min(total_grain_fd$MWST), 
                                          max(total_grain_fd$MWST), 
                                          length.out = 100))


g_mwst_total_pfd1 <- constructCIRibbon(g_mwst_total_pfd, grain_mwst1)
g_mwst_total_pfd1$type <- c("Straw C")
g_mwst_total_pfd1$TT <- "MWST"


# grain_mwst_h 
grain_mwst_h_fd <- Diff_crop_biomass_c %>% 
  filter(CT == 'Grain' & LT == "Mid-high latitude")
grain_mwst_h_model <- lme(fixed = `Crop biomass carbon` ~ MWST, 
                          data = grain_mwst_h_fd,
                          random = list(Climate = ~ MWST, Crop = ~ MWST),
                          na.action = "na.omit",
                          method = "REML", 
                          weights = varIdent(form =  ~ MWST | Climate),
                          control = lmeControl(opt = "optim", 
                                               msMaxIter = 500, 
                                               msTol = 1e-4))
summary(grain_mwst_h_model)
anova(grain_mwst_h_model)
AIC(logLik( grain_mwst_h_model))

g_mwst_h_pfd <- data.frame(MWST = seq(min(grain_mwst_h_fd$MWST), 
                                      max(grain_mwst_h_fd$MWST), 
                                      length.out = 100))
g_mwst_h_pfd1 <- constructCIRibbon(g_mwst_h_pfd, grain_mwst_h_model)
g_mwst_h_pfd1$type <- c("Grain C")
g_mwst_h_pfd1$TT <- "MWST"
g_mwst_h_pfd1$LT <- "Mid-high Latitude"

# grain_mwst_l
grain_mwst_l_fd <- Diff_crop_biomass_c %>% 
  filter(CT == 'Grain' & LT == "Low latitude")
grain_mwst_l_model <- lme(fixed = `Crop biomass carbon` ~ MWST, 
                          data = grain_mwst_l_fd,
                          random = list(Climate = ~ MWST, Crop = ~ MWST),
                          na.action = "na.omit",
                          method = "REML", 
                          weights = varIdent(form =  ~ MWST | Climate),
                          control = lmeControl(opt = "optim", 
                                               msMaxIter = 500, 
                                               msTol = 1e-4))
summary(grain_mwst_l_model)
anova(grain_mwst_l_model)
AIC(logLik( grain_mwst_l_model))

g_mwst_l_pfd <- data.frame(MWST = seq(min(grain_mwst_l_fd$MWST), 
                                      max(grain_mwst_l_fd$MWST), 
                                      length.out = 100))
g_mwst_l_pfd1 <- constructCIRibbon(g_mwst_l_pfd, grain_mwst_l_model)
g_mwst_l_pfd1$type <- c("Grain C")
g_mwst_l_pfd1$TT <- "MWST"
g_mwst_l_pfd1$LT <- "Low Latitude"

# Fixed effects and random effects adaptation
# plot the results
cal_lme_obs <- function(model, data, formula) {
  # cal the fixed effect
  modmat <- model.matrix(formula, data = data)
  fixed_effects <- modmat %*% fixef(model)

  # set up the initial vector
  random_effects <- rep(0, nrow(data))
  
  # extract random effect
  random_effects_climate <- ranef(model)$Climate
  random_effects_crop <- ranef(model)$Crop
  
  for(i in 1:nrow(data)) {
    climate_name <- data$Climate[i]
    climate_effects <- random_effects_climate[climate_name, ]
    crop_name <- data$Crop[i]
    crop_effects <- random_effects_crop[paste(climate_name, crop_name, sep = "/"), ]
    mwst_value <- data$MWST[i] 
    random_slope <- as.numeric(climate_effects['MWST'] * mwst_value)
    random_effects[i] <- as.numeric(climate_effects['(Intercept)']) + 
      as.numeric(crop_effects['(Intercept)']) + random_slope
  }
  
  data$renew_obs <- fixed_effects + random_effects
  return(data)
}

straw_mwst_total_1 <- cal_lme_obs(model = straw_mwst1, 
                                  data = total_straw_fd, 
                                  formula = ~ MWST)

straw_total <- ggplot()+
  geom_ribbon(data = s_mwst_total_pfd1,
              aes(x = MWST, ymin = lowCI, ymax = HiCI), 
              alpha = 0.1, 
              fill = "#C2A179")+
  geom_point(data = straw_mwst_total_1 %>% filter(CT == "StrawC"),
             aes(x = MWST, y = renew_obs), 
             alpha = 0.5,
             shape = 21,  
             size = 5, 
             fill = "#C2A179", 
             stroke = 0.2)+
  geom_line(data = s_mwst_total_pfd1,
            aes(x = MWST, y = Predict), 
            size = 1.5, 
            color = "#8C6A40")+
  scale_y_continuous(limits = c(-100, 190), breaks = seq(-90, 180, 90))+
  scale_x_continuous(expand = c(0, 1))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 14, vjust = 1),
    axis.text.y = element_text(size = 14, hjust = 1),
    axis.title.y = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.margin = margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = expression(paste(Delta, "Winter soil temperature (", degree, "C)")), 
       y = expression(paste(Delta, "Straw C (g kg"^-1,")")),
       title = "Global") ;straw_total

# distribution of crop biomass carbon
data <- read.xlsx("data.xlsx")
str(data)
dim(data)
# mid-high straw  
result.s.md <- t.test(data[data$Latitude.type == "Mid-high latitude" & 
                               data$Crop.Component == "Straw", ]$Crop.biomass.carbon)
result.s.md
ci.s.md <- result.s.md$conf.int; ci.s.md

# mid-high grain
result.g.md <- t.test(data[data$Latitude.type == "Mid-high latitude" & 
                               data$Crop.Component == "Grain", ]$Crop.biomass.carbon)
result.g.md
ci.g.md <- result.g.md$conf.int; ci.g.md

# low straw  
result.s.l <- t.test(data[data$Latitude.type == "Low latitude" & 
                              data$Crop.Component == "Straw", ]$Crop.biomass.carbon)
result.s.l
ci.s.l <- result.s.l$conf.int;ci.s.l

# low grain
result.g.l <- t.test(data[data$Latitude.type == "Low latitude" & 
                              data$Crop.Component == "Grain", ]$Crop.biomass.carbon)
result.g.l
ci.g.l <- result.g.l$conf.int;ci.g.l


total.meta.n.m <- data %>% 
  dplyr::group_by (Latitude.type, Crop.Component) %>%  
  dplyr::summarise(mean = mean(Crop.biomass.carbon));total.meta.n.m

total.meta.n.m$Latitude.type <- factor(total.meta.n.m$Latitude.type, 
                                       levels = c("Mid-high latitude", "Low latitude"))

total.meta.n.m.total <- data %>% 
  dplyr::group_by (Crop.Component) %>%  
  dplyr::summarise(mean = mean(Crop.biomass.carbon));total.meta.n.m.total

strawc_total <- 
  ggplot(data %>% filter(Crop.Component == "Straw") %>% 
           mutate(type = "Global"), aes(x = Crop.Component, y = Crop.biomass.carbon))+
  stat_boxplot(geom = "errorbar", width = 0.2)+
  geom_boxplot(notch = F, notchwidth = 0.7, 
               width = 0.35, alpha = 0.9, 
               outlier.colour = NA)+
  geom_point(data = total.meta.n.m.total %>% 
               filter(Crop.Component == "Straw") %>% 
               mutate(type = "Global"),
             mapping = aes(x = Crop.Component, y = mean), 
             fill = "#C2A179",
             size = 3, 
             shape = 21)+
  scale_color_manual(values = c("#C2A179", "#ECA8A9"))+
  scale_fill_manual(values = c("#C2A179", "#ECA8A9"))+
  scale_y_continuous(limits = c(290, 520), breaks = seq(300, 500, 100))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_text(size = 12, vjust = 0.5),
    axis.ticks.x = element_blank(),
    axis.title.y = element_text(size = 12, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 12, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.background = element_blank(),
    plot.margin = margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = "", y = expression(paste("Straw C (g kg "^-1,")")));strawc_total

straw_linear_box <- straw_total1 + annotation_custom(ggplotGrob(strawc_total),
                                                     xmin = 2, xmax = 10.8,
                                                     ymin = 0, ymax = 180)
straw_linear_box



grain_mwst_total_1 <- cal_lme_obs(model = grain_mwst1, 
                                  data = total_grain_fd, 
                                  formula = ~ MWST)

nrow(grain_mwst_total_1)

grain_total1 <- ggplot()+
  geom_point(data = grain_mwst_total_1,
             aes(x = MWST, y = renew_obs ),
             shape = 21,
             size = 5,
             color = "#C2A179",
             alpha = 0.3,
             stroke = 0.55)+
  geom_line(data = g_mwst_total_pfd1 ,
            aes(x = MWST, y = Predict),
            size = 1.5,
            color = "#8C6A40")+
  scale_y_continuous(limits = c(-300, 380), breaks = seq(-300, 300, 150))+
  scale_x_continuous(expand = c(0, 1.5))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 14, vjust = 1),
    axis.text.y = element_text(size = 14, hjust = 1),
    axis.title.y = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 15, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", size = 15, hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.margin = margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = expression(paste(Delta, "Winter soil temperature (", degree, "C)")), 
       y = expression(paste(Delta, "Grain C (g kg"^-1,")")),
       title = "Global") ;grain_total1


grainc_total <- 
  ggplot(metadf %>% filter(Crop.Component == "Grain") %>% 
           mutate(type = "Global"), aes(x = type, y = Crop.biomass.carbon))+
  stat_boxplot(geom = "errorbar", width = 0.2)+
  geom_boxplot(notch = F, notchwidth = 0.7, 
               width = 0.35, alpha = 0.9, 
               outlier.colour = NA)+
  geom_point(data = total.meta.n.m.total %>% 
               filter(Crop.Component == "Grain") %>% mutate(type = "Global"),
             mapping = aes(x = type, y = mean),
             color = "#8C6A40", 
             size = 3, 
             shape = 21, 
             stroke = 0.55)+
  scale_color_manual(values = c("#C2A179", "#ECA8A9"))+
  scale_fill_manual(values = c("#C2A179", "#ECA8A9"))+
  scale_y_continuous(limits = c(290, 525), breaks = seq(300, 500, 100))+
  theme_bw()+
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_blank(),
    axis.text.y = element_text(size = 12, vjust = 0.5),
    axis.ticks.x = element_blank(),
    axis.title.y = element_text(size = 12, vjust = 0.5, hjust = 0.5),
    axis.title.x = element_text(size = 12, vjust = 0.5, hjust = 0.5),
    plot.title = element_text(face = "bold", hjust = 0.5, vjust = 0.5),
    legend.position = "none",
    axis.ticks.length.y = unit(0.1, 'cm'),
    axis.ticks.length.x = unit(0.1, 'cm'),
    plot.background = element_blank(),
    plot.margin = margin(t = 0, r = 0.5, b = 1, l = 0.5, unit = "lines")
  )+
  labs(x = "", y = expression(paste("Grain C (g kg "^-1,")")));grainc_total


grain_linear_box <- grain_total1 + annotation_custom(ggplotGrob(grainc_total),
                                                     xmin = 4.5, xmax = 18,
                                                     ymin = -50, ymax = 370)
grain_linear_box


# random forest
data <- read.xlsx("data.xlsx", 1)
Strawdata <- data[data$Crop.Component == "Straw", ]
str(Strawdata)

Graindata <- data[data$Crop.Component == "Grain", ]
str(Graindata)

# Straw
set.seed(111)
Strawdata_rf <- randomForest(Crop.biomass.carbon ~ ., 
                             data = Strawdata[, c(-1:-2, -4, -5, -7:-10, -16, -23, -38)], 
                             importance = TRUE, proximity = TRUE, na.action = na.omit)
Strawdata_rf


set.seed(111)
Strawdata_rf_im_var <- rfPermute(Crop.biomass.carbon ~ ., 
                                 data = Strawdata[, c(-1:-2, -4, -5, -7:-10, -16, -23, -38)], 
                                 ntree = 500, na.action = na.omit, nrep = 100, num.cores = 14)
Strawdata_rf_im_var1 <- rfPermute::importance(Strawdata_rf_im_var, sort.by = NULL, decreasing = TRUE)
Strawdata_rf_im_var1 <- as.data.frame(Strawdata_rf_im_var1)
Strawdata_rf_im_var1$name <- rownames(Strawdata_rf_im_var1)
Strawdata_rf_im_var1
nrow(Strawdata_rf_im_var1)

Strawdata_rf_im_var1 %>% arrange(desc(IncNodePurity)) %>% 
  select(name, IncNodePurity, IncNodePurity.pval, everything())

Straw_rf <- rbind(as.data.frame(Strawdata_rf_im_var1[,c("IncNodePurity","IncNodePurity.pval")] %>%
                                  as_tibble(rownames = "names") %>%
                                  mutate(label = if_else(`IncNodePurity.pval`<0.001,"***",
                                                         if_else(`IncNodePurity.pval`<0.01,"**",
                                                                 if_else(`IncNodePurity.pval`<0.05,"*","ns")))) %>%
                                  arrange(`IncNodePurity`) %>%
                                  mutate(group = if_else(label=="ns","In_sig","Sig"),
                                         names = forcats::fct_inorder(names),
                                         crop = "Grain")))
Straw_rf$rel = round(Straw_rf$IncNodePurity / max(Straw_rf$IncNodePurity) * 100, 2)

# Straw #### 
straw.rf <- ggplot(Straw_rf, aes(
  x = names,            
  y = `rel`,  
))+
  geom_bar(stat = 'identity', position = position_dodge(0.9), color = "black", fill = "#EF853D")+
  geom_text(aes(x = names, y = `rel` + 0.5, label = label,
              hjust = -0.5), size = 5,
          position = position_dodge(0.9))+
  scale_y_continuous(limits = c(0, 105), expand = c(0, 0), breaks = seq(0, 100, 20))+
  theme_bw()+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(color = "black", size = 12, vjust = 0.5, hjust = 0.5),
        axis.text.y = element_text(color = "black", size = 12),
        legend.position = "none",
        axis.title = element_text(color = "black", size = 12),
        axis.ticks.length.y = unit(0.15, 'cm'),
        axis.ticks.length.x = unit(0.15, 'cm'),
        strip.text = element_blank(),
  )+
  coord_flip()+
  labs(x = "", y = "Relative importance (%)",
       title = "Straw purity");straw.rf


# grain
set.seed(111)
Graindata_rf <- randomForest(Crop.biomass.carbon ~ ., 
                             data = Graindata[, c(-1:-2, -4, -5, -7:-10, -16, -23, -38)], 
                             importance = TRUE, proximity = TRUE, na.action = na.omit)
Graindata_rf

set.seed(111)
Graindata_rf_im_var <- rfPermute(Crop.biomass.carbon ~ ., 
                                 data = Graindata[, c(-1:-2, -4, -5, -7:-10, -16, -23, -38)], 
                                 ntree = 500, na.action = na.omit, nrep = 100, num.cores = 14)
Graindata_rf_im_var1 <- rfPermute::importance(Graindata_rf_im_var, sort.by = NULL, decreasing = TRUE)
Graindata_rf_im_var1 <- as.data.frame(Graindata_rf_im_var1)
Graindata_rf_im_var1$name <- rownames(Graindata_rf_im_var1)
Graindata_rf_im_var1
nrow(Graindata_rf_im_var1)

Graindata_rf_im_var1 %>% arrange(desc(IncNodePurity)) %>% 
  select(name, IncNodePurity, IncNodePurity.pval, everything())

Grain_rf <- rbind(as.data.frame(Graindata_rf_im_var1[,c("IncNodePurity","IncNodePurity.pval")] %>%
                                  as_tibble(rownames = "names") %>%
                                  mutate(label = if_else(`IncNodePurity.pval`<0.001,"***",
                                                         if_else(`IncNodePurity.pval`<0.01,"**",
                                                                 if_else(`IncNodePurity.pval`<0.05,"*","ns")))) %>%
                                  arrange(`IncNodePurity`) %>%
                                  mutate(group = if_else(label=="ns","In_sig","Sig"),
                                         names = forcats::fct_inorder(names),
                                         crop = "Grain")))
Grain_rf$rel = round(Grain_rf$IncNodePurity / max(Grain_rf$IncNodePurity) * 100, 2)

# Straw #### 
grain.rf <- ggplot(Grain_rf, aes(
  x = names,      
  y = `rel`,
))+
  geom_bar(stat = 'identity', position = position_dodge(0.9), color = "black", fill = "#EF853D")+                       
  geom_text(aes(x = names, y = `rel` + 0.5, label = label,
              hjust = -0.5), size = 5,
          position = position_dodge(0.9))+
  scale_y_continuous(limits = c(0, 105), expand = c(0, 0), breaks = seq(0, 100, 20))+
  theme_bw()+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_text(color = "black", size = 12, vjust = 0.5, hjust = 0.5),
        axis.text.y = element_text(color = "black", size = 12),
        legend.position = "none",
        axis.title = element_text(color = "black", size = 12),
        axis.ticks.length.y = unit(0.15, 'cm'),
        axis.ticks.length.x = unit(0.15, 'cm'),
        strip.text = element_blank(),
  )+
  coord_flip()+
  labs(x = "", y = "Relative importance (%)",
       title = "Grain purity");grain.rf


#### controlling for the MWST variable to investgate ####
#### the relationship between SOM and crop biomass carbon ####
straw_mwst_lme <- lme(CropC ~ MWST,
                      data = straw,
                      random = list(Climate = ~ MWST, Crop = ~ MWST),
                      na.action = "na.omit",
                      method = "ML",
                      weights = varIdent(form =  ~ MWST | Climate),
                      control = lmeControl(opt = "optim",
                                           msMaxIter = 500, 
                                           msTol = 1e-4))
# resid(straw_mwst_lme)
summary(straw_mwst_lme)

straw$straw_mwst_resid = residuals(straw_mwst_lme)

som_mwst_lme = lme(straw_mwst_resid ~ SOM,
                   data = straw,
                   random = list(Climate = ~ MWST, Crop = ~ MWST),
                   na.action = "na.omit",
                   method = "ML",
                   weights = varIdent(form =  ~ MWST | Climate),
                   control = lmeControl(opt = "optim",
                                        msMaxIter = 500, 
                                        msTol = 1e-4))
summary(som_mwst_lme)


grain_mwst_lme <- lme(CropC ~ MWST,
                      data = grain,
                      random = list(Climate = ~ MWST, Crop = ~ MWST),
                      na.action = "na.omit",
                      method = "ML",
                      weights = varIdent(form =  ~ MWST | Climate),
                      control = lmeControl(opt = "optim",
                                           msMaxIter = 500, 
                                           msTol = 1e-4))
summary(grain_mwst_lme)

grain$grain_mwst_resid = residuals(grain_mwst_lme)

som_mwst_lme = lme(grain_mwst_resid ~ SOM,
                   data = grain,
                   random = list(Climate = ~ MWST, Crop = ~ MWST),
                   na.action = "na.omit",
                   method = "ML",
                   weights = varIdent(form =  ~ MWST | Climate),
                   control = lmeControl(opt = "optim",
                                        msMaxIter = 500, 
                                        msTol = 1e-4))
summary(som_mwst_lme)



####Original data####
Supdata = read.xlsx("data.xlsx", 1)

####grain####
grain = Supdata |> subset(Crop.Component == "Grain")

grain_lm = lm(Crop.biomass.carbon ~ MWST, grain)
grain_climate = lmer(Crop.biomass.carbon ~ MWST + (1 | Climate.type), data = grain)
grain_crop = lmer(Crop.biomass.carbon ~ MWST + (1 | Crop.type), data = grain)
grain_total = lmer(Crop.biomass.carbon ~ MWST + (1 | Crop.type) + (1 | Climate.type), 
                   data = grain)

AIC(grain_lm);AIC(grain_climate);AIC(grain_crop);AIC(grain_total)
BIC(grain_lm);BIC(grain_climate);BIC(grain_crop);BIC(grain_total)

grain_mwst_lme = lmer(Crop.biomass.carbon ~ MWST + (1 | Crop.type) + (1 | Climate.type), 
                      data = grain)
summary(grain_mwst_lme)

grain$grain_mwst_resid = residuals(grain_mwst_lme)
grain_som_mwstresid = lmer(grain_mwst_resid ~ SOM + (1 | Crop.type) + (1 | Climate.type), 
                           data = grain)
summary(grain_som_mwstresid)

####grain partial plot####
grain_mwst_som_lme = lmer(Crop.biomass.carbon ~ MWST + SOM + 
                            (1 | Crop.type) + (1 | Climate.type), 
                          data = grain)
visreg(grain_mwst_som_lme, "SOM", partial = TRUE, scale = "response")
plot(effect("SOM", grain_mwst_som_lme))

####straw####
straw = Supdata |> subset(Crop.Component == "Straw")

straw_lm = lm(Crop.biomass.carbon ~ MWST, straw)
straw_climate = lmer(Crop.biomass.carbon ~ MWST + (1 | Climate.type), data = straw)
straw_crop = lmer(Crop.biomass.carbon ~ MWST + (1 | Crop.type), data = straw)
straw_total = lmer(Crop.biomass.carbon ~ MWST + (1 | Crop.type) + (1 | Climate.type), 
                   data = straw)

AIC(straw_lm);AIC(straw_climate);AIC(straw_crop);AIC(straw_total)
BIC(straw_lm);BIC(straw_climate);BIC(straw_crop);BIC(straw_total)

straw_mwst_lme = lmer(Crop.biomass.carbon ~ MWST + (1 | Crop.type) + (1 | Climate.type), 
                      data = straw)
summary(straw_mwst_lme)

straw$straw_mwst_resid = residuals(straw_mwst_lme)
straw_som_mwstresid = lmer(straw_mwst_resid ~ SOM + (1 | Crop.type) + (1 | Climate.type), 
                           data = straw)
summary(straw_som_mwstresid)

####staw partial plot####
straw_mwst_som_lme = lmer(Crop.biomass.carbon ~ MWST + SOM + 
                            (1 | Crop.type) + (1 | Climate.type), 
                          data = straw)
visreg(straw_mwst_som_lme, "SOM", partial = TRUE, scale = "response")
plot(effect("SOM", straw_mwst_som_lme))