combine_sim_plots_highdim.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/20231016", "outputs/20231019", "outputs/20231030", 6, 'TRUE', 'TRUE', 'TRUE') # 40 and 100 covars
output.path <- as.character(args[1]) # args <- c("outputs/20231109", "outputs/20231110", "outputs/20231114", 6, 'TRUE', 'FALSE', 'TRUE') # just 100 covars
output.path2 <- as.character(args[2])
output.path3 <- as.character(args[3])
J <- as.numeric(args[4])
use.SL <- as.character(args[5])
covars40 <- as.character(args[6])
covars100 <- as.character(args[7])

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

overlap.setting <- c("adequate","inadequate","rct")
if(covars100){
  gamma.setting <- c("zero","low") 
}else{
  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_100_"),filenames)] # keep only highdem runs (R=100)

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))
  print(ncovars)
  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)
  }
}

# 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" 

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,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, 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"))  # 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")  

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)) %>% 
  group_by(Estimator,overlap.setting,gamma.setting,ncovars) %>% 
  mutate(CP = mean(Coverage)) %>%
  dplyr::slice(n()) %>%
  select(c(CP,Estimator,overlap.setting,gamma.setting,ncovars))

if(covars100){
  # 100 covariates
  sim.results.cp.avg.100 <- ggplot(data=coverage.df[coverage.df$ncovars==100,],
                                   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_discrete(name = "") +
    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.100)
  
  # 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,"_covars40_FALSE","_covars100_",covars100,"_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,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.cp.avg,scale=2)
  }
}else{
  
  # 40 covariates
  sim.results.cp.avg.40 <- ggplot(data=coverage.df[coverage.df$ncovars==40,],
                                  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_discrete(name = "") +
    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.40)
  
  # 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,"_covars40_",covars40,"_covars100_FALSE","_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,"_covars40_",covars40,"_covars100_FALSE","_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)) %>% 
  group_by(Estimator,overlap.setting,gamma.setting,ncovars) %>% 
  mutate(bias = mean(abs.bias)) %>%
  dplyr::slice(n()) %>%
  select(c(bias,Estimator,overlap.setting,gamma.setting,ncovars))

if(covars40){
  # 40 covars
  sim.results.bias.avg.40 <- ggplot(data=bias.df[bias.df$ncovars==40,],
                                    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_discrete(name = "") +
    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.40)
  
  # 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,"_covars40_",covars40,"_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,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.bias.avg,scale=2)
  }
}else{
  # 100 covars
  sim.results.bias.avg.100 <- ggplot(data=bias.df[bias.df$ncovars==100,],
                                     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_discrete(name = "") +
    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.100)
  
  # 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,"_covars40_FALSE","_covars100_",covars100,"_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, "_covars40_FALSE","_covars100_",covars100,"_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)) %>% 
  group_by(Estimator,overlap.setting,gamma.setting,ncovars) %>% 
  mutate(CIW = mean(CIW)) %>%
  dplyr::slice(n()) %>%
  select(c(CIW,Estimator,overlap.setting,gamma.setting,ncovars))

if(covars40){
  # 40 covars
  sim.results.CIW.avg.40 <- ggplot(data=CIW.df[CIW.df$ncovars==40,],
                                   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_discrete(name = "") +
    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.40)
  
  # 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,"_covars40_",covars40,"_covars100_FALSE_","_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,"_covars40_",covars40,"_covars100_FALSE_","_R_",results[[1]][[1]]$R,".png"),plot = z.CIW.avg,scale=2)
  }
}else{
  #100 covars
  sim.results.CIW.avg.100 <- ggplot(data=CIW.df[CIW.df$ncovars==100,],
                                    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_discrete(name = "") +
    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.100)
  
  # 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,"_covars40_FALSE","_covars100_",covars100,"_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,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.CIW.avg,scale=2)
  }
}
if(use.SL){
  
  yinitial_tmle <- NULL
  Ahat_tmle <-NULL
  yhat_tmle <- NULL
  ess_tmle <- NULL
  
  yinitial_tmle_bin <- NULL
  Ahat_tmle_bin <- NULL
  yhat_tmle_bin <- NULL
  ess_tmle_bin <- NULL
  
  yinitial_tmle_glm <- NULL
  Ahat_tmle_glm <-NULL
  yhat_tmle_glm <- NULL
  ess_tmle_glm <- NULL
  
  yinitial_tmle_glm_bin <- NULL
  Ahat_tmle_glm_bin <- NULL
  yhat_tmle_glm_bin <- NULL
  ess_tmle_glm_bin <- NULL
  
  ## plot estimated and observed Ys and As, initial estimates, and ESS (GLM + SL)
  
  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("R_100_"),filenames)] # keep only highdem runs (R=100)
  
  if(any( duplicated(substring(filenames, 18)))){
    print("removing duplicate filenames")
    filenames <- filenames[-which(duplicated(substring(filenames, 18)))]
  }
  
  results.AY <- list() # structure is: [[filename]][[metric]]
  for(f in filenames){
    print(f)
    result.AY.matrix <- readRDS(f)
    R <- ncol(result.AY.matrix )
    ncovars <- ifelse(grep("covars40_TRUE",f), 40, ifelse(grep("covars100_TRUE",f), 100, 6))
    if(isTRUE(grep("use_SL_FALSE",f)==1)){
      estimator <- estimators[grep("glm", estimators)]
      rownames(result.AY.matrix)[grep("tmle",rownames(result.AY.matrix))] <- gsub("tmle","tmle_glm", rownames(result.AY.matrix)[grep("tmle",rownames(result.AY.matrix))])
      rownames(result.AY.matrix)[grep("gcomp",rownames(result.AY.matrix))] <- gsub("gcomp","gcomp_glm", rownames(result.AY.matrix)[grep("gcomp",rownames(result.AY.matrix))])
      rownames(result.AY.matrix)[grep("_iptw",rownames(result.AY.matrix))] <- gsub("_iptw","_iptw_glm", rownames(result.AY.matrix)[grep("_iptw",rownames(result.AY.matrix))])
      rownames(result.AY.matrix)[grep("_aiptw",rownames(result.AY.matrix))] <- gsub("_aiptw","_aiptw_glm", rownames(result.AY.matrix)[grep("_aiptw",rownames(result.AY.matrix))])
    }else{
      estimator <- setdiff(estimators,estimators[grep("glm", estimators)])
    }
    
    obs_outcome <- t(sapply(result.AY.matrix["obs_outcome",],cbind))
    obs_treatment <- t(sapply(result.AY.matrix["obs_treatment",],cbind))
    obs_covariates <- t(sapply(result.AY.matrix["obs_covariates",],cbind))
    if(estimand=="ate"){
      true_taus <- t(sapply(result.AY.matrix["trueATE",],cbind))
    } else if(estimand=="att"){
      true_taus <- t(sapply(result.AY.matrix["trueATT",],cbind))
    }
    
    if(isTRUE(grep("use_SL_FALSE",f)==1)){
      yinitial_tmle_glm <- t(sapply(lapply(result.AY.matrix["yinitial_tmle_glm",],colMeans),cbind))
      Ahat_tmle_glm <-t(sapply(lapply(1:R,function(i) colMeans(result.AY.matrix["Ahat_tmle_glm",][[i]])),cbind))
      yhat_tmle_glm <- t(sapply(result.AY.matrix["yhat_tmle_glm",],cbind))
      ess_tmle_glm <- t(sapply(lapply(1:R,function(i) result.AY.matrix["ess_ate_tmle_glm",][[i]]),cbind))
      yinitial_tmle_glm_bin <- yinitial_tmle_glm # same
      Ahat_tmle_glm_bin <- t(sapply(lapply(result.AY.matrix["Ahat_tmle_glm_bin",],colMeans),cbind))
      yhat_tmle_glm_bin <- t(sapply(result.AY.matrix["yhat_tmle_glm_bin",],cbind))
      ess_tmle_glm_bin <- t(sapply(lapply(1:R,function(i) result.AY.matrix["ess_ate_tmle_glm_bin",][[i]]),cbind))
    }else{
      yinitial_tmle <- t(sapply(lapply(result.AY.matrix["yinitial_tmle",],colMeans),cbind))
      Ahat_tmle <-t(sapply(lapply(1:R,function(i) colMeans(result.AY.matrix["Ahat_tmle",][[i]])),cbind))
      yhat_tmle <- t(sapply(result.AY.matrix["yhat_tmle",],cbind))
      ess_tmle <- t(sapply(lapply(1:R,function(i) result.AY.matrix["ess_ate_tmle",][[i]]),cbind))
      yinitial_tmle_bin <- yinitial_tmle # same
      Ahat_tmle_bin <- t(sapply(lapply(result.AY.matrix["Ahat_tmle_bin",],colMeans),cbind))
      yhat_tmle_bin <- t(sapply(result.AY.matrix["yhat_tmle_bin",],cbind))
      ess_tmle_bin <- t(sapply(lapply(1:R,function(i) result.AY.matrix["ess_ate_tmle_bin",][[i]]),cbind))
    }
    
    results.AY[[f]] <- list("obs_outcome"=obs_outcome,"obs_treatment"=obs_treatment,"obs_covariates"=obs_covariates,"true_taus"=true_taus,
                            "yinitial_tmle"=yinitial_tmle,"Ahat_tmle"=Ahat_tmle,"yhat_tmle"= yhat_tmle,"ess_tmle"=ess_tmle,
                            "Ahat_tmle_bin"=Ahat_tmle_bin,"yhat_tmle_bin"= yhat_tmle_bin,"yinitial_tmle_bin"= yinitial_tmle_bin,"ess_tmle_bin"=ess_tmle_bin,
                            "yinitial_tmle_glm"=yinitial_tmle_glm,"Ahat_tmle_glm"=Ahat_tmle_glm,"yhat_tmle_glm"= yhat_tmle_glm,"ess_tmle_glm"=ess_tmle_glm,
                            "Ahat_tmle_glm_bin"=Ahat_tmle_glm_bin,"yhat_tmle_glm_bin"= yhat_tmle_glm_bin,"yinitial_tmle_glm_bin"= yinitial_tmle_glm_bin,"ess_tmle_glm_bin"=ess_tmle_glm_bin,
                            "R"=R)
  }
  
  # Create New lists
  # structure is: [[filename]]
  
  obs_outcome <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$obs_outcome)
  obs_treatment <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$obs_treatment)
  obs_treatment_n <- sapply(1:length(obs_treatment), function(i) obs_treatment[[i]]*n)
  obs_covariates <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$obs_covariates)
  true_taus <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$true_taus)
  yinitial_tmle <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$yinitial_tmle)
  ess_tmle <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$ess_tmle)
  Ahat_tmle <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$Ahat_tmle)
  yinitial_tmle_bin <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$yinitial_tmle_bin)
  ess_tmle_bin <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$ess_tmle_bin)
  Ahat_tmle_bin <- lapply(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), function(f) results.AY[[f]]$Ahat_tmle_bin)
  yinitial_tmle_glm <- lapply(setdiff(filenames,filenames[grep("use_SL_TRUE",filenames)]), function(f) results.AY[[f]]$yinitial_tmle_glm)
  ess_tmle_glm <- lapply(setdiff(filenames,filenames[grep("use_SL_TRUE",filenames)]), function(f) results.AY[[f]]$ess_tmle_glm)
  Ahat_tmle_glm <- lapply(setdiff(filenames,filenames[grep("use_SL_TRUE",filenames)]), function(f) results.AY[[f]]$Ahat_tmle_glm)
  yinitial_tmle_glm_bin <- lapply(setdiff(filenames,filenames[grep("use_SL_TRUE",filenames)]), function(f) results.AY[[f]]$yinitial_tmle_glm_bin)
  ess_tmle_glm_bin <- lapply(setdiff(filenames,filenames[grep("use_SL_TRUE",filenames)]), function(f) results.AY[[f]]$ess_tmle_glm_bin)
  Ahat_tmle_glm_bin <- lapply(setdiff(filenames,filenames[grep("use_SL_TRUE",filenames)]), function(f) results.AY[[f]]$Ahat_tmle_glm_bin)
  
  # Create dataframe for plot
  results.AY.df <- data.frame("Obs.A"=c(unlist(obs_treatment)),
                              "Obs.Y"=c(unlist(obs_outcome)),
                              "filename"=rep(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]), each=J*R)) # 54000
  
  results.AY.df$ncovars <- "6"
  results.AY.df$ncovars[grep("covars40_TRUE",results.AY.df$filename)] <- "40" 
  results.AY.df$ncovars[grep("covars100_TRUE",results.AY.df$filename)] <- "100" 
  
  results.AY.preds.df <- data.frame("Y.initial"=c(unlist(yinitial_tmle),unlist(yinitial_tmle_bin),unlist(yinitial_tmle_glm),unlist(yinitial_tmle_glm_bin)),
                                    "Y.initial.diff"=c(abs(unlist(yinitial_tmle)-unlist(obs_outcome)),abs(unlist(yinitial_tmle_bin)-unlist(obs_outcome)),abs(unlist(yinitial_tmle_glm)-unlist(obs_outcome)),abs(unlist(yinitial_tmle_glm_bin)-unlist(obs_outcome))),
                                    "A.mean"= c(unlist(Ahat_tmle),unlist(Ahat_tmle_bin),unlist(Ahat_tmle_glm),unlist(Ahat_tmle_glm_bin)),
                                    "A.mean.diff"= c(abs(unlist(Ahat_tmle)-unlist(obs_treatment)),abs(unlist(Ahat_tmle_bin)-unlist(obs_treatment)),abs(unlist(Ahat_tmle_glm)-unlist(obs_treatment)),abs(unlist(Ahat_tmle_glm_bin)-unlist(obs_treatment))),
                                    "ESS"= c(unlist(ess_tmle),unlist(ess_tmle_bin),unlist(ess_tmle_glm),unlist(ess_tmle_glm_bin)),
                                    "ESS_ratio"= c(unlist(ess_tmle)/unlist(obs_treatment_n),unlist(ess_tmle_bin)/unlist(obs_treatment_n),unlist(ess_tmle_glm)/unlist(obs_treatment_n),unlist(ess_tmle_glm_bin)/unlist(obs_treatment_n)),
                                    "Estimator"=c(rep("Multinomial (SL)",length.out=length(unlist(Ahat_tmle))), rep("Binomial (SL)",length.out=length(unlist(Ahat_tmle_bin))), rep("Multinomial (GLM)",length.out=length(unlist(Ahat_tmle_glm))), rep("Binomial (GLM)",length.out=length(unlist(Ahat_tmle_glm_bin)))),
                                    filename=c(unlist(sapply(1:length(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)])), function(i) rep(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)])[i], length.out=R*J))))) # (100*9*6) *5 = (54000)*5 = 22800
  
  results.AY.preds.df$ncovars <- "6"
  results.AY.preds.df$ncovars[grep("covars40_TRUE",results.AY.preds.df$filename)] <- "40" 
  results.AY.preds.df$ncovars[grep("covars100_TRUE",results.AY.preds.df$filename)] <- "100" 
  
  results.AY.df$J <- ifelse(J==3,3,6)
  results.AY.preds.df$J <- ifelse(J==3,3,6)
  
  results.AY.df$Treatment <- rep(rep(1:J, each=R), length(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)])))
  results.AY.preds.df$Treatment <- rep(rep(rep(1:J, each=R), length(setdiff(filenames,filenames[grep("use_SL_FALSE",filenames)]))), 2)
  
  results.AY.df$overlap.setting <- NA
  results.AY.preds.df$overlap.setting <- NA
  
  results.AY.df$gamma.setting <- NA
  results.AY.preds.df$gamma.setting <- NA
  
  for(s in overlap.setting){
    print(s)
    results.AY.df[grep(s, results.AY.df$filename),]$overlap.setting <- s
    results.AY.preds.df[grep(s, results.AY.preds.df$filename),]$overlap.setting <- s
  }
  
  for(s in gamma.setting){
    print(s)
    results.AY.df[grep(s, results.AY.df$filename),]$gamma.setting <- s
    results.AY.preds.df[grep(s, results.AY.preds.df$filename),]$gamma.setting <- s
  }
  
  # reshape and plot
  results.AY.df$id <- with(results.AY.df, paste(overlap.setting, gamma.setting,ncovars, J, sep = "_"))
  results_AY_long <- reshape2::melt(results.AY.df[!colnames(results.AY.df) %in% c("J","id","filename")], id.vars=c("Treatment","overlap.setting","gamma.setting","ncovars"))  # convert to long format
  
  results.AY.preds.df$id <- with(results.AY.preds.df, paste(overlap.setting, gamma.setting, ncovars, J, sep = "_"))
  results_AY_preds_long <- reshape2::melt(results.AY.preds.df[!colnames(results.AY.preds.df) %in% c("J","id","filename")], id.vars=c("Treatment","Estimator","overlap.setting","gamma.setting","ncovars"))  # convert to long format
  
  if(J==6){
    x.labels.AY <- c('1'=parse(text = TeX(paste0('$1$'))),
                     '2'=parse(text = TeX(paste0('$2$'))),
                     '3'=parse(text = TeX(paste0('$3$'))),
                     '4'=parse(text = TeX(paste0('$4$'))),
                     '5'=parse(text = TeX(paste0('$5$'))),
                     '6'=parse(text = TeX(paste0('$6$'))))
  } else if(J==3){
    x.labels.AY <- c('1'=parse(text = TeX(paste0('$1$'))),
                     '2'=parse(text = TeX(paste0('$2$'))),
                     '3'=parse(text = TeX(paste0('$3$'))))
  }
  
  if(covars40){
    # Observed treatment probabilities (40 covariates)
    sim.results.A.40 <- ggplot(data=results_AY_long[results_AY_long$variable=="Obs.A" & results_AY_long$ncovars==40,],
                               aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(variable)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Simulated treatment probability") + #ggtitle(paste0("Observed treatment probabilities (J=",J,", n=", n,")")) +
      scale_fill_manual(name = "Variable:", values = c("orange")) + 
      scale_x_discrete(labels=x.labels.AY,
                       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=16)) +
      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"))
    
    # Get the ggplot grob
    z.A <- ggplotGrob(sim.results.A.40)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.A$layout, grepl("strip-r", 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.A$widths[max(posR$r)]    # width of current right strips
    
    z.A <- gtable_add_cols(z.A, width, max(posR$r))  
    
    # Position the grobs in the gtable
    z.A <- gtable_add_grob(z.A, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.A <- gtable_add_cols(z.A, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.A)
    
    ggsave(paste0("sim_results/static_simulation_obs_A_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.A,scale=1.75)
  }else{
    # Observed treatment probabilities (100 covariates)
    sim.results.A.100 <- ggplot(data=results_AY_long[results_AY_long$variable=="Obs.A" & results_AY_long$ncovars==100,],
                                aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(variable)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Simulated treatment probability") + #ggtitle(paste0("Observed treatment probabilities (J=",J,", n=", n,")")) +
      scale_fill_manual(name = "Variable:", values = c("orange")) + 
      scale_x_discrete(labels=x.labels.AY,
                       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=16)) +
      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"))
    
    # Get the ggplot grob
    z.A <- ggplotGrob(sim.results.A.100)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.A$layout, grepl("strip-r", 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.A$widths[max(posR$r)]    # width of current right strips
    
    z.A <- gtable_add_cols(z.A, width, max(posR$r))  
    
    # Position the grobs in the gtable
    z.A <- gtable_add_grob(z.A, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.A <- gtable_add_cols(z.A, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.A)
    
    ggsave(paste0("sim_results/static_simulation_obs_A_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.A,scale=1.75)
  }
  if(covars40){
    # observed outcomes (40 covariates)
    sim.results.Y.40 <- ggplot(data=results_AY_long[results_AY_long$variable=="Obs.Y" & results_AY_long$ncovars==40,],
                               aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(variable)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~  gamma.setting, scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Simulated outcome under each treatment level") + #ggtitle(paste0("Observed outcomes under each treatment level (J=",J,", n=", n,")")) +
      scale_fill_manual(name = "Variable:", values = c("orange")) + 
      scale_x_discrete(labels=x.labels.AY,
                       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=16)) +
      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"))
    
    # Get the ggplot grob
    z.Y <- ggplotGrob(sim.results.Y.40)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posT <- subset(z.Y$layout, grepl("strip-t", name), select = t:r)
    posR <- subset(z.Y$layout, grepl("strip-r", 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
    height <- z.Y$heights[min(posT$t)]  # height of current top strips
    width <- z.Y$widths[max(posR$r)]    # width of current right strips
    
    z.Y <- gtable_add_cols(z.Y, width, max(posR$r))  
    z.Y <- gtable_add_rows(z.Y, height, min(posT$t)-1)
    
    # Position the grobs in the gtable
    z.Y <- gtable_add_grob(z.Y, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
    z.Y <- gtable_add_grob(z.Y, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")
    
    # Add small gaps between strips
    z.Y <- gtable_add_cols(z.Y, unit(1/5, "line"), max(posR$r))
    z.Y <- gtable_add_rows(z.Y, unit(1/5, "line"), min(posT$t))
    
    # Draw it
    grid.newpage()
    grid.draw(z.Y)
    
    ggsave(paste0("sim_results/static_simulation_Y_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.Y,scale=1.75)
  }else{
    # observed outcomes (100 covariates)
    sim.results.Y.100 <- ggplot(data=results_AY_long[results_AY_long$variable=="Obs.Y" & results_AY_long$ncovars==100,],
                                aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(variable)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~  gamma.setting, scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Simulated outcome under each treatment level") + #ggtitle(paste0("Observed outcomes under each treatment level (J=",J,", n=", n,")")) +
      scale_fill_manual(name = "Variable:", values = c("orange")) + 
      scale_x_discrete(labels=x.labels.AY,
                       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=16)) +
      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"))
    
    # Get the ggplot grob
    z.Y <- ggplotGrob(sim.results.Y.100)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posT <- subset(z.Y$layout, grepl("strip-t", name), select = t:r)
    posR <- subset(z.Y$layout, grepl("strip-r", 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
    height <- z.Y$heights[min(posT$t)]  # height of current top strips
    width <- z.Y$widths[max(posR$r)]    # width of current right strips
    
    z.Y <- gtable_add_cols(z.Y, width, max(posR$r))  
    z.Y <- gtable_add_rows(z.Y, height, min(posT$t)-1)
    
    # Position the grobs in the gtable
    z.Y <- gtable_add_grob(z.Y, stripR, t = min(posR$t)+1, l = max(posR$r) + 1, b = max(posR$b)+1, name = "strip-right")
    z.Y <- gtable_add_grob(z.Y, stripT, t = min(posT$t), l = min(posT$l), r = max(posT$r), name = "strip-top")
    
    # Add small gaps between strips
    z.Y <- gtable_add_cols(z.Y, unit(1/5, "line"), max(posR$r))
    z.Y <- gtable_add_rows(z.Y, unit(1/5, "line"), min(posT$t))
    
    # Draw it
    grid.newpage()
    grid.draw(z.Y)
    
    ggsave(paste0("sim_results/static_simulation_Y_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.Y,scale=1.75)
  }
  
  if(covars40){
    # treatment estimates (40 covariates)
    sim.results.A.mean.40 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="A.mean" & results_AY_preds_long$ncovars==40,],  
                                    aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Estimated treatment probability") + #ggtitle(paste0("Estimated treatment probabilities (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Treatment model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       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(1, 'cm'), legend.spacing.y = unit(1, '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=16)) +
      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"))
    
    # Get the ggplot grob
    z.A.mean <- ggplotGrob(sim.results.A.mean.40)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.A.mean$layout, grepl("strip-r", 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.A.mean$widths[max(posR$r)]    # width of current right strips
    
    z.A.mean <- gtable_add_cols(z.A.mean, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.A.mean <- gtable_add_grob(z.A.mean, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.A.mean <- gtable_add_cols(z.A.mean, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.A.mean)
    
    ggsave(paste0("sim_results/static_simulation_est_A_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.A.mean,scale=1.75)
  }else{
    # treatment estimates (100 covariates)
    sim.results.A.mean.100 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="A.mean" & results_AY_preds_long$ncovars==100,],  
                                     aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Estimated treatment probability") + #ggtitle(paste0("Estimated treatment probabilities (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Treatment model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       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(1, 'cm'), legend.spacing.y = unit(1, '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=16)) +
      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"))
    
    # Get the ggplot grob
    z.A.mean <- ggplotGrob(sim.results.A.mean.100)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.A.mean$layout, grepl("strip-r", 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.A.mean$widths[max(posR$r)]    # width of current right strips
    
    z.A.mean <- gtable_add_cols(z.A.mean, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.A.mean <- gtable_add_grob(z.A.mean, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.A.mean <- gtable_add_cols(z.A.mean, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.A.mean)
    
    ggsave(paste0("sim_results/static_simulation_est_A_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.A.mean,scale=1.75)
  }
  if(covars40){
    # treatment estimates -diff (40 covariates)
    sim.results.A.mean.diff.40 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="A.mean.diff" & results_AY_preds_long$ncovars==40,],  
                                         aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "free", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Absolute difference between estimated and observed treatment probabilities") + #ggtitle(paste0("Accuracy of treatment estimation (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Treatment model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       limits = rev) +
      #coord_cartesian(ylim=c(0,0.0075)) +
      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(1, 'cm'), legend.spacing.y = unit(1, '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=16)) +
      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"))
    
    # Get the ggplot grob
    z.A.mean.diff <- ggplotGrob(sim.results.A.mean.diff.40)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.A.mean.diff$layout, grepl("strip-r", 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.A.mean.diff$widths[max(posR$r)]    # width of current right strips
    
    z.A.mean.diff <- gtable_add_cols(z.A.mean.diff, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.A.mean.diff <- gtable_add_grob(z.A.mean.diff, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.A.mean.diff <- gtable_add_cols(z.A.mean.diff, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.A.mean.diff)
    
    ggsave(paste0("sim_results/static_simulation_est_A_diff_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.A.mean.diff,scale=1.75)
  }else{
    # treatment estimates -diff (100 covariates)
    sim.results.A.mean.diff.100 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="A.mean.diff" & results_AY_preds_long$ncovars==100,],  
                                          aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "free", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Absolute difference between estimated and observed treatment probabilities") + #ggtitle(paste0("Accuracy of treatment estimation (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Treatment model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       limits = rev) +
      #coord_cartesian(ylim=c(0,0.0075)) +
      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(1, 'cm'), legend.spacing.y = unit(1, '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=16)) +
      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"))
    
    # Get the ggplot grob
    z.A.mean.diff <- ggplotGrob(sim.results.A.mean.diff.100)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.A.mean.diff$layout, grepl("strip-r", 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.A.mean.diff$widths[max(posR$r)]    # width of current right strips
    
    z.A.mean.diff <- gtable_add_cols(z.A.mean.diff, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.A.mean.diff <- gtable_add_grob(z.A.mean.diff, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.A.mean.diff <- gtable_add_cols(z.A.mean.diff, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.A.mean.diff)
    
    ggsave(paste0("sim_results/static_simulation_est_A_diff_estimand_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.A.mean.diff,scale=1.75)
  }
  
  if(covars40){
    # initial outcome estimates (40 covariates)
    sim.results.Y.initial.40 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="Y.initial" & results_AY_preds_long$ncovars==40,], 
                                       aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ gamma.setting, scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Initial outcome estimate") + #ggtitle(paste0("Initial outcome estimate (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Outcome model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       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(1, 'cm'), legend.spacing.y = unit(1, 'cm')) +    theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
      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=16)) +
      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"))
    
    # Get the ggplot grob
    z.Y.initial <- ggplotGrob(sim.results.Y.initial.40)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.Y.initial$layout, grepl("strip-r", 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.Y.initial$widths[max(posR$r)]    # width of current right strips
    
    z.Y.initial <- gtable_add_cols(z.Y.initial, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.Y.initial <- gtable_add_grob(z.Y.initial, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.Y.initial <- gtable_add_cols(z.Y.initial, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.Y.initial)
    
    ggsave(paste0("sim_results/static_simulation_est_Y_initial_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.Y.initial,scale=1.75)
  }else{
    # initial outcome estimates (100 covariates)
    sim.results.Y.initial.100 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="Y.initial" & results_AY_preds_long$ncovars==100,], 
                                        aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ gamma.setting, scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Initial outcome estimate") + #ggtitle(paste0("Initial outcome estimate (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Outcome model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       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(1, 'cm'), legend.spacing.y = unit(1, 'cm')) +    theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
      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=16)) +
      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"))
    
    # Get the ggplot grob
    z.Y.initial <- ggplotGrob(sim.results.Y.initial.100)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.Y.initial$layout, grepl("strip-r", 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.Y.initial$widths[max(posR$r)]    # width of current right strips
    
    z.Y.initial <- gtable_add_cols(z.Y.initial, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.Y.initial <- gtable_add_grob(z.Y.initial, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.Y.initial <- gtable_add_cols(z.Y.initial, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.Y.initial)
    
    ggsave(paste0("sim_results/static_simulation_est_Y_initial_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.Y.initial,scale=1.75)
  }
  if(covars40){
    # initial outcome estimates - diff (40 covars)
    sim.results.Y.initial.diff.40 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="Y.initial.diff" & results_AY_preds_long$ncovars==40,], 
                                            aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ gamma.setting, scales = "free", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Absolute difference between initial outcome estimate and observed outcome") + #ggtitle(paste0("Accuracy of outcome estimation (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Outcome model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       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(1, 'cm'), legend.spacing.y = unit(1, 'cm')) +    theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
      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=16)) +
      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"))
    
    # Get the ggplot grob
    z.Y.initial.diff <- ggplotGrob(sim.results.Y.initial.diff.40)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.Y.initial.diff$layout, grepl("strip-r", 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.Y.initial.diff$widths[max(posR$r)]    # width of current right strips
    
    z.Y.initial.diff <- gtable_add_cols(z.Y.initial.diff, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.Y.initial.diff <- gtable_add_grob(z.Y.initial.diff, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.Y.initial.diff <- gtable_add_cols(z.Y.initial.diff, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.Y.initial.diff)
    
    ggsave(paste0("sim_results/static_simulation_est_Y_initial_diff_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.Y.initial.diff,scale=1.75)
  }else{
    # initial outcome estimates - diff (100 covars)
    sim.results.Y.initial.diff.100 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="Y.initial.diff" & results_AY_preds_long$ncovars==100,], 
                                             aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ gamma.setting, scales = "free", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Absolute difference between initial outcome estimate and observed outcome") + #ggtitle(paste0("Accuracy of outcome estimation (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Outcome model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       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(1, 'cm'), legend.spacing.y = unit(1, 'cm')) +    theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
      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=16)) +
      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"))
    
    # Get the ggplot grob
    z.Y.initial.diff <- ggplotGrob(sim.results.Y.initial.diff.100)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.Y.initial.diff$layout, grepl("strip-r", 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.Y.initial.diff$widths[max(posR$r)]    # width of current right strips
    
    z.Y.initial.diff <- gtable_add_cols(z.Y.initial.diff, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.Y.initial.diff <- gtable_add_grob(z.Y.initial.diff, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.Y.initial.diff <- gtable_add_cols(z.Y.initial.diff, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.Y.initial.diff)
    
    ggsave(paste0("sim_results/static_simulation_est_Y_initial_diff_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.Y.initial.diff,scale=1.75)
  }
  if(covars40){
    # ESS (40 covars)
    sim.results.ESS.40 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="ESS" & results_AY_preds_long$ncovars==40,],
                                 aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Estimated effective sample size (ESS)") + #ggtitle(paste0("Effective sample size (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Treatment model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       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(1, 'cm'), legend.spacing.y = unit(1, 'cm')) +    theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
      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=16)) +
      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"))
    
    # Get the ggplot grob
    z.ESS <- ggplotGrob(sim.results.ESS.40)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.ESS$layout, grepl("strip-r", 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.ESS$widths[max(posR$r)]    # width of current right strips
    
    z.ESS <- gtable_add_cols(z.ESS, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.ESS <- gtable_add_grob(z.ESS, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.ESS <- gtable_add_cols(z.ESS, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.ESS)
    
    ggsave(paste0("sim_results/static_simulation_est_ESS_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.ESS,scale=1.75)
  }else{
    # ESS (100 covars)
    sim.results.ESS.100 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="ESS" & results_AY_preds_long$ncovars==100,],
                                  aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab("Estimated effective sample size (ESS)") + #ggtitle(paste0("Effective sample size (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Treatment model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       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(1, 'cm'), legend.spacing.y = unit(1, 'cm')) +    theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
      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=16)) +
      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"))
    
    # Get the ggplot grob
    z.ESS <- ggplotGrob(sim.results.ESS.100)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.ESS$layout, grepl("strip-r", 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.ESS$widths[max(posR$r)]    # width of current right strips
    
    z.ESS <- gtable_add_cols(z.ESS, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.ESS <- gtable_add_grob(z.ESS, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.ESS <- gtable_add_cols(z.ESS, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.ESS)
    
    ggsave(paste0("sim_results/static_simulation_est_ESS_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_FALSE","_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.ESS,scale=1.75)
  }
  if(covars40){
    # ESS_ratio (40 covariates)
    sim.results.ESS_ratio.40 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="ESS_ratio" & results_AY_preds_long$ncovars==40,],
                                       aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab(TeX(paste0('Estimated effective sample size (ESS) ratio, $ESS_j/n_j$'))) + #ggtitle(paste0("Effective sample size ratio (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Treatment model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       limits = rev) +
      coord_cartesian(ylim=c(0,1)) +
      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(1, 'cm'), legend.spacing.y = unit(1, 'cm')) +    theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
      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=16)) +
      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"))
    
    # Get the ggplot grob
    z.ESS_ratio <- ggplotGrob(sim.results.ESS_ratio.40)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.ESS_ratio$layout, grepl("strip-r", 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.ESS_ratio$widths[max(posR$r)]    # width of current right strips
    
    z.ESS_ratio <- gtable_add_cols(z.ESS_ratio, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.ESS_ratio <- gtable_add_grob(z.ESS_ratio, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.ESS_ratio <- gtable_add_cols(z.ESS_ratio, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.ESS_ratio)
    
    ggsave(paste0("sim_results/static_simulation_est_ESS_ratio_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_",covars40,"_covars100_FALSE","_R_",results[[1]][[1]]$R,".png"),plot = z.ESS_ratio,scale=1.75)
  }else{
    # ESS_ratio (100 covariates)
    sim.results.ESS_ratio.100 <- ggplot(data=results_AY_preds_long[results_AY_preds_long$variable=="ESS_ratio" & results_AY_preds_long$ncovars==100,],
                                        aes(x=factor(Treatment), y=value,fill=forcats::fct_rev(Estimator)))  + geom_boxplot(outlier.alpha = 0.3,outlier.size = 1, outlier.stroke = 0.1, lwd=0.25) +
      facet_grid(overlap.setting ~ ., scales = "fixed", labeller=labeller3)  +  xlab("Treatment level (j)") + ylab(TeX(paste0('Estimated effective sample size (ESS) ratio, $ESS_j/n_j$'))) + #ggtitle(paste0("Effective sample size ratio (J=",J,", n=", n,")")) +
      #scale_fill_discrete(name = "Treatment model: ") +
      scale_fill_manual(values= c("#F8766D", "#FF6C90", "#A3A500", "#D89000")) +
      scale_x_discrete(labels=x.labels.AY,
                       limits = rev) +
      coord_cartesian(ylim=c(0,1)) +
      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(1, 'cm'), legend.spacing.y = unit(1, 'cm')) +    theme(plot.title = element_text(hjust = 0.5, family="serif", size=16)) +
      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=16)) +
      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"))
    
    # Get the ggplot grob
    z.ESS_ratio <- ggplotGrob(sim.results.ESS_ratio.100)
    
    # Get the positions of the strips in the gtable: t = top, l = left, ...
    posR <- subset(z.ESS_ratio$layout, grepl("strip-r", 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.ESS_ratio$widths[max(posR$r)]    # width of current right strips
    
    z.ESS_ratio <- gtable_add_cols(z.ESS_ratio, width, max(posR$r))
    
    # Position the grobs in the gtable
    z.ESS_ratio <- gtable_add_grob(z.ESS_ratio, stripR, t = min(posR$t), l = max(posR$r) + 1, b = max(posR$b), name = "strip-right")
    
    # Add small gaps between strips
    z.ESS_ratio <- gtable_add_cols(z.ESS_ratio, unit(1/5, "line"), max(posR$r))
    
    # Draw it
    grid.newpage()
    grid.draw(z.ESS_ratio)
    
    ggsave(paste0("sim_results/static_simulation_est_ESS_ratio_",estimand,"_J_",J,"_n_",n,"_outcome_",outcome.type,"_use_SL_",use.SL,"_covars40_FALSE",covars40,"_covars100_",covars100,"_R_",results[[1]][[1]]$R,".png"),plot = z.ESS_ratio,scale=1.75)
  }
}