Performance indicators for the point estimators include the absolute and relative bias as well as the root mean squared error. For the interval estimators, we use the respective coverage rate.

Monte Carlo errors are calculated using the jackknife method.

# Clear workspace
rm(list = ls())

# Load packages
suppressPackageStartupMessages({
  library("here")
  library("parallel")
  library("pbapply")
})

# Load simulated data
simdat <- readRDS(here::here("data", "simdat.rds"))

1 Basic functions

#
# Calculate performance indicators for the simulation
#
# @param data   A data.frame with the results of the 
#                 simulation as returned by sim_run()
# @return       A data.frame with the estimated indicators
#
calc_indicators <- function(data) {
  
  out <- c()
  for (m in paste0("masem_", letters[1:6])) {

    o <- data.frame(
      
      # Condition
      cond = substr(m, 7, 7),
      
      # Percentage of replications with errors in drawing samples
      err_samp = mean(data$iter_samp > 0, na.rm = TRUE) * 100,
      
      # Percentage of replications with errors in MASEM
      err_masem1a = mean(data[["masem1a.status"]] > 1, na.rm = TRUE) * 100,
      err_masem2a = mean(data[["masem2a.status"]] > 1, na.rm = TRUE) * 100,
      err_masem1b = mean(data[["masem1b.status"]] > 1, na.rm = TRUE) * 100,
      err_masem2b = mean(data[["masem2b.status"]] > 1, na.rm = TRUE) * 100,
      err_masem1c = mean(data[["masem1c.status"]] > 1, na.rm = TRUE) * 100,
      err_masem2c = mean(data[["masem2c.status"]] > 1, na.rm = TRUE) * 100,
      err_masem0d = mean(data[["masem0d.status"]] > 1, na.rm = TRUE) * 100,
      err_masem0e = mean(data[["masem0e.status"]] > 1, na.rm = TRUE) * 100,
      err_masem0f = mean(data[["masem0f.status"]] > 1, na.rm = TRUE) * 100
    )
    
    # Summary statistics for each effect
    for (b in c("MONX", "WONX", "YONM", "YONW", 
                "WONM", "YONX", "IND", "IND1", "IND2")) {
      
      # Variable name containing effect
      v <- paste0(m, ".", b, "_est")
      if (!(v %in% names(data))) next
      
      # Raw bias
      o[[paste0("bias_", b)]] <- 
        mean(data[[v]], na.rm = TRUE) - mean(data[[b]])
      
      # Percent bias
      o[[paste0("pbias_", b)]] <- 
        o[[paste0("bias_", b)]] / mean(data[[b]]) * 100
        
      # Root mean squared error
      o[[paste0("rmse_", b)]] <- 
        mean((data[[v]] - data[[b]])^2, na.rm = TRUE)^.5
      
      # Coverage rate for confidence interval
      o[[paste0("co_", b)]] <- 
        mean(data[[b]] > data[[paste0(m, ".", b, "_lb")]] & 
             data[[b]] < data[[paste0(m, ".", b, "_ub")]], na.rm = TRUE) * 100

    }
    
    o[sapply(o, is.nan)] <- o[sapply(o, is.infinite)] <- NA
    
    if (length(o) > 4)
      out <- rbind(out, o)
  }
  rownames(out) <- out$cond
  
  return(out)
  
}


#
# Estimate standard errors of the performance
#   indicators using jackknife 
#
# @param data   A data.frame with the results of the 
#                 simulation as returned by sim_run()
# @param cores  The number of cores to use for parallel 
#                 processing (optional)
# @return       A data.frame with standard errors
# @require      parallel, pbapply
#
jk <- function(data, 
               cores = parallel::detectCores()) {

  # Calculate performance indicators repeatedly
  #  while excluding a single replication
  cl <-  parallel::makeCluster(cores)
  parallel::clusterExport(cl, c("calc_indicators", "data"))
  x_list <- 
    pbapply::pbsapply(
      X = seq_len(nrow(data)),
      FUN = \(i) { calc_indicators(data[-i, ])[, -1] },
      cl = cl,
      simplify = FALSE
    ) |>
    t() |>
    as.data.frame()
  parallel::stopCluster(cl)
  
  # Calculate standard errors
  n <- nrow(data)
  x <- array(unlist(x_list), dim = c(dim(x_list[[1]][[1]]), n))
  y <- array(t(apply(x, 1, rowMeans)), dim = dim(x))
  se <- t((apply((x - y)^2, 1, rowMeans) * (n - 1))^.5)
  
  # Set names
  p <- calc_indicators(data)
  se <- data.frame(cond = p$con, se)
  colnames(se) <- colnames(p)
  rownames(se) <- se$cond
  
  return(se)
}

2 Calculate indicators

# Performance indicators
pind <- lapply(simdat, calc_indicators)

# Monte Carlo errors
mce <- lapply(simdat, jk)

# Restructure results
pind_list <- list()
for (m in names(pind)) {
  
  pind_list[[m]] <- list(
    err = subset(pind[[m]], select = grepl("^err_", names(pind[[m]])))
  )
  
  for (i in c("bias", "pbias", "rmse", "co")) {
    
    # Estimate
    pind_list[[m]][[i]] <- subset(
      pind[[m]], 
      select = grepl(paste0("^", i, "_"), names(pind[[m]]))
    )
    
    # Error
    pind_list[[m]][[paste0(i, "_err")]] <- subset(
      mce[[m]],
      select = grepl(paste0("^", i, "_"), names(mce[[m]]))
    )
    
    # Confidence interval
    pind_list[[m]][[paste0(i, "_lb")]] <- 
      pind_list[[m]][[i]] - qnorm(.975) * pind_list[[m]][[paste0(i, "_err")]]
    pind_list[[m]][[paste0(i, "_ub")]] <- 
      pind_list[[m]][[i]] + qnorm(.975) * pind_list[[m]][[paste0(i, "_err")]]
    
  }
}
pind <- pind_list
rm(m, i, pind_list)

3 Save results

saveRDS(pind, file = here::here("data", "simind.rds"))
saveRDS(mce, file = here::here("data", "simind_err.rds"))