###################################################
# Plotting performance measure: bias, MSE. CP     #
# Input: simulation data, empirical data          #
# Output: Figure 3: simulation results plot       #
# Original: July 18, 2022                         #
# Last modified: September 20, 2022               #
###################################################

rm(list = ls())
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
library(data.table);library(reshape2);library(tidyverse);
library(magrittr); library(gridExtra); library(ggpubr)

n = 1000

file_name <- c("CLINWTB_EI","MEANWST_EI","HR_EI","SBP_EI","DBP_EI","VO2_EI",
               "CLINWTB_sod","MEANWST_sod","HR_sod","SBP_sod","DBP_sod","VO2_sod",
               "CLINWTB_pot","MEANWST_pot","HR_pot","SBP_pot","DBP_pot","VO2_pot",
               "CLINWTB_pro","MEANWST_pro","HR_pro","SBP_pro","DBP_pro","VO2_pro")

NI = c(rep("EI",6),rep("sodium",6),rep("potassium",6),rep("protein",6))
Intake = c(rep("Energy",18),rep("Sodium",18),rep("Potassium",18),rep("Protein",18))
HO = c("BW (kg)","WC (cm)","HR (beat/min)","SBP (mmHg)","DBP (mmHg)","VO2 max (L/min)")
Outcome = rep(HO,4)

folderName <- "simulation/"
Num = 1

TrueBetas <- read.csv("Empirical_data/true_betas.csv")
TrueBetas$X <- rep(HO,4)

betas <- list()
Mean_bias <- list()
CPs <- list()
mse1<-c(); mse2<-c()
for(i in 1:24){
  if(i %in% 7:18){factor = 1000}else{factor=1}
  
  beta1 = TrueBetas[i,2]*factor
  
  df <- read.csv(paste0(folderName,file_name[i],".csv"),header=FALSE)
  colnames(df)<-c("SR0","SR1","SRlow","SRup","G0","G1","Glow","Gup","BIO0","BIO1","BIOlow","BIOup" )
  
  beta_hat_SR = df[,2]*factor 
  beta_hat_G = df[,6] *factor
  beta_hat_BIO = df[,10]*factor
  
  # BIO
  ##### Bias #####
  Mean_bias_BIO = mean(beta_hat_BIO,na.rm = T)
  (bias_BIO = Mean_bias_BIO-beta1)
  
  ##### MSE #####
  mse11 <- 1/n*(sum((beta_hat_BIO-mean(beta_hat_BIO,na.rm = T))^2,na.rm=T))
  mse21 <- bias_BIO^2
  (MSE_BIO = 1/n*(sum((beta_hat_BIO-mean(beta_hat_BIO,na.rm = T))^2,na.rm=T))+bias_BIO^2)
  
  # SR
  ##### Bias #####
  Mean_bias_SR = mean(beta_hat_SR,na.rm = T)
  (bias_SR = Mean_bias_SR - beta1)
  
  ##### MSE #####
  mse12 <- 1/n*(sum((beta_hat_SR-mean(beta_hat_SR,na.rm = T))^2,na.rm=T))
  mse22 <- bias_SR^2
  (MSE_SR = 1/n*(sum((beta_hat_SR-mean(beta_hat_SR,na.rm = T))^2,na.rm=T))+bias_SR^2)
  
  # G
  ##### Bias #####
  Mean_bias_G = mean(beta_hat_G,na.rm = T)
  (bias_G = Mean_bias_G-beta1)
  
  ##### MSE #####
  mse13 <- 1/n*(sum((beta_hat_BIO-mean(beta_hat_BIO,na.rm = T))^2,na.rm=T))
  mse23 <- bias_G^2
  (MSE_G = 1/n*(sum((beta_hat_BIO-mean(beta_hat_BIO,na.rm = T))^2,na.rm=T))+bias_G^2)
  
  Mean_bias[[i]] <- data.frame(BIO = Mean_bias_BIO, SR = Mean_bias_SR, G = Mean_bias_G, T_beta = beta1)
  
  betas[[i]] <- data.frame(NI = rep(NI[i],3),
                           group = rep(Outcome[i],3),
                           Bias = c(bias_BIO,bias_SR,bias_G),
                           MSE = c(MSE_BIO,MSE_SR,MSE_G),
                           metric = c("BIO","SR","G"))
  
  mse1 <- c(mse1,mse11,mse12,mse13)
  mse2 <- c(mse2,mse21,mse22,mse23)
  ##### CP #####
  df$BIO = ifelse(test = (df$BIOlow*factor<beta1)&(beta1<df$BIOup*factor), yes="Included",no="Excluded")
  df$SR = ifelse(test = (df$SRlow*factor<beta1)&(beta1<df$SRup*factor), yes="Included",no="Excluded")
  df$G = ifelse(test = (df$Glow*factor<beta1)&(beta1<df$Gup*factor), yes="Included",no="Excluded")
  table(df$BIO);table(df$SR);table(df$G)
  
  CP = c(table(df$BIO)[2]/10,table(df$SR)[2]/10,table(df$G)[2]/10)
  CP[is.na(CP)]<-0
  CPs[[i]] <- data.frame(Nutrition = rep(file_name[i],3),
                         parameter = c("BIO", "SR", "G"),
                         true = c(beta1,beta1,beta1),
                         CP =CP)
}

Betas <- do.call("rbind",betas)
Mean4bias <- do.call("rbind",Mean_bias)

Coverage <- do.call("rbind",CPs)

PM <- cbind(Betas,Coverage[,c(3,4)])
write.csv(PM,paste0("PerformanceMeasure.csv"))


##############################################
# Use the following for the result section   #
Bias <- Betas[,c(3,5,6)]
G <- Bias[Bias$metric =="G",]
SR <- Bias[Bias$metric =="SR",]
pm <- abs(G$Bias)<abs(SR$Bias)
SR$pm <- pm
##############################################

############################
##### Plot the results #####
############################
Font1 = 40; Font2 = 38; Font3 = 36; SIZE=12; SIZE2=10; ALPHA = 0.7

theme_box = function() {
  theme_bw()+
    theme(#aspect.ratio=1,
      axis.text.x = element_text(size = Font3),
      axis.text.y = element_text(size = Font3,face=c("plain", "plain", "plain", "bold", "plain", "plain")),
      plot.title = element_text(face = "bold",size=Font1),
      plot.subtitle = element_text(face = "bold",size=Font2),
      axis.text = element_text(size =Font3),
      axis.title = element_text(size =Font2, face = "bold"),
      legend.title=element_blank(),
      strip.text.x = element_text(size = Font2,face = "bold"),
      legend.text = element_text(size=Font3))
}

EI <- PM[1:18,]
Sod <- PM[19:36,]
Pot <- PM[37:54,]
Pro <- PM[55:72,]

dummy <- data.frame(metric = c("Bias","MSE","CP"), group = c("Bias","MSE","CP"), 
                    variable =c("Bias","MSE","CP"), value = c(0,0,95))

### EI 
Bias_EI <- EI[,c(2,3,4,5,7)]
data2plot <- melt(Bias_EI)
data2plot$group <- factor(data2plot$group,levels=HO,ordered=TRUE)
data2plot$variable <- factor(data2plot$variable,levels=c("Bias","MSE","CP"),ordered=TRUE)

p1 <- ggplot(data2plot,aes(x=value, y=group, shape=metric, color=metric))+
  geom_vline(data = dummy, aes(xintercept=value), size=1,color="grey")+#linetype = "dashed")+
  geom_point(size = SIZE,position = position_jitter(height = 0.1),alpha=ALPHA)+
  scale_shape_manual(values=c(15,16,17))+
  labs(title = "A: Energy Intake (1000kcal per day)", x = "Value", y = "Health Outcomes") +    
  scale_y_discrete(limits=rev) +
  theme_box()+
  facet_wrap(~factor(variable), scales = "free_x")

### Sod 
Bias_sod <- Sod[,c(2,3,4,5,7)]
data2plot2 <- melt(Bias_sod)
data2plot2$group <- factor(data2plot2$group,levels=HO,ordered=TRUE)
data2plot2$variable <- factor(data2plot2$variable,levels=c("Bias","MSE","CP"),ordered=TRUE)

p2 <- ggplot(data2plot2,aes(x=value, y=group, shape=metric, color=metric))+
  geom_vline(data = dummy, aes(xintercept=value), size=1,color="grey")+#linetype = "dashed")+
  geom_point(size = SIZE,position = position_jitter(height = 0.1),alpha=ALPHA)+
  scale_shape_manual(values=c(15,16,17))+
  labs(title = "B: Sodium Intake (g per day)", x = "Value", y = "Health Outcomes") +    
  scale_y_discrete(limits=rev) +
  theme_box()+
  facet_wrap(~factor(variable), scales = "free_x")


### Pot 
Bias_pot <- Pot[,c(2,3,4,5,7)]
data2plot3 <- melt(Bias_pot)
data2plot3$group <- factor(data2plot3$group,levels=HO,ordered=TRUE)
data2plot3$variable <- factor(data2plot3$variable,levels=c("Bias","MSE","CP"),ordered=TRUE)

p3 <- ggplot(data2plot3,aes(x=value, y=group, shape=metric, color=metric))+
  geom_vline(data = dummy, aes(xintercept=value), size=1,color="grey")+#linetype = "dashed")+
  geom_point(size = SIZE,position = position_jitter(height = 0.1),alpha=ALPHA)+
  scale_shape_manual(values=c(15,16,17))+
  labs(title = "C: Potassium Intake (g per day)", x = "Value", y = "Health Outcomes") +    
  scale_y_discrete(limits=rev) +
  theme_box()+
  facet_wrap(~factor(variable), scales = "free_x")

### pro
Bias_pro <- Pro[,c(2,3,4,5,7)]
data2plot4 <- melt(Bias_pro)
data2plot4$group <- factor(data2plot4$group,levels=HO,ordered=TRUE)
data2plot4$variable <- factor(data2plot4$variable,levels=c("Bias","MSE","CP"),ordered=TRUE)

p4 <- ggplot(data2plot4,aes(x=value, y=group, shape=metric, color=metric))+
  geom_vline(data = dummy, aes(xintercept=value), size=1,color="grey")+#linetype = "dashed")+
  geom_point(size = SIZE,position = position_jitter(height = 0.1),alpha=ALPHA)+
  scale_shape_manual(values=c(15,16,17))+
  labs(title = "D: Protein Intake (g per day)",  x = "Value", y = "Health Outcomes") +    
  scale_y_discrete(limits=rev) +
  theme_box()+
  facet_wrap(~factor(variable), scales = "free_x")

g <- ggarrange(p1, p2, p3, p4, ncol=1, nrow=4, common.legend = TRUE, legend="bottom")
ggsave(filename = paste0("Figures/Figure3.png"), g, dpi = 400,width =900,height = 900,units = c("mm"))