diff --git a/.gitignore b/.gitignore
index 8e6d851a2..f8bf65c30 100644
--- a/.gitignore
+++ b/.gitignore
@@ -37,7 +37,6 @@ src/*.so
 src/*.dll
 testthat-problems.rds
 README.md
-**/sascfg_personal.py
 .Rprofile
 /doc/
 /Meta/
diff --git a/simulations/missing-data-benchmarks/.gitignore b/simulations/missing-data-benchmarks/.gitignore
index e69a4c2c1..914c4293f 100644
--- a/simulations/missing-data-benchmarks/.gitignore
+++ b/simulations/missing-data-benchmarks/.gitignore
@@ -1,2 +1,3 @@
 **/fit_results.[Rr]ds
 slurm.out
+.future
diff --git a/simulations/missing-data-benchmarks/R/customizations/experiment.R b/simulations/missing-data-benchmarks/R/customizations/experiment.R
new file mode 100644
index 000000000..2ac4eb5e6
--- /dev/null
+++ b/simulations/missing-data-benchmarks/R/customizations/experiment.R
@@ -0,0 +1,108 @@
+stopifnot(any(grepl("simChef", search())))
+
+# Custom new method which only kicks off the future jobs.
+# This is modified from the fit() method, see
+# https://github.com/Yu-Group/simChef/blob/main/R/experiment.R#L976C3-L980C62
+Experiment$set(
+  "public",
+  "onlystart",
+  function(n_reps = 1,
+           future.globals = NULL,
+           future.packages = NULL,
+           future.seed = TRUE,
+           verbose = 1, ...) {
+    parallel_strategy <- "reps"
+
+    dgp_list <- private$.get_obj_list("dgp")
+    method_list <- private$.get_obj_list("method")
+
+    if (length(dgp_list) == 0) {
+      private$.throw_empty_list_error("dgp", "generate data from")
+    }
+
+    if (length(method_list) == 0) {
+      private$.throw_empty_list_error("method", "fit methods in")
+    }
+
+    private$.update_fit_params()
+
+    checkpoint <- FALSE
+    n_reps_cached <- 0
+    n_reps_total <- n_reps
+    fit_results <- data.frame()
+
+    if (verbose >= 1) {
+      inform(sprintf("Fitting %s...", self$name))
+      start_time <- Sys.time()
+    }
+
+    if (is.null(future.packages)) {
+      future.packages <- private$.future.packages
+    }
+
+    if (is.null(future.globals)) {
+      future.globals <- private$.future.globals
+    }
+
+    dgp_params_list <- private$.combine_vary_params("dgp")
+    method_params_list <- private$.combine_vary_params("method")
+
+    # if new_fit_params is not NULL after the if statement below, then not all
+    # combos of (dgp_params_list, method_params_list) need to be rerun so need
+    # to check cache ids when fitting
+    new_fit_params <- NULL
+
+    duplicate_param_names <- private$.get_duplicate_param_names()
+
+    # simulations
+    n_reps <- min(n_reps, n_reps_total - n_reps_cached)
+
+    new_fit_results <- local({
+
+      # create an env with objs/funcs that the future workers need
+      workenv <- rlang::new_environment(
+        data = list(
+          verbose = verbose,
+          dgp_list = dgp_list,
+          method_list = method_list,
+          new_fit_params = new_fit_params,
+          dgp_params_list = dgp_params_list,
+          method_params_list = method_params_list,
+          duplicate_param_names = duplicate_param_names,
+          do_call_wrapper = function(name,
+                                     fun,
+                                     params,
+                                     verbose,
+                                     call) {
+            tryCatch(
+              do_call_handler(
+                name, fun, params, verbose, call
+              ),
+              error = identity
+            )
+          }
+        ),
+        parent = rlang::ns_env()
+      )
+
+      # get the experiment compute fun
+      compute_fun <- compute_rep
+
+      environment(compute_fun) <- workenv
+
+      # compute the experiment
+      compute_fun(n_reps,
+                  future.globals,
+                  future.packages,
+                  future.seed)
+    })
+
+    gc()
+
+    new_fit_results %>%
+      dplyr::mutate(
+        .rep = as.character(as.numeric(.rep) + n_reps_cached)
+      ) %>%
+      simplify_tibble()
+  }
+)
diff --git a/simulations/missing-data-benchmarks/R/customizations/format-fit-results.R b/simulations/missing-data-benchmarks/R/customizations/format-fit-results.R
new file mode 100644
index 000000000..5907da73b
--- /dev/null
+++ b/simulations/missing-data-benchmarks/R/customizations/format-fit-results.R
@@ -0,0 +1,44 @@
+# format the fit results
+format_fit_results <- function(fit_results) {
+  fit_results %>%
+    transmute(
+      effect_size = str_extract(.dgp_name, "^([^_])+"),
+      rep = .rep,
+      dgp_name = .dgp_name,
+      method_name = .method_name,
+      converged = pmap(
+        .l = list(fit, .method_name, converged),
+        .f = function(f, method_name, conv_status) {
+          get_convergence(method_name, f, conv_status)
+        }
+      ),
+      converged = unlist(converged),
+      fit_time = fit_time,
+      emmeans_output = pmap(
+        .l = list(fit, .dgp_name, .method_name, data, converged),
+        .f = function(f, dgp_name, method_name, dt, conv_status) {
+          get_emmeans_output(method_name, f, dt, conv_status)
+        }
+      ),
+      covmat_estimates = pmap(
+        .l = list(fit, .method_name, converged),
+        .f = function(f, method_name, conv_status) {
+          get_cov_mat_estimate(method_name, f, conv_status)
+        }
+      )
+    )
+}
+
+# The whole process
+format_fit_and_save <- function(experiment) {
+  gc()
+  fit_results <- experiment$get_cached_results("fit")
+  formatted_fit_results <- format_fit_results(fit_results)
+  formatted_fit_results_file <- file.path(experiment$get_save_dir(), "formatted_fit_results.rds")
+  saveRDS(formatted_fit_results, file = formatted_fit_results_file)
+  formatted_fit_results
+}
+
+
+
+
diff --git a/simulations/missing-data-benchmarks/R/eval/helpers.R b/simulations/missing-data-benchmarks/R/eval/helpers.R
index f3e874733..17311f2f4 100644
--- a/simulations/missing-data-benchmarks/R/eval/helpers.R
+++ b/simulations/missing-data-benchmarks/R/eval/helpers.R
@@ -242,19 +242,26 @@ format_emmeans_df <- function(emmeans_df) {
 
 # extract model fits output at each visit from mmrm fit
 get_mmrm_emmeans_output <- function(fit) {
-  # extract emmeans output
-  marginal_means <- emmeans(
-    fit,
-    spec = trt.vs.ctrl ~ trt | visit_num,
-    weights = "proportional"
-  )
-  emmeans_df <- as.data.frame(marginal_means$contrasts)
+  result <- try({
+    # extract emmeans output
+    marginal_means <- emmeans(
+      fit,
+      spec = trt.vs.ctrl ~ trt | visit_num,
+      weights = "proportional"
+    )
+    emmeans_df <- as.data.frame(marginal_means$contrasts)
 
-  # compute lower and upper 95% CI
-  emmeans_df <- get_95_ci(emmeans_df)
+    # compute lower and upper 95% CI
+    emmeans_df <- get_95_ci(emmeans_df)
 
-  # format to resemble SAS output
-  format_emmeans_df(emmeans_df)
+    # format to resemble SAS output
+    format_emmeans_df(emmeans_df)
+  })
+  if (!inherits(result, "try-error")) {
+    result
+  } else {
+    NA
+  }
 }
 
 # extract model fit outputs at each visit from glmmTMB fit
@@ -323,3 +330,49 @@ get_emmeans_output <- function(method, fit, dt, converged) {
     NA
   }
 }
+
+# get_cov_mat_estimate ----
+
+# extract covariance matrix estimate from mmrm fit
+get_mmrm_cov_mat_estimate <- function(fit) {
+  mmrm::VarCorr(fit)
+}
+
+# extract covariance matrix estimate from glmmTMB fit
+get_glmm_cov_mat_estimate <- function(fit) {
+  mat_with_attrs <- glmmTMB::VarCorr(fit)[[c("cond", "participant")]]
+  dim <- nrow(mat_with_attrs)
+  ind <- seq_len(dim)
+  mat_with_attrs[ind, ind]
+}
+
+# extract covariance matrix estimate from nlme fit
+get_nlme_cov_mat_estimate <- function(fit) {
+  mat_with_attrs <- nlme::getVarCov(fit)
+  dim <- nrow(mat_with_attrs)
+  ind <- seq_len(dim)
+  mat_with_attrs[ind, ind]
+}
+
+# extract covariance matrix estimate from proc mixed fit
+get_mixed_cov_mat_estimate <- function(fit) {
+  attr(fit, "covmat")
+}
+
+# general function for covariance matrix estimate
+get_cov_mat_estimate <- function(method, fit, converged) {
+  if (get_convergence(method, fit, converged)) {
+    if (str_detect(method, "mmrm")) {
+      get_mmrm_cov_mat_estimate(fit)
+    } else if (str_detect(method, "glmmtmb")) {
+      get_glmm_cov_mat_estimate(fit)
+    } else if (str_detect(method, "nlme")) {
+      get_nlme_cov_mat_estimate(fit)
+    } else if (str_detect(method, "proc_mixed")) {
+      get_mixed_cov_mat_estimate(fit)
+    }
+  } else {
+    NA
+  }
+}
+
diff --git a/simulations/missing-data-benchmarks/R/format-replicate-results/format-results.R b/simulations/missing-data-benchmarks/R/format-replicate-results/format-results.R
index 759352bb9..07b62248e 100644
--- a/simulations/missing-data-benchmarks/R/format-replicate-results/format-results.R
+++ b/simulations/missing-data-benchmarks/R/format-replicate-results/format-results.R
@@ -9,12 +9,26 @@ library(tidyr)
 library(emmeans)
 
 # load the required helper functions
-source("R/format-replicate-results/helpers.R")
+# This is just an example how to use these functions.
+
+source("R/customizations/format-fit-results.R")
 source("R/eval/helpers.R")
 
 # load a fit_results.rds file
-fit_results <- readRDS("R/format-replicate-results/fit_results_100.rds")
+fit_results <- readRDS("R/format-replicate-results/fit_results_test.rds")
 
 # format the results
-formatted_fit_results_df <- fit_results %>%
-  format_fit_results(missingness = "none", sample_size = 200)
+formatted_fit_results_df <- format_fit_results(fit_results)
+saveRDS(formatted_fit_results_df, file = "R/format-replicate-results/df_test.rds")
+
+# formatted_fit_results_df <- readRDS("df_test.rds")
+head(formatted_fit_results_df)
+
+# to obtain the estimated covariance matrices:
+formatted_fit_results_df$covmat_estimates[1:2]
+
+# to obtain the emmeans by visit:
+formatted_fit_results_df |>
+  select(- covmat_estimates, - fit_time, - converged) |>
+  unnest(cols = emmeans_output)
+
diff --git a/simulations/missing-data-benchmarks/R/format-replicate-results/helpers.R b/simulations/missing-data-benchmarks/R/format-replicate-results/helpers.R
deleted file mode 100644
index 1db13109d..000000000
--- a/simulations/missing-data-benchmarks/R/format-replicate-results/helpers.R
+++ /dev/null
@@ -1,134 +0,0 @@
-# compute the 95% CI from emmeans output
-get_95_ci <- function(emmeans_df) {
-  emmeans_df %>%
-    mutate(
-      lower = estimate - 1.96 * SE,
-      upper = estimate + 1.96 * SE,
-    )
-}
-
-format_emmeans_df <- function(emmeans_df) {
-  emmeans_df %>%
-    transmute(
-      visit_num = visit_num,
-      estimate = estimate,
-      stderr = SE,
-      df = df,
-      tvalue = t.ratio,
-      pvalue = p.value,
-      lower = lower,
-      upper = upper
-    )
-}
-
-# extract model fits output at each visit from mmrm fit
-get_mmrm_emmeans_output <- function(fit) {
-  # extract emmeans output
-  marginal_means <- emmeans(
-    fit,
-    spec = trt.vs.ctrl ~ trt | visit_num,
-    weights = "proportional"
-  )
-  emmeans_df <- as.data.frame(marginal_means$contrasts)
-
-  # compute lower and upper 95% CI
-  emmeans_df <- get_95_ci(emmeans_df)
-
-  # format to resemble SAS output
-  format_emmeans_df(emmeans_df)
-}
-
-# extract model fit outputs at each visit from glmmTMB fit
-get_glmm_emmeans_output <- function(fit) {
-  marginal_means <- emmeans(
-    fit,
-    spec = trt.vs.ctrl ~ trt | visit_num,
-    weights = "proportional"
-  )
-  emmeans_df <- as.data.frame(marginal_means$contrasts)
-
-  # compute lower and upper 95% CI
-  emmeans_df <- get_95_ci(emmeans_df)
-
-  # format to resemble SAS output
-  format_emmeans_df(emmeans_df)
-}
-
-# extract model fit outputs estimates at each visit from gls fit
-get_nlme_emmeans_output <- function(fit, data) {
-  marginal_means <- emmeans(
-    fit,
-    spec = trt.vs.ctrl ~ trt | visit_num,
-    weights = "proportional",
-    data = data,
-    mode = "df.error"
-  )
-  emmeans_df <- as.data.frame(marginal_means$contrasts)
-
-  # compute lower and upper 95% CI
-  emmeans_df <- get_95_ci(emmeans_df)
-
-  # format to resemble SAS output
-  format_emmeans_df(emmeans_df)
-}
-
-# extract model fit outputs estimates at each visit from PROC MIXED fit
-get_mixed_emmeans_like_output <- function(fit) {
-  fit %>% transmute(
-    visit_num = str_extract(contrast, "^([^:])+"),
-    estimate = Estimate,
-    stderr = StdErr,
-    df = DF,
-    tvalue = tValue,
-    pvalue = map2_dbl(tValue, DF, function(tvalue, df) {
-      2 * min(c(pt(tvalue, df), pt(tvalue, df, lower.tail = FALSE)))
-    }),
-    lower = Lower,
-    upper = Upper
-  )
-}
-
-# general function for emmeans like information
-get_emmeans_output <- function(method, fit, dt, converged) {
-  if (get_convergence(method, fit, converged)) {
-    if (str_detect(method, "mmrm")) {
-      get_mmrm_emmeans_output(fit)
-    } else if (str_detect(method, "glmmtmb")) {
-      get_glmm_emmeans_output(fit)
-    } else if (str_detect(method, "nlme")) {
-      get_nlme_emmeans_output(fit, dt)
-    } else if (str_detect(method, "proc_mixed")) {
-      get_mixed_emmeans_like_output(fit)
-    }
-  } else {
-    NA
-  }
-}
-
-# format the fit results
-format_fit_results <- function(fit_results, missingness, sample_size) {
-  fit_results %>%
-    transmute(
-      missingness = missingness,
-      sample_size = sample_size,
-      effect_size = str_extract(.dgp_name, "^([^_])+"),
-      rep = .rep,
-      dgp_name = .dgp_name,
-      method_name = .method_name,
-      converged = pmap(
-        .l = list(fit, .method_name, converged),
-        .f = function(f, method_name, conv_status) {
-          get_convergence(method_name, f, conv_status)
-        }
-      ),
-      converged = unlist(converged),
-      fit_time = fit_time,
-      emmeans_output = pmap(
-        .l = list(fit, .dgp_name, .method_name, data, converged),
-        .f = function(f, dgp_name, method_name, dt, conv_status) {
-          get_emmeans_output(method_name, f, dt, conv_status)
-        }
-      )
-    ) %>%
-    unnest(cols = emmeans_output)
-}
diff --git a/simulations/missing-data-benchmarks/R/meal.R b/simulations/missing-data-benchmarks/R/meal.R
index 7a0420dfb..1b80516ab 100644
--- a/simulations/missing-data-benchmarks/R/meal.R
+++ b/simulations/missing-data-benchmarks/R/meal.R
@@ -131,7 +131,7 @@ type_2_error_rate_eval <- create_evaluator(
   .eval_fun = type_2_error_rate_fun, true_params = true_params
 )
 
-# specify the resul summarizers
+# specify the result "summarizers"
 
 # create the experiment
 experiment <- create_experiment(
diff --git a/simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness-R.R b/simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness-R.R
new file mode 100644
index 000000000..892e95f36
--- /dev/null
+++ b/simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness-R.R
@@ -0,0 +1,268 @@
+################################################################################
+# Custom modification of simChef for
+# simulation with Extreme Levels of Missingness
+################################################################################
+
+# load required libraries
+library(simChef)
+library(mmrm)
+library(glmmTMB)
+library(nlme)
+library(sasr)
+library(stringr)
+library(dplyr)
+library(purrr)
+library(tidyr)
+library(emmeans)
+library(clusterGeneration)
+library(ssh)
+library(microbenchmark)
+library(future)
+library(future.batchtools)
+
+# source the R scripts
+sim_functions_files <- list.files(
+  c("R/dgp", "R/method", "R/eval", "R/viz", "R/customizations"),
+  pattern = "*.R$", full.names = TRUE, ignore.case = TRUE
+)
+sapply(sim_functions_files, source)
+
+# generate the possible covariance matrices
+us_cov_mat <- compute_unstructured_matrix()
+csh_cov_mat <- compute_csh_matrix()
+toep_cov_mat <- toeplitz(c(1, 0.5, 0.25, 0.125, rep(0, 6)))
+
+# set the sample size
+n_obs <- 600
+
+# dgps with no treatment effect
+no_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  missingness = "extreme"
+)
+no_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  missingness = "extreme"
+)
+no_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  missingness = "extreme"
+)
+
+# dgps with small treatment effect
+small_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "extreme"
+)
+small_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "extreme"
+)
+small_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "extreme"
+)
+
+# dgps with moderate treatment effect
+mod_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "extreme"
+)
+mod_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "extreme"
+)
+mod_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "extreme"
+)
+
+# specify the methods
+mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
+mmrm_csh_meth <- create_method(
+  .method_fun = mmrm_wrapper_fun, covar_type = "csh"
+)
+mmrm_toeph_meth <- create_method(
+  .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
+)
+glmmtmb_us_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
+)
+glmmtmb_csh_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
+)
+glmmtmb_toeph_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
+)
+nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
+nlme_csh_meth <- create_method(
+  .method_fun = nlme_wrapper_fun, covar_type = "csh"
+)
+# proc_mixed_us_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "us"
+# )
+# proc_mixed_csh_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "csh"
+# )
+# proc_mixed_toeph_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "toeph"
+# )
+
+# specify the evaluation metrics
+mean_time_eval <- create_evaluator(.eval_fun = mean_time_fun)
+true_params <- list(
+  "no_effect_us" = rep(0, 10),
+  "no_effect_csh" = rep(0, 10),
+  "no_effect_toeph" = rep(0, 10),
+  "small_effect_us" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "small_effect_csh" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "small_effect_toeph" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "mod_effect_us" = seq(from = 0.5, by = 0.5, length.out = 10),
+  "mod_effect_csh" = seq(from = 0.5, by = 0.5, length.out = 10),
+  "mod_effect_toeph" = seq(from = 0.5, by = 0.5, length.out = 10)
+)
+bias_eval <- create_evaluator(.eval_fun = bias_fun, true_params = true_params)
+variance_eval <- create_evaluator(.eval_fun = variance_fun)
+convergence_rate_eval <- create_evaluator(.eval_fun = convergence_rate_fun)
+coverage_eval <- create_evaluator(
+  .eval_fun = coverage_fun, true_params = true_params
+)
+type_1_error_rate_eval <- create_evaluator(
+  .eval_fun = type_1_error_rate_fun, true_params = true_params
+)
+type_2_error_rate_eval <- create_evaluator(
+  .eval_fun = type_2_error_rate_fun, true_params = true_params
+)
+
+# fixInNamespace
+fixInNamespace("readLog", "batchtools")
+
+# Define parallelization
+future_globals <- ls()
+
+options(future.globals.onReference = NULL)
+# options(future.globals.onReference = "error")
+
+plan(
+  batchtools_lsf,
+  workers = 1000,
+  resources = list(memory = 8000, walltime = 3 * 60 * 60)
+)
+
+# n_workers <- parallelly::availableCores() - 1L
+# plan(future::multicore, workers = n_workers)
+
+# create the experiment
+experiment <- create_experiment(
+  name = "mmrm-benchmark-extreme-missingness-n-600-R",
+  save_dir = "results/extreme-miss/n-600-R",
+  future.globals = future_globals
+) %>%
+  add_dgp(no_effect_us_dgp, name = "no_effect_us") %>%
+  add_dgp(no_effect_csh_dgp, name = "no_effect_csh") %>%
+  add_dgp(no_effect_toeph_dgp, name = "no_effect_toeph") %>%
+  add_dgp(small_effect_us_dgp, name = "small_effect_us") %>%
+  add_dgp(small_effect_csh_dgp, name = "small_effect_csh") %>%
+  add_dgp(small_effect_toeph_dgp, name = "small_effect_toeph") %>%
+  add_dgp(mod_effect_us_dgp, name = "mod_effect_us") %>%
+  add_dgp(mod_effect_csh_dgp, name = "mod_effect_csh") %>%
+  add_dgp(mod_effect_toeph_dgp, name = "mod_effect_toeph") %>%
+  add_method(mmrm_us_meth, name = "mmrm_us") %>%
+  add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
+  add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
+  add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
+  add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
+  add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
+  add_method(nlme_us_meth, name = "nlme_us") %>%
+  add_method(nlme_csh_meth, name = "nlme_csh") %>%
+  # We don't run the SAS routines here.
+  # add_method(proc_mixed_us_meth, name = "proc_mixed_us") %>%
+  # add_method(proc_mixed_csh_meth, name = "proc_mixed_csh") %>%
+  # add_method(proc_mixed_toeph_meth, name = "proc_mixed_toeph") %>%
+  add_evaluator(mean_time_eval, name = "mean_fit_time") %>%
+  add_evaluator(convergence_rate_eval, name = "convergence_rate") %>%
+  add_evaluator(bias_eval, name = "bias") %>%
+  add_evaluator(variance_eval, name = "variance") %>%
+  add_evaluator(coverage_eval, name = "coverage_rate") %>%
+  add_evaluator(type_1_error_rate_eval, name = "type_1_error_rate") %>%
+  add_evaluator(type_2_error_rate_eval, name = "type_2_error_rate")
+
+# No need for SAS here
+
+# # Determine host name
+# hostname <- "sabanesd-nwpfwh-eu.ocean"
+#
+# # Do this to register host
+# session <- ssh::ssh_connect(hostname)
+# ssh_disconnect(session)
+#
+# # We can also be specific if there are multiple SAS containers we want to connect to.
+# cfg_name <- "~/sascfg_personal.py"
+# sasr.roche::setup_sasr_ocean(hostname, sascfg = cfg_name)
+# options(sascfg = cfg_name)
+#
+# # Test if it works.
+# df <- data.frame(a = 1, b = 2)
+# sasr::df2sd(df)
+
+
+# run the experiment
+# we break it down into 10 batches so that we stay within the memory limits
+
+exp_dir <- experiment$get_save_dir()
+
+for (i in 7:10) {
+  # run this batch
+  set.seed(23412 + i * 10)
+  results <- experiment$onlystart(
+    n_reps = 100,
+    verbose = 2
+  )
+
+  # format results
+  formatted_fit_results <- format_fit_results(results)
+
+  # correct repetition information
+  formatted_fit_results$rep <- as.integer(formatted_fit_results$rep) + 100L * (i - 1L)
+
+  # save this batch
+  file_name <- paste0("formatted_fit_results_", i, ".rds")
+  file_name <- file.path(exp_dir, file_name)
+  saveRDS(formatted_fit_results, file = file_name)
+
+  # clean memory
+  rm(formatted_fit_results)
+  rm(results)
+  gc()
+}
+
+# Finally aggregate all tibbles together:
+all_files <- grep(pattern = "^formatted_fit_results_", x = dir(exp_dir), value = TRUE)
+all_files <- file.path(exp_dir, all_files)
+all_results <- lapply(all_files, readRDS)
+result <- do.call(rbind, all_results)
+saveRDS(result, file = file.path(exp_dir, "result.rds"))
diff --git a/simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness.R b/simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness-SAS.R
similarity index 71%
rename from simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness.R
rename to simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness-SAS.R
index 9f96d044f..9f8aa87a0 100644
--- a/simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness.R
+++ b/simulations/missing-data-benchmarks/R/meals/meal-extreme-missingness-SAS.R
@@ -13,6 +13,9 @@ library(dplyr)
 library(purrr)
 library(tidyr)
 library(emmeans)
+library(clusterGeneration)
+library(ssh)
+library(microbenchmark)
 
 # source the R scripts
 sim_functions_files <- list.files(
@@ -96,26 +99,26 @@ mod_effect_toeph_dgp <- create_dgp(
 )
 
 # specify the methods
-mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
-mmrm_csh_meth <- create_method(
-  .method_fun = mmrm_wrapper_fun, covar_type = "csh"
-)
-mmrm_toeph_meth <- create_method(
-  .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
-)
-glmmtmb_us_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
-)
-glmmtmb_csh_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
-)
-glmmtmb_toeph_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
-)
-nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
-nlme_csh_meth <- create_method(
-  .method_fun = nlme_wrapper_fun, covar_type = "csh"
-)
+# mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
+# mmrm_csh_meth <- create_method(
+#   .method_fun = mmrm_wrapper_fun, covar_type = "csh"
+# )
+# mmrm_toeph_meth <- create_method(
+#   .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
+# )
+# glmmtmb_us_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
+# )
+# glmmtmb_csh_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
+# )
+# glmmtmb_toeph_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
+# )
+# nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
+# nlme_csh_meth <- create_method(
+#   .method_fun = nlme_wrapper_fun, covar_type = "csh"
+# )
 proc_mixed_us_meth <- create_method(
   .method_fun = proc_mixed_wrapper_fun, covar_type = "us"
 )
@@ -154,8 +157,8 @@ type_2_error_rate_eval <- create_evaluator(
 
 # create the experiment
 experiment <- create_experiment(
-  name = "mmrm-benchmark-extreme-missingness-n-600",
-  save_dir = "results/extreme-miss/n-600"
+  name = "mmrm-benchmark-extreme-missingness-n-600-SAS",
+  save_dir = "results/extreme-miss/n-600-SAS"
 ) %>%
   add_dgp(no_effect_us_dgp, name = "no_effect_us") %>%
   add_dgp(no_effect_csh_dgp, name = "no_effect_csh") %>%
@@ -166,14 +169,14 @@ experiment <- create_experiment(
   add_dgp(mod_effect_us_dgp, name = "mod_effect_us") %>%
   add_dgp(mod_effect_csh_dgp, name = "mod_effect_csh") %>%
   add_dgp(mod_effect_toeph_dgp, name = "mod_effect_toeph") %>%
-  add_method(mmrm_us_meth, name = "mmrm_us") %>%
-  add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
-  add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
-  add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
-  add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
-  add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
-  add_method(nlme_us_meth, name = "nlme_us") %>%
-  add_method(nlme_csh_meth, name = "nlme_csh") %>%
+  # add_method(mmrm_us_meth, name = "mmrm_us") %>%
+  # add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
+  # add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
+  # add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
+  # add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
+  # add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
+  # add_method(nlme_us_meth, name = "nlme_us") %>%
+  # add_method(nlme_csh_meth, name = "nlme_csh") %>%
   add_method(proc_mixed_us_meth, name = "proc_mixed_us") %>%
   add_method(proc_mixed_csh_meth, name = "proc_mixed_csh") %>%
   add_method(proc_mixed_toeph_meth, name = "proc_mixed_toeph") %>%
@@ -185,11 +188,32 @@ experiment <- create_experiment(
   add_evaluator(type_1_error_rate_eval, name = "type_1_error_rate") %>%
   add_evaluator(type_2_error_rate_eval, name = "type_2_error_rate")
 
+# # Determine host name
+# hostname <- "sabanesd-nwpfwh-eu.ocean"
+#
+# # Do this to register host
+# session <- ssh::ssh_connect(hostname)
+# ssh_disconnect(session)
+#
+# # We can also be specific if there are multiple SAS containers we want to connect to.
+# cfg_name <- "~/sascfg_personal.py"
+# sasr.roche::setup_sasr_ocean(hostname, sascfg = cfg_name)
+# options(sascfg = cfg_name)
+#
+# # Test if it works.
+# df <- data.frame(a = 1, b = 2)
+# sasr::df2sd(df)
+
 # run the experiment
-set.seed(713452)
-results <- experiment$run(
-  n_reps = 100,
-  save = TRUE,
-  verbose = 2,
-  checkpoint_n_reps = 25
-)
+# set.seed(713452)
+# results <- experiment$run(
+#   n_reps = 1000,
+#   save = TRUE,
+#   verbose = 2,
+#   checkpoint_n_reps = 10
+# )
+
+source("R/format-replicate-results/helpers.R")
+
+format_fit_and_save(experiment)
+
diff --git a/simulations/missing-data-benchmarks/R/meals/meal-high-missingness-R.R b/simulations/missing-data-benchmarks/R/meals/meal-high-missingness-R.R
new file mode 100644
index 000000000..e3cd60483
--- /dev/null
+++ b/simulations/missing-data-benchmarks/R/meals/meal-high-missingness-R.R
@@ -0,0 +1,266 @@
+################################################################################
+# Simulation with High Levels of Missingness
+################################################################################
+
+# load required libraries
+library(simChef)
+library(mmrm)
+library(glmmTMB)
+library(nlme)
+library(sasr)
+library(stringr)
+library(dplyr)
+library(purrr)
+library(tidyr)
+library(emmeans)
+library(clusterGeneration)
+library(ssh)
+library(microbenchmark)
+library(future)
+library(future.batchtools)
+
+# source the R scripts
+sim_functions_files <- list.files(
+  c("R/dgp", "R/method", "R/eval", "R/viz", "R/customizations"),
+  pattern = "*.R$", full.names = TRUE, ignore.case = TRUE
+)
+sapply(sim_functions_files, source)
+
+# generate the possible covariance matrices
+us_cov_mat <- compute_unstructured_matrix()
+csh_cov_mat <- compute_csh_matrix()
+toep_cov_mat <- toeplitz(c(1, 0.5, 0.25, 0.125, rep(0, 6)))
+
+# set the sample size
+n_obs <- 600
+
+# dgps with no treatment effect
+no_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  missingness = "high"
+)
+no_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  missingness = "high"
+)
+no_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  missingness = "high"
+)
+
+# dgps with small treatment effect
+small_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "high"
+)
+small_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "high"
+)
+small_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "high"
+)
+
+# dgps with moderate treatment effect
+mod_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "high"
+)
+mod_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "high"
+)
+mod_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "high"
+)
+
+# specify the methods
+mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
+mmrm_csh_meth <- create_method(
+  .method_fun = mmrm_wrapper_fun, covar_type = "csh"
+)
+mmrm_toeph_meth <- create_method(
+  .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
+)
+glmmtmb_us_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
+)
+glmmtmb_csh_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
+)
+glmmtmb_toeph_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
+)
+nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
+nlme_csh_meth <- create_method(
+  .method_fun = nlme_wrapper_fun, covar_type = "csh"
+)
+# proc_mixed_us_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "us"
+# )
+# proc_mixed_csh_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "csh"
+# )
+# proc_mixed_toeph_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "toeph"
+# )
+
+# specify the evaluation metrics
+mean_time_eval <- create_evaluator(.eval_fun = mean_time_fun)
+true_params <- list(
+  "no_effect_us" = rep(0, 10),
+  "no_effect_csh" = rep(0, 10),
+  "no_effect_toeph" = rep(0, 10),
+  "small_effect_us" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "small_effect_csh" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "small_effect_toeph" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "mod_effect_us" = seq(from = 0.5, by = 0.5, length.out = 10),
+  "mod_effect_csh" = seq(from = 0.5, by = 0.5, length.out = 10),
+  "mod_effect_toeph" = seq(from = 0.5, by = 0.5, length.out = 10)
+)
+bias_eval <- create_evaluator(.eval_fun = bias_fun, true_params = true_params)
+variance_eval <- create_evaluator(.eval_fun = variance_fun)
+convergence_rate_eval <- create_evaluator(.eval_fun = convergence_rate_fun)
+coverage_eval <- create_evaluator(
+  .eval_fun = coverage_fun, true_params = true_params
+)
+type_1_error_rate_eval <- create_evaluator(
+  .eval_fun = type_1_error_rate_fun, true_params = true_params
+)
+type_2_error_rate_eval <- create_evaluator(
+  .eval_fun = type_2_error_rate_fun, true_params = true_params
+)
+
+# fixInNamespace
+fixInNamespace("readLog", "batchtools")
+
+# Define parallelization
+future_globals <- ls()
+
+options(future.globals.onReference = NULL)
+# options(future.globals.onReference = "error")
+
+plan(
+  batchtools_lsf,
+  workers = 100,
+  resources = list(memory = 10000, walltime = 3 * 60 * 60) # queue = "short",
+)
+
+# create the experiment
+experiment <- create_experiment(
+  name = "mmrm-benchmark-high-missingness-n-600-R",
+  save_dir = "results/high-miss/n-600-R",
+  future.globals = future_globals
+) %>%
+  add_dgp(no_effect_us_dgp, name = "no_effect_us") %>%
+  add_dgp(no_effect_csh_dgp, name = "no_effect_csh") %>%
+  add_dgp(no_effect_toeph_dgp, name = "no_effect_toeph") %>%
+  add_dgp(small_effect_us_dgp, name = "small_effect_us") %>%
+  add_dgp(small_effect_csh_dgp, name = "small_effect_csh") %>%
+  add_dgp(small_effect_toeph_dgp, name = "small_effect_toeph") %>%
+  add_dgp(mod_effect_us_dgp, name = "mod_effect_us") %>%
+  add_dgp(mod_effect_csh_dgp, name = "mod_effect_csh") %>%
+  add_dgp(mod_effect_toeph_dgp, name = "mod_effect_toeph") %>%
+  add_method(mmrm_us_meth, name = "mmrm_us") %>%
+  add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
+  add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
+  add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
+  add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
+  add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
+  add_method(nlme_us_meth, name = "nlme_us") %>%
+  add_method(nlme_csh_meth, name = "nlme_csh") %>%
+  # We don't run the SAS routines here.
+  # add_method(proc_mixed_us_meth, name = "proc_mixed_us") %>%
+  # add_method(proc_mixed_csh_meth, name = "proc_mixed_csh") %>%
+  # add_method(proc_mixed_toeph_meth, name = "proc_mixed_toeph") %>%
+  add_evaluator(mean_time_eval, name = "mean_fit_time") %>%
+  add_evaluator(bias_eval, name = "bias") %>%
+  add_evaluator(variance_eval, name = "variance") %>%
+  add_evaluator(convergence_rate_eval, name = "convergence_rate") %>%
+  add_evaluator(coverage_eval, name = "coverage_rate") %>%
+  add_evaluator(type_1_error_rate_eval, name = "type_1_error_rate") %>%
+  add_evaluator(type_2_error_rate_eval, name = "type_2_error_rate")
+
+# No need for SAS here
+
+# # Input here the name of the SAS container.
+# hostname <- "sabanesd-kjii4i-eu.ocean"
+# # ssh::ssh_connect(hostname)
+#
+# # We can also be specific if there are multiple SAS containers we want to connect to.
+# cfg_name <- "sascfg_personal_2.py"
+# sasr.roche::setup_sasr_ocean(hostname, sascfg = cfg_name)
+# options(sascfg = cfg_name)
+#
+# df2sd(mtcars, "mt")
+#
+# result <- run_sas("
+#   proc freq data = mt;
+#   run;
+# ")
+
+# run the experiment
+# we break it down into 10 batches so that we stay within the memory limits
+
+exp_dir <- experiment$get_save_dir()
+if (!dir.exists(exp_dir)) {
+  dir.create(exp_dir)
+}
+
+for (i in 9:10) {
+  # run this batch
+  set.seed(2144 + i * 10)
+  results <- experiment$onlystart(
+    n_reps = 100,
+    verbose = 2
+  )
+
+  # format results
+  formatted_fit_results <- format_fit_results(results)
+
+  # correct repetition information
+  formatted_fit_results$rep <- as.integer(formatted_fit_results$rep) + 100L * (i - 1L)
+
+  # save this batch
+  file_name <- paste0("formatted_fit_results_", i, ".rds")
+  file_name <- file.path(exp_dir, file_name)
+  saveRDS(formatted_fit_results, file = file_name)
+
+  # clean memory
+  rm(formatted_fit_results)
+  rm(results)
+  gc()
+}
+
+# Finally aggregate all tibbles together:
+all_files <- grep(pattern = "^formatted_fit_results_", x = dir(exp_dir), value = TRUE)
+all_files <- file.path(exp_dir, all_files)
+all_results <- lapply(all_files, readRDS)
+result <- do.call(rbind, all_results)
+saveRDS(result, file = file.path(exp_dir, "result.rds"))
diff --git a/simulations/missing-data-benchmarks/R/meals/meal-high-missingness.R b/simulations/missing-data-benchmarks/R/meals/meal-high-missingness-SAS.R
similarity index 71%
rename from simulations/missing-data-benchmarks/R/meals/meal-high-missingness.R
rename to simulations/missing-data-benchmarks/R/meals/meal-high-missingness-SAS.R
index fed10058f..9596b2974 100644
--- a/simulations/missing-data-benchmarks/R/meals/meal-high-missingness.R
+++ b/simulations/missing-data-benchmarks/R/meals/meal-high-missingness-SAS.R
@@ -13,6 +13,9 @@ library(dplyr)
 library(purrr)
 library(tidyr)
 library(emmeans)
+library(clusterGeneration)
+library(ssh)
+library(microbenchmark)
 
 # source the R scripts
 sim_functions_files <- list.files(
@@ -96,26 +99,26 @@ mod_effect_toeph_dgp <- create_dgp(
 )
 
 # specify the methods
-mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
-mmrm_csh_meth <- create_method(
-  .method_fun = mmrm_wrapper_fun, covar_type = "csh"
-)
-mmrm_toeph_meth <- create_method(
-  .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
-)
-glmmtmb_us_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
-)
-glmmtmb_csh_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
-)
-glmmtmb_toeph_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
-)
-nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
-nlme_csh_meth <- create_method(
-  .method_fun = nlme_wrapper_fun, covar_type = "csh"
-)
+# mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
+# mmrm_csh_meth <- create_method(
+#   .method_fun = mmrm_wrapper_fun, covar_type = "csh"
+# )
+# mmrm_toeph_meth <- create_method(
+#   .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
+# )
+# glmmtmb_us_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
+# )
+# glmmtmb_csh_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
+# )
+# glmmtmb_toeph_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
+# )
+# nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
+# nlme_csh_meth <- create_method(
+#   .method_fun = nlme_wrapper_fun, covar_type = "csh"
+# )
 proc_mixed_us_meth <- create_method(
   .method_fun = proc_mixed_wrapper_fun, covar_type = "us"
 )
@@ -154,8 +157,8 @@ type_2_error_rate_eval <- create_evaluator(
 
 # create the experiment
 experiment <- create_experiment(
-  name = "mmrm-benchmark-high-missingness-n-600",
-  save_dir = "results/high-miss/n-600"
+  name = "mmrm-benchmark-high-missingness-n-600-SAS",
+  save_dir = "results/high-miss/n-600-SAS"
 ) %>%
   add_dgp(no_effect_us_dgp, name = "no_effect_us") %>%
   add_dgp(no_effect_csh_dgp, name = "no_effect_csh") %>%
@@ -166,14 +169,14 @@ experiment <- create_experiment(
   add_dgp(mod_effect_us_dgp, name = "mod_effect_us") %>%
   add_dgp(mod_effect_csh_dgp, name = "mod_effect_csh") %>%
   add_dgp(mod_effect_toeph_dgp, name = "mod_effect_toeph") %>%
-  add_method(mmrm_us_meth, name = "mmrm_us") %>%
-  add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
-  add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
-  add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
-  add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
-  add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
-  add_method(nlme_us_meth, name = "nlme_us") %>%
-  add_method(nlme_csh_meth, name = "nlme_csh") %>%
+  # add_method(mmrm_us_meth, name = "mmrm_us") %>%
+  # add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
+  # add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
+  # add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
+  # add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
+  # add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
+  # add_method(nlme_us_meth, name = "nlme_us") %>%
+  # add_method(nlme_csh_meth, name = "nlme_csh") %>%
   add_method(proc_mixed_us_meth, name = "proc_mixed_us") %>%
   add_method(proc_mixed_csh_meth, name = "proc_mixed_csh") %>%
   add_method(proc_mixed_toeph_meth, name = "proc_mixed_toeph") %>%
@@ -185,11 +188,31 @@ experiment <- create_experiment(
   add_evaluator(type_1_error_rate_eval, name = "type_1_error_rate") %>%
   add_evaluator(type_2_error_rate_eval, name = "type_2_error_rate")
 
+# # Input here the name of the SAS container.
+# hostname <- "sabanesd-k2n72p-eu.ocean"
+#
+# # Do this to register host
+# session <- ssh::ssh_connect(hostname)
+# ssh_disconnect(session)
+#
+# # We can also be specific if there are multiple SAS containers we want to connect to.
+# cfg_name <- "sascfg_personal_2.py"
+# sasr.roche::setup_sasr_ocean(hostname, sascfg = cfg_name)
+# options(sascfg = cfg_name)
+#
+# # Test if it works.
+# df <- data.frame(a = 1, b = 2)
+# sasr::df2sd(df)
+
 # run the experiment
-set.seed(713452)
-results <- experiment$run(
-  n_reps = 100,
-  save = TRUE,
-  verbose = 2,
-  checkpoint_n_reps = 25
-)
+# set.seed(713452)
+# results <- experiment$run(
+#   n_reps = 1000,
+#   save = TRUE,
+#   verbose = 2,
+#   checkpoint_n_reps = 10
+# )
+
+source("R/format-replicate-results/helpers.R")
+
+format_fit_and_save(experiment)
diff --git a/simulations/missing-data-benchmarks/R/meals/meal-low-missingness-R.R b/simulations/missing-data-benchmarks/R/meals/meal-low-missingness-R.R
new file mode 100644
index 000000000..36fc7c94a
--- /dev/null
+++ b/simulations/missing-data-benchmarks/R/meals/meal-low-missingness-R.R
@@ -0,0 +1,258 @@
+################################################################################
+# Simulation with Low Levels of Missingness
+################################################################################
+
+# load required libraries
+library(simChef)
+library(mmrm)
+library(glmmTMB)
+library(nlme)
+library(sasr)
+library(stringr)
+library(dplyr)
+library(purrr)
+library(tidyr)
+library(emmeans)
+library(clusterGeneration)
+library(ssh)
+library(microbenchmark)
+library(future.batchtools)
+
+# source the R scripts
+sim_functions_files <- list.files(
+  c("R/dgp", "R/method", "R/eval", "R/viz", "R/customizations"),
+  pattern = "*.R$", full.names = TRUE, ignore.case = TRUE
+)
+sapply(sim_functions_files, source)
+
+# generate the possible covariance matrices
+us_cov_mat <- compute_unstructured_matrix()
+csh_cov_mat <- compute_csh_matrix()
+toep_cov_mat <- toeplitz(c(1, 0.5, 0.25, 0.125, rep(0, 6)))
+
+# set the sample size
+n_obs <- 600
+
+# dgps with no treatment effect
+no_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  missingness = "low"
+)
+no_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  missingness = "low"
+)
+no_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  missingness = "low"
+)
+
+# dgps with small treatment effect
+small_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "low"
+)
+small_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "low"
+)
+small_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  trt_visit_coef = 0.25,
+  missingness = "low"
+)
+
+# dgps with moderate treatment effect
+mod_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "low"
+)
+mod_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "low"
+)
+mod_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat,
+  trt_visit_coef = 0.5,
+  missingness = "low"
+)
+
+# specify the methods
+mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
+mmrm_csh_meth <- create_method(
+  .method_fun = mmrm_wrapper_fun, covar_type = "csh"
+)
+mmrm_toeph_meth <- create_method(
+  .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
+)
+glmmtmb_us_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
+)
+glmmtmb_csh_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
+)
+glmmtmb_toeph_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
+)
+nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
+nlme_csh_meth <- create_method(
+  .method_fun = nlme_wrapper_fun, covar_type = "csh"
+)
+# proc_mixed_us_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "us"
+# )
+# proc_mixed_csh_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "csh"
+# )
+# proc_mixed_toeph_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "toeph"
+# )
+
+# specify the evaluation metrics
+mean_time_eval <- create_evaluator(.eval_fun = mean_time_fun)
+true_params <- list(
+  "no_effect_us" = rep(0, 10),
+  "no_effect_csh" = rep(0, 10),
+  "no_effect_toeph" = rep(0, 10),
+  "small_effect_us" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "small_effect_csh" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "small_effect_toeph" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "mod_effect_us" = seq(from = 0.5, by = 0.5, length.out = 10),
+  "mod_effect_csh" = seq(from = 0.5, by = 0.5, length.out = 10),
+  "mod_effect_toeph" = seq(from = 0.5, by = 0.5, length.out = 10)
+)
+bias_eval <- create_evaluator(.eval_fun = bias_fun, true_params = true_params)
+variance_eval <- create_evaluator(.eval_fun = variance_fun)
+convergence_rate_eval <- create_evaluator(.eval_fun = convergence_rate_fun)
+coverage_eval <- create_evaluator(
+  .eval_fun = coverage_fun, true_params = true_params
+)
+type_1_error_rate_eval <- create_evaluator(
+  .eval_fun = type_1_error_rate_fun, true_params = true_params
+)
+type_2_error_rate_eval <- create_evaluator(
+  .eval_fun = type_2_error_rate_fun, true_params = true_params
+)
+
+# fixInNamespace
+fixInNamespace("readLog", "batchtools")
+
+# Define parallelization
+future_globals <- ls()
+
+options(future.globals.onReference = NULL)
+# options(future.globals.onReference = "error")
+
+plan(
+  batchtools_lsf,
+  workers = 100,
+  resources = list(memory = 10000, walltime = 5 * 60 * 60) # queue = "short",
+)
+
+# create the experiment
+experiment <- create_experiment(
+  name = "mmrm-benchmark-low-missingness-n-600-R",
+  save_dir = "results/low-miss/n-600-R",
+  future.globals = future_globals
+) %>%
+  add_dgp(no_effect_us_dgp, name = "no_effect_us") %>%
+  add_dgp(no_effect_csh_dgp, name = "no_effect_csh") %>%
+  add_dgp(no_effect_toeph_dgp, name = "no_effect_toeph") %>%
+  add_dgp(small_effect_us_dgp, name = "small_effect_us") %>%
+  add_dgp(small_effect_csh_dgp, name = "small_effect_csh") %>%
+  add_dgp(small_effect_toeph_dgp, name = "small_effect_toeph") %>%
+  add_dgp(mod_effect_us_dgp, name = "mod_effect_us") %>%
+  add_dgp(mod_effect_csh_dgp, name = "mod_effect_csh") %>%
+  add_dgp(mod_effect_toeph_dgp, name = "mod_effect_toeph") %>%
+  add_method(mmrm_us_meth, name = "mmrm_us") %>%
+  add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
+  add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
+  add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
+  add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
+  add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
+  add_method(nlme_us_meth, name = "nlme_us") %>%
+  add_method(nlme_csh_meth, name = "nlme_csh") %>%
+  # We don't run the SAS routines here.
+  # add_method(proc_mixed_us_meth, name = "proc_mixed_us") %>%
+  # add_method(proc_mixed_csh_meth, name = "proc_mixed_csh") %>%
+  # add_method(proc_mixed_toeph_meth, name = "proc_mixed_toeph") %>%
+  add_evaluator(mean_time_eval, name = "mean_fit_time") %>%
+  add_evaluator(bias_eval, name = "bias") %>%
+  add_evaluator(variance_eval, name = "variance") %>%
+  add_evaluator(convergence_rate_eval, name = "convergence_rate") %>%
+  add_evaluator(coverage_eval, name = "coverage_rate") %>%
+  add_evaluator(type_1_error_rate_eval, name = "type_1_error_rate") %>%
+  add_evaluator(type_2_error_rate_eval, name = "type_2_error_rate")
+
+# No need for SAS here
+
+# # Input here the name of the SAS container.
+# hostname <- "sabanesd-3qamsv-eu.ocean"
+# ssh::ssh_connect(hostname)
+#
+# # We can also be specific if there are multiple SAS containers we want to connect to.
+# cfg_name <- "sascfg_personal_3.py"
+# sasr.roche::setup_sasr_ocean(hostname, sascfg = cfg_name)
+# options(sascfg = cfg_name)
+
+# run the experiment
+# we break it down into 10 batches so that we stay within the memory limits
+
+exp_dir <- experiment$get_save_dir()
+if (!dir.exists(exp_dir)) {
+  dir.create(exp_dir)
+}
+
+for (i in 1:10) {
+  # run this batch
+  set.seed(1341 + i * 10)
+  results <- experiment$onlystart(
+    n_reps = 100,
+    verbose = 2
+  )
+
+  # format results
+  formatted_fit_results <- format_fit_results(results)
+
+  # correct repetition information
+  formatted_fit_results$rep <- as.integer(formatted_fit_results$rep) + 100L * (i - 1L)
+
+  # save this batch
+  file_name <- paste0("formatted_fit_results_", i, ".rds")
+  file_name <- file.path(exp_dir, file_name)
+  saveRDS(formatted_fit_results, file = file_name)
+
+  # clean memory
+  rm(formatted_fit_results)
+  rm(results)
+  gc()
+}
+
+# Finally aggregate all tibbles together:
+all_files <- grep(pattern = "^formatted_fit_results_", x = dir(exp_dir), value = TRUE)
+all_files <- file.path(exp_dir, all_files)
+all_results <- lapply(all_files, readRDS)
+result <- do.call(rbind, all_results)
+saveRDS(result, file = file.path(exp_dir, "result.rds"))
diff --git a/simulations/missing-data-benchmarks/R/meals/meal-low-missingness.R b/simulations/missing-data-benchmarks/R/meals/meal-low-missingness-SAS.R
similarity index 71%
rename from simulations/missing-data-benchmarks/R/meals/meal-low-missingness.R
rename to simulations/missing-data-benchmarks/R/meals/meal-low-missingness-SAS.R
index 6ddd383db..2ac74258a 100644
--- a/simulations/missing-data-benchmarks/R/meals/meal-low-missingness.R
+++ b/simulations/missing-data-benchmarks/R/meals/meal-low-missingness-SAS.R
@@ -13,6 +13,9 @@ library(dplyr)
 library(purrr)
 library(tidyr)
 library(emmeans)
+library(clusterGeneration)
+library(ssh)
+library(microbenchmark)
 
 # source the R scripts
 sim_functions_files <- list.files(
@@ -96,26 +99,26 @@ mod_effect_toeph_dgp <- create_dgp(
 )
 
 # specify the methods
-mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
-mmrm_csh_meth <- create_method(
-  .method_fun = mmrm_wrapper_fun, covar_type = "csh"
-)
-mmrm_toeph_meth <- create_method(
-  .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
-)
-glmmtmb_us_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
-)
-glmmtmb_csh_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
-)
-glmmtmb_toeph_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
-)
-nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
-nlme_csh_meth <- create_method(
-  .method_fun = nlme_wrapper_fun, covar_type = "csh"
-)
+# mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
+# mmrm_csh_meth <- create_method(
+#   .method_fun = mmrm_wrapper_fun, covar_type = "csh"
+# )
+# mmrm_toeph_meth <- create_method(
+#   .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
+# )
+# glmmtmb_us_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
+# )
+# glmmtmb_csh_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
+# )
+# glmmtmb_toeph_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
+# )
+# nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
+# nlme_csh_meth <- create_method(
+#   .method_fun = nlme_wrapper_fun, covar_type = "csh"
+# )
 proc_mixed_us_meth <- create_method(
   .method_fun = proc_mixed_wrapper_fun, covar_type = "us"
 )
@@ -154,8 +157,8 @@ type_2_error_rate_eval <- create_evaluator(
 
 # create the experiment
 experiment <- create_experiment(
-  name = "mmrm-benchmark-low-missingness-n-600",
-  save_dir = "results/low-miss/n-600"
+  name = "mmrm-benchmark-low-missingness-n-600-SAS",
+  save_dir = "results/low-miss/n-600-SAS"
 ) %>%
   add_dgp(no_effect_us_dgp, name = "no_effect_us") %>%
   add_dgp(no_effect_csh_dgp, name = "no_effect_csh") %>%
@@ -166,14 +169,14 @@ experiment <- create_experiment(
   add_dgp(mod_effect_us_dgp, name = "mod_effect_us") %>%
   add_dgp(mod_effect_csh_dgp, name = "mod_effect_csh") %>%
   add_dgp(mod_effect_toeph_dgp, name = "mod_effect_toeph") %>%
-  add_method(mmrm_us_meth, name = "mmrm_us") %>%
-  add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
-  add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
-  add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
-  add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
-  add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
-  add_method(nlme_us_meth, name = "nlme_us") %>%
-  add_method(nlme_csh_meth, name = "nlme_csh") %>%
+  # add_method(mmrm_us_meth, name = "mmrm_us") %>%
+  # add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
+  # add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
+  # add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
+  # add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
+  # add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
+  # add_method(nlme_us_meth, name = "nlme_us") %>%
+  # add_method(nlme_csh_meth, name = "nlme_csh") %>%
   add_method(proc_mixed_us_meth, name = "proc_mixed_us") %>%
   add_method(proc_mixed_csh_meth, name = "proc_mixed_csh") %>%
   add_method(proc_mixed_toeph_meth, name = "proc_mixed_toeph") %>%
@@ -185,11 +188,31 @@ experiment <- create_experiment(
   add_evaluator(type_1_error_rate_eval, name = "type_1_error_rate") %>%
   add_evaluator(type_2_error_rate_eval, name = "type_2_error_rate")
 
+# # Input here the name of the SAS container.
+# hostname <- "sabanesd-byknxg-eu.ocean"
+#
+# # Do this to register host
+# session <- ssh::ssh_connect(hostname)
+# ssh_disconnect(session)
+#
+# # We can also be specific if there are multiple SAS containers we want to connect to.
+# cfg_name <- "sascfg_personal_3.py"
+# sasr.roche::setup_sasr_ocean(hostname, sascfg = cfg_name)
+# options(sascfg = cfg_name)
+#
+# # Test if it works.
+# df <- data.frame(a = 1, b = 2)
+# sasr::df2sd(df)
+
 # run the experiment
-set.seed(56129)
-results <- experiment$run(
-  n_reps = 100,
-  save = TRUE,
-  verbose = 2,
-  checkpoint_n_reps = 25
-)
+# set.seed(56129)
+# results <- experiment$run(
+#   n_reps = 1000,
+#   save = TRUE,
+#   verbose = 2,
+#   checkpoint_n_reps = 10
+# )
+
+source("R/format-replicate-results/helpers.R")
+
+format_fit_and_save(experiment)
diff --git a/simulations/missing-data-benchmarks/R/meals/meal-no-missingness-R.R b/simulations/missing-data-benchmarks/R/meals/meal-no-missingness-R.R
new file mode 100644
index 000000000..ee072081a
--- /dev/null
+++ b/simulations/missing-data-benchmarks/R/meals/meal-no-missingness-R.R
@@ -0,0 +1,244 @@
+################################################################################
+# Simulation with No Missingness
+################################################################################
+
+# load required libraries
+library(simChef)
+library(mmrm)
+library(glmmTMB)
+library(nlme)
+library(sasr)
+library(stringr)
+library(dplyr)
+library(purrr)
+library(tidyr)
+library(emmeans)
+library(clusterGeneration)
+library(ssh)
+library(microbenchmark)
+library(future)
+library(future.batchtools)
+
+# source the R scripts
+sim_functions_files <- list.files(
+  c("R/dgp", "R/method", "R/eval", "R/viz", "R/customizations"),
+  pattern = "*.R$", full.names = TRUE, ignore.case = TRUE
+)
+sapply(sim_functions_files, source)
+
+# generate the possible covariance matrices
+us_cov_mat <- compute_unstructured_matrix()
+csh_cov_mat <- compute_csh_matrix()
+toep_cov_mat <- toeplitz(c(1, 0.5, 0.25, 0.125, rep(0, 6)))
+
+# set the sample size
+n_obs <- 600
+
+# dgps with no treatment effect
+no_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = us_cov_mat
+)
+no_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = csh_cov_mat
+)
+no_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  n_obs = n_obs,
+  outcome_covar_mat = toep_cov_mat
+)
+
+# dgps with small treatment effect
+small_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  outcome_covar_mat = us_cov_mat,
+  n_obs = n_obs,
+  trt_visit_coef = 0.25
+)
+small_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  outcome_covar_mat = csh_cov_mat,
+  n_obs = n_obs,
+  trt_visit_coef = 0.25
+)
+small_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  outcome_covar_mat = toep_cov_mat,
+  n_obs = n_obs,
+  trt_visit_coef = 0.25
+)
+
+# dgps with moderate treatment effect
+mod_effect_us_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  outcome_covar_mat = us_cov_mat,
+  n_obs = n_obs,
+  trt_visit_coef = 0.5
+)
+mod_effect_csh_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  outcome_covar_mat = csh_cov_mat,
+  n_obs = n_obs,
+  trt_visit_coef = 0.5
+)
+mod_effect_toeph_dgp <- create_dgp(
+  .dgp_fun = rct_dgp_fun,
+  outcome_covar_mat = toep_cov_mat,
+  n_obs = n_obs,
+  trt_visit_coef = 0.5
+)
+
+# specify the methods
+mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
+mmrm_csh_meth <- create_method(
+  .method_fun = mmrm_wrapper_fun, covar_type = "csh"
+)
+mmrm_toeph_meth <- create_method(
+  .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
+)
+glmmtmb_us_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
+)
+glmmtmb_csh_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
+)
+glmmtmb_toeph_meth <- create_method(
+  .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
+)
+nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
+nlme_csh_meth <- create_method(
+  .method_fun = nlme_wrapper_fun, covar_type = "csh"
+)
+# proc_mixed_us_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "us"
+# )
+# proc_mixed_csh_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "csh"
+# )
+# proc_mixed_toeph_meth <- create_method(
+#   .method_fun = proc_mixed_wrapper_fun, covar_type = "toeph"
+# )
+
+# specify the evaluation metrics
+mean_time_eval <- create_evaluator(.eval_fun = mean_time_fun)
+true_params <- list(
+  "no_effect_us" = rep(0, 10),
+  "no_effect_csh" = rep(0, 10),
+  "no_effect_toeph" = rep(0, 10),
+  "small_effect_us" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "small_effect_csh" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "small_effect_toeph" = seq(from = 0.25, by = 0.25, length.out = 10),
+  "mod_effect_us" = seq(from = 0.5, by = 0.5, length.out = 10),
+  "mod_effect_csh" = seq(from = 0.5, by = 0.5, length.out = 10),
+  "mod_effect_toeph" = seq(from = 0.5, by = 0.5, length.out = 10)
+)
+bias_eval <- create_evaluator(.eval_fun = bias_fun, true_params = true_params)
+variance_eval <- create_evaluator(.eval_fun = variance_fun)
+convergence_rate_eval <- create_evaluator(.eval_fun = convergence_rate_fun)
+coverage_eval <- create_evaluator(
+  .eval_fun = coverage_fun, true_params = true_params
+)
+type_1_error_rate_eval <- create_evaluator(
+  .eval_fun = type_1_error_rate_fun, true_params = true_params
+)
+type_2_error_rate_eval <- create_evaluator(
+  .eval_fun = type_2_error_rate_fun, true_params = true_params
+)
+
+# fixInNamespace
+fixInNamespace("readLog", "batchtools")
+
+# Define parallelization
+future_globals <- ls()
+
+options(future.globals.onReference = NULL)
+# options(future.globals.onReference = "error")
+
+plan(
+  batchtools_lsf,
+  workers = 100,
+  resources = list(memory = 10000, walltime = 5 * 60 * 60) # queue = "short",
+)
+
+# create the experiment
+experiment <- create_experiment(
+  name = "mmrm-benchmark-no-missingness-n-600-R",
+  save_dir = "results/no-miss/n-600-R",
+  future.globals = future_globals
+) %>%
+  add_dgp(no_effect_us_dgp, name = "no_effect_us") %>%
+  add_dgp(no_effect_csh_dgp, name = "no_effect_csh") %>%
+  add_dgp(no_effect_toeph_dgp, name = "no_effect_toeph") %>%
+  add_dgp(small_effect_us_dgp, name = "small_effect_us") %>%
+  add_dgp(small_effect_csh_dgp, name = "small_effect_csh") %>%
+  add_dgp(small_effect_toeph_dgp, name = "small_effect_toeph") %>%
+  add_dgp(mod_effect_us_dgp, name = "mod_effect_us") %>%
+  add_dgp(mod_effect_csh_dgp, name = "mod_effect_csh") %>%
+  add_dgp(mod_effect_toeph_dgp, name = "mod_effect_toeph") %>%
+  add_method(mmrm_us_meth, name = "mmrm_us") %>%
+  add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
+  add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
+  add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
+  add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
+  add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
+  add_method(nlme_us_meth, name = "nlme_us") %>%
+  add_method(nlme_csh_meth, name = "nlme_csh") %>%
+  # add_method(proc_mixed_us_meth, name = "proc_mixed_us") %>%
+  # add_method(proc_mixed_csh_meth, name = "proc_mixed_csh") %>%
+  # add_method(proc_mixed_toeph_meth, name = "proc_mixed_toeph") %>%
+  add_evaluator(mean_time_eval, name = "mean_fit_time") %>%
+  add_evaluator(bias_eval, name = "bias") %>%
+  add_evaluator(variance_eval, name = "variance") %>%
+  add_evaluator(convergence_rate_eval, name = "convergence_rate") %>%
+  add_evaluator(coverage_eval, name = "coverage_rate") %>%
+  add_evaluator(type_1_error_rate_eval, name = "type_1_error_rate") %>%
+  add_evaluator(type_2_error_rate_eval, name = "type_2_error_rate")
+
+# # Input here the name of the SAS container.
+# hostname <- "sabanesd-trvhpl-eu.ocean"
+# ssh::ssh_connect(hostname)
+#
+# # We can also be specific if there are multiple SAS containers we want to connect to.
+# cfg_name <- "sascfg_personal_4.py"
+# sasr.roche::setup_sasr_ocean(hostname, sascfg = cfg_name)
+# options(sascfg = cfg_name)
+
+exp_dir <- experiment$get_save_dir()
+if (!dir.exists(exp_dir)) {
+  dir.create(exp_dir)
+}
+
+for (i in 1:10) {
+  # run this batch
+  set.seed(1341 + i * 10)
+  results <- experiment$onlystart(
+    n_reps = 100,
+    verbose = 2
+  )
+
+  # format results
+  formatted_fit_results <- format_fit_results(results)
+
+  # correct repetition information
+  formatted_fit_results$rep <- as.integer(formatted_fit_results$rep) + 100L * (i - 1L)
+
+  # save this batch
+  file_name <- paste0("formatted_fit_results_", i, ".rds")
+  file_name <- file.path(exp_dir, file_name)
+  saveRDS(formatted_fit_results, file = file_name)
+
+  # clean memory
+  rm(formatted_fit_results)
+  rm(results)
+  gc()
+}
+
+# Finally aggregate all tibbles together:
+all_files <- grep(pattern = "^formatted_fit_results_", x = dir(exp_dir), value = TRUE)
+all_files <- file.path(exp_dir, all_files)
+all_results <- lapply(all_files, readRDS)
+result <- do.call(rbind, all_results)
+saveRDS(result, file = file.path(exp_dir, "result.rds"))
diff --git a/simulations/missing-data-benchmarks/R/meals/meal-no-missingness.R b/simulations/missing-data-benchmarks/R/meals/meal-no-missingness-SAS.R
similarity index 69%
rename from simulations/missing-data-benchmarks/R/meals/meal-no-missingness.R
rename to simulations/missing-data-benchmarks/R/meals/meal-no-missingness-SAS.R
index 4621bdee3..f625c24dc 100644
--- a/simulations/missing-data-benchmarks/R/meals/meal-no-missingness.R
+++ b/simulations/missing-data-benchmarks/R/meals/meal-no-missingness-SAS.R
@@ -13,6 +13,9 @@ library(dplyr)
 library(purrr)
 library(tidyr)
 library(emmeans)
+library(clusterGeneration)
+library(ssh)
+library(microbenchmark)
 
 # source the R scripts
 sim_functions_files <- list.files(
@@ -87,26 +90,26 @@ mod_effect_toeph_dgp <- create_dgp(
 )
 
 # specify the methods
-mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
-mmrm_csh_meth <- create_method(
-  .method_fun = mmrm_wrapper_fun, covar_type = "csh"
-)
-mmrm_toeph_meth <- create_method(
-  .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
-)
-glmmtmb_us_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
-)
-glmmtmb_csh_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
-)
-glmmtmb_toeph_meth <- create_method(
-  .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
-)
-nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
-nlme_csh_meth <- create_method(
-  .method_fun = nlme_wrapper_fun, covar_type = "csh"
-)
+# mmrm_us_meth <- create_method(.method_fun = mmrm_wrapper_fun, covar_type = "us")
+# mmrm_csh_meth <- create_method(
+#   .method_fun = mmrm_wrapper_fun, covar_type = "csh"
+# )
+# mmrm_toeph_meth <- create_method(
+#   .method_fun = mmrm_wrapper_fun, covar_type = "toeph"
+# )
+# glmmtmb_us_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "us"
+# )
+# glmmtmb_csh_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "csh"
+# )
+# glmmtmb_toeph_meth <- create_method(
+#   .method_fun = glmmtmb_wrapper_fun, covar_type = "toeph"
+# )
+# nlme_us_meth <- create_method(.method_fun = nlme_wrapper_fun, covar_type = "us")
+# nlme_csh_meth <- create_method(
+#   .method_fun = nlme_wrapper_fun, covar_type = "csh"
+# )
 proc_mixed_us_meth <- create_method(
   .method_fun = proc_mixed_wrapper_fun, covar_type = "us"
 )
@@ -145,8 +148,8 @@ type_2_error_rate_eval <- create_evaluator(
 
 # create the experiment
 experiment <- create_experiment(
-  name = "mmrm-benchmark-no-missingness-n-600",
-  save_dir = "results/no-miss/n-600"
+  name = "mmrm-benchmark-no-missingness-n-600-SAS",
+  save_dir = "results/no-miss/n-600-SAS"
 ) %>%
   add_dgp(no_effect_us_dgp, name = "no_effect_us") %>%
   add_dgp(no_effect_csh_dgp, name = "no_effect_csh") %>%
@@ -157,14 +160,14 @@ experiment <- create_experiment(
   add_dgp(mod_effect_us_dgp, name = "mod_effect_us") %>%
   add_dgp(mod_effect_csh_dgp, name = "mod_effect_csh") %>%
   add_dgp(mod_effect_toeph_dgp, name = "mod_effect_toeph") %>%
-  add_method(mmrm_us_meth, name = "mmrm_us") %>%
-  add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
-  add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
-  add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
-  add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
-  add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
-  add_method(nlme_us_meth, name = "nlme_us") %>%
-  add_method(nlme_csh_meth, name = "nlme_csh") %>%
+  # add_method(mmrm_us_meth, name = "mmrm_us") %>%
+  # add_method(mmrm_csh_meth, name = "mmrm_csh") %>%
+  # add_method(mmrm_toeph_meth, name = "mmrm_toeph") %>%
+  # add_method(glmmtmb_us_meth, name = "glmmtmb_us") %>%
+  # add_method(glmmtmb_csh_meth, name = "glmmtmb_csh") %>%
+  # add_method(glmmtmb_toeph_meth, name = "glmmtmb_toeph") %>%
+  # add_method(nlme_us_meth, name = "nlme_us") %>%
+  # add_method(nlme_csh_meth, name = "nlme_csh") %>%
   add_method(proc_mixed_us_meth, name = "proc_mixed_us") %>%
   add_method(proc_mixed_csh_meth, name = "proc_mixed_csh") %>%
   add_method(proc_mixed_toeph_meth, name = "proc_mixed_toeph") %>%
@@ -175,10 +178,33 @@ experiment <- create_experiment(
   add_evaluator(coverage_eval, name = "coverage_rate") %>%
   add_evaluator(type_1_error_rate_eval, name = "type_1_error_rate") %>%
   add_evaluator(type_2_error_rate_eval, name = "type_2_error_rate")
+#
+# # Input here the name of the SAS container.
+# hostname <- "sabanesd-ddi1zu-eu.ocean"
+#
+# # Do this to register host
+# session <- ssh::ssh_connect(hostname)
+# ssh_disconnect(session)
+#
+# # We can also be specific if there are multiple SAS containers we want to connect to.
+# cfg_name <- "sascfg_personal_4.py"
+# sasr.roche::setup_sasr_ocean(hostname, sascfg = cfg_name)
+# options(sascfg = cfg_name)
+#
+# # Test if it works.
+# df <- data.frame(a = 1, b = 2)
+# sasr::df2sd(df)
+#
+# # run the experiment
+# set.seed(62342)
+# results <- experiment$run(
+#   n_reps = 1000,
+#   save = TRUE,
+#   verbose = 2,
+#   checkpoint_n_reps = 10
+# )
+
+source("R/format-replicate-results/helpers.R")
+
+format_fit_and_save(experiment)
 
-# run the experiment
-set.seed(62342)
-results <- experiment$run(
-  n_reps = 100,
-  save = TRUE
-)
diff --git a/simulations/missing-data-benchmarks/R/method/proc-mixed-wrapper.R b/simulations/missing-data-benchmarks/R/method/proc-mixed-wrapper.R
index 0e4d690ae..42c4efd23 100644
--- a/simulations/missing-data-benchmarks/R/method/proc-mixed-wrapper.R
+++ b/simulations/missing-data-benchmarks/R/method/proc-mixed-wrapper.R
@@ -46,7 +46,7 @@ proc_mixed_wrapper_fun <- function(
   ## specify the SAS code: only return covariance matrix estimates for
   ## repeated measures
   if (covar_type == "csh") {
-    sas_code <- "ODS OUTPUT LSMEANS = lsmeans_out DIFFS = diffs_out
+    sas_code <- "ODS OUTPUT LSMEANS = lsmeans_out DIFFS = diffs_out R = covmat_out
       ModelInfo = model_info ConvergenceStatus = conv_status;
     PROC MIXED DATA = sas_df METHOD = reml;
       CLASS trt(ref = '0') visit_num strata participant;
@@ -55,7 +55,7 @@ proc_mixed_wrapper_fun <- function(
       LSMEANS visit_num*trt / pdiff slice=visit_num cl alpha = 0.05 OBSMARGINS;
     RUN;"
   } else if (covar_type == "toeph") {
-    sas_code <- "ODS OUTPUT LSMEANS = lsmeans_out DIFFS = diffs_out
+    sas_code <- "ODS OUTPUT LSMEANS = lsmeans_out DIFFS = diffs_out R = covmat_out
       ModelInfo = model_info ConvergenceStatus = conv_status;
     PROC MIXED DATA = sas_df METHOD = reml;
       CLASS trt(ref = '0') visit_num strata participant;
@@ -64,7 +64,7 @@ proc_mixed_wrapper_fun <- function(
       LSMEANS visit_num*trt / pdiff slice=visit_num cl alpha = 0.05 OBSMARGINS;
     RUN;"
   } else if (covar_type == "us") {
-    sas_code <- "ODS OUTPUT LSMEANS = lsmeans_out DIFFS = diffs_out
+    sas_code <- "ODS OUTPUT LSMEANS = lsmeans_out DIFFS = diffs_out R = covmat_out
       ModelInfo = model_info ConvergenceStatus = conv_status;
     PROC MIXED DATA = sas_df METHOD = reml;
       CLASS trt(ref = '0') visit_num strata participant;
@@ -96,6 +96,12 @@ proc_mixed_wrapper_fun <- function(
       ## extract the model information
       mod_inf_df <- sasr::sd2df("model_info")
 
+      ## extract the estimated covariance matrix
+      covmat_df <- sasr::sd2df("covmat_out")
+      ## remove the first two columns which are just "Index" and "Row",
+      ## the remaining part is the matrix we need
+      covmat <- as.matrix(covmat_df[, -c(1, 2)])
+
       ## extract the ATEs across visits
       ates_df <- sasr::sd2df("diffs_out") %>%
         dplyr::filter(visit_num == `_visit_num`) %>%
@@ -107,11 +113,13 @@ proc_mixed_wrapper_fun <- function(
       mod_inf_df <- data.frame(NA)
       lsmeans_df <- data.frame(NA)
       ates_df <- data.frame(NA)
+      covmat <- matrix(NA)
     }
   } else {
     mod_inf_df <- data.frame(NA)
     lsmeans_df <- data.frame(NA)
     ates_df <- data.frame(NA)
+    covmat <- matrix(NA)
     converged <- FALSE
   }
 
@@ -120,13 +128,16 @@ proc_mixed_wrapper_fun <- function(
     stringr::str_extract("(?<=user cpu time)\\s*[0-9.]+") %>%
     as.numeric()
 
-  return(list(
-    "fit" = ates_df,
-    "fit_time" = fit_time,
-    "output" = sas_result$result$LST,
-    "lsmeans_df" = lsmeans_df,
-    "converged" = converged,
-    "mod_inf_df" = mod_inf_df,
-    "log" = sas_result$result$LOG
-  ))
+  list(
+    fit = structure(ates_df, covmat = covmat),
+    # Note: this is not really used below, but otherwise the downstream
+    # evaluation might fail if there is no `data` column.
+    data = df,
+    fit_time = fit_time,
+    output = sas_result$result$LST,
+    lsmeans_df = lsmeans_df,
+    converged = converged,
+    mod_inf_df = mod_inf_df,
+    log = sas_result$result$LOG
+  )
 }
diff --git a/simulations/missing-data-benchmarks/README.md b/simulations/missing-data-benchmarks/README.md
index 0f89c4564..2d74244a0 100644
--- a/simulations/missing-data-benchmarks/README.md
+++ b/simulations/missing-data-benchmarks/README.md
@@ -81,7 +81,7 @@ This directory contains the following folders:
   - `dgp/`: Code relating to DGP definition.
   - `eval/`: Functions for measuring empirical metrics and operating
      characteristics.
-  - `format-replicate-results/`: Functions to post-process existing tibble of fitted model objects to extract the `emmeans` results.
+  - `format-replicate-results/`: Functions to post-process existing tibble of fitted model objects to extract the `emmeans` results and fitted covariance matrices.
   - `meal.R`: Main R script to test the simulation study, only runs 2 repetitions.
   - `meals/`: Main R scripts for running the simulation studies. 
      Specify here the number of repetitions via `n_reps` at the end of each script.
diff --git a/simulations/missing-data-benchmarks/batchtools.lsf.tmpl.doesnotwork b/simulations/missing-data-benchmarks/batchtools.lsf.tmpl.doesnotwork
new file mode 100644
index 000000000..dac2ee765
--- /dev/null
+++ b/simulations/missing-data-benchmarks/batchtools.lsf.tmpl.doesnotwork
@@ -0,0 +1,6 @@
+#BSUB -J <%= if (exists("job.name", mode = "character")) job.name else job.hash %>
+#BSUB -q short
+#BSUB -n 1
+#BSUB -o <%= log.file %>
+singularity exec https://ross.science.roche.com/singularity-cbs/coderoche/r_container/bee_4.2.1.sif Rscript --vanilla 'batchtools::doJobCollection("<%= uri %>")'
+
diff --git a/simulations/missing-data-benchmarks/results/combine-df.R b/simulations/missing-data-benchmarks/results/combine-df.R
new file mode 100644
index 000000000..754eae0d0
--- /dev/null
+++ b/simulations/missing-data-benchmarks/results/combine-df.R
@@ -0,0 +1,18 @@
+# Let's finally combine all the results together.
+# copy from laptop to Ocean drive
+# combine all
+# share with Gonzalo
+
+ocean_path <- "/ocean/harbour/CDTpractice/PR002/demo-biostats-workflow/dev/data/other/sabanesd"
+sub_dirs <- c("extreme-miss", "high-miss", "low-miss", "no-miss")
+r_sas_dirs <- c("n-600-R", "n-600-SAS")
+
+combined_dirs <- expand.grid(sub_dirs, r_sas_dirs)
+combined_dirs <- apply(combined_dirs, 1L, paste, collapse = "/")
+combined_dirs <- file.path(ocean_path, combined_dirs, "formatted_fit_results.rds")
+
+df_list <- lapply(combined_dirs, readRDS)
+str(df_list, 1)
+
+length(unique(df_list[[7]]$rep))
+combined_dirs[7]
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/eval_results.rds b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/eval_results.rds
deleted file mode 100644
index 17319aad2..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/eval_results.rds and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/experiment.rds b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/experiment.rds
deleted file mode 100644
index 085d88a33..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/experiment.rds and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/experiment_cached_params.rds b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/experiment_cached_params.rds
deleted file mode 100644
index dc963a9ea..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/experiment_cached_params.rds and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_bias.jpeg b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_bias.jpeg
deleted file mode 100644
index 900e2fa1f..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_bias.jpeg and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_convergence_rate.jpeg b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_convergence_rate.jpeg
deleted file mode 100644
index 3c6ee62ed..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_convergence_rate.jpeg and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_coverage_rate.jpeg b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_coverage_rate.jpeg
deleted file mode 100644
index db7e7af14..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_coverage_rate.jpeg and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_mean_fit_time.jpeg b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_mean_fit_time.jpeg
deleted file mode 100644
index 025765138..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_mean_fit_time.jpeg and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_type_1_error_rate.jpeg b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_type_1_error_rate.jpeg
deleted file mode 100644
index 27b61b32e..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_type_1_error_rate.jpeg and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_type_2_error_rate.jpeg b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_type_2_error_rate.jpeg
deleted file mode 100644
index 23d813c8d..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_type_2_error_rate.jpeg and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_variance.jpeg b/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_variance.jpeg
deleted file mode 100644
index 9cfd826d3..000000000
Binary files a/simulations/missing-data-benchmarks/results/extreme-miss/n-600/result-plots/empirical_variance.jpeg and /dev/null differ
diff --git a/simulations/missing-data-benchmarks/sascfg_personal.py b/simulations/missing-data-benchmarks/sascfg_personal.py
new file mode 100644
index 000000000..b30e9eba8
--- /dev/null
+++ b/simulations/missing-data-benchmarks/sascfg_personal.py
@@ -0,0 +1,2 @@
+SAS_config_names=['default']
+default={'host': 'sabanesd-nwpfwh-eu.ocean', 'saspath': '/opt/sas/SASHome/SASFoundation/9.4/sas', 'ssh': '/usr/bin/ssh', 'encoding': 'latin1', 'options': ['-fullstimer'], 'tunnel': 9991, 'rtunnel': 9992}
diff --git a/simulations/missing-data-benchmarks/sascfg_personal_2.py b/simulations/missing-data-benchmarks/sascfg_personal_2.py
new file mode 100644
index 000000000..04307d0ca
--- /dev/null
+++ b/simulations/missing-data-benchmarks/sascfg_personal_2.py
@@ -0,0 +1,2 @@
+SAS_config_names=['default']
+default={'host': 'sabanesd-k2n72p-eu.ocean', 'saspath': '/opt/sas/SASHome/SASFoundation/9.4/sas', 'ssh': '/usr/bin/ssh', 'encoding': 'latin1', 'options': ['-fullstimer'], 'tunnel': 9991, 'rtunnel': 9992}
diff --git a/simulations/missing-data-benchmarks/sascfg_personal_3.py b/simulations/missing-data-benchmarks/sascfg_personal_3.py
new file mode 100644
index 000000000..1552b1c65
--- /dev/null
+++ b/simulations/missing-data-benchmarks/sascfg_personal_3.py
@@ -0,0 +1,2 @@
+SAS_config_names=['default']
+default={'host': 'sabanesd-byknxg-eu.ocean', 'saspath': '/opt/sas/SASHome/SASFoundation/9.4/sas', 'ssh': '/usr/bin/ssh', 'encoding': 'latin1', 'options': ['-fullstimer'], 'tunnel': 9991, 'rtunnel': 9992}
diff --git a/simulations/missing-data-benchmarks/sascfg_personal_4.py b/simulations/missing-data-benchmarks/sascfg_personal_4.py
new file mode 100644
index 000000000..299883fb2
--- /dev/null
+++ b/simulations/missing-data-benchmarks/sascfg_personal_4.py
@@ -0,0 +1,2 @@
+SAS_config_names=['default']
+default={'host': 'sabanesd-ddi1zu-eu.ocean', 'saspath': '/opt/sas/SASHome/SASFoundation/9.4/sas', 'ssh': '/usr/bin/ssh', 'encoding': 'latin1', 'options': ['-fullstimer'], 'tunnel': 9991, 'rtunnel': 9992}
diff --git a/simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-helpers.R b/simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-helpers.R
index 58c27ebc0..e2131c333 100644
--- a/simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-helpers.R
+++ b/simulations/missing-data-benchmarks/tests/testthat/eval-tests/test-helpers.R
@@ -507,3 +507,107 @@ test_that("get_mixed_emmeans_like_output extracts emmeans contrasts", {
   expected_cols <- c("visit_num", "estimate", "stderr", "df", "tvalue", "pvalue", "lower", "upper")
   expect_named(result, expected_cols)
 })
+
+# get_cov_mat_estimate ----
+
+test_that("get_mmrm_cov_mat_estimate works as expected", {
+  # generate data
+  set.seed(51235)
+  dgp_output <- rct_dgp_fun(
+    n_obs = 100,
+    outcome_covar_mat = diag(3),
+    trt_visit_coef = 0.5
+  )
+
+  # fit the mmrm with mmrm
+  mmrm_fit <- mmrm_wrapper_fun(
+    participant = dgp_output$participant,
+    visit_num = dgp_output$visit_num,
+    base_bcva = dgp_output$base_bcva,
+    strata = dgp_output$strata,
+    trt = dgp_output$trt,
+    bcva_change = dgp_output$bcva_change,
+    covar_type = "us"
+  )
+
+  result <- expect_silent(get_mmrm_cov_mat_estimate(mmrm_fit$fit))
+  expect_matrix(result, nrows = 3, ncols = 3)
+  expect_silent(chol(result))
+})
+
+test_that("get_glmm_cov_mat_estimate works as expected", {
+  # generate data
+  set.seed(51235)
+  dgp_output <- rct_dgp_fun(
+    n_obs = 100,
+    outcome_covar_mat = diag(3),
+    trt_visit_coef = 0.5
+  )
+
+  # fit the mmrm with mmrm
+  glmmtmb_fit <- glmmtmb_wrapper_fun(
+    participant = dgp_output$participant,
+    visit_num = dgp_output$visit_num,
+    base_bcva = dgp_output$base_bcva,
+    strata = dgp_output$strata,
+    trt = dgp_output$trt,
+    bcva_change = dgp_output$bcva_change,
+    covar_type = "us"
+  )
+
+  result <- expect_silent(get_glmm_cov_mat_estimate(glmmtmb_fit$fit))
+  expect_matrix(result, nrows = 3, ncols = 3)
+  expect_silent(chol(result))
+})
+
+test_that("get_nlme_cov_mat_estimate works as expected", {
+  # generate data
+  set.seed(51235)
+  dgp_output <- rct_dgp_fun(
+    n_obs = 100,
+    outcome_covar_mat = diag(3),
+    trt_visit_coef = 0.5
+  )
+
+  # fit the mmrm with mmrm
+  nlme_fit <- nlme_wrapper_fun(
+    participant = dgp_output$participant,
+    visit_num = dgp_output$visit_num,
+    base_bcva = dgp_output$base_bcva,
+    strata = dgp_output$strata,
+    trt = dgp_output$trt,
+    bcva_change = dgp_output$bcva_change,
+    covar_type = "us"
+  )
+
+  result <- expect_silent(get_nlme_cov_mat_estimate(nlme_fit$fit))
+  expect_matrix(result, nrows = 3, ncols = 3)
+  expect_silent(chol(result))
+})
+
+test_that("get_mixed_cov_mat_estimate works as expected", {
+  skip_if_not(run_sas_tests)
+  # generate data
+  set.seed(51235)
+  dgp_output <- rct_dgp_fun(
+    n_obs = 100,
+    outcome_covar_mat = diag(3),
+    trt_visit_coef = 0.5
+  )
+
+  # fit the mmrm with mmrm
+  proc_mixed_fit <- proc_mixed_wrapper_fun(
+    participant = dgp_output$participant,
+    visit_num = dgp_output$visit_num,
+    base_bcva = dgp_output$base_bcva,
+    strata = dgp_output$strata,
+    trt = dgp_output$trt,
+    bcva_change = dgp_output$bcva_change,
+    covar_type = "us"
+  )
+
+  # extract estimates
+  result <- expect_silent(get_mixed_cov_mat_estimate(proc_mixed_fit$fit))
+  expect_matrix(result, nrows = 3, ncols = 3)
+  expect_silent(chol(result))
+})
diff --git a/simulations/missing-data-benchmarks/tests/testthat/method-tests/test-proc-mixed-wrapper.R b/simulations/missing-data-benchmarks/tests/testthat/method-tests/test-proc-mixed-wrapper.R
index 9575faac6..eb40b0c31 100644
--- a/simulations/missing-data-benchmarks/tests/testthat/method-tests/test-proc-mixed-wrapper.R
+++ b/simulations/missing-data-benchmarks/tests/testthat/method-tests/test-proc-mixed-wrapper.R
@@ -19,6 +19,7 @@ test_that(paste(
     covar_type = "us"
   )
 
+  expect_true("data" %in% names(proc_mixed_fit))
   trt_effect_ests <- proc_mixed_fit$fit$Estimate
   expect_equal(trt_effect_ests, c(0, 0, 0), tolerance = 0.05)
 })