library(metafor)
library(ggplot2)
library(glmulti)

random1<-rma(yi,vi, data=d11, method="REML")
summary(random1)
###Biochar application rates(t/ha):
r2<-rma(yi,vi, data=subset(d11,Dosage=="M<=10"), method="REML")
r3<-rma(yi,vi, data=subset(d11,Dosage=="10<M<=40"), method="REML")
r4<-rma(yi,vi, data=subset(d11,Dosage=="40<M<=80"), method="REML")
r5<-rma(yi,vi, data=subset(d11,Dosage=="M>80"), method="REML")

summary(r2)
summary(r3)
summary(r4)
summary(r5)


confint(r2)
confint(r3)
confint(r4)
confint(r5)

estimates2 <- c(coef(r2), coef(r3),coef(r4), coef(r5))
variances2 <- c(vcov(r2), vcov(r3),vcov(r4), vcov(r5))
labels2 <- c("M<=10", "10<M<=40","40<M<=80","M>80")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 7),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.111,0.305), 5.5, c("Application rates(t/ha)", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.04), 5.5, c("(P<0.0001)"),cex=1.3,font=7)

r41<-rma(yi,vi, mods=~Dosage,data=d11, method="REML")
summary(r41)
r432<-rma(yi,vi, mods=~Dosage-1,data=d11, method="REML")
summary(r432)
anova(r432, L=c(1,-1,0,0))
anova(r432,L=c(1,0,-1,0))
anova(r432,L=c(0,0,1,-1))
anova(r432,L=c(1,0,0,-1))
anova(r432,L=c(0,1,-1,0))
anova(r432,L=c(0,1,0,-1))



####Pyrolysis temperature:
r10<-rma(yi,vi, data=subset(d11,temperature=="<=400"), method="REML")
r11<-rma(yi,vi, data=subset(d11,temperature=="401~500"), method="REML")
r12<-rma(yi,vi, data=subset(d11,temperature=="501~600"), method="REML")
r13<-rma(yi,vi, data=subset(d11,temperature==">600"), method="REML")

summary(r10)
summary(r11)
summary(r12)
summary(r13)


confint(r10)
confint(r11)
confint(r12)
confint(r13)


estimates2 <- c(coef(r10), coef(r11),coef(r12), coef(r13))
variances2 <- c(vcov(r10), vcov(r11),vcov(r12), vcov(r13) )
labels2 <- c("<=400", "401~500","501~600",">600")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 7),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.020,0.152), 5.4, c("Pyrolysis temperature", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.04), 5.4, c("(P=0.8093)"),cex=1.3,font=7)
r42<-rma(yi,vi, mods=~temperature,data=d11, method="REML")
summary(r42)


r28<-rma(yi,vi, mods=~temperature-1,data=d11, method="REML")
summary(r28)
anova(r28, L=c(1,-1,0,0))
anova(r28,L=c(1,0,-1,0))
anova(r28,L=c(0,0,1,-1))
anova(r28,L=c(1,0,0,-1))
anova(r28,L=c(0,1,-1,0))
anova(r28,L=c(0,1,0,-1))







####Biochar.species:
r101<-rma(yi,vi, data=subset(d11,Biochar.species=="Crop straw"), method="REML")
r111<-rma(yi,vi, data=subset(d11,Biochar.species=="Wood group"), method="REML")
r121<-rma(yi,vi, data=subset(d11,Biochar.species=="Shell-slag group"), method="REML")
r131<-rma(yi,vi, data=subset(d11,Biochar.species=="Animal waste"), method="REML")

summary(r101)
summary(r111)
summary(r121)
summary(r131)


confint(r101)
confint(r111)
confint(r121)
confint(r131)

estimates2 <- c(coef(r101), coef(r111),coef(r121), coef(r131))
variances2 <- c(vcov(r101), vcov(r111),vcov(r121), vcov(r131) )
labels2 <- c("Straw crop", "Wood group","Shell-slag group","Animal waste")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 6.5),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.04,0.178), 4.8, c("Biochar type", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.01), 4.8, c("(p=0.0005)"),cex=1.3,font=7)

r43<-rma(yi,vi, mods=~Biochar.species,data=d11, method="REML")
summary(r43)

r431<-rma(yi,vi, mods=~Biochar.species-1,data=d11, method="REML")
summary(r431)
anova(r431, L=c(0,1,-1,0,0))
anova(r431,L=c(0,1,0,-1,0))
anova(r431,L=c(0,0,0,1,-1))
anova(r431,L=c(0,1,0,0,-1))
anova(r431,L=c(0,0,1,-1,0))
anova(r431,L=c(0,0,1,0,-1))








####Biochar pH:
r14<-rma(yi,vi, data=subset(d11,amendments.pH.1=="pH<=8"), method="REML")
r15<-rma(yi,vi, data=subset(d11,amendments.pH.1=="8<pH<=10"), method="REML")
r16<-rma(yi,vi, data=subset(d11,amendments.pH.1=="pH>10"), method="REML")

summary(r14)
summary(r15)
summary(r16)



confint(r14)
confint(r15)
confint(r16)

estimates2 <- c(coef(r14), coef(r15),coef(r16))
variances2 <- c(vcov(r14), vcov(r15),vcov(r16))
labels2 <- c("pH<=8", "8<pH<=10","pH>10")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 5.5),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.07,0.21), 3.8, c("Biochar pH", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.-0.011), 3.8, c("(p=0.0006)"),cex=1.3,font=7)

r431<-rma(yi,vi, mods=~amendments.pH.1-1,data=d11, method="REML")
summary(r431)

r43<-rma(yi,vi, mods=~amendments.pH.1,data=d11, method="REML")
summary(r43)

anova(r431, L=c(1,-1,0,0))
anova(r431,L=c(1,0,-1,0))
anova(r431,L=c(0,0,1,-1))
anova(r431,L=c(1,0,0,-1))
anova(r431,L=c(0,1,-1,0))
anova(r431,L=c(0,1,0,-1))









#####Initial soil pH:
r14<-rma(yi,vi, data=subset(d11,Soil.pH.1=="pH<=4.5"), method="REML")
r15<-rma(yi,vi, data=subset(d11,Soil.pH.1=="4.5<pH<=5.5"), method="REML")
r16<-rma(yi,vi, data=subset(d11,Soil.pH.1=="5.5<pH<=6.5"), method="REML")


summary(r14)
summary(r15)
summary(r16)



confint(r14)
confint(r15)
confint(r16)


estimates2 <- c(coef(r14), coef(r15),coef(r16))
variances2 <- c(vcov(r14), vcov(r15),vcov(r16))
labels2 <- c("pH<=4.5", "4.5<pH<=5.5","5.5<pH<=6.5")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 5.5),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.012,0.179), 3.8, c("Initial soil pH", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.034), 3.8, c("(p=0.5969)"),cex=1.3,font=7)
r7<-rma(yi,vi, mods=~Soil.pH.1,data=d11, method="REML")
summary(r7)
r8<-rma(yi,vi, mods=~Soil.pH.1-1,data=d11, method="REML")
summary(r8)

anova(r8, L=c(0,1,-1,0))
anova(r8, L=c(0,1,0,-1))
anova(r8, L=c(0,0,1,-1))












####Experiment type:
r17<-rma(yi,vi, data=subset(d11,experiment.type=="field experiment"), method="REML")
r18<-rma(yi,vi, data=subset(d11,experiment.type=="pot experiment"), method="REML")
r19<-rma(yi,vi, data=subset(d11,experiment.type=="Incubation experiment"), method="REML")


summary(r17)
summary(r18)
summary(r19)



confint(r17)
confint(r18)
confint(r19)


estimates2 <- c(coef(r17), coef(r18),coef(r19))
variances2 <- c(vcov(r17), vcov(r18),vcov(r19))
labels2 <- c("field", "pot","incubation")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 5.5),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.02,0.162), 3.8, c("Experiment type", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.03), 3.8, c("(p<0.0001)"),cex=1.3,font=7)


####soil pH (lime):
r14<-rma(yi,vi, data=subset(d31,Soil.pH.1=="pH<=4.5"), method="REML")
r15<-rma(yi,vi, data=subset(d31,Soil.pH.1=="4.5<pH<=5.5"), method="REML")



summary(r14)
summary(r15)




confint(r14)
confint(r15)




estimates2 <- c(coef(r14), coef(r15))
variances2 <- c(vcov(r14), vcov(r15))
labels2 <- c("pH<=4.5", "4.5<pH<=5.5")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 4.5),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(0.0483,0.268), 2.8, c("Initial soil pH", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.10), 2.8, c("(p<0.0001)"),cex=1.3,font=7)

r7<-rma(yi,vi, mods=~Soil.pH.1,data=d31, method="REML")
summary(r7)
r8<-rma(yi,vi, mods=~Soil.pH.1-1,data=d31, method="REML")
summary(r8)

anova(r8,L=c(1,-1))
















####Amendments typle:
r29<-rma(yi,vi, data=subset(d21,amendments=="Biochar"), method="REML")
r39<-rma(yi,vi, data=subset(d21,amendments=="Lime"), method="REML")
r49<-rma(yi,vi, data=subset(d21,amendments=="Crop straw"), method="REML")
r59<-rma(yi,vi, data=subset(d21,amendments=="Rapeseed cake"), method="REML")
r69<-rma(yi,vi, data=subset(d21,amendments=="Animal waste"), method="REML")

summary(r29)
summary(r39)
summary(r49)
summary(r59)
summary(r69)

confint(r29)
confint(r39)
confint(r49)
confint(r59)
confint(r69)
estimates2 <- c(coef(r29), coef(r39),coef(r49), coef(r59),coef(r69))
variances2 <- c(vcov(r29), vcov(r39),vcov(r49), vcov(r59),vcov(r69))
labels2 <- c("biochar", "lime","crop straw","rapeseed cake","animal waste")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 7.3),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.16827,0.3015), 5.6, c("Amendments typle", "Mean [95% CI]"),cex=1.3,font=7)
text(c(-0.03), 5.6, c("(P<0.0001)"),cex=1.3,font=7)

r41<-rma(yi,vi, mods=~amendments,data=d21, method="REML")
summary(r41)
r432<-rma(yi,vi, mods=~amendments-1,data=d21, method="REML")
summary(r432)
anova(r432, L=c(1,0,0,-1,0))
anova(r432,L=c(1,0,-1,0,0))
anova(r432,L=c(0,0,0,1,-1))
anova(r432,L=c(1,0,0,0,-1))
anova(r432,L=c(0,1,-1,0,0))
anova(r432,L=c(0,1,0,0,-1))





####Lime application rates:
r2<-rma(yi,vi, data=subset(d31,Dosage=="M<=1.5"), method="REML")
r3<-rma(yi,vi, data=subset(d31,Dosage=="1.5<M<=4.5"), method="REML")
r4<-rma(yi,vi, data=subset(d31,Dosage=="M>4.5"), method="REML")

summary(r2)
summary(r3)
summary(r4)



confint(r2)
confint(r3)
confint(r4)


estimates2 <- c(coef(r2), coef(r3),coef(r4))
variances2 <- c(vcov(r2), vcov(r3),vcov(r4))
labels2 <- c("M<=1.5", "1.5<M<=4.5","M>4.5")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 5.6),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.111,0.423), 4, c("Application rates(t/ha)", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.08), 4, c("(P<0.0001)"),cex=1.3,font=7)

r41<-rma(yi,vi, mods=~Dosage,data=d31, method="REML")
summary(r41)
r432<-rma(yi,vi, mods=~Dosage-1,data=d31, method="REML")
summary(r432)
anova(r432, L=c(1,-1,0))
anova(r432,L=c(1,0,-1))
anova(r432,L=c(0,1,-1))





#####Experiment type(lime)：
r17<-rma(yi,vi, data=subset(d31,experiment.type=="field experiment"), method="REML")

r19<-rma(yi,vi, data=subset(d31,experiment.type=="Incubation experiment"), method="REML")


summary(r17)

summary(r19)



confint(r17)

confint(r19)


estimates2 <- c(coef(r17),coef(r19))
variances2 <- c(vcov(r17),vcov(r19))
labels2 <- c("field", "incubation")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 4.5),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.05,0.392), 2.8, c("Experiment type", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.07), 2.8, c("(p=0.0038)"),cex=1.3,font=7)

r46<-rma(yi,vi, mods=~experiment.type,data=d31, method="REML")
summary(r46)




####Application rates(crop straw)
r2<-rma(yi,vi, data=subset(d41,Dosage=="M<=20"), method="REML")
r3<-rma(yi,vi, data=subset(d41,Dosage=="20<M<=40"), method="REML")
r4<-rma(yi,vi, data=subset(d41,Dosage=="M>40"), method="REML")

summary(r2)
summary(r3)
summary(r4)



confint(r2)
confint(r3)
confint(r4)

estimates2 <- c(coef(r2), coef(r3),coef(r4))
variances2 <- c(vcov(r2), vcov(r3),vcov(r4))
labels2 <- c("M<=20", "20<M<=40","M>40")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 5.6),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.031,0.125), 4, c("Application rates(t/ha)", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.023), 4, c("(P=0.1746)"),cex=1.3,font=7)

r41<-rma(yi,vi, mods=~Dosage,data=d41, method="REML")
summary(r41)
r432<-rma(yi,vi, mods=~Dosage-1,data=d41, method="REML")
summary(r432)
anova(r432, L=c(1,-1,0))
anova(r432,L=c(1,0,-1))
anova(r432,L=c(0,1,-1))



###Application rates(rapeseed cake)
r2<-rma(yi,vi, data=subset(d51,Dosage=="M<=4"), method="REML")
r3<-rma(yi,vi, data=subset(d51,Dosage=="4<M<=10"), method="REML")
r4<-rma(yi,vi, data=subset(d51,Dosage=="M>10"), method="REML")

summary(r2)
summary(r3)
summary(r4)



confint(r2)
confint(r3)
confint(r4)


estimates2 <- c(coef(r2), coef(r3),coef(r4))
variances2 <- c(vcov(r2), vcov(r3),vcov(r4))
labels2 <- c("M<=4", "4<M<=10","M>10")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 5.6),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.159,0.1051), 4, c("Application rates(t/ha)", "Mean [95% CI]"),cex=1.3,font=7)
text(c(-0.066), 4, c("(P=0.2347)"),cex=1.3,font=7)

r41<-rma(yi,vi, mods=~Dosage,data=d51, method="REML")
summary(r41)
r432<-rma(yi,vi, mods=~Dosage-1,data=d51, method="REML")
summary(r432)
anova(r432, L=c(1,-1,0))
anova(r432,L=c(1,0,-1))
anova(r432,L=c(0,1,-1))





####Animal excrement application rates
r2<-rma(yi,vi, data=subset(d61,Dosage=="M<=10"), method="REML")
r3<-rma(yi,vi, data=subset(d61,Dosage=="10<M<=30"), method="REML")
r4<-rma(yi,vi, data=subset(d61,Dosage=="M>30"), method="REML")

summary(r2)
summary(r3)
summary(r4)



confint(r2)
confint(r3)
confint(r4)


estimates2 <- c(coef(r2), coef(r3),coef(r4))
variances2 <- c(vcov(r2), vcov(r3),vcov(r4))
labels2 <- c("M<=10", "10<M<=30","M>30")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 5.6),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.030,0.1851), 4, c("Application rates(t/ha)", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.045), 4, c("(P=0.9412)"),cex=1.3,font=7)

r41<-rma(yi,vi, mods=~Dosage,data=d61, method="REML")
summary(r41)
r432<-rma(yi,vi, mods=~Dosage-1,data=d61, method="REML")
summary(r432)
anova(r432, L=c(1,-1,0))
anova(r432,L=c(1,0,-1))
anova(r432,L=c(0,1,-1))



####animal waste pH:

r15<-rma(yi,vi, data=subset(d61,amendments.Ph.1=="6<pH<=8"), method="REML")
r16<-rma(yi,vi, data=subset(d61,amendments.Ph.1=="pH>8"), method="REML")



summary(r15)
summary(r16)




confint(r15)
confint(r16)


estimates2 <- c( coef(r15),coef(r16))
variances2 <- c( vcov(r15),vcov(r16))
labels2 <- c("6<pH<=8","pH>8")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 4.5),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.03,0.245), 2.8, c("Animal waste pH", "Mean [95% CI]"),cex=1.3,font=7)
text(c(0.0465), 2.8, c("(p=0.0531)"),cex=1.3,font=7)

r431<-rma(yi,vi, mods=~amendments.Ph.1-1,data=d61, method="REML")
summary(r431)

r43<-rma(yi,vi, mods=~amendments.Ph.1,data=d61, method="REML")
summary(r43)

anova(r431, L=c(1,-1,0,0))
anova(r431,L=c(1,0,-1,0))
anova(r431,L=c(0,0,1,-1))
anova(r431,L=c(1,0,0,-1))
anova(r431,L=c(0,1,-1,0))
anova(r431,L=c(0,1,0,-1))


####Initial soil pH (animal waste)
r14<-rma(yi,vi, data=subset(d61,Soil.pH.1=="pH<=4.5"), method="REML")
r15<-rma(yi,vi, data=subset(d61,Soil.pH.1=="4.5<pH<=5.5"), method="REML")
r16<-rma(yi,vi, data=subset(d61,Soil.pH.1=="5.5<pH<=6.5"), method="REML")

summary(r14)
summary(r15)
summary(r16)



confint(r14)
confint(r15)
confint(r16)


estimates2 <- c(coef(r14), coef(r15),coef(r16))
variances2 <- c(vcov(r14), vcov(r15),vcov(r16))
labels2 <- c("pH<=4.5", "4.5<pH<=5.5","5.5<pH<=6.5")
forest(estimates2, variances2,slab=labels2,ylim=c(0.8, 5.5),font=7,cex=1,psize=1,xlab="Effect size (log response ratio)")
text(c(-0.262,0.503), 3.8, c("Initial soil pH", "Mean [95% CI]"),cex=1.3,font=7)
text(c(-0.073), 3.8, c("(p=0.0069)"),cex=1.3,font=7)


r7<-rma(yi,vi, mods=~Soil.pH.1,data=d61, method="REML")
summary(r7)
r8<-rma(yi,vi, mods=~Soil.pH.1-1,data=d61, method="REML")
summary(r8)

anova(r8, L=c(0,1,-1,0))
anova(r8, L=c(0,1,0,-1))
anova(r8, L=c(0,0,1,-1))