# This file is the R codes for Figures and Tables  in this manuscript 
###### install packages
install.packages (c("metafor","VIM","mice","lmodel2","smatr","ggplot2", "maps","maptools","scales"))
library(VIM)
library(mice) 
library(metafor)
library(lmodel2)
library(smatr)
rm(list=ls())
library(ggplot2)
library(splines)
library(grid)

##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
# Figure 1: plotting sites in maps
library(mapproj)
map(database= "world", col="grey80", fill=TRUE, projection="gilbert")
site<-read.csv(file.choose())
coord <- mapproject(site$longitude, site$Latitude, proj="gilbert")  #convert points to projected lat/long
points(coord, pch=20, cex=1.2, col="red")  #plot converted points
map("world", fill = TRUE, col = rainbow(200), ylim = c(-60, 90), mar = c(0, 0, 0, 0))
points(site$longitude, site$Latitude, col="grey80", pch=18)  
map.axes() 
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
##### Figure 2-4a): effects between plant life forms 
### read in data
mydata<-read.csv(file.choose())
### explore data
head(mydata)
dir()
ls()
names(mydata)
dim(mydata)  
## subset data
growthform<-mydata[,-(1:8)]; growthform<-growthform[,-(2:7)]; growthform<-growthform[,-(23:41)]  
forb<-subset(growthform, growthform$PFT=="01Forb")
grass<-subset(growthform, growthform$PFT=="02Grass")
conifer<-subset(growthform, growthform$PFT=="03Conifer")
evergreen<-subset(growthform, growthform$PFT=="04Evergreen")
deciduous<-subset(growthform, growthform$PFT=="05Deciduous")
mixed<-subset(growthform, growthform$PFT=="13Mixed")
fgcedm<-rbind(forb,grass, conifer, evergreen,deciduous, mixed) 
overall<-growthform 
overall$PFT<-"12all"
data<-rbind(fgcedm,overall) 
##### 1.1 [N]green
plantN<-data[,1:24]; plantN<-plantN[,-(5:10)]; plantN<-plantN[,-(11:17)]
greenleafN<-subset(plantN, plantN$Tissue2=="01greenleaf")
liverootN<-subset(plantN, plantN$Tissue2=="03liveroot")
greenN<-rbind(greenleafN,liverootN)
### re-name variables
names(greenN)[5]<-"Xc" 
names(greenN)[6]<-"Xe" 
names(greenN)[7]<-"Sc" 
names(greenN)[8]<-"Se" 
names(greenN)[9]<-"Nc" 
names(greenN)[10]<-"Ne"
greenN<-na.omit(greenN)
##### 1.2 [P]green
plantP<-data[,1:24]; plantP<-plantP[,-(5:16)]; plantP<-plantP[,-(11)]
greenleafP<-subset(plantP, plantP$Tissue2=="01greenleaf")
liverootP<-subset(plantP, plantP$Tissue2=="03liveroot")
greenP<-rbind(greenleafP,liverootP)
### re-name variables
names(greenP)[5]<-"Xc" 
names(greenP)[6]<-"Xe" 
names(greenP)[7]<-"Sc" 
names(greenP)[8]<-"Se" 
names(greenP)[9]<-"Nc" 
names(greenP)[10]<-"Ne"
greenP<-na.omit(greenP)
##### 1.3 NRP
deadN<-data[,1:24]; deadN<-deadN[,-(5:10)]; deadN<-deadN[,-(11:17)]
deadleafN<-subset(deadN, deadN$Tissue2=="02deadleaf")
deadrootN<-subset(deadN, deadN$Tissue2=="04deadroot")
NRP<-rbind(deadleafN,deadrootN)
### re-name variables
names(NRP)[5]<-"Xc" 
names(NRP)[6]<-"Xe" 
names(NRP)[7]<-"Sc" 
names(NRP)[8]<-"Se" 
names(NRP)[9]<-"Nc" 
names(NRP)[10]<-"Ne"
greenP<-na.omit(greenP)
##### 1.4 PRP
deadP<-data[,1:24]; deadP<-deadP[,-(5:16)]; deadP<-deadP[,-(11)]
deadleafP<-subset(deadP, deadP$Tissue2=="02deadleaf")
deadrootP<-subset(deadP, deadP$Tissue2=="04deadroot")
PRP<-rbind(deadleafP,deadrootP)
### re-name variables
names(PRP)[5]<-"Xc" 
names(PRP)[6]<-"Xe" 
names(PRP)[7]<-"Sc" 
names(PRP)[8]<-"Se" 
names(PRP)[9]<-"Nc" 
names(PRP)[10]<-"Ne"
PRP<-na.omit(PRP)
##### 1.5 NRE
NRE<-data[,1:30]; NRE<-NRE[,-(5:23)]
NRE1 <- NRE[,-(5)]; NRE2<- NRE[,1:5]; NRE3 <- cbind(NRE1, NRE2); NRE<-NRE3[,-(11:14)]
### re-name variables
names(NRE)[5]<-"Xc" 
names(NRE)[6]<-"Xe" 
names(NRE)[7]<-"Sc" 
names(NRE)[8]<-"Se" 
names(NRE)[9]<-"Nc" 
names(NRE)[10]<-"Ne"
NRE<-na.omit(NRE)
##### 1.6 PRE
PRE<-data[,1:36]; PRE<-PRE[,-(5:23)]; PRE<-PRE[,-(6:11)]
PRE1 <- PRE[,-(5)]; PRE2<- PRE[,1:5]; PRE3 <- cbind(PRE1, PRE2); PRE<-PRE3[,-(11:14)]
### re-name variables
names(PRE)[5]<-"Xc" 
names(PRE)[6]<-"Xe" 
names(PRE)[7]<-"Sc" 
names(PRE)[8]<-"Se" 
names(PRE)[9]<-"Nc" 
names(PRE)[10]<-"Ne"
PRE<-na.omit(PRE)
### subset data
## grreen-leaf [N] subset
NfertgreenN<-subset(greenN, greenN$treatment2=="05N")
PfertgreenN<-subset(greenN, greenN$treatment2=="06P")
NPfertgreenN<-subset(greenN, greenN$treatment2=="21NP")
## grreen-leaf [P] subset
NfertgreenP<-subset(greenP, greenP$treatment2=="05N")
PfertgreenP<-subset(greenP, greenP$treatment2=="06P")
NPfertgreenP<-subset(greenP, greenP$treatment2=="21NP")
## senesced-leaf [N] subset
NfertNRP<-subset(NRP, NRP$treatment2=="05N")
PfertNRP<-subset(NRP, NRP$treatment2=="06P")
NPfertNRP<-subset(NRP, NRP$treatment2=="21NP")
## senesced-leaf [P] subset
NfertPRP<-subset(PRP, PRP$treatment2=="05N")
PfertPRP<-subset(PRP, PRP$treatment2=="06P")
NPfertPRP<-subset(PRP, PRP$treatment2=="21NP")
## Reesorption of N subset
NfertNRE<-subset(NRE, NRE$treatment2=="05N")
PfertNRE<-subset(NRE, NRE$treatment2=="06P")
NPfertNRE<-subset(NRE, NRE$treatment2=="21NP")
## Reesorption of P subset
NfertPRE<-subset(PRE, PRE$treatment2=="05N")
PfertPRE<-subset(PRE, PRE$treatment2=="06P")
NPfertPRE<-subset(PRE, PRE$treatment2=="21NP")
###### Run R and get the response ratios 
#### Figure 2a): N fert on green [N]
NfertgreenNdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertgreenN, append=TRUE)
head(NfertgreenNdat)
yi <- with(NfertgreenN, log(Xe/Xc))
vi <- with(NfertgreenN, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertgreenN<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NfertgreenNdat, digits=8)  
summary(rNfertgreenN)
#### Figure 2a): P fert on green [N]
PfertgreenNdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertgreenN, append=TRUE)
head(PfertgreenNdat)
yi <- with(PfertgreenN, log(Xe/Xc))
vi <- with(PfertgreenN, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertgreenN<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=PfertgreenNdat, digits=8)  
summary(rPfertgreenN)
#### Figure 2a): N+P fert on green [N]
NPfertgreenNdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertgreenN, append=TRUE)
head(NPfertgreenNdat)
yi <- with(NPfertgreenN, log(Xe/Xc))
vi <- with(NPfertgreenN, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertgreenN<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NPfertgreenNdat, digits=8)  
summary(rNPfertgreenN)
#### Figure 2a): N fert on green [P]
NfertgreenPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertgreenP, append=TRUE)
head(NfertgreenPdat)
yi <- with(NfertgreenP, log(Xe/Xc))
vi <- with(NfertgreenP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertgreenP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NfertgreenPdat, digits=8)  
summary(rNfertgreenP)
#### Figure 2a): P fert on green [P]
PfertgreenPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertgreenP, append=TRUE)
head(PfertgreenPdat)
yi <- with(PfertgreenP, log(Xe/Xc))
vi <- with(PfertgreenP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertgreenP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=PfertgreenPdat, digits=8)  
summary(rPfertgreenP)
#### Figure 2a): N+P fert on green [P]
NPfertgreenPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertgreenP, append=TRUE)
head(NPfertgreenPdat)
yi <- with(NPfertgreenP, log(Xe/Xc))
vi <- with(NPfertgreenP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertgreenP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NPfertgreenPdat, digits=8)  
summary(rNPfertgreenP)
#### Figure 3a): N fert on senesced-leaf [N]
NfertNRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertNRP, append=TRUE)
head(NfertNRPdat)
yi <- with(NfertNRP, log(Xe/Xc))
vi <- with(NfertNRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertNRP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NfertNRPdat, digits=8)  
summary(rNfertNRP)
#### Figure 3a): P fert on senesced-leaf [N]
PfertNRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertNRP, append=TRUE)
head(PfertNRPdat)
yi <- with(PfertNRP, log(Xe/Xc))
vi <- with(PfertNRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertNRP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=PfertNRPdat, digits=8)  
summary(rPfertNRP)
#### Figure 3a): N+P fert on senesced-leaf [N]
NPfertNRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertNRP, append=TRUE)
head(NPfertNRPdat)
yi <- with(NPfertNRP, log(Xe/Xc))
vi <- with(NPfertNRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertNRP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NPfertNRPdat, digits=8)  
summary(rNPfertNRP)

#### Figure 3a): N fert on senesced-leaf [P]
NfertPRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertPRP, append=TRUE)
head(NfertPRPdat)
yi <- with(NfertPRP, log(Xe/Xc))
vi <- with(NfertPRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertPRP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NfertPRPdat, digits=8)  
summary(rNfertPRP)
#### Figure 3a): P fert on senesced-leaf [P]
PfertPRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertPRP, append=TRUE)
head(PfertPRPdat)
yi <- with(PfertPRP, log(Xe/Xc))
vi <- with(PfertPRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertPRP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=PfertPRPdat, digits=8)  
summary(rPfertPRP)
 #### Figure 3a): N+P fert on senesced-leaf [P]
NPfertPRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertPRP, append=TRUE)
head(NPfertPRPdat)
yi <- with(NPfertPRP, log(Xe/Xc))
vi <- with(NPfertPRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertPRP<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NPfertPRPdat, digits=8)  
summary(rNPfertPRP)

#### Figure 4a): N fert on NRE
NfertNREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertNRE, append=TRUE)
head(NfertNREdat)
yi <- with(NfertNRE, log(Xe/Xc))
vi <- with(NfertNRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertNRE<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NfertNREdat, digits=8)  
summary(rNfertNRE)
#### Figure 4a): P fert on NRE
PfertNREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertNRE, append=TRUE)
head(PfertNREdat)
yi <- with(PfertNRE, log(Xe/Xc))
vi <- with(PfertNRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertNRE<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=PfertNREdat, digits=8)  
summary(rPfertNRE)
#### Figure 4a): N+P fert on NRE
NPfertNREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertNRE, append=TRUE)
head(NPfertNREdat)
yi <- with(NPfertNRE, log(Xe/Xc))
vi <- with(NPfertNRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertNRE<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NPfertNREdat, digits=8)  
summary(rNPfertNRE)

#### Figure 4a): N fert on PRE
NfertPREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertPRE, append=TRUE)
head(NfertPREdat)
yi <- with(NfertPRE, log(Xe/Xc))
vi <- with(NfertPRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertPRE<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NfertPREdat, digits=8)  
summary(rNfertPRE)
#### Figure 4a): P fert on PRE
PfertPREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertPRE, append=TRUE)
head(PfertPREdat)
yi <- with(PfertPRE, log(Xe/Xc))
vi <- with(PfertPRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertPRE<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=PfertPREdat, digits=8)  
summary(rPfertPRE)
#### Figure 4a): N+P fert on PRE
NPfertPREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertPRE, append=TRUE)
head(NPfertPREdat)
yi <- with(NPfertPRE, log(Xe/Xc))
vi <- with(NPfertPRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertPRE<-rma.uni(yi,vi,mods=~factor(PFT)-1, method="REML", data=NPfertPREdat, digits=8)  
summary(rNPfertPRE)


####### Figure 2-4b): effects between ecosystems
## read in data
mydata<-read.csv(file.choose())
## subset data
bf<-subset(mydata, mydata$BIOME=="01boreal")
temf<-subset(mydata, mydata$BIOME=="02temF")
tropf<-subset(mydata, mydata$BIOME=="03tropF")
temg<-subset(mydata, mydata$BIOME=="04temG")
trog<-subset(mydata, mydata$BIOME=="05tropG")
tundra<-subset(mydata, mydata$BIOME=="06tundra")
wetland<-subset(mydata, mydata$BIOME=="07wetland")
desert<-subset(mydata, mydata$BIOME=="08desert")
crop<-subset(mydata, mydata$BIOME=="09crop")
overall<-mydata 
overall$BIOME<-"12all"
dataecoall<-rbind(bf,temf,tropf,temg,trog,tundra, wetland, desert, crop, overall) 
dataeco<-dataecoall[,-(1:15)]; dataeco<-dataeco[,-(4:8)]; dataeco<-dataeco[,-(20:35)] 
##### [N]green by ecosystems
plantNeco<-dataeco[,1:10]
greenleafNeco<-subset(plantNeco, plantNeco$Tissue2=="01greenleaf")
liverootNeco<-subset(plantNeco, plantNeco$Tissue2=="03liveroot")
greenNeco<-rbind(greenleafNeco,liverootNeco)
### re-name variables
names(greenNeco)[5]<-"Xc" 
names(greenNeco)[6]<-"Xe" 
names(greenNeco)[7]<-"Sc" 
names(greenNeco)[8]<-"Se" 
names(greenNeco)[9]<-"Nc" 
names(greenNeco)[10]<-"Ne"
greenNeco<-na.omit(greenNeco)
#####[P]green by ecosystems
plantPeco<-dataeco[,1:16]
greenleafPeco<-subset(plantPeco, plantPeco$Tissue2=="01greenleaf")
liverootPeco<-subset(plantPeco, plantPeco$Tissue2=="03liveroot")
greenPeco<-rbind(greenleafPeco,liverootPeco)
greenPeco<-greenPeco[,-(5:10)]
### re-name variables
names(greenPeco)[5]<-"Xc" 
names(greenPeco)[6]<-"Xe" 
names(greenPeco)[7]<-"Sc" 
names(greenPeco)[8]<-"Se" 
names(greenPeco)[9]<-"Nc" 
names(greenPeco)[10]<-"Ne"
greenPeco<-na.omit(greenPeco)
###NRP by biome
deadNeco<-dataeco[,1:10]
deadleafNeco<-subset(deadNeco, deadNeco$Tissue2=="02deadleaf")
deadrootNeco<-subset(deadNeco, deadNeco$Tissue2=="04deadroot")    
NRPeco<-rbind(deadleafNeco,deadrootNeco)
### re-name variables
names(NRPeco)[5]<-"Xc" 
names(NRPeco)[6]<-"Xe" 
names(NRPeco)[7]<-"Sc" 
names(NRPeco)[8]<-"Se" 
names(NRPeco)[9]<-"Nc" 
names(NRPeco)[10]<-"Ne"
NRPeco<-na.omit(NRPeco)
#####PRP
deadPeco<-dataeco[,1:16]; deadPeco<-deadPeco[,-(5:10)]
deadleafPeco<-subset(deadPeco, deadPeco$Tissue2=="02deadleaf")
deadrootPeco<-subset(deadPeco, deadPeco$Tissue2=="04deadroot")    
PRPeco<-rbind(deadleafPeco,deadrootPeco)
### re-name variables
names(PRPeco)[5]<-"Xc" 
names(PRPeco)[6]<-"Xe" 
names(PRPeco)[7]<-"Sc" 
names(PRPeco)[8]<-"Se" 
names(PRPeco)[9]<-"Nc" 
names(PRPeco)[10]<-"Ne"
PRPeco<-na.omit(PRPeco)
#####NRE by ecosystems
NREeco<-dataeco[,1:27]; NREeco<-NREeco[,-(5:21)]
### re-name variables
names(NREeco)[5]<-"Xc" 
names(NREeco)[6]<-"Xe" 
names(NREeco)[7]<-"Sc" 
names(NREeco)[8]<-"Se" 
names(NREeco)[9]<-"Nc" 
names(NREeco)[10]<-"Ne"
NREeco<-na.omit(NREeco)
#####PRE by ecosystems
PREeco<-dataeco[,-(5:27)]
### re-name variables
names(PREeco)[5]<-"Xc" 
names(PREeco)[6]<-"Xe" 
names(PREeco)[7]<-"Sc" 
names(PREeco)[8]<-"Se" 
names(PREeco)[9]<-"Nc" 
names(PREeco)[10]<-"Ne"
PREeco<-na.omit(PREeco)
### subset data 
## green-leaf [N]
NfertgreenNeco<-subset(greenNeco, greenNeco$treatment2=="05N")
PfertgreenNeco<-subset(greenNeco, greenNeco$treatment2=="06P")
NPfertgreenNeco<-subset(greenNeco, greenNeco$treatment2=="21NP")
## green-leaf [P]
NfertgreenPeco<-subset(greenPeco, greenPeco$treatment2=="05N")
PfertgreenPeco<-subset(greenPeco, greenPeco$treatment2=="06P")
NPfertgreenPeco<-subset(greenPeco, greenPeco$treatment2=="21NP")
## senesced-leaf [N]
NfertNRPeco<-subset(NRPeco, NRPeco$treatment2=="05N")
PfertNRPeco<-subset(NRPeco, NRPeco$treatment2=="06P")
NPfertNRPeco<-subset(NRPeco, NRPeco$treatment2=="21NP")
## senesced-leaf [P]
NfertPRPeco<-subset(PRPeco, PRPeco$treatment2=="05N")
PfertPRPeco<-subset(PRPeco, PRPeco$treatment2=="06P")
NPfertPRPeco<-subset(PRPeco, PRPeco$treatment2=="21NP")
## NRE
NfertNREeco<-subset(NREeco, NREeco$treatment2=="05N")
PfertNREeco<-subset(NREeco, NREeco$treatment2=="06P")
NPfertNREeco<-subset(NREeco, NREeco$treatment2=="21NP")
## PRE
NfertPREeco<-subset(PREeco, PREeco$treatment2=="05N")
PfertPREeco<-subset(PREeco, PREeco$treatment2=="06P")
NPfertPREeco<-subset(PREeco, PREeco$treatment2=="21NP")
###### Run R codes to get Figs 2-4b) data
## Fig. 2b): N fert on green-leaf [N]
NfertgreenNecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertgreenNeco, append=TRUE)
head(NfertgreenNecodat)
yi <- with(NfertgreenNeco, log(Xe/Xc))
vi <- with(NfertgreenNeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertgreenNeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NfertgreenNecodat, digits=8)  
summary(rNfertgreenNeco)
## Fig. 2b): P fert on green-leaf [N]
PfertgreenNecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertgreenNeco, append=TRUE)
head(PfertgreenNecodat)
yi <- with(PfertgreenNeco, log(Xe/Xc))
vi <- with(PfertgreenNeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertgreenNeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=PfertgreenNecodat, digits=8)  
summary(rPfertgreenNeco)
## Fig. 2b): N+P fert on green-leaf [N]
NPfertgreenNecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertgreenNeco, append=TRUE)
head(NPfertgreenNecodat)
yi <- with(NPfertgreenNeco, log(Xe/Xc))
vi <- with(NPfertgreenNeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertgreenNeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NPfertgreenNecodat, digits=8)  
summary(rNPfertgreenNeco)
## Fig. 2b): N fert on green-leaf [P]
NfertgreenPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertgreenPeco, append=TRUE)
head(NfertgreenPecodat)
yi <- with(NfertgreenPeco, log(Xe/Xc))
vi <- with(NfertgreenPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertgreenPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NfertgreenPecodat, digits=8)  
summary(rNfertgreenPeco)
## Fig. 2b): P fert on green-leaf [P]
PfertgreenPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertgreenPeco, append=TRUE)
head(PfertgreenPecodat)
yi <- with(PfertgreenPeco, log(Xe/Xc))
vi <- with(PfertgreenPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertgreenPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=PfertgreenPecodat, digits=8)  
summary(rPfertgreenPeco)
## Fig. 2b): N+P fert on green-leaf [P]
NPfertgreenPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertgreenPeco, append=TRUE)
head(NPfertgreenPecodat)
yi <- with(NPfertgreenPeco, log(Xe/Xc))
vi <- with(NPfertgreenPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertgreenPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NPfertgreenPecodat, digits=8)  
summary(rNPfertgreenPeco)
## Fig. 3b): N fert on senesced-leaf [N]
NfertNRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertNRPeco, append=TRUE)
head(NfertNRPecodat)
yi <- with(NfertNRPeco, log(Xe/Xc))
vi <- with(NfertNRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertNRPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NfertNRPecodat, digits=8)  
summary(rNfertNRPeco)
## Fig. 3b): P fert on senesced-leaf [N]
PfertNRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertNRPeco, append=TRUE)
head(PfertNRPecodat)
yi <- with(PfertNRPeco, log(Xe/Xc))
vi <- with(PfertNRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertNRPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=PfertNRPecodat, digits=8)  
summary(rPfertNRPeco)
## Fig. 3b): N+P fert on senesced-leaf [N]
NPfertNRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertNRPeco, append=TRUE)
head(NPfertNRPecodat)
yi <- with(NPfertNRPeco, log(Xe/Xc))
vi <- with(NPfertNRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertNRPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NPfertNRPecodat, digits=8)  
summary(rNPfertNRPeco)
## Fig. 3b): N fert on senesced-leaf [P]
NfertPRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertPRPeco, append=TRUE)
head(NfertPRPecodat)
yi <- with(NfertPRPeco, log(Xe/Xc))
vi <- with(NfertPRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertPRPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NfertPRPecodat, digits=8)  
summary(rNfertPRPeco)
## Fig. 3b): P fert on senesced-leaf [P]
PfertPRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertPRPeco, append=TRUE)
head(PfertPRPecodat)
yi <- with(PfertPRPeco, log(Xe/Xc))
vi <- with(PfertPRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertPRPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=PfertPRPecodat, digits=8)  
summary(rPfertPRPeco)
## Fig. 3b): N+P fert on senesced-leaf [P]
NPfertPRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertPRPeco, append=TRUE)
head(NPfertPRPecodat)
yi <- with(NPfertPRPeco, log(Xe/Xc))
vi <- with(NPfertPRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertPRPeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NPfertPRPecodat, digits=8)  
summary(rNPfertPRPeco)
## Fig. 4b): N fert on NRE
NfertNREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertNREeco, append=TRUE)
head(NfertNREecodat)
yi <- with(NfertNREeco, log(Xe/Xc))
vi <- with(NfertNREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertNREeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NfertNREecodat, digits=8)  
summary(rNfertNREeco)
## Fig. 4b): P fert on NRE
PfertNREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertNREeco, append=TRUE)
head(PfertNREecodat)
yi <- with(PfertNREeco, log(Xe/Xc))
vi <- with(PfertNREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertNREeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=PfertNREecodat, digits=8)  
summary(rPfertNREeco)
## Fig. 4b): N+P fert on NRE
NPfertNREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertNREeco, append=TRUE)
head(NPfertNREecodat)
yi <- with(NPfertNREeco, log(Xe/Xc))
vi <- with(NPfertNREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertNREeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NPfertNREecodat, digits=8)  
summary(rNPfertNREeco)
## Fig. 4b): N fert on PRE
NfertPREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertPREeco, append=TRUE)
head(NfertPREecodat)
yi <- with(NfertPREeco, log(Xe/Xc))
vi <- with(NfertPREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNfertPREeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NfertPREecodat, digits=8)  
summary(rNfertPREeco)
## Fig. 4b): P fert on PRE
PfertPREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertPREeco, append=TRUE)
head(PfertPREecodat)
yi <- with(PfertPREeco, log(Xe/Xc))
vi <- with(PfertPREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rPfertPREeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=PfertPREecodat, digits=8)  
summary(rPfertPREeco)
## Fig. 4b): N+P fert on PRE
NPfertPREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertPREeco, append=TRUE)
head(NPfertPREecodat)
yi <- with(NPfertPREeco, log(Xe/Xc))
vi <- with(NPfertPREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rNPfertPREeco<-rma.uni(yi,vi,mods=~factor(BIOME)-1, method="REML", data=NPfertPREecodat, digits=8)  
summary(rNPfertPREeco)

##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
##### Fig. 5 ggplot responses to N&P rates
### read in data
ggplotdataNPrate<-read.csv(file.choose())
## subset data
NrateGreenN<-subset(ggplotdataNPrate, ggplotdataNPrate$TYPE=="01N-greenN")
NrateNRP<-subset(ggplotdataNPrate, ggplotdataNPrate$TYPE=="03N-NRP")
NrateNRE<-subset(ggplotdataNPrate, ggplotdataNPrate$TYPE=="05N-NRE")
PrateGreenP<-subset(ggplotdataNPrate, ggplotdataNPrate$TYPE=="02P-greenP")
PratePRP<-subset(ggplotdataNPrate, ggplotdataNPrate$TYPE=="04P-PRP")
PratePRE<-subset(ggplotdataNPrate, ggplotdataNPrate$TYPE=="06P-PRE")
NPrateRR<-rbind(NrateGreenN, NrateNRP, NrateNRE, PrateGreenP, PratePRP, PratePRE)
### explore data
head(NPrateRR)
dir()
ls()
names(NPrateRR)
dim(NPrateRR)  
qplot(rate, rate, data = NPrateRR)
qplot(rate, response, data = NPrateRR, colour = TYPE)
qplot(rate, response, data=NPrateRR, facets = . ~ TYPE) + geom_smooth()
## Nfert-[N]green
pNNg<-ggplot(NrateGreenN, aes(rate, response))
pNNg+ geom_point(colour="black", size=3.5)+
      geom_point(fill="#009e23", size=3, shape=21)+
      theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),panel.background=element_blank())+
      xlim(-10,800)+
      ylim(-1,2)+
     stat_smooth(fill="#19a781", colour="red", size=1, alpha = 0.2)+
     geom_smooth(method = 'lm', colour="blue", size=1, aes(fill = 'confidence'), alpha = 0.5)
## Nfert-[N]dead
pNNd<-ggplot(NrateNRP, aes(rate, response))
pNNd+ geom_point(colour="black", size=3.5)+
      geom_point(fill="#009e23", size=3, shape=21)+
      theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),panel.background=element_blank())+
      xlim(-10,1000)+
      ylim(-1,2)+
     stat_smooth(fill="#19a781", colour="red", size=1, alpha = 0.2)+
     geom_smooth(method = 'lm', size=1, aes(fill = 'confidence'), alpha = 0.5)
## Nfert-NRE
pNNRE<-ggplot(NrateNRE, aes(rate, response))
pNNRE+ geom_point(colour="black", size=3.5)+
      geom_point(fill="#009e23", size=3, shape=21)+
      theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),panel.background=element_blank())+
      xlim(-10,800)+
      ylim(-2.5,1)+
     stat_smooth(fill="#19a781", colour="red", size=1, alpha = 0.2)+
     geom_smooth(method = 'lm',size=1, aes(fill = 'confidence'), alpha = 0.5)
## Pfert-[P]green
pPPg<-ggplot(PrateGreenP, aes(rate, response))
pPPg+ geom_point(colour="black", size=3.5)+
      geom_point(fill="#009e23", size=3, shape=21)+
      theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),panel.background=element_blank())+
      xlim(-10,600)+
      ylim(-1.5,3)+
     stat_smooth(fill="#19a781", colour="red", size=1, alpha = 0.2)+
     geom_smooth(method = 'lm', size=1,aes(fill = 'confidence'), alpha = 0.5)
## Pfert-[P]dead
pPPd<-ggplot(PratePRP, aes(rate, response))
pPPd+ geom_point(colour="black", size=3.5)+
      geom_point(fill="#009e23", size=3, shape=21)+
      theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),panel.background=element_blank())+
      xlim(-10,150)+
      ylim(-0.5,3)+
     geom_smooth(method = 'lm', size=1,aes(fill = 'confidence'), alpha = 0.5)+
     stat_smooth(fill="#19a781", colour="red", size=1, alpha = 0.2)
## Pfert-PRE
pPPRE<-ggplot(PratePRE, aes(rate, response))
pPPRE+ geom_point(colour="black", size=3.5)+
      geom_point(fill="#009e23", size=3, shape=21)+
      theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),panel.background=element_blank())+
      xlim(-10,350)+
      ylim(-1,2)+
     geom_smooth(method = 'lm', size=1,aes(fill = 'confidence'), alpha = 0.5)+
     stat_smooth(fill="#19a781", colour="red", size=1, alpha = 0.2)
### combining figures
p<-ggplot(NPrateRR, aes(x=rate, y=response, group=TYPE, fill=TYPE)) 
p+geom_point(colour="black", size=3, shape=21)+ 
scale_fill_manual( values = c("#009e23", "#009e23", "#009e23", 
                                 "#009e23", "#009e23", "#009e23"))+
theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank(),panel.background=element_blank())+
facet_wrap(~TYPE, ncol=2, scales="free") + 
xlim(-10,150)+
ylim(-0.5,3)+
stat_smooth(fill="#19a781", colour="red", size=0.5, alpha = 0.2)

##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
### Table 1: Compare difference beteen plant groups and between biome
## between plant functional groups
####  N fert on green [N]
NfertgreenNdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertgreenN, append=TRUE)
yi <- with(NfertgreenN, log(Xe/Xc))
vi <- with(NfertgreenN, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NfertgreenNecodat)
####  P fert on green [N]
PfertgreenNdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertgreenN, append=TRUE)
yi <- with(PfertgreenN, log(Xe/Xc))
vi <- with(PfertgreenN, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = PfertgreenNecodat)
####  N+P fert on green [N]
NPfertgreenNdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertgreenN, append=TRUE)
yi <- with(NPfertgreenN, log(Xe/Xc))
vi <- with(NPfertgreenN, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NPfertgreenNecodat)
####  N fert on green [P]
NfertgreenPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertgreenP, append=TRUE)
yi <- with(NfertgreenP, log(Xe/Xc))
vi <- with(NfertgreenP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NfertgreenPecodat)
####  P fert on green [P]
PfertgreenPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertgreenP, append=TRUE)
yi <- with(PfertgreenP, log(Xe/Xc))
vi <- with(PfertgreenP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = PfertgreenPecodat)
####  N+P fert on green [P]
NPfertgreenPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertgreenP, append=TRUE)
yi <- with(NPfertgreenP, log(Xe/Xc))
vi <- with(NPfertgreenP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NPfertgreenPecodat)
####  N fert on senesced-leaf [N]
NfertNRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertNRP, append=TRUE)
yi <- with(NfertNRP, log(Xe/Xc))
vi <- with(NfertNRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NfertNRPecodat)
####  P fert on senesced-leaf [N]
PfertNRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertNRP, append=TRUE)
yi <- with(PfertNRP, log(Xe/Xc))
vi <- with(PfertNRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = PfertNRPecodat)
####  N+P fert on senesced-leaf [N]
NPfertNRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertNRP, append=TRUE)
yi <- with(NPfertNRP, log(Xe/Xc))
vi <- with(NPfertNRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NPfertNRPecodat)
####  N fert on senesced-leaf [P]
NfertPRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertPRP, append=TRUE)
yi <- with(NfertPRP, log(Xe/Xc))
vi <- with(NfertPRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NfertPRPecodat)
####  P fert on senesced-leaf [P]
PfertPRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertPRP, append=TRUE)
yi <- with(PfertPRP, log(Xe/Xc))
vi <- with(PfertPRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = PfertPRPecodat)
 ####  N+P fert on senesced-leaf [P]
NPfertPRPdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertPRP, append=TRUE)
yi <- with(NPfertPRP, log(Xe/Xc))
vi <- with(NPfertPRP, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NPfertPRPecodat)

####  N fert on NRE
NfertNREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertNRE, append=TRUE)
yi <- with(NfertNRE, log(Xe/Xc))
vi <- with(NfertNRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NfertNREecodat)
####  P fert on NRE
PfertNREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertNRE, append=TRUE)
yi <- with(PfertNRE, log(Xe/Xc))
vi <- with(PfertNRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = PfertNREecodat)
####  N+P fert on NRE
NPfertNREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertNRE, append=TRUE)
yi <- with(NPfertNRE, log(Xe/Xc))
vi <- with(NPfertNRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NPfertNREecodat)

####  N fert on PRE
NfertPREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertPRE, append=TRUE)
yi <- with(NfertPRE, log(Xe/Xc))
vi <- with(NfertPRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NfertPREecodat)
####  P fert on PRE
PfertPREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertPRE, append=TRUE)
yi <- with(PfertPRE, log(Xe/Xc))
vi <- with(PfertPRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = PfertPREecodat)
####  N+P fert on PRE
NPfertPREdat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertPRE, append=TRUE)
yi <- with(NPfertPRE, log(Xe/Xc))
vi <- with(NPfertPRE, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(PFT), data = NPfertPREecodat)

#### between ecosystems
##  N fert on green-leaf [N]
NfertgreenNecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertgreenNeco, append=TRUE)
yi <- with(NfertgreenNeco, log(Xe/Xc))
vi <- with(NfertgreenNeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NfertgreenNecodat)
##  P fert on green-leaf [N]
PfertgreenNecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertgreenNeco, append=TRUE)
yi <- with(PfertgreenNeco, log(Xe/Xc))
vi <- with(PfertgreenNeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = PfertgreenNecodat)
##  N+P fert on green-leaf [N]
NPfertgreenNecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertgreenNeco, append=TRUE)
yi <- with(NPfertgreenNeco, log(Xe/Xc))
vi <- with(NPfertgreenNeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NPfertgreenNecodat)
##  N fert on green-leaf [P]
NfertgreenPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertgreenPeco, append=TRUE)
yi <- with(NfertgreenPeco, log(Xe/Xc))
vi <- with(NfertgreenPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NfertgreenPecodat)
##  P fert on green-leaf [P]
PfertgreenPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertgreenPeco, append=TRUE)
yi <- with(PfertgreenPeco, log(Xe/Xc))
vi <- with(PfertgreenPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = PfertgreenPecodat)
##  N+P fert on green-leaf [P]
NPfertgreenPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertgreenPeco, append=TRUE)
yi <- with(NPfertgreenPeco, log(Xe/Xc))
vi <- with(NPfertgreenPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NPfertgreenPecodat)
##  N fert on senesced-leaf [N]
NfertNRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertNRPeco, append=TRUE)
yi <- with(NfertNRPeco, log(Xe/Xc))
vi <- with(NfertNRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NfertNRPecodat)
##  P fert on senesced-leaf [N]
PfertNRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertNRPeco, append=TRUE)
yi <- with(PfertNRPeco, log(Xe/Xc))
vi <- with(PfertNRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = PfertNRPecodat)
##  N+P fert on senesced-leaf [N]
NPfertNRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertNRPeco, append=TRUE)
yi <- with(NPfertNRPeco, log(Xe/Xc))
vi <- with(NPfertNRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NPfertNRPecodat)
##  N fert on senesced-leaf [P]
NfertPRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertPRPeco, append=TRUE)
yi <- with(NfertPRPeco, log(Xe/Xc))
vi <- with(NfertPRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NfertPRPecodat)
##  P fert on senesced-leaf [P]
PfertPRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertPRPeco, append=TRUE)
yi <- with(PfertPRPeco, log(Xe/Xc))
vi <- with(PfertPRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = PfertPRPecodat)
##  N+P fert on senesced-leaf [P]
NPfertPRPecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertPRPeco, append=TRUE)
yi <- with(NPfertPRPeco, log(Xe/Xc))
vi <- with(NPfertPRPeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NPfertPRPecodat)
##  N fert on NRE
NfertNREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertNREeco, append=TRUE)
yi <- with(NfertNREeco, log(Xe/Xc))
vi <- with(NfertNREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NfertNREecodat)
##  P fert on NRE
PfertNREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertNREeco, append=TRUE)
yi <- with(PfertNREeco, log(Xe/Xc))
vi <- with(PfertNREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = PfertNREecodat)
##  N+P fert on NRE
NPfertNREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertNREeco, append=TRUE)
yi <- with(NPfertNREeco, log(Xe/Xc))
vi <- with(NPfertNREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NPfertNREecodat)
##  N fert on PRE
NfertPREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NfertPREeco, append=TRUE)
yi <- with(NfertPREeco, log(Xe/Xc))
vi <- with(NfertPREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NfertPREecodat)
##  P fert on PRE
PfertPREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=PfertPREeco, append=TRUE)
yi <- with(PfertPREeco, log(Xe/Xc))
vi <- with(PfertPREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = PfertPREecodat)
##  N+P fert on PRE
NPfertPREecodat <- escalc(m1i=Xe, sd1i=Se, n1i=Ne, m2i=Xc, sd2i=Sc, n2i=Nc, measure="ROM", data=NPfertPREeco, append=TRUE)
yi <- with(NPfertPREeco, log(Xe/Xc))
vi <- with(NPfertPREeco, (1/Ne)*(Se/Xe)^2 + (1/Nc)*(Sc/Xc)^2)
rma(yi,vi, mods = ~ factor(BIOME), data = NPfertPREecodat)
##%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The End