##########################################################
# Input: beta-jack-summary.csv, pearson-jack-summary.csv #
# Output: Figure2.tif                                    #
# Modification: June 20, 2022                            #
# Last modified: September 20, 2022                      #
##########################################################

rm(list = ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

library(data.table)

tiff("Figures/Figure2.tif", res = 1000, width = 6, height = 6, units = 'in')
layout(matrix(c(1:3), 1, 3, byrow = F))
#### Read data ####
data = read.csv("Empirical_data/beta-jack-summary.csv",header=F)[,12:23]*100
data = data[rev(1:24),]

Categories = c("Energy, Weight",
               "Energy, Waist circumference",
               "Energy, HR post fitness test",
               "Energy, Resting SBP",
               "Energy, Resting DBP",
               "Energy, VO2 max",
               "Sodium, Weight",
               "Sodium, Waist circumference",
               "Sodium, HR post fitness test",
               "Sodium, Resting SBP",
               "Sodium, Resting DBP",
               "Sodium, VO2 max",
               "Potassium, Weight",
               "Potassium, Waist circumference",
               "Potassium, HR post fitness test",
               "Potassium, Resting SBP",
               "Potassium, Resting DBP",
               "Potassium, VO2 max",
               "Protein, Weight",
               "Protein, Waist circumference",
               "Protein, HR post fitness test",
               "Protein, Resting SBP",
               "Protein, Resting DBP",
               "Protein, VO2 max")

# Percent bias #
# (SR-DLW)/DLW
par(mgp=c(1.5,1.0,0.2), mai=c(0.4,0.9,0.1,0.0), ps=8)

# Plot the left label
plot(0, 1, xlab="",ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, type="n")
mtext(text=rev(Categories), 
      font=rev(c(c(4,4,1,1,1,1),c(4,4,1,1,1,1),c(1,1,1,1,1,4),c(4,4,1,1,1,1))), 
      side=2,line=6.5,at=1:24,las=1,adj=0)

par(mgp=c(1.5,1.0,0.2), mai=c(0.4,0.25,0.1,0.1), ps=12)
plot(0, 1, xlab="Percent bias (%)",
     ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, type="n")
for(i in 1:length(data[,1])){
  par(new=T)
  if(data[i,1] < 200 & data[i,1] > -200){
    plot(data[i,1], i, xlab="",ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, pch=0, main="")
  }else if(data[i,1] > 200){
    plot(215, i, xlab="",ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, pch=15, main="")
  }else if(data[i,1] < -200){
    plot(-215, i, xlab="",ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, pch=15, main="")
  }
}
for(i in 1:length(data[,1])){
  if(data[i,2] < -200 & data[i,3] < 200){
    arrows(-200, i, data[i,3], i, angle = 90, length = 0.02)
  }else if(data[i,2] > -200 & data[i,3] < 200){
    arrows(data[i,2], i, data[i,3], i, angle = 90, length = 0.02)
    arrows(data[i,3], i, data[i,2], i, angle = 90, length = 0.02)
  }else if(data[i,2] > -200 & data[i,3] > 200){
    arrows(200, i, data[i,2], i, angle = 90, length = 0.02)
  }else if(data[i,2] < -200 & data[i,3] > 200){
    arrows(-200, i, 200, i, angle = 90, length = 0.0)
  }
}
abline(v=0, lty=2)
axis(1, pos=0.5, at=seq(-200,200,100), labels=seq(-200,200,100), adj = 0) # x axis
axis(2, pos=280, at=1:24, labels=rev(c("*","*","","","","", # Energy
                                       "","","","","","", # Sodium
                                       "","","","","","", # Potassium
                                       "","","","","","")), adj = 0, tick = TRUE, lwd = 0, cex.axis=1.4, las=1) # Protein

# Percent remaining bias #
# (G-DLW)/DLW
par(mgp=c(1.5,1.0,0.2), mai=c(0.4,0.1,0.1,0.1), ps=12)
plot(0, 1, xlab="Percent remaining bias (%)",
     ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, type="n")
for(i in 1:length(data[,1])){
  par(new=T)
  if(data[i,10] < 200 & data[i,10] > -200){
    plot(data[i,10], i, xlab="",ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, pch=0, main="")
  }else if(data[i,10] > 200){
    plot(215, i, xlab="",ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, pch=15, main="")
  }else if(data[i,10] < -200){
    plot(-215, i, xlab="",ylab="",xlim=c(-215,215),ylim=c(1,24), axes=F, pch=15, main="")
  }
}
for(i in 1:length(data[,10])){
  if(data[i,11] < -200 & data[i,12] < 200){
    arrows(-200, i, data[i,12], i, angle = 90, length = 0.02)
  }else if(data[i,11] > -200 & data[i,12] < 200){
    arrows(data[i,11], i, data[i,12], i, angle = 90, length = 0.02)
    arrows(data[i,12], i, data[i,11], i, angle = 90, length = 0.02)
  }else if(data[i,11] > -200 & data[i,12] > 200){
    arrows(200, i, data[i,11], i, angle = 90, length = 0.02)
  }else if(data[i,11] < -200 & data[i,12] > 200){
    arrows(-200, i, 200, i, angle = 90, length = 0.0)
  }
}
axis(1, pos=0.5, at=seq(-200,200,100), labels=seq(-200,200,100), adj = 0) # x axis
axis(2, pos=280, at=1:24, labels=rev(c("*","*","","","","", # Energy
                                       "","","","","","", # Sodium
                                       "","","","","","", # Potassium
                                       "","","","","","")), adj = 0, tick = TRUE, lwd = 0, cex.axis=1.4, las=1) # Protein

abline(v=0, lty=2)
dev.off()