
##Test2: RER IN SITU RHIZOTRONS_Figure5

library (doBy) # to use the function "summarizeBy()"
library (reshape)

addr_input = "D://"
addr_output="D://"
##
morpDF = read.table(file = paste (addr_input, "Data_Figures 5&6.r.csv", sep = ""), 
                    header = T,  sep = ";", dec=".")

morp2DF = morpDF[!is.na(morpDF$CumLength_cm), ]
morp3DF = morp2DF[!is.na(morp2DF$RER), ]

summary (morp3DF$Method)

morp3DF$MethodName = factor (ifelse (morp3DF$Method == 1, "Flatbed",
                                    ifelse (morp3DF$Method == 2, "Handheld",
                                            ifelse (morp3DF$Method == 3, "Smartphone",
                                                    ifelse (morp3DF$Method == 4, "Transparentsheet",  "Error")))))

morp3DF$Split = with (morp3DF, paste(RootName, "_", Sequence, sep = "" ) )
morp3DF$Split = with (morp3DF, factor(paste(RootName, "_", Sequence, sep = "" )) )
morp3DF$Split = with (morp3DF, factor(paste( "M",Method,"_R", RootName, "_S", Sequence, sep = "" )) )
morp3DF$Split = with (morp3DF, factor(paste("R", RootName, "_S", Sequence, sep = "" )) )


# --------------------------------
# COMPUTE THE RELATIVE RER

morp3DF <- morp3DF[!(morp3DF$RER == 0 & morp3DF$Method == 1),] # Remove the data where method = 1 and RER = 0
totalLevel = levels(factor(morp3DF$Split))
for (level in totalLevel) {
  if(length(morp3DF$RER[morp3DF$Split == level & morp3DF$Method == 1]) > 0){
    morp3DF$RER_rel[morp3DF$Split == level] <- morp3DF$RER[morp3DF$Split == level] / 
                                                    morp3DF$RER[morp3DF$Split == level & morp3DF$Method == 1]
  }
}
morp3DF <- morp3DF[!is.na(morp3DF$RER_rel),]

# --------------------------------

nameFig = "boxplot_RER_GL.Fig5.pdf"

pdf (file = paste(addr_output, nameFig, sep =""), width = 7, height = 6 )
  x0 = morp3DF$MethodName
  y0 = morp3DF$RER_rel
  
  boxplot ( y0 ~ x0, type ="n", xlab = NA, ylab = NA,
            names = c ("Flatbed", "Handheld", "Smartphone", "TransparentSheet"), 
            outline=T, boxwex=0.5)
  
  totalLevel = levels(factor(morp3DF$Split))
  for (level in totalLevel) {
    # level = "R0_S2"
    dataI = morp3DF[morp3DF$Split == level,]
    points (x = 1:4,
            y = dataI$RER_rel,
            col="#91765C60", pch = 21, bg = "#FF0D1760", cex = 1)
  }
  
  boxplot ( y0 ~ x0, type ="n", xlab = NA, ylab = NA, col="#FFFFFF50",
            names=c ("Flatbed", "Handheld", "Smartphone", "TransparentSheet"), 
            outline=T, boxwex=0.5, add=T)
  
  abline(h=c(1.1, 0.9), lty=2)
  
  mtext (side = 1, "Methods", cex = 1.2, line = 3)
  mtext (side = 2, "Relative root elongation rate (mm/day)", cex = 1.2, line = 2.5)
  
dev.off()


## ANOVA AND NORMALITY TEST


lm1 = lm (formula = RER_rel ~ MethodName, data = morp3DF)
summary(lm1)
summary (aov (lm1))
lm2 = lm (formula = RER_rel ~ MethodName + as.factor(Sequence), data = morp3DF)
summary (lm2)
summary (aov (lm2))
lm3 = lm (formula = RER_rel ~ MethodName * as.factor(Sequence), data = morp3DF)
summary (lm3)
summary (aov (lm3))

shapiro.test(morp3DF$RER_rel)
shapiro.test(morp3DF$RER_rel ^ 2)
shapiro.test(morp3DF$RER_rel ^ 0.5)
plot (density(morp3DF$RER_rel), col = 3 )
hist (density(morp3DF$RER_rel), col = 3 )
hist (morp3DF$RER_rel )
par (mfrow = c(2,2) )
plot (density(morp3DF$RER_rel), col = 3 )
plot (density(morp3DF$RER_rel ^ 0.5), col = 3 )
plot (density(log(morp3DF$RER_rel + 1) ), col = 3 )


lm4 = lm (formula = RER_rel ~ MethodName + as.factor(RootName), data = morp3DF)
summary (lm4)
lm4 = lm (formula = RER_rel ~ MethodName * as.factor(RootName), data = morp3DF)
summary (lm4)


