### needs Script_upload_1200_fixed to be run before

all_D$col<-ifelse(all_D$group=="fibrous","black","darkgrey") #for coloring symbols
all_D$overlap<-all_D$overlap*100 #transform to %
##### Errors diam ####


tiff(file="Errors_diam.tif",height=13,width=19.5, unit='in',res=300)
par(mfrow= c(2,3), 
    bty="l", 
    oma=c(5,5,0,0),# outer margins c(bottom, left, top, right)
    mar=c(1.5,1.5,0,0.1), #c(bottom, left, top, right)
    mgp=c(4, 1, 0), #axes,labels,ticks
    xpd = NA)
range(all_D$d_d_WR)
range(all_D$d_d_IJC)
##### rel WRhizo ~ diam ####
all_D<-all_D[order(-all_D$diam) , ]
plot(abs(all_D$d_d_WR*100)~all_D$diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="WinRhizo",
     #xlim=Xlim,
     ylim=c(0,16),
     tck=0.02,
     xlab="",
     xaxt="n",
     ylab=expression(Relative~error~diameter~WR~~"("~"%"~")")
)  
Axis(side = 1, labels = FALSE, tck = 0.02)
mod1<-lm(abs(all_D$d_d_WR*100)~all_D$diam)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$diam), max(all_D$diam), length.out=length(all_D$diam))
  preds <- predict(mod1, newdata = data.frame(diam=newx), 
                   interval = 'confidence')
  lines(all_D$diam, preds[ ,3], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  lines(all_D$diam, preds[ ,2], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  polygon(c(rev(all_D$diam), all_D$diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topleft",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                              paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                              ""))
}
legend(x="topleft", 
       c("fibrous","taproot"), 
       pch=c(21,21), lty=c(2,1),pt.cex=3,
       pt.bg = c("black","lightgrey"),
       col=c("white","white"), 
       bty="n",cex=1.5)

####rel WRhizo~overlap ####
all_D<-all_D[order(all_D$overlap) , ]
plot(abs(all_D$d_d_WR*100)~all_D$overlap,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="WinRhizo",
     #xlim=Xlim,
     ylim=c(0,16),
     tck=0.02,
     yaxt="n",
     xaxt="n",
     xlab="",
     ylab=""
)  
Axis(side = 1, labels = FALSE, tck = 0.02)
Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_d_WR*100)~all_D$overlap)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$overlap), max(all_D$overlap), length.out=length(all_D$overlap))
  preds <- predict(mod1, newdata = data.frame(overlap=newx), 
                   interval = 'confidence')
  lines(all_D$overlap, preds[ ,3], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  lines(all_D$overlap, preds[ ,2], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  polygon(c(rev(all_D$overlap), all_D$overlap), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$overlap, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topleft",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                              paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                              ""))
}

####rel WRhizo~CV_D ####
all_D<-all_D[order(all_D$CV_diam) , ]
plot(abs(all_D$d_d_WR*100)~all_D$CV_diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="WinRhizo",
     #xlim=Xlim,
     ylim=c(0,16),
     tck=0.02,
     yaxt="n",
     xaxt="n",
     xlab="",
     ylab=""
)  
Axis(side = 1, labels = FALSE, tck = 0.02)
Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_d_WR*100)~all_D$CV_diam)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$CV_diam,na.rm=T), max(all_D$CV_diam,na.rm=T), length.out=length(all_D$overlap))
  preds <- predict(mod1, newdata = data.frame(CV_diam=newx), 
                   interval = 'confidence')
  lines(all_D$CV_diam, preds[ ,3], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  lines(all_D$CV_diam, preds[ ,2], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  polygon(c(rev(all_D$CV_diam), all_D$CV_diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$CV_diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topleft",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                              paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                              ""))
}
#### rel IJ ~ diam ####
all_D<-all_D[order(-all_D$diam) , ]
plot(abs(all_D$d_d_IJC*100)~all_D$diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,16),
     tck=0.02,
     xlab=expression(Root~diameter~~"("~mm~")"),
     ylab=expression(Relative~error~diameter~IJ_Rhizo~~"("~"%"~")")
)  
Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_d_IJC*100)~all_D$diam)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$diam), max(all_D$diam), length.out=length(all_D$diam))
  preds <- predict(mod1, newdata = data.frame(diam=newx), 
                   interval = 'confidence')
  lines(all_D$diam, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$diam, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$diam), all_D$diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topleft",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                              paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                              ""))
}

##### rel IJ~ overlap ####
all_D<-all_D[order(all_D$overlap) , ]
plot(abs(all_D$d_d_IJC*100)~all_D$overlap,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="WinRhizo",
     #xlim=Xlim,
     ylim=c(0,16),
     tck=0.02,
     yaxt="n",
     xlab=expression(Root~overlap~~"("~"%"~")"),
     ylab=""
)  

Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_d_IJC*100)~all_D$overlap)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$overlap), max(all_D$overlap), length.out=length(all_D$overlap))
  preds <- predict(mod1, newdata = data.frame(overlap=newx), 
                   interval = 'confidence')
  lines(all_D$overlap, preds[ ,3], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  lines(all_D$overlap, preds[ ,2], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  polygon(c(rev(all_D$overlap), all_D$overlap), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$overlap, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topleft",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                              paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                              ""))
}

##### rel IJ~ CV_D ####
all_D<-all_D[order(all_D$CV_diam) , ]
plot(abs(all_D$d_d_IJC*100)~all_D$CV_diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="WinRhizo",
     #xlim=Xlim,
     ylim=c(0,16),
     tck=0.02,
     yaxt="n",
     xlab=expression(CV~root~diameter~~"("~"%"~")"),
     ylab=""
)  

Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_d_IJC*100)~all_D$CV_diam)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$CV_diam,na.rm=T), max(all_D$CV_diam,na.rm=T), length.out=length(all_D$CV_diam))
  preds <- predict(mod1, newdata = data.frame(overlap=newx), 
                   interval = 'confidence')
  lines(all_D$CV_diam, preds[ ,3], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  lines(all_D$CV_diam, preds[ ,2], lty = 'dashed',lwd=2, col = 'darkgrey') # intervals
  polygon(c(rev(all_D$CV_diam), all_D$CV_diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$CV_diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topleft",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                          paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                          ""))
}

dev.off()

#### Errors length ####
tiff(file="Errors_length.tif",height=13,width=19.5, unit='in',res=300)

par(mfrow= c(2,3), 
    bty="l", 
    oma=c(5,5,0,0),# outer margins c(bottom, left, top, right)
    mar=c(1.5,1.5,0,0.1), #c(bottom, left, top, right)
    mgp=c(4, 1, 0), #axes,labels,ticks
    xpd = NA)

#### WR ####
all_D<-all_D[order(all_D$Sum_length) , ]
plot(abs(all_D$d_l_WR*100)~all_D$Sum_length,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(-0,32),
     tck=0.02,
     xaxt="n",
     xlab="",
     ylab=expression(Relative~error~length~WR~~"("~"%"~")")
)  
Axis(side = 1, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_l_WR*100)~(all_D$Sum_length))
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$Sum_length), max(all_D$Sum_length), length.out=length(all_D$Sum_length))
  preds <- predict(mod1, newdata = data.frame(Sum_length=newx), 
                   interval = 'confidence')
  lines(all_D$Sum_length, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$Sum_length, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$Sum_length), all_D$Sum_length), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$Sum_length, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}
legend(x="topleft", 
       c("fibrous","taproot"), 
       pch=c(21,21), lty=c(2,1),pt.cex=3,
       pt.bg = c("black","lightgrey"),
       col=c("white","white"), 
       bty="n",cex=1.5)

all_D<-all_D[order(all_D$density) , ]
#rel WR
plot(abs(all_D$d_l_WR*100)~all_D$density,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,32),
     tck=0.02,
     xaxt="n",
     yaxt="n",
     xlab="",
     ylab=""
)  
Axis(side = 1, labels = FALSE, tck = 0.02)
Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_l_WR*100)~all_D$density)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$density), max(all_D$density), length.out=length(all_D$density))
  preds <- predict(mod1, newdata = data.frame(density=newx), 
                   interval = 'confidence')
  lines(all_D$density, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$density, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$density), all_D$density), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$density, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}

all_D<-all_D[order(all_D$overlap) , ]
#rel WR
plot(abs(all_D$d_l_WR*100)~all_D$overlap,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,32),
     tck=0.02,
     xaxt="n",
     yaxt="n",
     xlab="",
     ylab=""
)  
Axis(side = 1, labels = FALSE, tck = 0.02)
Axis(side = 2, labels = FALSE, tck = 0.02)


mod1<-lm(abs(all_D$d_l_WR*100)~all_D$overlap)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$overlap), max(all_D$overlap), length.out=length(all_D$density))
  preds <- predict(mod1, newdata = data.frame(overlap=newx), 
                   interval = 'confidence')
  lines(all_D$overlap, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$overlap, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$overlap), all_D$overlap), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$overlap, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}
#### IJ ####
all_D<-all_D[order(all_D$Sum_length) , ]
plot(abs(all_D$d_l_IJ*100)~all_D$Sum_length,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(-0,32),
     tck=0.02,
     xlab=expression(Root~length~~"("~cm~")"),
     ylab=expression(Relative~error~length~IJ_Rhizo~~"("~"%"~")")
)  


mod1<-lm(abs(all_D$d_l_IJ*100)~(all_D$Sum_length))
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$Sum_length), max(all_D$Sum_length), length.out=length(all_D$Sum_length))
  preds <- predict(mod1, newdata = data.frame(Sum_length=newx), 
                   interval = 'confidence')
  lines(all_D$Sum_length, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$Sum_length, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$Sum_length), all_D$Sum_length), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$Sum_length, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}

#rel IJ
all_D<-all_D[order(all_D$density) , ]
plot(abs(all_D$d_l_IJ*100)~all_D$density,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,32),
     tck=0.02,
     yaxt="n",
     ylab="",
     xlab=expression(Root~length~density~~"("~cm~cm^{-2}~")")
)  
Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_l_IJ*100)~all_D$density)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$density), max(all_D$density), length.out=length(all_D$density))
  preds <- predict(mod1, newdata = data.frame(density=newx), 
                   interval = 'confidence')
  lines(all_D$density, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$density, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$density), all_D$density), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$density, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                              paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                              ""))
}

all_D<-all_D[order(all_D$overlap) , ]
#rel WR
plot(abs(all_D$d_l_IJ*100)~all_D$overlap,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,32),
     tck=0.02,
     yaxt="n",
     xlab=expression(Root~overlap~~"("~"%"~")"),
     ylab=""
) 

Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_l_IJ*100)~all_D$overlap)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$overlap), max(all_D$overlap), length.out=length(all_D$density))
  preds <- predict(mod1, newdata = data.frame(overlap=newx), 
                   interval = 'confidence')
  lines(all_D$overlap, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$overlap, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$overlap), all_D$overlap), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$overlap, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}

dev.off()

#### Errors volume ####
# log relationship for Vsum ~ V

tiff(file="Errors_volume.tif",height=19.5,width=19.5, unit='in',res=300)
par(mfrow= c(3,3), 
    bty="l", 
    oma=c(5,5,0,0),# outer margins c(bottom, left, top, right)
    mar=c(1.5,1.5,0,0.1), #c(bottom, left, top, right)
    mgp=c(4, 1, 0), #axes,labels,ticks
    xpd = NA)
#### IJ Rhizo ####
all_D<-all_D[order(all_D$vol) , ]
plot(abs(all_D$d_V_IJ*100)~all_D$vol,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     xaxt="n",
     xlab="",
     ylab=expression(Relative~error~volume~IJ_Rhizo~~"("~"%"~")")
)  

mod1<-lm(abs(all_D$d_V_IJ*100)~all_D$vol)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$vol), max(all_D$vol), length.out=length(all_D$vol))
  preds <- predict(mod1, newdata = data.frame(vol=newx), 
                   interval = 'confidence')
  lines(all_D$vol, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$vol, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$vol), all_D$vol), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$vol, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}

legend(x="topleft", 
       c("fibrous","taproot"), 
       pch=c(21,21), lty=c(2,1),pt.cex=3,
       pt.bg = c("black","lightgrey"),
       col=c("white","white"), 
       bty="n",cex=1.5)

all_D<-all_D[order(all_D$diam) , ]
plot(abs(all_D$d_V_IJ*100)~all_D$diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     yaxt="n",
     xaxt="n",
     xlab="",
     ylab=""
)  
Axis(side = 2, labels = FALSE, tck = 0.02)
Axis(side = 1, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_V_IJ*100)~log(all_D$diam))
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$diam), max(all_D$diam), length.out=length(all_D$diam))
  preds <- predict(mod1, newdata = data.frame(diam=newx), 
                   interval = 'confidence')
  lines(all_D$diam, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$diam, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$diam), all_D$diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}

#rel WR CVRD
all_D<-all_D[order(all_D$CV_diam) , ]
plot(abs(all_D$d_V_IJ*100)~all_D$CV_diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     yaxt="n",
     xaxt="n",
     xlab="",
     ylab=""
)  
Axis(side = 2, labels = FALSE, tck = 0.02)
Axis(side = 1, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_V_IJ*100)~all_D$CV_diam)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$CV_diam, na.rm=T), max(all_D$CV_diam, na.rm=T), length.out=length(all_D$CV_diam))
  preds <- predict(mod1, newdata = data.frame(CV_diam=newx), 
                   interval = 'confidence')
  lines(all_D$CV_diam, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$CV_diam, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$CV_diam), all_D$CV_diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$CV_diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                              paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                              ""))
}

#### WR ####
all_D<-all_D[order(all_D$vol) , ]
plot(abs(all_D$d_V_WR*100)~all_D$vol,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     xaxt="n",
     xlab="",
     ylab=expression(Relative~error~volume~WR~"("~"%"~")")
) 

Axis(side = 1, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_V_WR*100)~all_D$vol)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$vol), max(all_D$vol), length.out=length(all_D$vol))
  preds <- predict(mod1, newdata = data.frame(vol=newx), 
                   interval = 'confidence')
  lines(all_D$vol, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$vol, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$vol), all_D$vol), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$vol, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}

all_D<-all_D[order(all_D$diam) , ]
plot(abs(all_D$d_V_WR*100)~all_D$diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     yaxt="n",
     yaxt="n",
     xaxt="n",
     xlab="",
     ylab=""
)  
Axis(side = 2, labels = FALSE, tck = 0.02)
Axis(side = 1, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_V_WR*100)~log(all_D$diam))
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$diam), max(all_D$diam), length.out=length(all_D$diam))
  preds <- predict(mod1, newdata = data.frame(diam=newx), 
                   interval = 'confidence')
  lines(all_D$diam, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$diam, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$diam), all_D$diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}

#rel WR CVRD
all_D<-all_D[order(all_D$CV_diam) , ]
plot(abs(all_D$d_V_WR*100)~all_D$CV_diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     yaxt="n",
     xaxt="n",
     xlab="",
     ylab=""
)  
Axis(side = 2, labels = FALSE, tck = 0.02)
Axis(side = 1, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_V_WR*100)~all_D$CV_diam)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$CV_diam,na.rm=T), max(all_D$CV_diam,na.rm=T), length.out=length(all_D$CV_diam))
  preds <- predict(mod1, newdata = data.frame(CV_diam=newx), 
                   interval = 'confidence')
  lines(all_D$CV_diam, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$CV_diam, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$CV_diam), all_D$CV_diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$CV_diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                              paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                              ""))
}

#### WR sum ####
#rel WR RD
all_D<-all_D[order(all_D$vol) , ]
plot(abs(all_D$d_Vc_WR*100)~all_D$vol,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     xlab=expression(Root~volume~~"("~cm^3~")"),
     ylab=expression(Relative~error~volume~WR[SUM]~~"("~"%"~")")
)  

mod1<-lm(abs(all_D$d_Vc_WR*100)~log(all_D$vol))
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$vol), max(all_D$vol), length.out=length(all_D$vol))
  preds <- predict(mod1, newdata = data.frame(vol=newx), 
                   interval = 'confidence')
  lines(all_D$vol, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$vol, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$vol), all_D$vol), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$vol, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}


all_D<-all_D[order(all_D$diam) , ]
plot(abs(all_D$d_Vc_WR*100)~all_D$diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     yaxt="n",
     xlab=expression(Root~diameter~~"("~mm~")"),
     ylab=""
)  
Axis(side = 2, labels = FALSE, tck = 0.02)

mod1<-lm(abs(all_D$d_Vc_WR*100)~log(all_D$diam))
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$diam), max(all_D$diam), length.out=length(all_D$diam))
  preds <- predict(mod1, newdata = data.frame(diam=newx), 
                   interval = 'confidence')
  lines(all_D$diam, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$diam, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$diam), all_D$diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("topright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}

#rel WR CVRD
all_D<-all_D[order(all_D$CV_diam) , ]
plot(abs(all_D$d_Vc_WR*100)~all_D$CV_diam,
     pch=21,
     col='white',
     bg=all_D$col,
     cex=1.5,
     cex.main=3, cex.lab=2.5, cex.axis=1.5,
     #main="IJ_Rhizo",
     #xlim=Xlim,
     ylim=c(0,40),
     tck=0.02,
     yaxt="n",
     xlab=expression(CV~root~diameter),
     ylab=""
)  

mod1<-lm(abs(all_D$d_Vc_WR*100)~all_D$CV_diam)
if(summary(mod1)$coefficients[2,4]<0.05) {
  newx <- seq(min(all_D$CV_diam,na.rm=T), max(all_D$CV_diam,na.rm=T), length.out=length(all_D$CV_diam))
  preds <- predict(mod1, newdata = data.frame(CV_diam=newx), 
                   interval = 'confidence')
  lines(all_D$CV_diam, preds[ ,3], lty = 'dashed', col = 'darkgrey') # intervals
  lines(all_D$CV_diam, preds[ ,2], lty = 'dashed', col = 'darkgrey') # intervals
  polygon(c(rev(all_D$CV_diam), all_D$CV_diam), c(rev(preds[ ,3]), preds[ ,2]), col = 'lightgrey', border = NA)
  lines(all_D$CV_diam, preds[ ,1], lty = 'solid',lwd=2, col = 'black')# model
  
  legend("bottomright",bty="n",cex=2,legend=c(paste("Rē =",format(summary(mod1)$r.squared, digits=2)),
                                           paste("P =",format(summary(mod1)$coefficients[2,4], digits=1)),
                                           ""))
}


dev.off()


##### boxplot ####
wilcox.test(d_d_WR~group,all_D)
wilcox.test(d_d_IJC~group,all_D) #
wilcox.test(d_l_WR~group,all_D) #
wilcox.test(d_l_IJ~group,all_D) #
wilcox.test(d_V_WR~group,all_D) #

mycol<-c("lightgrey","darkgrey")
tiff(file="boxplot.tif",height=6.5,width=16.75, unit='in',res=300)
par(mfcol=c(1,5), 
    mar=c(3,4,1,0.1),  #c(bottom, left, top, right)
    oma=c(0,1,1.5,0), #c(bottom, left, top, right)
    mgp=c(3,1,0), # distance of label, numbers and ticks to axis
    bty="l")
boxplot(d_d_WR*100~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Relative~error~diameter~WR~~"("~"%"~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,outer=F) 
mtext("***", outer=F, at=c(1.5),cex=1.5)
boxplot(d_d_IJC*100~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Relative~error~diameter~IJ_Rhizo~~"("~"%"~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,outer=F) 
mtext("ns", outer=F, at=c(1.5),cex=1.5)
boxplot(d_l_WR*100~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Relative~error~length~WR~~"("~"%"~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,outer=F) 
mtext("**", outer=F, at=c(1.5),cex=1.5)
boxplot(d_l_IJ*100~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Relative~error~length~IJ_Rhizo~~"("~"%"~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,outer=F) 
mtext("ns", outer=F, at=c(1.5),cex=1.5)
boxplot(d_V_WR*100~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Relative~error~volume~~"("~"%"~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,outer=F) 
mtext("***", outer=F, at=c(1.5),cex=1.5)

dev.off()

#### boxplot overview ####
wilcox.test(Sum_length~group,all_D) #highly sign
wilcox.test(diam~group,all_D) #highly sign
wilcox.test(vol~group,all_D) #highly sign
wilcox.test(CV_diam~group,all_D) #
wilcox.test(overlap~group,all_D) #

mycol<-c("lightgrey","darkgrey")
tiff(file="boxplot_over.tif",height=6.5,width=16.75, unit='in',res=300)
par(mfcol=c(1,5), 
    mar=c(3,4,1,0.1),  #c(bottom, left, top, right)
    oma=c(0,1,1.5,0), #c(bottom, left, top, right)
    mgp=c(3,1,0), # distance of label, numbers and ticks to axis
    bty="l")
boxplot(diam~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Diameter~~"("~mm~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,outer=F) 
mtext("***", outer=F, at=c(1.5),cex=1.5)
boxplot(Sum_length~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Length~~"("~cm~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,font=2,outer=F) 
mtext("***", outer=F, at=c(1.5),cex=1.5)
boxplot(vol~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Volume~~"("~cm^3~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,font=2,outer=F) 
mtext("***", outer=F, at=c(1.5),cex=1.5)
boxplot(CV_diam~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(CV~root~diameter), side=2, 
      line=2.5,cex=1.5,adj=0.55,font=2,outer=F) 
mtext("***", outer=F, at=c(1.5),cex=1.5)
boxplot(overlap~group,all_D, cex.lab=2,cex.axis=1.5,col=mycol,boxwex = 0.5,notch=T)
mtext(text=expression(Overlap~~"("~"%"~")"), side=2, 
      line=2.5,cex=1.5,adj=0.55,font=2,outer=F) 
mtext("***", outer=F, at=c(1.5),cex=1.5)
dev.off()