########################################################################
##############Created by Shawn (Shaofei) Jin-PeerJ######################
########################################################################

#Library needed
#
library(metafor)
library(ggplot2)


# Creat Figure 1====
obs.rmd=c("0.667","0.855","0.939","0.979","1")
obs.depth=c("50","100","150","200","240")
obs=data.frame(as.numeric(obs.rmd),as.numeric(obs.depth))
a=7.188*10^4
b=1.716
c=45
d=1836.97
intergrand= function(x){ a/((x^b)*(exp(c/x)+d))}


inter.rmdend=integrate(intergrand,lower=0,upper=300)
inter.rmd=function(y) {
  result=integrate(intergrand,lower=1,upper=y,
                   stop.on.error = FALSE)$value/(inter.rmdend$value)
  return(as.numeric(result))
}
#inter.rmd1=function(x)inter.rmd(x)$value

vi=Vectorize(inter.rmd,SIMPLIFY = T)
plot(vi ,from=0,to=300,lwd=2,
     xlab="Soil depth(cm)",ylab="Accumulative ratio")
points(obs$as.numeric.obs.depth.,obs$as.numeric.obs.rmd.,type="p",
       xlim=c(0,300))
lines( curve(vi ,from=0,to=300,lwd=2,
             xlab="Soil depth(cm)",ylab="Accumulative ratio"))  

abline(h=c(0.9,0.97),v=c(80,180))
#curve(intergrand,from=0,to=300)
lines(obs.depth,obs.rmd,type="p")

root=read.csv("root.csv",header = T)
ggplot(data.frame(x = c(0, 300)), aes(x = x)) +
  stat_function(fun = vi)+
  geom_point(data=root,aes(x=depth,y=ratio))


p=ggplot(root,aes(x=depth,y=ratio))+
  #geom_line()+
  stat_function(fun = vi,size=1.2,colour="red")+
  geom_point(size=1.5)+
  geom_hline(yintercept=0.9,size=1.1,linetype=2)+
  geom_vline(xintercept=82,size=1.1,linetype=2)+
  theme_bw()+
  theme(text=element_text(family="sans",size=14))+
  xlab("Depth (cm)")+
  ylab("Accumulative ratio of root mass")+
  coord_flip()+
  scale_x_reverse()
ggsave("figure1peerJ.pdf",width=12,heigh=12,units="cm")


nmeta=read.csv("nitrogendata.csv",
               header = T)

#Creat the data====
nes=escalc(measure = "ROM",
           m1i=mean2014,
           sd1i=sd2014,
           n1i=n2014,
           m2i=mean2002,
           sd2i=sd2002/sqrt(6),
           n2i=n2002,
           data=nmeta)

nes$Depth=factor(nes$Depth,order=TRUE,
                 levels=c("180-200","160-180","140-160",
                          "120-140","100-120","80-100","60-80",
                          "40-60","20-40","0-20")) 
ggplot(nes,
       aes(x=Depth,y=yi,
           ymax=yi+sqrt(vi),
           ymin=yi-sqrt(vi)))+
  geom_pointrange()+
  coord_flip()+
  geom_hline(yintercept=0,linetype=2)+
  facet_grid(N~Carbon,scales = "free")

nfit=rma.uni(yi, vi, 
             data = nes[which(nes$Carbon=="SIC" & nes$bio=="root"),])

nfit=rma.uni(yi, vi, 
             data = nes[which(nes$Carbon=="SIC" ),],
             subset=(bio=="root"),
             mods=cbind(Depth))
nfit$subset
example(forest)

forest(nfit)
funnel(nfit)
print.rma.uni(nfit)

trial=matrix(NA,nrow=18,ncol =9)
n=1
for (i in levels(nes$N)){
  for (j in levels(nes$Carbon)){
    for (k in levels(nes$bio)){
      nesdata=subset(nes,N==i &
                       Carbon==j &
                       bio==k)
      nesfit=rma.uni(yi,vi,data=nesdata)
      trial[n,]=c(i,j,k,nesfit$b[1],
                  nesfit$se,
                  nesfit$ci.lb,
                  nesfit$ci.ub,
                  nesfit$b[1]-nesfit$se*0.675, 
                  nesfit$b[1]+nesfit$se*0.675)
      n=n+1
    }
  }
}
colnames(trial)=c("nitrogen","carbon","root",
                  "rr","rrse","rrl","rrup","rr5l","rr5up")
trial1=matrix(NA,nrow=9,ncol =9)
n=1
for (i in levels(nes$N)){
  for (j in levels(nes$Carbon)){
    
    nesdata=subset(nes,N==i &
                     Carbon==j )
    nesfit=rma.uni(yi,vi,data=nesdata)
    trial1[n,]=c(i,j,"whole",nesfit$b[1],
                 nesfit$se,
                 nesfit$ci.lb,
                 nesfit$ci.ub,
                 nesfit$b[1]-nesfit$se*0.675, 
                 nesfit$b[1]+nesfit$se*0.675)
    n=n+1
  }
}
colnames(trial1)=c("nitrogen","carbon","root",
                   "rr","rrse","rrl","rrup","rr5l","rr5up")
trial=rbind(trial,trial1)
write.csv(trial,"result2.csv")


# Creat Figure 2====
rr=read.csv("result2.csv",header = T)
rr$root=factor(rr$root,order=TRUE,
               levels=c("Root","Non-root","Whole"))
ggplot(rr,aes(x=root,y=exp(rr)))+
  
  coord_flip()+
  geom_errorbar(aes(ymin=exp(rrl),ymax=exp(rrup)),
                position=position_dodge(0.9),width=0.5,stat="identity")+
  geom_errorbar(aes(ymin=exp(rr5l),ymax=exp(rr5up),size=0.5),colour="red",
                position=position_dodge(0.9),width=0.0,stat="identity")+
  geom_point()+
  labs(y="Effect size",x="",size=12)+
  geom_hline(yintercept=1,linetype=2)+
  scale_x_discrete(limits=rev(c("Root","Non-root","Whole")))+
  theme_bw()+
  theme(text=element_text(family="sans",size=12))+
  theme(strip.text=element_text(family="sans",size=12))+
  theme(strip.background = element_rect(fill="white"))+
  theme(strip.text=element_text(size=10))+
  theme(panel.grid.major=element_blank(),
        panel.grid.minor=element_blank())+
  theme(legend.position="none")+
  facet_grid(nitrogen~carbon,scales = "free") 
ggsave("figure2pperJ.pdf",width=16,heigh=12,units="cm")



# Creat Figure 3====
sic=read.csv("carbonchange.csv",header = T)

ggplot(sic,aes(x=X,y=D))+
  geom_bar(stat="identity", position="identity", 
           colour="grey", size=0.1,width = 0.5)+
  geom_errorbar(aes(ymin=low,ymax=up),
                position=position_dodge(0.9),width=0.2,stat="identity")+
  geom_hline(yintercept = 0)+
  xlab("")+
  ylab("Changes in soil carbon (PgC)") +
  theme_bw()+
  theme(text=element_text(family="sans",size=12))+
  facet_grid(carbon~.,scales = "free")

ggsave("####",width=10,heigh=12,units="cm",dpi=800)