combine_sim_plots_misspec.R

############################################################################################
# Combine plots from static setting simulations, comparing 6 vs 40 vs 100 covariates       #
############################################################################################

library(tidyverse)
library(ggplot2)
library(data.table)
library(latex2exp)
library(dplyr)
library(grid)
library(gtable)

source("./src/misc_fns.R")

results_dir <- './sim_results/'
if(!dir.exists(results_dir)){
  print(paste0('create folder for results plot at: ', results_dir))
  dir.create(results_dir)
}

# Define parameters

options(echo=TRUE)
args <- commandArgs(trailingOnly = TRUE) # args <- c("outputs/20231109", "outputs/20231110", "outputs/20231114", 6, 'TRUE') 
output.path <- as.character(args[1]) 
output.path2 <- as.character(args[2])
output.path3 <- as.character(args[3])
J <- as.numeric(args[4])
use.SL <- as.character(args[5])

estimand <- "ate"
m <- choose(J, 2)
n <- ifelse(J==6, 10000, 5000)
outcome.type <- 'binomial'

overlap.setting <- c("adequate","inadequate","rct")
gamma.setting <- c("zero","low","yang") 

estimators <- c("tmle","tmle_bin", "tmle_glm","tmle_glm_bin", "gcomp", "gcomp_glm", "iptw", "iptw_bin",  "iptw_glm", "iptw_glm_bin","aiptw", "aiptw_bin","aiptw_glm", "aiptw_glm_bin")
if(use.SL){
  estimators <- estimators[grep("glm", estimators, invert=TRUE)]
  estimator.color <- c("#F8766D", "#A3A500", "#00BF7D", "#00B0F6", "#E76BF3")
}else{
  estimators <- estimators[grep("glm", estimators)]
  estimator.color <- c("#FF6C90", "#D89000", "#39B600", "#00BFC4", "#FF62BC")
}

n.estimators <- as.numeric(length(estimators))

filenames <- c(list.files(path=output.path, pattern = ".rds", full.names = TRUE), 
               list.files(path=output.path2, pattern = ".rds", full.names = TRUE),
               list.files(path=output.path3, pattern = ".rds", full.names = TRUE))

filenames <- filenames[grep(paste0("J_",J),filenames)]
filenames <- filenames[grep(paste0("n_",n,"_"),filenames)]
filenames <- filenames[grep(paste0("outcome_type_",outcome.type),filenames)]
filenames <- filenames[grep(paste0("use_SL_",use.SL),filenames)]
filenames <- filenames[grep(paste0("use_SL_",use.SL),filenames)]
filenames <- filenames[grep(paste0("R_1000_"),filenames)] # keep only misspc runs (R=1000)

if(any( duplicated(substring(filenames, 18)))){
  print("removing duplicate filenames")
  filenames <- filenames[-which(duplicated(substring(filenames, 18)))]
}

results <- list() # structure is: [[filename]][[comparison]][[metric]]
for(f in filenames){
  print(f)
  result.matrix <- readRDS(f)
  ncovars <- ifelse(any(grep("covars40_TRUE",f)), 40, ifelse(any(grep("covars100_TRUE",f)), 100, 6))
  misspec <- ifelse(any(grep("misTreat_TRUE",f)), "misTreat", ifelse(any(grep("misOut_TRUE",f)), "misOut", "misBoth"))
  R <- ncol(result.matrix )
  if(isTRUE(grep("use_SL_FALSE",f)==1)){
    estimator <- estimators[grep("glm", estimators)]
    rownames(result.matrix)[grep("tmle",rownames(result.matrix))] <- gsub("tmle","tmle_glm", rownames(result.matrix)[grep("tmle",rownames(result.matrix))])
    rownames(result.matrix)[grep("gcomp",rownames(result.matrix))] <- gsub("gcomp","gcomp_glm", rownames(result.matrix)[grep("gcomp",rownames(result.matrix))])
    rownames(result.matrix)[grep("_iptw",rownames(result.matrix))] <- gsub("_iptw","_iptw_glm", rownames(result.matrix)[grep("_iptw",rownames(result.matrix))])
    rownames(result.matrix)[grep("_aiptw",rownames(result.matrix))] <- gsub("_aiptw","_aiptw_glm", rownames(result.matrix)[grep("_aiptw",rownames(result.matrix))])
  }else{
    estimator <- setdiff(estimators,estimators[grep("glm", estimators)])
  }
  bias <- matrix(NA, R, n.estimators)
  colnames(bias) <- paste0("bias_",estimand,"_",estimators)
  CP <- matrix(NA, R, n.estimators)
  colnames(CP) <- paste0("CP_",estimand,"_",estimators)
  CIW <- matrix(NA, R,  n.estimators)
  colnames(CIW) <- paste0("CIW_",estimand,"_",estimators)
  for(j in 1:m){
    for(i in paste0("bias_",estimand,"_",estimator)){
      bias[,which(i==colnames(bias))] <- unlist(lapply(result.matrix[i,], "[", j))[1:R]
    }
    for(i in paste0("CP_",estimand,"_",estimator)){
      CP[,which(i==colnames(CP))] <- unlist(lapply(result.matrix[i,], "[", j))[1:R]
    }
    for(i in paste0("CIW_",estimand,"_",estimator)){
      CIW[,which(i==colnames(CIW))] <- unlist(lapply(result.matrix[i,], "[", j))[1:R]
    }
    results[[f]][[j]] <- list("bias"=bias,"CP"=CP,"CIW"=CIW,"R"=R,"ncovars"=ncovars,"misspec"=misspec)
  }
}

# Create New lists
# structure is: [[estimator]][[filename]][[comparison]]

bias <- list()
for(i in 1:length(paste0("bias_",estimand,"_",estimators))){
  bias[[i]] <- lapply(filenames, function(f) lapply(1:m, function(j) results[[f]][[j]]$bias[,i]))
}

CP <- list()
for(i in 1:length(paste0("CP_",estimand,"_",estimators))){
  CP[[i]] <- lapply(filenames, function(f) lapply(1:m, function(j) results[[f]][[j]]$CP[,i]))
}

CIW <- list()
for(i in 1:length(paste0("CIW_",estimand,"_",estimators))){
  CIW[[i]] <- lapply(filenames, function(f) lapply(1:m, function(j) results[[f]][[j]]$CIW[,i]))
}

# Create dataframe for plot
results.df <- data.frame("abs.bias"=abs(unlist(bias)),
                         "Coverage"=unlist(CP),
                         "CIW"=unlist(CIW),
                         "filename"=c(rep(unlist(sapply(1:length(filenames), function(i) rep(filenames[i], length.out=results[[i]][[1]]$R*m))), n.estimators)))

results.df$ncovars <- "6"
results.df$ncovars[grep("covars40_TRUE",results.df$filename)] <- "40" 
results.df$ncovars[grep("covars100_TRUE",results.df$filename)] <- "100" 

results.df$misspec <- "misBoth"
results.df$misspec[grep("misTreat_TRUE",results.df$filename)] <- "misTreat" 
results.df$misspec[grep("misOut_TRUE",results.df$filename)] <- "misOut" 

if(use.SL){
  results.df$Estimator <- c(rep("TMLE-multi. (SL)",length.out=length(c(unlist(CP[[1]])))), 
                            rep("TMLE-bin. (SL)",length.out=length(c(unlist(CP[[2]])))),
                            rep("G-comp. (SL)",length.out=length(c(unlist(CP[[3]])))),
                            rep("IPTW-multi. (SL)",length.out=length(c(unlist(CP[[4]])))),
                            rep("IPTW-bin. (SL)",length.out=length(c(unlist(CP[[5]])))),
                            rep("AIPTW-multi. (SL)",length.out=length(c(unlist(CP[[6]])))),
                            rep("AIPTW-bin. (SL)",length.out=length(c(unlist(CP[[7]])))))
}else{
  results.df$Estimator <- c(rep("TMLE-multi. (GLM)",length.out=length(c(unlist(CP[[1]])))), 
                            rep("TMLE-bin. (GLM)",length.out=length(c(unlist(CP[[2]])))),
                            rep("G-comp. (GLM)",length.out=length(c(unlist(CP[[3]])))),
                            rep("IPTW-multi. (GLM)",length.out=length(c(unlist(CP[[4]])))),
                            rep("IPTW-bin. (GLM)",length.out=length(c(unlist(CP[[5]])))),
                            rep("AIPTW-multi. (GLM)",length.out=length(c(unlist(CP[[6]])))),
                            rep("AIPTW-bin. (GLM)",length.out=length(c(unlist(CP[[7]])))))
}


if(J==3){
  results.df$comparison <- factor(rep(unlist(sapply(1:length(filenames), function (i) c(rep("32",length.out=results[[i]][[1]]$R), rep("31",length.out=results[[i]][[1]]$R), rep("21",length.out=results[[i]][[1]]$R)))), n.estimators))                        
}else if(J==6){
  results.df$comparison <- factor(rep(unlist(sapply(1:length(filenames), function (i) c(rep("65",length.out=results[[i]][[1]]$R),
                                                                                        rep("64",length.out=results[[i]][[1]]$R),
                                                                                        rep("63",length.out=results[[i]][[1]]$R),
                                                                                        rep("62",length.out=results[[i]][[1]]$R),
                                                                                        rep("61",length.out=results[[i]][[1]]$R),
                                                                                        rep("54",length.out=results[[i]][[1]]$R),
                                                                                        rep("53",length.out=results[[i]][[1]]$R),
                                                                                        rep("52",length.out=results[[i]][[1]]$R),
                                                                                        rep("51",length.out=results[[i]][[1]]$R),
                                                                                        rep("43",length.out=results[[i]][[1]]$R),
                                                                                        rep("42",length.out=results[[i]][[1]]$R),
                                                                                        rep("41",length.out=results[[i]][[1]]$R),
                                                                                        rep("32",length.out=results[[i]][[1]]$R),
                                                                                        rep("31",length.out=results[[i]][[1]]$R),
                                                                                        rep("21",length.out=results[[i]][[1]]$R)))), n.estimators))                        
}
results.df$J <- ifelse(J==3,3,6)
results.df$overlap.setting <- NA
results.df$gamma.setting <- NA

for(s in overlap.setting){
  print(s)
  results.df[grep(s, results.df$filename),]$overlap.setting <- s
}

for(s in gamma.setting){
  print(s)
  results.df[grep(s, results.df$filename),]$gamma.setting <- s
}

results.df <- results.df[rowSums(is.na(results.df[,1:3])) ==0,]

# create coverage rate variable and filter Estimators

results.df <- results.df %>%
  group_by(Estimator,comparison,overlap.setting,gamma.setting,ncovars,misspec,J) %>% 
  mutate(CP = mean(Coverage)) 

if(use.SL){
  results.df <- results.df %>%
    filter(Estimator %in% c("TMLE-multi. (SL)",
                            "TMLE-bin. (SL)",
                            "G-comp. (SL)",
                            "IPTW-multi. (SL)",
                            "IPTW-bin. (SL)"))
}else{
  results.df <- results.df %>%
    filter(Estimator %in% c("TMLE-multi. (GLM)",
                            "TMLE-bin. (GLM)",
                            "G-comp. (GLM)",
                            "IPTW-multi. (GLM)",
                            "IPTW-bin. (GLM)"))
}
n.estimators <- length(unique(results.df$Estimator))

# reshape and plot
results.df$id <- with(results.df, paste(overlap.setting, gamma.setting, ncovars, misspec, J, sep = "_"))
results_long <- reshape2::melt(results.df[!colnames(results.df) %in% c("J","id","filename")], id.vars=c("Estimator","comparison","overlap.setting","gamma.setting","ncovars","misspec"))  # convert to long format

variable_names3 <- list(
  'low'="Low event rate",
  'yang'="Moderate event rate",
  'zero'= "No treatment effect",
  'adequate'= "Adequate overlap",
  'inadequate'= "Inadequate overlap",
  "rct"= "RCT",
  "40" = "40 covariates",
  "100" = "100 covariates",
  "misBoth" = "Outcome and treatment models misspecified",
  "misTreat" = "Only treatment model misspecified",
  "misOut" = "Only outcome model misspecified")  

labeller3 <- function(variable,value){
  return(variable_names3[value])
}

## coverage bar plot (avg. across comparisons)

coverage.df <- data.frame(results.df) %>% 
  select(c(Coverage,Estimator,overlap.setting,gamma.setting,ncovars,misspec)) %>% 
  group_by(Estimator,overlap.setting,gamma.setting,ncovars,misspec) %>% 
  mutate(CP = mean(Coverage)) %>%
  dplyr::slice(n()) %>%
  select(c(CP,Estimator,overlap.setting,gamma.setting,ncovars,misspec))

#misBoth
sim.results.cp.avg.misBoth <- ggplot(data=coverage.df[coverage.df$misspec=="misBoth",],
                                     aes(x=Estimator, y=CP, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~ gamma.setting, scales = "fixed", labeller=labeller3)   + ylab("Average coverage probability over all pairwise comparisons") + 
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  geom_hline(yintercept = 0.95, linetype="dotted")+
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.cp.avg <- ggplotGrob(sim.results.cp.avg.misBoth)

# Labels 
labelR <- "Treatment setting" #TeX("$\\beta_j$ treatment model coefficients")
labelT <- "Outcome setting" #TeX("$\\gamma_j$ outcome model coefficients")

# Construct the new strip grobs
stripR <- gTree(name = "Strip_right", children = gList(
  rectGrob(gp = gpar(col = NA, fill = "grey85")),
  textGrob(labelR, rot = -90, gp = gpar(fontsize=22, col = "grey10"))))

stripT <- gTree(name = "Strip_top", children = gList(
  rectGrob(gp = gpar(col = NA, fill = "grey85")),
  textGrob(labelT, gp = gpar(fontsize=22, col = "grey10"))))

# Get the positions of the strips in the gtable: t = top, l = left, ...
posR <- subset(z.cp.avg$layout, grepl("strip-r", name), select = t:r)
posT <- subset(z.cp.avg$layout, grepl("strip-t", name), select = t:r)

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.cp.avg$widths[max(posR$r)]    # width of current right strips
height <- z.cp.avg$heights[min(posT$t)]  # height of current top strips

z.cp.avg <- gtable_add_cols(z.cp.avg, width, max(posR$r))  
z.cp.avg <- gtable_add_rows(z.cp.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.cp.avg <- gtable_add_grob(z.cp.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.cp.avg <- gtable_add_grob(z.cp.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.cp.avg <- gtable_add_cols(z.cp.avg, unit(1/5, "line"), max(posR$r))
z.cp.avg <- gtable_add_rows(z.cp.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.cp.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_cp_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misBoth","_R_",results[[1]][[1]]$R,".png"),plot = z.cp.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_cp_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misBoth","_R_",results[[1]][[1]]$R,".png"),plot = z.cp.avg,scale=2)
}

#misTreat
sim.results.cp.avg.misTreat <- ggplot(data=coverage.df[coverage.df$misspec=="misTreat",],
                                     aes(x=Estimator, y=CP, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~ gamma.setting, scales = "fixed", labeller=labeller3)   + ylab("Average coverage probability over all pairwise comparisons") + 
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  geom_hline(yintercept = 0.95, linetype="dotted")+
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.cp.avg <- ggplotGrob(sim.results.cp.avg.misTreat)

# Labels 
labelR <- "Treatment setting" #TeX("$\\beta_j$ treatment model coefficients")
labelT <- "Outcome setting" #TeX("$\\gamma_j$ outcome model coefficients")

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.cp.avg$widths[max(posR$r)]    # width of current right strips
height <- z.cp.avg$heights[min(posT$t)]  # height of current top strips

z.cp.avg <- gtable_add_cols(z.cp.avg, width, max(posR$r))  
z.cp.avg <- gtable_add_rows(z.cp.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.cp.avg <- gtable_add_grob(z.cp.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.cp.avg <- gtable_add_grob(z.cp.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.cp.avg <- gtable_add_cols(z.cp.avg, unit(1/5, "line"), max(posR$r))
z.cp.avg <- gtable_add_rows(z.cp.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.cp.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_cp_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misTreat","_R_",results[[1]][[1]]$R,".png"),plot = z.cp.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_cp_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misTreat","_R_",results[[1]][[1]]$R,".png"),plot = z.cp.avg,scale=2)
}

#misOut
sim.results.cp.avg.misOut <- ggplot(data=coverage.df[coverage.df$misspec=="misOut",],
                                     aes(x=Estimator, y=CP, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~ gamma.setting, scales = "fixed", labeller=labeller3)   + ylab("Average coverage probability over all pairwise comparisons") + 
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  geom_hline(yintercept = 0.95, linetype="dotted")+
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.cp.avg <- ggplotGrob(sim.results.cp.avg.misOut)

# Labels 
labelR <- "Treatment setting" #TeX("$\\beta_j$ treatment model coefficients")
labelT <- "Outcome setting" #TeX("$\\gamma_j$ outcome model coefficients")

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.cp.avg$widths[max(posR$r)]    # width of current right strips
height <- z.cp.avg$heights[min(posT$t)]  # height of current top strips

z.cp.avg <- gtable_add_cols(z.cp.avg, width, max(posR$r))  
z.cp.avg <- gtable_add_rows(z.cp.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.cp.avg <- gtable_add_grob(z.cp.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.cp.avg <- gtable_add_grob(z.cp.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.cp.avg <- gtable_add_cols(z.cp.avg, unit(1/5, "line"), max(posR$r))
z.cp.avg <- gtable_add_rows(z.cp.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.cp.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_cp_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misOut","_R_",results[[1]][[1]]$R,".png"),plot = z.cp.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_cp_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misOut","_R_",results[[1]][[1]]$R,".png"),plot = z.cp.avg,scale=2)
}

## bias bar plot (avg. across comparisons)

bias.df <- data.frame(results.df) %>% 
  select(c(abs.bias,Estimator,overlap.setting,gamma.setting,ncovars,misspec)) %>% 
  group_by(Estimator,overlap.setting,gamma.setting,ncovars,misspec) %>% 
  mutate(bias = mean(abs.bias)) %>%
  dplyr::slice(n()) %>%
  select(c(bias,Estimator,overlap.setting,gamma.setting,ncovars,misspec))

# misBoth
sim.results.bias.avg.misBoth <- ggplot(data=bias.df[bias.df$misspec=="misBoth",],
                                       aes(x=Estimator, y=bias, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~  gamma.setting, scales = "free", labeller=labeller3)   + ylab("Average bias over all pairwise comparisons") +  
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.bias.avg <- ggplotGrob(sim.results.bias.avg.misBoth)

# Get the positions of the strips in the gtable: t = top, l = left, ...
posR <- subset(z.bias.avg$layout, grepl("strip-r", name), select = t:r)
posT <- subset(z.bias.avg$layout, grepl("strip-t", name), select = t:r)

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.bias.avg$widths[max(posR$r)]    # width of current right strips
height <- z.bias.avg$heights[min(posT$t)]  # height of current top strips

z.bias.avg <- gtable_add_cols(z.bias.avg, width, max(posR$r))  
z.bias.avg <- gtable_add_rows(z.bias.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.bias.avg <- gtable_add_grob(z.bias.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.bias.avg <- gtable_add_grob(z.bias.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.bias.avg <- gtable_add_cols(z.bias.avg, unit(1/5, "line"), max(posR$r))
z.bias.avg <- gtable_add_rows(z.bias.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.bias.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_bias_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misBoth","_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.bias.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_bias_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misBoth","_R_",results[[1]][[1]]$R,".png"),plot = z.bias.avg,scale=2)
}

# misOut
sim.results.bias.avg.misOut <- ggplot(data=bias.df[bias.df$misspec=="misOut",],
                                       aes(x=Estimator, y=bias, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~  gamma.setting, scales = "free", labeller=labeller3)   + ylab("Average bias over all pairwise comparisons") +  
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.bias.avg <- ggplotGrob(sim.results.bias.avg.misOut)

# Get the positions of the strips in the gtable: t = top, l = left, ...
posR <- subset(z.bias.avg$layout, grepl("strip-r", name), select = t:r)
posT <- subset(z.bias.avg$layout, grepl("strip-t", name), select = t:r)

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.bias.avg$widths[max(posR$r)]    # width of current right strips
height <- z.bias.avg$heights[min(posT$t)]  # height of current top strips

z.bias.avg <- gtable_add_cols(z.bias.avg, width, max(posR$r))  
z.bias.avg <- gtable_add_rows(z.bias.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.bias.avg <- gtable_add_grob(z.bias.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.bias.avg <- gtable_add_grob(z.bias.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.bias.avg <- gtable_add_cols(z.bias.avg, unit(1/5, "line"), max(posR$r))
z.bias.avg <- gtable_add_rows(z.bias.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.bias.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_bias_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misOut","_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.bias.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_bias_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misOut","_R_",results[[1]][[1]]$R,".png"),plot = z.bias.avg,scale=2)
}

# misTreat
sim.results.bias.avg.misTreat <- ggplot(data=bias.df[bias.df$misspec=="misTreat",],
                                       aes(x=Estimator, y=bias, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~  gamma.setting, scales = "free", labeller=labeller3)   + ylab("Average bias over all pairwise comparisons") +  
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.bias.avg <- ggplotGrob(sim.results.bias.avg.misTreat)

# Get the positions of the strips in the gtable: t = top, l = left, ...
posR <- subset(z.bias.avg$layout, grepl("strip-r", name), select = t:r)
posT <- subset(z.bias.avg$layout, grepl("strip-t", name), select = t:r)

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.bias.avg$widths[max(posR$r)]    # width of current right strips
height <- z.bias.avg$heights[min(posT$t)]  # height of current top strips

z.bias.avg <- gtable_add_cols(z.bias.avg, width, max(posR$r))  
z.bias.avg <- gtable_add_rows(z.bias.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.bias.avg <- gtable_add_grob(z.bias.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.bias.avg <- gtable_add_grob(z.bias.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.bias.avg <- gtable_add_cols(z.bias.avg, unit(1/5, "line"), max(posR$r))
z.bias.avg <- gtable_add_rows(z.bias.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.bias.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_bias_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misTreat","_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.bias.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_bias_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misTreat","_R_",results[[1]][[1]]$R,".png"),plot = z.bias.avg,scale=2)
}

## CIW bar plot (avg. across comparisons)

CIW.df <- data.frame(results.df) %>% 
  select(c(CIW,Estimator,overlap.setting,gamma.setting,ncovars,misspec)) %>% 
  group_by(Estimator,overlap.setting,gamma.setting,ncovars,misspec) %>% 
  mutate(CIW = mean(CIW)) %>%
  dplyr::slice(n()) %>%
  select(c(CIW,Estimator,overlap.setting,gamma.setting,ncovars,misspec))

# misBoth
sim.results.CIW.avg.misBoth <- ggplot(data=CIW.df[CIW.df$misspec=="misBoth",],
                                 aes(x=Estimator, y=CIW, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~  gamma.setting, scales = "free", labeller=labeller3)   + ylab("Average confidence interval width over all pairwise comparisons") +  
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.CIW.avg <- ggplotGrob(sim.results.CIW.avg.misBoth)

# Get the positions of the strips in the gtable: t = top, l = left, ...
posR <- subset(z.CIW.avg$layout, grepl("strip-r", name), select = t:r)
posT <- subset(z.CIW.avg$layout, grepl("strip-t", name), select = t:r)

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.CIW.avg$widths[max(posR$r)]    # width of current right strips
height <- z.CIW.avg$heights[min(posT$t)]  # height of current top strips

z.CIW.avg <- gtable_add_cols(z.CIW.avg, width, max(posR$r))  
z.CIW.avg <- gtable_add_rows(z.CIW.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.CIW.avg <- gtable_add_grob(z.CIW.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.CIW.avg <- gtable_add_grob(z.CIW.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.CIW.avg <- gtable_add_cols(z.CIW.avg, unit(1/5, "line"), max(posR$r))
z.CIW.avg <- gtable_add_rows(z.CIW.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.CIW.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_CIW_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misBoth","_R_",results[[1]][[1]]$R,".png"),plot = z.CIW.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_CIW_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misBoth","_R_",results[[1]][[1]]$R,".png"),plot = z.CIW.avg,scale=2)
}

# misTreat
sim.results.CIW.avg.misTreat <- ggplot(data=CIW.df[CIW.df$misspec=="misTreat",],
                                      aes(x=Estimator, y=CIW, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~  gamma.setting, scales = "free", labeller=labeller3)   + ylab("Average confidence interval width over all pairwise comparisons") +  
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.CIW.avg <- ggplotGrob(sim.results.CIW.avg.misTreat)

# Get the positions of the strips in the gtable: t = top, l = left, ...
posR <- subset(z.CIW.avg$layout, grepl("strip-r", name), select = t:r)
posT <- subset(z.CIW.avg$layout, grepl("strip-t", name), select = t:r)

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.CIW.avg$widths[max(posR$r)]    # width of current right strips
height <- z.CIW.avg$heights[min(posT$t)]  # height of current top strips

z.CIW.avg <- gtable_add_cols(z.CIW.avg, width, max(posR$r))  
z.CIW.avg <- gtable_add_rows(z.CIW.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.CIW.avg <- gtable_add_grob(z.CIW.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.CIW.avg <- gtable_add_grob(z.CIW.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.CIW.avg <- gtable_add_cols(z.CIW.avg, unit(1/5, "line"), max(posR$r))
z.CIW.avg <- gtable_add_rows(z.CIW.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.CIW.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_CIW_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misTreat","_R_",results[[1]][[1]]$R,".png"),plot = z.CIW.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_CIW_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misTreat","_R_",results[[1]][[1]]$R,".png"),plot = z.CIW.avg,scale=2)
}

# misOut
sim.results.CIW.avg.misOut <- ggplot(data=CIW.df[CIW.df$misspec=="misOut",],
                                      aes(x=Estimator, y=CIW, fill=forcats::fct_rev(Estimator)))  + geom_col() +
  facet_grid(overlap.setting ~  gamma.setting, scales = "free", labeller=labeller3)   + ylab("Average confidence interval width over all pairwise comparisons") +  
  scale_fill_manual(values= estimator.color) +
  scale_x_discrete(labels=NULL, limits = rev) +
  theme(legend.position = "none",legend.margin=margin(1,5,5,5), legend.justification="center",
        legend.box.margin=margin(0,0,0,0),legend.text=element_text(size=14), legend.key.width = unit(0.75, "cm"),legend.spacing.x = unit(0.75, 'cm'), legend.spacing.y = unit(0.75, 'cm')) +  
  theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
  theme(axis.title=element_text(family="serif", size=16)) +
  theme(axis.text.y=element_text(family="serif", size=16)) +
  theme(axis.text.x=element_text(family="serif", size=14, angle = 0, vjust = 0.5, hjust=0.25)) +
  theme(legend.text=element_text(family="serif", size=14)) +
  theme(legend.title=element_text(family="serif", size=14)) +
  theme(strip.text.x = element_text(family="serif", size=16)) +
  theme(strip.text.y = element_text(family="serif", size=16)) +
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l =0))) +
  theme(axis.title.x = element_text(margin = margin(t = 20, r = 0, b = 0, l =0))) +
  theme(panel.spacing = unit(1, "lines")) +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

# Get the ggplot grob
z.CIW.avg <- ggplotGrob(sim.results.CIW.avg.misOut)

# Get the positions of the strips in the gtable: t = top, l = left, ...
posR <- subset(z.CIW.avg$layout, grepl("strip-r", name), select = t:r)
posT <- subset(z.CIW.avg$layout, grepl("strip-t", name), select = t:r)

# Add a new column to the right of current right strips, 
# and a new row on top of current top strips
width <- z.CIW.avg$widths[max(posR$r)]    # width of current right strips
height <- z.CIW.avg$heights[min(posT$t)]  # height of current top strips

z.CIW.avg <- gtable_add_cols(z.CIW.avg, width, max(posR$r))  
z.CIW.avg <- gtable_add_rows(z.CIW.avg, height, min(posT$t)-1)

# Position the grobs in the gtable
z.CIW.avg <- gtable_add_grob(z.CIW.avg, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
z.CIW.avg <- gtable_add_grob(z.CIW.avg, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")

# Add small gaps between strips
z.CIW.avg <- gtable_add_cols(z.CIW.avg, unit(1/5, "line"), max(posR$r))
z.CIW.avg <- gtable_add_rows(z.CIW.avg, unit(1/5, "line"), min(posT$t))

# Draw it
grid.newpage()
grid.draw(z.CIW.avg)

if(use.SL){
  ggsave(paste0("sim_results/static_simulation_CIW_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misOut","_R_",results[[1]][[1]]$R,".png"),plot = z.CIW.avg,scale=1.75)
}else{
  ggsave(paste0("sim_results/static_simulation_CIW_avg_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_misspec_misOut","_R_",results[[1]][[1]]$R,".png"),plot = z.CIW.avg,scale=2)
}