library(readxl)
library(metafor)
library(ggplot2)
library(glmulti)
production<-read_excel("H:/ZYY/论文/meta分析_氮沉降/投稿/Plant and Soil/resubmit/Appendices.xlsx")
View(production)

la50<-subset(production,abs(Latitude)<=50)
View(la50)
effect_prod50<-escalc(measure ="ROM",data=la50,m1i=Prod.mean.E,sd1i=Prod.sd.E,n1i=Prod.n.E,m2i=Prod.mean.C,sd2i=Prod.sd.C,n2i=Prod.n.C)
library(dplyr)
effect_prod50<- select(effect_prod50,-Prod.mean.E,-Prod.sd.E,-Prod.n.E,-Prod.mean.C,-Prod.sd.C,-Prod.n.C,-ln.prod.RR,-Var.Prod.RR)
View(effect_prod50)

random_prod50<-rma(yi,vi, data=effect_prod50, method="REML")
summary(random_prod50)
attach(effect_prod50)
effect_prod501<- effect_prod50[order(-vi),] 
attach(effect_prod501)
forest(yi, vi,annotate=FALSE, col="grey", slab=NA, xlab="Log response ratio", pch=20,efac=0,font=2,cex=1,psize=1,xlim=c(-2,2), lty=c("solid","blank"))
abline(v=coef(random_prod50), lty="dashed",col="black")  
abline(v=0, lty="dotted",col="red") 
title("Forest of effect size for production(latitude<50°)",cex.main=1.2,font.main=2)



r_Nload50<-rma(yi,vi, mods=~Nload,data=effect_prod50, method="REML")
summary(r_Nload50)
pred_Nload50 <- predict(r_Nload50, newmods=c(0:180))
wi<- 1/sqrt(effect_prod50$vi)
size<-wi*0.1+0.8
plot(effect_prod50$Nload, effect_prod50$yi, pch=19, cex=size, col="grey",ylab="Change in production (ln resposne ratio)",
     xlab="N application(kg N ha-2 yr-1)",
     las=1, bty="l", font=1,cex.lab=1.2)
lines(0:180, pred_Nload50$pred)
lines(0:180, pred_Nload50$ci.lb, lty="dashed")
lines(0:180,pred_Nload50$ci.ub, lty="dashed")
abline(h=0, lty="dotted",col="red",lwd=2)
abline(v=15.75, lty="dotted",col="blue",lwd=2)
text(c(140), -2.5, c("Qm=6.9678  p= 0.0083"),cex=1,font=1,col="red")
text(c(140), -3, c("y=-0.0059x+0.0929"),cex=1,font=1,col="red")
text(c(28), -3.2, c("15.75"),cex=1,font=1,col="blue")
title(main="Critical of N concentration(latitude<50°)", line=1.5, cex.lab=1.5)


r_totN50<-rma(yi,vi, mods=~tot_N,data=effect_prod50, method="REML")
summary(r_totN50)
pred_totN50 <- predict(r_totN50, newmods=c(0:190))
wi<- 1/sqrt(effect_prod50$vi)
size<-wi*0.1+0.8
plot(effect_prod50$tot_N, effect_prod50$yi, pch=19, cex=size, col="grey",ylab="Change in production (ln resposne ratio)",
     xlab="Total N (kg N ha-2 yr-1)",
     las=1, bty="l", font=1,cex.lab=1.2)
lines(0:190, pred_totN50$pred)
lines(0:190, pred_totN50$ci.lb, lty="dashed")
lines(0:190, pred_totN50$ci.ub, lty="dashed")
abline(h=0, lty="dotted",col="red",lwd=2)
abline(v=19.74, lty="dotted",col="blue",lwd=2)
text(c(140), -2.5, c("Qm=8.2251  p= 0.0041"),cex=1,font=1,col="red")
text(c(140), -3, c("y=-0.008x+0.1579"),cex=1,font=1,col="red")
text(c(33), -3.2, c("19.74"),cex=1,font=1,col="blue")
title(main="Critical of N concentration(latitude<50°)", line=1.5, cex.lab=1.5)


r_tempmean50<-rma(yi,vi, mods=~tempmean,data=effect_prod50, method="REML")
summary(r_tempmean50)
r_Julytem50<-rma(yi,vi, mods=~Julytem,data=effect_prod50, method="REML")
summary(r_Julytem50)
r_prectot50<-rma(yi,vi, mods=~prectot,data=effect_prod50, method="REML")
summary(r_prectot50)
rprod_micro50<-rma(yi,vi, mods=~Hummock,data=effect_prod50, method="REML")
summary(rprod_micro50)
rprod_hummock50<-rma(yi,vi, data=subset(effect_prod50,Hummock=="1"), method="REML")
summary(rprod_hummock50)
rprod_hollow50<-rma(yi,vi, data=subset(effect_prod50,Hummock=="0"), method="REML")
summary(rprod_hollow50)


r_species50<-rma(yi,vi, mods=~Species,data=effect_prod50, method="REML")
summary(r_species50)
r_species501<-rma(yi,vi, mods=~Species-1,data=effect_prod50, method="REML")
summary(r_species501)
coef(summary(r_species501))
estimates_species501 <- coef(summary(r_species501))[,1]
se_species501 <- coef(summary(r_species501))[,2]
variances_species501<- se_species501^2
labels_species501 <- c( "angustifolium", "capillifolium", "fallax", "fuscum", "magellanicum",
                        "palustre","papillosum","papillosum*magellanicum")
forest(estimates_species501, variances_species501, slab=labels_species501, font=1,cex=1,psize=1,
       xlim=c(-1.2,0.6),
       xlab="Effect size (log response ratio)")
text(c(-1.4,0.85), 9.5, c("Species", "Mean [95% CI]"),cex=1,font=1)
anova(r_species501, L=c(0,1,0,0,-1,0,0,0))
anova(r_species501, L=c(0,0,1,0,0,0,0,-1))



total50<-rma(yi,vi,mods=~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory,data=effect_prod50, method="REML")
summary(total50)
h_total50<-rma.mv(yi,vi,mods = ~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory,data=effect_prod50,random =list(~1|IDnr,~1|ID),method = "ML")
summary(h_total50)
rma.glmulti <- function(formula, data, ...) { rma(formula, vi, data=data, method="ML", ...)} 
modelselection50<- glmulti(yi ~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory, data=effect_prod50, level=1, fitfunction=rma.glmulti, crit="aicc", confsetsize=1000)
plot(modelselection50)


modelweights50<- weightable(modelselection50)
modelweights50
summary(modelselection50@objects[[1]])
bestmodelweights50<- modelweights50[modelweights50$aicc <= min(modelweights50$aicc) + 2,]


bestmodelweights50
plot(modelselection50, type="s")


setOldClass("rma.uni")
setMethod('getfit', 'rma.uni', function(object, ...) {
  if (object$test=="z") {
    cbind(estimate=coef(object), se=sqrt(diag(vcov(object))), df=Inf)
  } else {
    cbind(estimate=coef(object), se=sqrt(diag(vcov(object))), df=object$k-object$p)
  }
})

round(coef(modelselection50), 4)


r11<-rma(yi,vi, data=effect_prod50, method="REML")
funnel(r11,main="For Production")
regtest(r11, model="rma") ##Egger's regression test for funnel plot asymmetry,for the meta-analytic models
summary(r11)
t<-trimfill(r11)
t
funnel(t)

fsn(yi,vi, data=effect_prod50)

inf<-influence(r11)
plot(inf,plotdfbs = T)




la502<-subset(production,abs(Latitude)<=60)
la50_60<-subset(la502,abs(Latitude)>50)
View(la50_60)

effect_prod50_60<-escalc(measure ="ROM",data=la50_60,m1i=Prod.mean.E,sd1i=Prod.sd.E,n1i=Prod.n.E,m2i=Prod.mean.C,sd2i=Prod.sd.C,n2i=Prod.n.C)
library(dplyr)
effect_prod50_60<- select(effect_prod50_60,-Prod.mean.E,-Prod.sd.E,-Prod.n.E,-Prod.mean.C,-Prod.sd.C,-Prod.n.C,-ln.prod.RR,-Var.Prod.RR)
View(effect_prod50_60)

random_prod50_60<-rma(yi,vi, data=effect_prod50_60, method="REML")
summary(random_prod50_60)
attach(effect_prod50_60)
effect_prod50601<- effect_prod50_60[order(-vi),] 
attach(effect_prod50601)
forest(yi, vi,annotate=FALSE, col="grey", slab=NA, xlab="Log response ratio", pch=20,efac=0,font=2,cex=1,psize=1,xlim=c(-2,2), lty=c("solid","blank"))
abline(v=coef(random_prod50), lty="dashed",col="black")  
abline(v=0, lty="dotted",col="red") 
title("Forest of effect size for production(50<latitude<60°)",cex.main=1.2,font.main=2)


r_Nload50_60<-rma(yi,vi, mods=~Nload,data=effect_prod50_60, method="REML")
summary(r_Nload50_60)
r_totN50_60<-rma(yi,vi, mods=~tot_N,data=effect_prod50_60, method="REML")
summary(r_totN50_60)
pred_totN50_60 <- predict(r_totN50_60, newmods=c(0:190))
wi<- 1/sqrt(effect_prod50_60$vi)
size<-wi*0.1+0.8
plot(effect_prod50_60$tot_N, effect_prod50_60$yi, pch=19, cex=size, col="grey",ylab="Change in production (ln resposne ratio)",
     xlab="Total N (kg N ha-2 yr-1)",
     las=1, bty="l", font=1,cex.lab=1.2)
lines(0:190, pred_totN50_60$pred)
lines(0:190, pred_totN50_60$ci.lb, lty="dashed")
lines(0:190, pred_totN50_60$ci.ub, lty="dashed")
abline(h=0, lty="dotted",col="red",lwd=2)
abline(v=10, lty="dotted",col="blue",lwd=2)
text(c(120), 1.2, c("Qm=8.3337  p= 0.0039"),cex=1,font=1,col="red")
text(c(120), -1.5, c("y=-0.0018x+0.018"),cex=1,font=1,col="red")
text(c(15), -2.2, c("10"),cex=1,font=1,col="blue")
title(main="Critical of N concentration(50°<latitude<60°)", line=1.5, cex.lab=1.5)


r_tempmean50_60<-rma(yi,vi, mods=~tempmean,data=effect_prod50_60, method="REML")
summary(r_tempmean50_60)
r_Julytem50_60<-rma(yi,vi, mods=~Julytem,data=effect_prod50_60, method="REML")
summary(r_Julytem50_60)
r_prectot50_60<-rma(yi,vi, mods=~prectot,data=effect_prod50_60, method="REML")
summary(r_prectot50_60)
rprod_micro50_60<-rma(yi,vi, mods=~Hummock,data=effect_prod50_60, method="REML")
summary(rprod_micro50_60)
rprod_hummock50_60<-rma(yi,vi, data=subset(effect_prod50_60,Hummock=="1"), method="REML")
summary(rprod_hummock50_60)
rprod_hollow50_60<-rma(yi,vi, data=subset(effect_prod50_60,Hummock=="0"), method="REML")
summary(rprod_hollow50_60)


r_species50_60<-rma(yi,vi, mods=~Species,data=effect_prod50_60, method="REML")
summary(r_species50_60)
r_species50_601<-rma(yi,vi, mods=~Species-1,data=effect_prod50_60, method="REML")
summary(r_species50_601)
coef(summary(r_species50_601))
estimates_species50_601 <- coef(summary(r_species50_601))[,1]
se_species50_601 <- coef(summary(r_species50_601))[,2]
variances_species50_601<- se_species50_601^2
labels_species50_601 <- c( "Green capillifolium","Red capillifolium","balticum*fuscum",
                           "capillifolium ","contortum","cuspidatum","fallax", 
                           "fuscum", "magellanicum","papillosum","rubellum","squarrosum")
forest(estimates_species50_601, variances_species50_601, slab=labels_species50_601, font=1,cex=1,psize=1,
       xlim=c(-0.7,0.7),
       xlab="Effect size (log response ratio)")
text(c(-0.65,0.62), 13.5, c("Species", "Mean [95% CI]"),cex=1,font=1)
anova(r_species50_601, L=c(1,0,0,0,-1,0,0,0,0,0,0,0))
anova(r_species50_601, L=c(0,0,1,0,0,-1,0,0,0,0,0,0))



total50_60<-rma(yi,vi,mods=~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory,data=effect_prod50_60, method="REML")
summary(total50_60)
h_total50_60<-rma.mv(yi,vi,mods = ~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory,data=effect_prod50_60,random =list(~1|IDnr,~1|ID),method = "ML")
summary(h_total50_60)
rma.glmulti <- function(formula, data, ...) { rma(formula, vi, data=data, method="ML", ...)} 
modelselection50_60<- glmulti(yi ~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory, data=effect_prod50_60, level=1, fitfunction=rma.glmulti, crit="aicc", confsetsize=1000)
plot(modelselection50_60)

modelweights50_60<- weightable(modelselection50_60)
modelweights50_60
summary(modelselection50_60@objects[[1]])

bestmodelweights50_60<- modelweights50_60[modelweights50_60$aicc <= min(modelweights50_60$aicc) + 2,]

bestmodelweights50_60
plot(modelselection50_60, type="s")

setOldClass("rma.uni")
setMethod('getfit', 'rma.uni', function(object, ...) {
  if (object$test=="z") {
    cbind(estimate=coef(object), se=sqrt(diag(vcov(object))), df=Inf)
  } else {
    cbind(estimate=coef(object), se=sqrt(diag(vcov(object))), df=object$k-object$p)
  }
})

round(coef(modelselection50_60), 4)
r11<-rma(yi,vi, data=effect_prod50_60, method="REML")
funnel(r11,main="For Production")
regtest(r11, model="rma") ##Egger's regression test for funnel plot asymmetry,for the meta-analytic models
summary(r11)
t<-trimfill(r11)
t
funnel(t)


fsn(yi,vi, data=effect_prod50_60)


inf<-influence(r11)
plot(inf,plotdfbs = T)





la60<-subset(production,abs(Latitude)>60)
View(la60)
effect_prod60<-escalc(measure ="ROM",data=la60,m1i=Prod.mean.E,sd1i=Prod.sd.E,n1i=Prod.n.E,m2i=Prod.mean.C,sd2i=Prod.sd.C,n2i=Prod.n.C)
library(dplyr)
effect_prod60<- select(effect_prod60,-Prod.mean.E,-Prod.sd.E,-Prod.n.E,-Prod.mean.C,-Prod.sd.C,-Prod.n.C,-ln.prod.RR,-Var.Prod.RR)
View(effect_prod60)


random_prod60<-rma(yi,vi, data=effect_prod60, method="REML")
summary(random_prod60)
attach(effect_prod60)
effect_prod601<- effect_prod60[order(-vi),] 
attach(effect_prod601)
forest(yi, vi,annotate=FALSE, col="grey", slab=NA, xlab="Log response ratio", pch=20,efac=0,font=2,cex=1,psize=1,xlim=c(-2,2), lty=c("solid","blank"))
abline(v=coef(random_prod50), lty="dashed",col="black")  
abline(v=0, lty="dotted",col="red") 
title("Forest of effect size for production(latitude>60°)",cex.main=1.2,font.main=2)


r_Nload60<-rma(yi,vi, mods=~Nload,data=effect_prod60, method="REML")
summary(r_Nload60)
pred_Nload60 <- predict(r_Nload60, newmods=c(0:100))
wi<- 1/sqrt(effect_prod60$vi)
size<-wi*0.1+0.8
plot(effect_prod60$Nload, effect_prod60$yi, pch=19, cex=size, col="grey",ylab="Change in production (ln resposne ratio)",
     xlab="N application(kg N ha-2 yr-1)",
     las=1, bty="l", font=1,cex.lab=1.2)
lines(0:100, pred_Nload60$pred)
lines(0:100, pred_Nload60$ci.lb, lty="dashed")
lines(0:100,pred_Nload60$ci.ub, lty="dashed")
abline(h=0, lty="dotted",col="red",lwd=2)
abline(v=27.21, lty="dotted",col="blue",lwd=2)
text(c(80), 1, c("Qm=7.4575  p= 0.0063"),cex=1,font=1,col="red")
text(c(80), -1, c("y=-0.0048x+0.1306"),cex=1,font=1,col="red")
text(c(35), -1.4, c("27.21"),cex=1,font=1,col="blue")
title(main="Critical of N concentration(latitude>60°)", line=1.5, cex.lab=1.5)


r_totN60<-rma(yi,vi, mods=~tot_N,data=effect_prod60, method="REML")
summary(r_totN60)
pred_totN60 <- predict(r_totN60, newmods=c(0:110))
wi<- 1/sqrt(effect_prod60$vi)
size<-wi*0.1+0.8
plot(effect_prod60$tot_N, effect_prod60$yi, pch=19, cex=size, col="grey",ylab="Change in production (ln resposne ratio)",
     xlab="Total N (kg N ha-2 yr-1)",
     las=1, bty="l", font=1,cex.lab=1.2)
lines(0:110, pred_totN60$pred)
lines(0:110, pred_totN60$ci.lb, lty="dashed")
lines(0:110, pred_totN60$ci.ub, lty="dashed")
abline(h=0, lty="dotted",col="red",lwd=2)
abline(v=29.23, lty="dotted",col="blue",lwd=2)
text(c(80), 1, c("Qm=7.2848  p= 0.0070"),cex=1,font=1,col="red")
text(c(80), -1, c("y=-0.0047x+0.1374"),cex=1,font=1,col="red")
text(c(37), -1.3, c("29.23"),cex=1,font=1,col="blue")
title(main="Critical of N concentration(latitude>60°)", line=1.5, cex.lab=1.5)


r_tempmean60<-rma(yi,vi, mods=~tempmean,data=effect_prod60, method="REML")
summary(r_tempmean60)
pred_tempmean60 <- predict(r_tempmean60, newmods=c(-3:9))
wi<- 1/sqrt(effect_prod60$vi)
size<-wi*0.1+0.8
plot(effect_prod60$tempmean, effect_prod60$yi, pch=19, cex=size, col="grey",ylab="Change in production (ln resposne ratio)",
     xlab="Annual temperature (℃)",
     las=1, bty="l", font=1,cex.lab=1.2)
lines(-3:9, pred_tempmean60$pred)
lines(-3:9, pred_tempmean60$ci.lb, lty="dashed")
lines(-3:9, pred_tempmean60$ci.ub, lty="dashed")
abline(h=0, lty="dotted",col="red",lwd=2)

r_Julytem60<-rma(yi,vi, mods=~Julytem,data=effect_prod60, method="REML")
summary(r_Julytem60)
r_prectot60<-rma(yi,vi, mods=~prectot,data=effect_prod60, method="REML")
summary(r_prectot60)
pred_prectot60 <- predict(r_prectot60, newmods=c(530:690))
wi<- 1/sqrt(effect_prod60$vi)
size<-wi*0.1+0.8
plot(effect_prod60$prectot, effect_prod60$yi, pch=19, cex=size, col="grey",ylab="Change in production (ln resposne ratio)",
     xlab="Annual precipitation (mm)",
     las=1, bty="l", font=1,cex.lab=1.2)
lines(530:690, pred_prectot60$pred)
lines(530:690, pred_prectot60$ci.lb, lty="dashed")
lines(530:690, pred_prectot60$ci.ub, lty="dashed")
abline(h=0, lty="dotted",col="red",lwd=2)

rprod_micro60<-rma(yi,vi, mods=~Hummock,data=effect_prod60, method="REML")
summary(rprod_micro60)
rprod_hummock60<-rma(yi,vi, data=subset(effect_prod60,Hummock=="1"), method="REML")
summary(rprod_hummock60)
rprod_hollow60<-rma(yi,vi, data=subset(effect_prod60,Hummock=="0"), method="REML")
summary(rprod_hollow60)

r_species60<-rma(yi,vi, mods=~Species,data=effect_prod60, method="REML")
summary(r_species60)
r_species601<-rma(yi,vi, mods=~Species-1,data=effect_prod60, method="REML")
summary(r_species601)
coef(summary(r_species601))
estimates_species601 <- coef(summary(r_species601))[,1]
se_species601 <- coef(summary(r_species601))[,2]
variances_species601<- se_species601^2
labels_species601 <- c( "angustifolium","balticum","balticum*fuscum","fallax", 
                        "fuscum", "magellanicum","papillosum","rubellum","warnstorfii")
forest(estimates_species601, variances_species601, slab=labels_species601, font=1,cex=1,psize=1,
       xlim=c(-0.7,0.7),
       xlab="Effect size (log response ratio)")
text(c(-0.9,0.88), 10.5, c("Species", "Mean [95% CI]"),cex=1,font=1)
anova(r_species601, L=c(1,0,-1,0,0,0,0,0,0))
anova(r_species601, L=c(0,0,0,1,-1,0,0,0,0))



total60<-rma(yi,vi,mods=~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory,data=effect_prod60, method="REML")
summary(total60)
h_total60<-rma.mv(yi,vi,mods = ~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory,data=effect_prod60,random =list(~1|IDnr,~1|ID),method = "ML")
summary(h_total60)
rma.glmulti <- function(formula, data, ...) { rma(formula, vi, data=data, method="ML", ...)} 
modelselection60<- glmulti(yi ~ Nload+Julytem+Hummock+prectot+tempmean+Nhistory, data=effect_prod60, level=1, fitfunction=rma.glmulti, crit="aicc", confsetsize=1000)
plot(modelselection60)


modelweights60<- weightable(modelselection60)
modelweights60
summary(modelselection60@objects[[1]])

bestmodelweights60<- modelweights60[modelweights60$aicc <= min(modelweights60$aicc) + 2,]


bestmodelweights60
plot(modelselection60, type="s")


setOldClass("rma.uni")
setMethod('getfit', 'rma.uni', function(object, ...) {
  if (object$test=="z") {
    cbind(estimate=coef(object), se=sqrt(diag(vcov(object))), df=Inf)
  } else {
    cbind(estimate=coef(object), se=sqrt(diag(vcov(object))), df=object$k-object$p)
  }
})


round(coef(modelselection60), 4)


r11<-rma(yi,vi, data=effect_prod60, method="REML")
funnel(r11,main="For Production")
regtest(r11, model="rma") ##Egger's regression test for funnel plot asymmetry,for the meta-analytic models
summary(r11)
t<-trimfill(r11)
t
funnel(t)


fsn(yi,vi, data=effect_prod60)


inf<-influence(r11)
plot(inf,plotdfbs = T)