tmle_MultinomialTrts.R

############################################################################################
# Static setting (T=1) simulations: Compare multinomial TMLE with binary TMLE             #
############################################################################################

# load utils
source('./src/data_generation.R')
source('./src/tmle_calculation.R')
source('./src/tmleContrast.R')
source('./src/misc_fns.R')

#######################
# Simulation function #
######################

staticSim <- function(r, J, n, gbound, ybound, n.folds, overlap.setting, gamma.setting, outcome.type, target.gwt, use.SL, scale.continuous, covars40, covars100, misTreatment, misOutcome, misBoth){
  
  library(purrr)
  library(origami)
  library(sl3)
  library(nnet)
  library(ranger)
  library(xgboost)
  library(glmnet)
  library(Rsolnp)
  library(MASS)
  library(VGAM)
  library(extraDistr)
  
  print(paste("This is simulation run number",r, "\n"))
  
  if(!J%in%c(3,6)){
    stop("J must be 3 or 6")
  }
  
  if(J==3 & (covars40 | covars100)){
    stop("J must be 6 if covars40 or covars100 is TRUE")
  }
  
  ## Generate data
  generated.data <- generateData(r, J, n, overlap.setting, gamma.setting, outcome.type, scale.continuous, covars40, covars100)
  true.ates <- generated.data$trueATE
  obs.treatment <- generated.data$observed.treatment
  
  C <- generated.data$C # contrast matrix
  
  Y <- generated.data$data$Y # extract outcome
  
  A <- as.factor(generated.data$data$Z) # extract treatment
  
  # define a dataframe of covariates
  L <- generated.data$data[!colnames(generated.data$data)%in%c("Y","Z")] 
  rm(generated.data)
  
  # store mean observed outcome under treatment, mean treatment, and true ATEs
  obs_outcome <- colSums(Y*obs.treatment)/colSums(obs.treatment)
  obs_treatment <- colMeans(obs.treatment)
  obs_covariates <- colMeans(L)
  
  ## Manual TMLE (ours)
  
  # stack learners into a model
  learner_stack_A <- make_learner_stack(list("Lrnr_xgboost",nrounds=20, objective="multi:softprob", eval_metric="mlogloss",num_class=J), list("Lrnr_ranger",num.trees=100),list("Lrnr_ranger",num.trees=500), list("Lrnr_glmnet",nfolds = n.folds,alpha = 1, family = "multinomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.25, family = "multinomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.5, family = "multinomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.75, family = "multinomial"))  
  learner_stack_A_bin <- make_learner_stack(list("Lrnr_xgboost",nrounds=20, objective = "reg:logistic"), list("Lrnr_ranger",num.trees=100),list("Lrnr_ranger",num.trees=500), list("Lrnr_glmnet",nfolds = n.folds,alpha = 1, family = "binomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.25, family = "binomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.5, family = "binomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.75, family = "binomial"))  
  
  # metalearner defaults 
  if(outcome.type=="continuous"){
    metalearner_Y <- make_learner(Lrnr_solnp,learner_function=metalearner_linear,eval_function=loss_squared_error) # nonlinear optimization via augmented Lagrange
    learner_stack_Y <- make_learner_stack(list("Lrnr_xgboost",nrounds=20, objective = "reg:squarederror"), list("Lrnr_ranger",num.trees=100), list("Lrnr_ranger",num.trees=500), list("Lrnr_glmnet",nfolds = n.folds,alpha = 1, family = "gaussian"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.25, family = "gaussian"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.5, family = "gaussian"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.75, family = "gaussian")) 
  }else if (outcome.type =="binomial"){
    metalearner_Y <- make_learner(Lrnr_solnp,learner_function=metalearner_logistic_binomial,eval_function=loss_loglik_binomial)
    learner_stack_Y <- make_learner_stack(list("Lrnr_xgboost",nrounds=20, objective = "reg:logistic"), list("Lrnr_ranger",num.trees=100), list("Lrnr_ranger",num.trees=500), list("Lrnr_glmnet",nfolds = n.folds,alpha = 1, family = "binomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.25, family = "binomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.5, family = "binomial"), list("Lrnr_glmnet",nfolds = n.folds,alpha = 0.75, family = "binomial")) 
  }
  metalearner_A <- make_learner(Lrnr_solnp,learner_function=metalearner_linear_multinomial, eval_function=loss_loglik_multinomial)
  metalearner_A_bin <- make_learner(Lrnr_solnp,learner_function=metalearner_logistic_binomial,eval_function=loss_loglik_binomial)
  
  ## fit initial outcome model
  
  if(misOutcome | misBoth){
    L.outcome <- L[,1:(ncol(L)-1)]
  }else{
    L.outcome <- L
  }
  
  if(use.SL){
    
    # define task and candidate learners
    initial_model_for_Y_task <- make_sl3_Task(cbind(Y,L.outcome,obs.treatment), covariates = c(colnames(L.outcome),colnames(obs.treatment)), outcome = "Y", outcome_type=outcome.type, 
                                              folds = origami::make_folds(cbind(Y,L.outcome,obs.treatment), fold_fun = folds_vfold, V = n.folds))
    
    # Super Learner algorithm
    
    initial_model_for_Y_sl <- make_learner(Lrnr_sl, # cross-validates base models
                                           learners = learner_stack_Y,
                                           metalearner = metalearner_Y,
                                           keep_extra=FALSE)
    initial_model_for_Y_sl_fit <- initial_model_for_Y_sl$train(initial_model_for_Y_task)
    initial_model_for_Y_preds <- initial_model_for_Y_sl_fit$predict(initial_model_for_Y_task) # predicted probs.
    
    for(j in 1:J){
      newdata <- cbind(Y,L.outcome,matrix(0, nrow=nrow(obs.treatment), ncol=ncol(obs.treatment), dimnames = dimnames(obs.treatment)))
      newdata[paste0("D",j)] <- 1
      assign(paste0("Q",j), initial_model_for_Y_sl_fit$predict(sl3_Task$new(newdata, covariates = c(colnames(L.outcome),colnames(obs.treatment)), outcome = "Y", outcome_type="binomial")))
    }
  } else{
    if(outcome.type=="continuous"){
      initial_model_for_Y <- glm(Y~.,data=cbind(Y,obs.treatment[,-2],L.outcome), family=gaussian(), control = glm.control(maxit = 100))
      initial_model_for_Y_preds <- predict(initial_model_for_Y)
    }else{
      initial_model_for_Y <- glm(Y~.,data=cbind(Y,obs.treatment[,-2],L.outcome), family=binomial(), control = glm.control(maxit = 100))
      initial_model_for_Y_preds <- predict(initial_model_for_Y,type="response") # predicted probs.
    }
    
    for(j in 1:J){
      newdata <- cbind(Y,L.outcome,matrix(0, nrow=nrow(obs.treatment), ncol=ncol(obs.treatment), dimnames = dimnames(obs.treatment)))
      newdata[paste0("D",j)] <- 1
      assign(paste0("Q",j), predict(initial_model_for_Y,newdata=newdata, type="response"))
    }
  }
  rm(newdata)
  
  Qs <- t(do.call(mapply, c(FUN = cbind, mget(paste0("Q", 1:J)))))
  colnames(Qs) <- paste0("Q", 1:J)
  
  QAW <- data.frame(apply(cbind(QA=initial_model_for_Y_preds,Qs), 2, boundProbs, bounds=ybound)) # bound predictions
  
  ##  fit initial treatment model
  
  if(misTreatment | misBoth){
    L.treatment <- L[,1:(ncol(L)-1)]
  }else{
    L.treatment <- L
  }
  
  if(use.SL){
    
    # multinomial
    
    initial_model_for_A_task <- make_sl3_Task(cbind(A,L.treatment), covariates = c(names(L.treatment)), outcome = "A", outcome_type="categorical", 
                                              folds = origami::make_folds(cbind(A,L.treatment), fold_fun = folds_vfold, V = n.folds))
    
    initial_model_for_A_sl <- make_learner(Lrnr_sl, # cross-validates base models
                                           learners = learner_stack_A,
                                           metalearner = metalearner_A,
                                           keep_extra=FALSE)
    initial_model_for_A_sl_fit <- initial_model_for_A_sl$train(initial_model_for_A_task)
    
    g_preds  <- initial_model_for_A_sl_fit$predict(initial_model_for_A_task)  # estimated propensity scores 
    
    g_preds  <-  data.frame(matrix(unlist(lapply(g_preds, unlist)), nrow=length(lapply(g_preds, unlist)), byrow=TRUE)) 
    
    colnames(g_preds) <- colnames(obs.treatment)
    
    # binomial
    
    initial_model_for_A_task_bin <- lapply(1:J, function(j) make_sl3_Task(cbind("A"=obs.treatment[,j],L.treatment), covariates = c(names(L.treatment)), outcome = "A", outcome_type="binomial",
                                                                          folds = origami::make_folds(cbind(A,L.treatment), fold_fun = folds_vfold, V = n.folds)))
    
    initial_model_for_A_sl_bin <- make_learner(Lrnr_sl, # cross-validates base models
                                               learners = learner_stack_A_bin,
                                               metalearner = metalearner_A_bin,
                                               keep_extra=FALSE)
    initial_model_for_A_sl_fit_bin <- lapply(1:J, function(j) initial_model_for_A_sl_bin$train(initial_model_for_A_task_bin[[j]]))
    
    g_preds_bin  <- sapply(1:J, function(j) initial_model_for_A_sl_fit_bin[[j]]$predict(initial_model_for_A_task_bin[[j]]))  # estimated propensity scores 
    
    colnames(g_preds_bin) <- colnames(obs.treatment)
    
  } else{
    
    # multinomial logistic regression
    initial_model_for_A <- nnet::multinom(formula=A ~ ., data=cbind(A, L.treatment), maxit = 500, trace = FALSE) 
    g_preds <- fitted(initial_model_for_A)
    colnames(g_preds) <- colnames(obs.treatment)
    
    # binomial logistic regression
    initial_model_for_A_bin <- lapply(1:J, function(j) glm(formula=treatment ~ ., data=cbind("treatment"=obs.treatment[,j], L.treatment), family="binomial"))
    g_preds_bin <- sapply(1:J, function(j) predict(initial_model_for_A_bin[[j]], type="response"))
    colnames(g_preds_bin) <- colnames(obs.treatment)
  }
  
  g_preds_bounded <- boundProbs(g_preds, gbound) # (trimmed) estimated propensity scores
  g_preds_bin_bounded <- boundProbs(g_preds_bin, gbound)
  
  # outcome estimates on response scale
  
  TMLE <- getTMLE(QAW, g_preds_bounded, obs.treatment, Y, outcome.type, iptw=FALSE, aiptw = FALSE, target.gwt)
  TMLE_bin <- getTMLE(QAW, g_preds_bin_bounded, obs.treatment, Y, outcome.type, iptw=FALSE, aiptw = FALSE, target.gwt)
  
  IPTW <- getTMLE(QAW=NULL, g_preds_bounded, obs.treatment, Y, outcome.type, iptw=TRUE, aiptw = FALSE, target.gwt)
  IPTW_bin <- getTMLE(QAW=NULL, g_preds_bin_bounded, obs.treatment, Y, outcome.type, iptw=TRUE, aiptw = FALSE, target.gwt)
  
  AIPTW <- getTMLE(QAW, g_preds_bounded, obs.treatment, Y, outcome.type, iptw=FALSE, aiptw = TRUE, target.gwt)
  AIPTW_bin <- getTMLE(QAW, g_preds_bin_bounded, obs.treatment, Y, outcome.type, iptw=FALSE, aiptw = TRUE, target.gwt)
  
  # calculate contrasts
  tmle_contrasts <- tmle_contrast(Qstar=TMLE, C, g_preds_bounded, obs.treatment, Y, outcome.type, alpha=0.05, gcomp = FALSE, iptw = FALSE, aiptw = FALSE, multi.adjust=FALSE)
  tmle_contrasts_bin <- tmle_contrast(Qstar=TMLE_bin, C, g_preds_bin_bounded, obs.treatment, Y, outcome.type, alpha=0.05, gcomp = FALSE, iptw = FALSE, aiptw = FALSE, multi.adjust=FALSE)
  
  gcomp_contrasts <- tmle_contrast(Qstar=NULL, C, g_preds_bounded=NULL, obs.treatment, Y, outcome.type, alpha=0.05, gcomp = TRUE, iptw = FALSE, aiptw = FALSE, multi.adjust=FALSE, QAW=QAW)
  
  iptw_contrasts <- tmle_contrast(Qstar=IPTW, C, g_preds_bounded, obs.treatment, Y, outcome.type, alpha=0.05, gcomp = FALSE, iptw = TRUE, aiptw = FALSE, multi.adjust=FALSE)
  iptw_contrasts_bin <- tmle_contrast(Qstar=IPTW_bin, C, g_preds_bin_bounded, obs.treatment, Y, outcome.type, alpha=0.05, gcomp = FALSE, iptw = TRUE, aiptw = FALSE, multi.adjust=FALSE)
  
  aiptw_contrasts <- tmle_contrast(Qstar=AIPTW, C, g_preds_bounded, obs.treatment, Y, outcome.type, alpha=0.05, gcomp = FALSE, iptw = FALSE, aiptw = TRUE, multi.adjust=FALSE)
  aiptw_contrasts_bin <- tmle_contrast(Qstar=AIPTW_bin, C, g_preds_bin_bounded, obs.treatment, Y, outcome.type, alpha=0.05, gcomp = FALSE, iptw = FALSE, aiptw = TRUE, multi.adjust=FALSE)
  
  # store TMLE diagnostics (estimated initial Y, treatment probs.)
  
  yinitial_tmle  <- QAW[-1] 
  Ahat_tmle  <- g_preds
  Ahat_tmle_bin  <- g_preds_bin
  
  yhat_tmle <- colMeans(tmle_contrasts$Qstar)
  yhat_tmle_bin <- colMeans(tmle_contrasts_bin$Qstar)
  
  ess_ate_tmle <- (colSums(obs.treatment*(1/g_preds))**2)/(colSums(obs.treatment*(1/g_preds)**2)) 
  ess_ate_tmle_bin <- (colSums(obs.treatment*(1/g_preds_bin))**2)/(colSums(obs.treatment*(1/g_preds_bin)**2)) 
  
  # point and variance estimates for ATE
  
  bias_ate_tmle  <- tmle_contrasts$taus$ATE- true.ates # tmle
  var_ate_tmle <- sapply(tmle_contrasts$var,"[[",1)
  CP_ate_tmle <- as.numeric((sapply(tmle_contrasts$CI,"[[",1)[1,] < true.ates) & (sapply(tmle_contrasts$CI,"[[",1)[2,] > true.ates))
  CIW_ate_tmle  <- sapply(tmle_contrasts$CI,"[[",1)[2,] -sapply(tmle_contrasts$CI,"[[",1)[1,] 
  
  bias_ate_tmle_bin  <- tmle_contrasts_bin$taus$ATE- true.ates
  var_ate_tmle_bin <- sapply(tmle_contrasts_bin$var,"[[",1)
  CP_ate_tmle_bin <- as.numeric((sapply(tmle_contrasts_bin$CI,"[[",1)[1,] < true.ates) & (sapply(tmle_contrasts_bin$CI,"[[",1)[2,] > true.ates))
  CIW_ate_tmle_bin  <- sapply(tmle_contrasts_bin$CI,"[[",1)[2,] -sapply(tmle_contrasts_bin$CI,"[[",1)[1,] 
  
  bias_ate_gcomp  <- gcomp_contrasts$taus$ATE- true.ates # gcomp
  var_ate_gcomp <- sapply(gcomp_contrasts$var,"[[",1)
  CP_ate_gcomp <- as.numeric((sapply(gcomp_contrasts$CI,"[[",1)[1,] < true.ates) & (sapply(gcomp_contrasts$CI,"[[",1)[2,] > true.ates))
  CIW_ate_gcomp  <- sapply(gcomp_contrasts$CI,"[[",1)[2,] -sapply(gcomp_contrasts$CI,"[[",1)[1,] 
  
  bias_ate_iptw  <- iptw_contrasts$taus$ATE- true.ates # iptw
  var_ate_iptw <- sapply(iptw_contrasts$var,"[[",1)
  CP_ate_iptw <- as.numeric((sapply(iptw_contrasts$CI,"[[",1)[1,] < true.ates) & (sapply(iptw_contrasts$CI,"[[",1)[2,] > true.ates))
  CIW_ate_iptw  <- sapply(iptw_contrasts$CI,"[[",1)[2,] -sapply(iptw_contrasts$CI,"[[",1)[1,] 
  
  bias_ate_iptw_bin  <- iptw_contrasts_bin$taus$ATE- true.ates
  var_ate_iptw_bin <- sapply(iptw_contrasts_bin$var,"[[",1)
  CP_ate_iptw_bin <- as.numeric((sapply(iptw_contrasts_bin$CI,"[[",1)[1,] < true.ates) & (sapply(iptw_contrasts_bin$CI,"[[",1)[2,] > true.ates))
  CIW_ate_iptw_bin  <- sapply(iptw_contrasts_bin$CI,"[[",1)[2,] -sapply(iptw_contrasts_bin$CI,"[[",1)[1,]
  
  bias_ate_aiptw  <- aiptw_contrasts$taus$ATE- true.ates # aiptw
  var_ate_aiptw <- sapply(aiptw_contrasts$var,"[[",1)
  CP_ate_aiptw <- as.numeric((sapply(aiptw_contrasts$CI,"[[",1)[1,] < true.ates) & (sapply(aiptw_contrasts$CI,"[[",1)[2,] > true.ates))
  CIW_ate_aiptw  <- sapply(aiptw_contrasts$CI,"[[",1)[2,] -sapply(aiptw_contrasts$CI,"[[",1)[1,] 
  
  bias_ate_aiptw_bin  <- aiptw_contrasts_bin$taus$ATE- true.ates
  var_ate_aiptw_bin <- sapply(aiptw_contrasts_bin$var,"[[",1)
  CP_ate_aiptw_bin <- as.numeric((sapply(aiptw_contrasts_bin$CI,"[[",1)[1,] < true.ates) & (sapply(aiptw_contrasts_bin$CI,"[[",1)[2,] > true.ates))
  CIW_ate_aiptw_bin  <- sapply(aiptw_contrasts_bin$CI,"[[",1)[2,] -sapply(aiptw_contrasts_bin$CI,"[[",1)[1,] 
  
  # return results
  return(list("obs_outcome"=obs_outcome,"obs_treatment"=obs_treatment,"obs_covariates"=obs_covariates,"trueATE"=true.ates,
              "yinitial_tmle"=yinitial_tmle,"Ahat_tmle"=Ahat_tmle,"yhat_tmle"= yhat_tmle, "ess_ate_tmle"=ess_ate_tmle,
              "Ahat_tmle_bin"=Ahat_tmle_bin,"yhat_tmle_bin"= yhat_tmle_bin, "ess_ate_tmle_bin"=ess_ate_tmle_bin,
              "bias_ate_tmle"= unlist(bias_ate_tmle),"var_ate_tmle"= unlist(var_ate_tmle),"CP_ate_tmle"=unlist(CP_ate_tmle),"CIW_ate_tmle"=unlist(CIW_ate_tmle),
              "bias_ate_tmle_bin"= unlist(bias_ate_tmle_bin),"var_ate_tmle_bin"= unlist(var_ate_tmle_bin),"CP_ate_tmle_bin"=unlist(CP_ate_tmle_bin),"CIW_ate_tmle_bin"=unlist(CIW_ate_tmle_bin),
              "bias_ate_gcomp"= unlist(bias_ate_gcomp),"var_ate_gcomp"= unlist(var_ate_gcomp),"CP_ate_gcomp"=unlist(CP_ate_gcomp),"CIW_ate_gcomp"=unlist(CIW_ate_gcomp),
              "bias_ate_iptw"= unlist(bias_ate_iptw),"var_ate_iptw"= unlist(var_ate_iptw),"CP_ate_iptw"=unlist(CP_ate_iptw),"CIW_ate_iptw"=unlist(CIW_ate_iptw),
              "bias_ate_iptw_bin"= unlist(bias_ate_iptw_bin),"var_ate_iptw_bin"= unlist(var_ate_iptw_bin),"CP_ate_iptw_bin"=unlist(CP_ate_iptw_bin),"CIW_ate_iptw_bin"=unlist(CIW_ate_iptw_bin),
              "bias_ate_aiptw"= unlist(bias_ate_aiptw),"var_ate_aiptw"= unlist(var_ate_aiptw),"CP_ate_aiptw"=unlist(CP_ate_aiptw),"CIW_ate_aiptw"=unlist(CIW_ate_aiptw),
              "bias_ate_aiptw_bin"= unlist(bias_ate_aiptw_bin),"var_ate_aiptw_bin"= unlist(var_ate_aiptw_bin),"CP_ate_aiptw_bin"=unlist(CP_ate_aiptw_bin),"CIW_ate_aiptw_bin"=unlist(CIW_ate_aiptw_bin)))
}

#####################
# Set parameters    #
#####################

# define settings for simulation
settings <- expand.grid("J"=c(3,6),
                        "n"=c(10000),
                        "overlap.setting"=c("adequate","inadequate","rct"),
                        "gamma.setting"=c("zero","yang","low"))
settings$n <- ifelse(settings$J==6, 10000, 5000)

options(echo=TRUE)
args <- commandArgs(trailingOnly = TRUE) # command line arguments # args <- c("8",'binomial', 'TRUE', 'FALSE', 'FALSE', 'TRUE', 'FALSE', 'FALSE', 'FALSE')
thisrun <- settings[as.numeric(args[1]),]
outcome.type <-  as.character(args[2]) # "continuous" or "binomial" 
use.SL <- as.logical(args[3])  # When TRUE, use Super Learner for initial Y model and treatment model estimation; if FALSE, use GLM
doMPI <- as.logical(args[4]) # When TRUE, use MPI parallel processing
covars40 <- as.logical(args[5]) # When TRUE, generate 40 covariates instead of 6
covars100 <- as.logical(args[6]) # When TRUE, generate 100 covariates instead of 6
misTreatment <- as.logical(args[7]) # When TRUE, misspecify treatment model
misOutcome <- as.logical(args[8]) # When TRUE, misspecify outcome model
misBoth <- as.logical(args[9]) # When TRUE, misspecify both treatment and outcome models

# define parameters

J <- as.numeric(thisrun[1]) # number of treatments
n <- as.numeric(thisrun[2]) # total sample size
overlap.setting <- thisrun[[3]] # "adequate"= adequate overlap scenario of Li and Li (2019); "inadequate" overlap scenario of Yang et al. (2016); "rct": kappa values all zero
gamma.setting <- thisrun[[4]] # zero"= gamma values are all zero, so that there is no treatment effect; "yang"= use values from Yang et al. (2016); "li" use values from Li and Li (2019)

if(covars40 | covars100){
  R <- 100 # number of simulation runs
}else{
  R <- 1000 # number of simulation runs
}

gbound <- c(0.001,0.999) # define bounds to be used for the propensity score

if(outcome.type=="continuous"){
  ybound <- c(0,1000)
}else{
  ybound <- c(0.0001, 0.9999) # define bounds to be used for the Y predictions
}

n.folds <- 5

target.gwt <- TRUE # When TRUE, moves propensity weights from denominator of clever covariate to regression weight when fitting updated model for Y

scale.continuous <- TRUE # standardize continuous covariates

output_dir <- './outputs/'
simulation_version <- paste0(format(Sys.time(), "%Y%m%d"),"/")
if(!dir.exists(output_dir)){
  print(paste0('create folder for outputs at: ', output_dir))
  dir.create(output_dir)
}
output_dir = paste0(output_dir, simulation_version)
if(!dir.exists(output_dir)){
  print(paste0('create folder for outputs at: ', output_dir))
  dir.create(output_dir)
}

filename <- paste0(output_dir, 
                   "_R_", R,
                   "_n_", n,
                   "_J_", J,
                   "_n_folds_",n.folds,
                   "_overlap_setting_", overlap.setting,
                   "_gamma_setting_", gamma.setting,
                   "_outcome_type_", outcome.type,
                   "_use_SL_", use.SL,
                   "_covars40_", covars40,
                   "_covars100_", covars100,
                   "_misTreat_", misTreatment,
                   "_misOut_", misOutcome,
                   "_misBoth_", misBoth,
                   "_scale_cont_", scale.continuous,
                   "_target_gwt_", target.gwt,".rds")

# Setup parallel processing
if(doMPI){
  library(doMPI)
  
  # Start cluster
  cl <- startMPIcluster(verbose=TRUE)
  
  # Register cluster
  registerDoMPI(cl)
  
  # Check cluster size
  print(paste0("cluster size: ", clusterSize(cl)))
  
} else{
  library(parallel)
  library(doParallel)
  library(foreach)
  
  cores <- parallel::detectCores()
  print(paste0("number of cores used: ", cores))
  
  cl <- parallel::makeCluster(cores, outfile="")
  
  doParallel::registerDoParallel(cl) # register cluster
}

#####################
# Run simulation #
#####################

sim.results <- foreach(r = 1:R, .combine='cbind', .verbose = FALSE) %dopar% {
  staticSim(r=r, J, n, gbound, ybound, n.folds, overlap.setting, gamma.setting, outcome.type, target.gwt, use.SL, scale.continuous, covars40, covars100, misTreatment,misOutcome,misBoth)
}
sim.results

saveRDS(sim.results, filename)

if(doMPI){
  closeCluster(cl) # close down MPIcluster
  mpi.finalize()
}else{
  stopCluster(cl)
}