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
<- readRDS(here::here("data", "simdat.rds")) simdat
#
# 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
#
<- function(data) {
calc_indicators
<- c()
out for (m in paste0("masem_", letters[1:6])) {
<- data.frame(
o
# 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
<- paste0(m, ".", b, "_est")
v if (!(v %in% names(data))) next
# Raw bias
paste0("bias_", b)]] <-
o[[mean(data[[v]], na.rm = TRUE) - mean(data[[b]])
# Percent bias
paste0("pbias_", b)]] <-
o[[paste0("bias_", b)]] / mean(data[[b]]) * 100
o[[
# Root mean squared error
paste0("rmse_", b)]] <-
o[[mean((data[[v]] - data[[b]])^2, na.rm = TRUE)^.5
# Coverage rate for confidence interval
paste0("co_", b)]] <-
o[[mean(data[[b]] > data[[paste0(m, ".", b, "_lb")]] &
< data[[paste0(m, ".", b, "_ub")]], na.rm = TRUE) * 100
data[[b]]
}
sapply(o, is.nan)] <- o[sapply(o, is.infinite)] <- NA
o[
if (length(o) > 4)
<- rbind(out, o)
out
}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
#
<- function(data,
jk cores = parallel::detectCores()) {
# Calculate performance indicators repeatedly
# while excluding a single replication
<- parallel::makeCluster(cores)
cl ::clusterExport(cl, c("calc_indicators", "data"))
parallel<-
x_list ::pbsapply(
pbapplyX = seq_len(nrow(data)),
FUN = \(i) { calc_indicators(data[-i, ])[, -1] },
cl = cl,
simplify = FALSE
|>
) t() |>
as.data.frame()
::stopCluster(cl)
parallel
# Calculate standard errors
<- nrow(data)
n <- array(unlist(x_list), dim = c(dim(x_list[[1]][[1]]), n))
x <- array(t(apply(x, 1, rowMeans)), dim = dim(x))
y <- t((apply((x - y)^2, 1, rowMeans) * (n - 1))^.5)
se
# Set names
<- calc_indicators(data)
p <- data.frame(cond = p$con, se)
se colnames(se) <- colnames(p)
rownames(se) <- se$cond
return(se)
}
# Performance indicators
<- lapply(simdat, calc_indicators)
pind
# Monte Carlo errors
<- lapply(simdat, jk)
mce
# Restructure results
<- list()
pind_list for (m in names(pind)) {
<- list(
pind_list[[m]] err = subset(pind[[m]], select = grepl("^err_", names(pind[[m]])))
)
for (i in c("bias", "pbias", "rmse", "co")) {
# Estimate
<- subset(
pind_list[[m]][[i]]
pind[[m]], select = grepl(paste0("^", i, "_"), names(pind[[m]]))
)
# Error
paste0(i, "_err")]] <- subset(
pind_list[[m]][[
mce[[m]],select = grepl(paste0("^", i, "_"), names(mce[[m]]))
)
# Confidence interval
paste0(i, "_lb")]] <-
pind_list[[m]][[- qnorm(.975) * pind_list[[m]][[paste0(i, "_err")]]
pind_list[[m]][[i]] paste0(i, "_ub")]] <-
pind_list[[m]][[+ qnorm(.975) * pind_list[[m]][[paste0(i, "_err")]]
pind_list[[m]][[i]]
}
}<- pind_list
pind rm(m, i, pind_list)
saveRDS(pind, file = here::here("data", "simind.rds"))
saveRDS(mce, file = here::here("data", "simind_err.rds"))