LME
library(readxl)
library(lme4)
library(nlme)
data0 <- read_excel("E:/SR/Mann.xlsx")
a<-lme(N_effect ~ (PWM)^2+ PWM, random = ~ 1 | sample, na.action=na.omit, data = data0)
summary(a)

b<-lme(N_effect ~ AP, random = ~ 1 | sample, na.action=na.omit, data = data0)
summary(b)
data0 <- read_excel("E:/SR/Mann.xlsx")
c<-lme(N_effect ~ (MTCQ)^2+ MTCQ, random = ~ 1 | sample, na.action=na.omit, data = data0)
summary(c)

d<-lme(N_effect ~ PWM, random = ~ 1 | sample, na.action=na.omit, data = data0)
summary(d)

c<-lm(N_effect ~ (MinTCM)^2+ MinTCM,  na.action=na.omit, data = data0)
summary(c)


piecewiseSEM

library(nlme)
library(lme4)
library(piecewiseSEM)
library(QuantPsyc)

Rt<-read_excel("E:/Q10/data/Rt.xlsx")


#composite predictor of climate#
model1 <- lm(N_effect ~  legacy17 + legacy19 + current18, Rt)
coefs(model1, standardize = 'scale')
beta_legacy17 <-  summary(model1)$coefficients[2, 1]
beta_legacy19 <- summary(model1)$coefficients[3, 1]
beta_current18 <-  summary(model1)$coefficients[4, 1]
climate <- beta_legacy17 * Rt$legacy17 + beta_legacy19 * Rt$legacy19 + beta_current18 * Rt$current18
Rt$climate <- climate
summary(lm(N_effect ~ climate, Rt))
coefs(lm(N_effect ~ climate, Rt))
#composite predictor of soil#
model2 <- lm(N_effect ~ Cropland + Forest + wetland + Grassland, Rt)
summary(model2)$coefficients
beta_Cropland <- summary(model2)$coefficients[2, 1]
beta_Forest <-  summary(model2)$coefficients[3, 1]
beta_wetland <- summary(model2)$coefficients[4, 1]
ecosystem <- beta_Cropland * Rt$Cropland + beta_Forest * Rt$Forest +  beta_wetland * Rt$wetland
Rt$ecosystem <- ecosystem
summary(lm(N_effect ~ ecosystem, Rt))
coefs(lm(N_effect ~ ecosystem, Rt))

#composite predictor of microbe#
model3 <- lm(N_effect ~ SOC + TN + C_N, Rt)
coefs(model3, standardize = 'scale')
beta_SOC <- summary(model3)$coefficients[2, 1]
beta_TN <-  summary(model3)$coefficients[3, 1]
beta_C_N <- summary(model3)$coefficients[4, 1]
soil <- beta_SOC * Rt$SOC + beta_TN * Rt$TN + beta_C_N * Rt$C_N
Rt$soil <- soil
summary(lm(N_effect ~ soil, Rt))
coefs(lm(N_effect ~ soil, Rt)) ##### stander effect ######
model0 <- lm(N_effect ~ Lon + Lat, Rt)
coefs(model0, standardize = 'scale')
beta_Lon <-  summary(model0)$coefficients[2, 1]
beta_Lat <- summary(model0)$coefficients[3, 1]
geography <- beta_Lon * Rt$Lon + beta_Lat * Rt$Lat
Rt$geography <- geography
summary(lm(N_effect ~ geography, Rt))
coefs(lm(N_effect ~ geography, Rt))


#multiple regression#
microbe.list <- list(
  lme(climate ~ geography , random = ~ 1 | site , na.action = na.omit,
      data = Rt),
  lme(soil ~ climate + Rate+ ecosystem, random = ~ 1 | site , na.action = na.omit,
      data = Rt),
  lme(N_effect ~   climate + ecosystem + soil, random = ~ 1 | site , na.action = na.omit,
      data = Rt)
)
microbe.psem <- as.psem(microbe.list)
(new.summary <- summary(microbe.psem, .progressBar = F))
plot(microbe.psem,return = FALSE,alpha = 0.05,show = "std")


Rh<-read_excel("E:/Q10/data/Rh.xlsx")
model1 <- lm(N_effect ~  legacy2 + legacy15 + legacy7, Rh)
coefs(model1, standardize = 'scale')
beta_legacy2 <-  summary(model1)$coefficients[2, 1]
beta_legacy15 <- summary(model1)$coefficients[3, 1]
beta_legacy7 <-  summary(model1)$coefficients[4, 1]
climate <- beta_legacy2 * Rh$legacy2 + beta_legacy15 * Rh$legacy15 + beta_legacy7 * Rh$legacy7
Rh$climate <- climate
summary(lm(N_effect ~ climate, Rh))
coefs(lm(N_effect ~ climate, Rh))

model2 <- lm(N_effect ~ Cropland + Forest + Grassland, Rh)
summary(model2)$coefficients
beta_Cropland <- summary(model2)$coefficients[2, 1]
beta_Forest <-  summary(model2)$coefficients[3, 1]
beta_Grassland <- summary(model2)$coefficients[4, 1]
ecosystem <- beta_Cropland * Rh$Cropland + beta_Forest * Rh$Forest 
Rh$ecosystem <- ecosystem
summary(lm(N_effect ~ ecosystem, Rh))
coefs(lm(N_effect ~ ecosystem, Rh))

model3 <- lm(N_effect ~ SOC + TN + C_N, Rh)
summary(model3)$coefficients
beta_SOC <- summary(model3)$coefficients[2, 1]
beta_TN <-  summary(model3)$coefficients[3, 1]
beta_C_N <- summary(model3)$coefficients[4, 1]
soil <- beta_SOC * Rh$SOC + beta_TN * Rh$TN + beta_C_N * Rh$C_N
Rh$soil <- soil
summary(lm(N_effect ~ soil, Rh))
coefs(lm(N_effect ~ soil, Rh))

model0 <- lm(N_effect ~ Lon + Lat, Rh)
coefs(model0, standardize = 'scale')
beta_Lon <-  summary(model0)$coefficients[2, 1]
beta_Lat <- summary(model0)$coefficients[3, 1]
geography <- beta_Lon * Rh$Lon + beta_Lat * Rh$Lat
Rh$geography <- geography
summary(lm(N_effect ~ geography, Rh))
coefs(lm(N_effect ~ geography, Rh))

microbe.list <- list(
  lme(climate ~ geography , random = ~ 1 | site , na.action = na.omit,
      data = Rh),
  lme(soil ~ climate + Rate+ ecosystem, random = ~ 1 | site , na.action = na.omit,
      data = Rh),
  lme(N_effect ~   climate + ecosystem + soil, random = ~ 1 | site , na.action = na.omit,
      data = Rh)
)
microbe.psem <- as.psem(microbe.list)
(new.summary <- summary(microbe.psem, .progressBar = F))
plot(microbe.psem,return = FALSE,alpha = 0.05,show = "std")



Ra<-read_excel("E:/Q10/data/Ra.xlsx")
model1 <- lm(N_effect ~  current18 + legacy12, Ra)
coefs(model1, standardize = 'scale')
beta_current18 <-  summary(model1)$coefficients[2, 1]
beta_legacy12 <- summary(model1)$coefficients[3, 1]
climate <- beta_current18 * Ra$current18 + beta_legacy12 * Ra$legacy12
Ra$climate <- climate
summary(lm(N_effect ~ climate, Ra))
coefs(lm(N_effect ~ climate, Ra))

model2 <- lm(N_effect ~ Cropland + Forest +  Grass, Ra)
summary(model2)$coefficients
beta_Cropland <- summary(model2)$coefficients[2, 1]
beta_Forest <-  summary(model2)$coefficients[3, 1]
beta_Grassland <- summary(model2)$coefficients[4, 1]
ecosystem <- beta_Cropland * Ra$Cropland + beta_Forest * Ra$Forest +  beta_Grassland * Ra$Grass
Ra$ecosystem <- ecosystem
summary(lm(N_effect ~ ecosystem, Ra))
coefs(lm(N_effect ~ ecosystem, Ra))

model3 <- lm(N_effect ~ SOC + TN + C_N, Ra)
coefs(model3, standardize = 'scale')
beta_SOC <- summary(model3)$coefficients[2, 1]
beta_TN <-  summary(model3)$coefficients[3, 1]
beta_C_N <- summary(model3)$coefficients[4, 1]
soil <- beta_SOC * Ra$SOC + beta_TN * Ra$TN + beta_C_N * Ra$C_N
Ra$soil <- soil
summary(lm(N_effect ~ soil, Ra))
coefs(lm(N_effect ~ soil, Ra))

model0 <- lm(N_effect ~ Lon + Lat, Ra)
coefs(model0, standardize = 'scale')
beta_Lon <-  summary(model0)$coefficients[2, 1]
beta_Lat <- summary(model0)$coefficients[3, 1]
geography <- beta_Lon * Ra$Lon + beta_Lat * Ra$Lat
Ra$geography <- geography
summary(lm(N_effect ~ geography, Ra))
coefs(lm(N_effect ~ geography, Ra))

microbe.list <- list(
  lme(climate ~ geography , random = ~ 1 | site , na.action = na.omit,
      data = Ra),
  lme(soil ~ climate + Rate+ ecosystem, random = ~ 1 | site , na.action = na.omit,
      data = Ra),
  lme(N_effect ~  climate + ecosystem + soil, random = ~ 1 | site , na.action = na.omit,
      data = Ra)
)
microbe.psem <- as.psem(microbe.list)
(new.summary <- summary(microbe.psem, .progressBar = F))
plot(microbe.psem,return = FALSE,alpha = 0.05,show = "std")


VP

library(readxl)
SR<-read_excel("E:/SR/SR2.xlsx")
Envir<-read_excel("E:/SR/Env2.xlsx")


library(vegan)
rda.result<-rda(SR,Envir)
rda.result
summary(rda.result)
vif.cca(rda.result)
rda.vp<-varpart(SR,Envir[c("Rate","Type","Fertilizer_Type")],Envir[c("legacy1", "legacy2","legacy3","legacy4","legacy5","legacy6","legacy7","legacy8","legacy9","legacy10","legacy11","legacy12","legacy13","legacy14","legacy15","legacy16","legacy17","legacy18","legacy19")],
                Envir[c("current1","current2","current3","current4","current5","current6","current7","current8","current9","current10","current11","current12","current13","current14","current15","current16","current17","current18","current19")])
rda.vp
plot(rda.vp, digits = 2, Xnames = c('N_addition', 'Contemporary_Climate',"Paleoclimate"), bg = c('blue','red',"green"))