| Title: | Simulation Infrastructure for Mediation Analysis |
|---|---|
| Description: | Provides standardized infrastructure for conducting Monte Carlo simulation studies in mediation analysis. Features include environment-aware execution (local vs HPC cluster), parallel processing with progress reporting, ground truth caching, automated result analysis, and publication-ready visualization. Designed to work seamlessly with the mediationverse ecosystem: 'medfit', 'probmed', 'RMediation', and 'medrobust'. |
| Authors: | Davood Tofighi [aut, cre] (ORCID: <https://orcid.org/0000-0001-8523-7776>) |
| Maintainer: | Davood Tofighi <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 0.6.0.9000 |
| Built: | 2026-10-25 06:25:13 UTC |
| Source: | https://github.com/Data-Wise/medsim |
Thin, documented amputer aligning mice::ampute semantics to the medsim DGM contract.
Returns the input frame with NAs inserted in target column(s); column names/order
are preserved (DGM data contract).
medsim_amputate( data, target, mechanism = c("MCAR", "MAR", "MNAR"), prop = 0.2, predictors = NULL, weights = NULL, type = c("RIGHT", "LEFT", "MID", "TAIL") )medsim_amputate( data, target, mechanism = c("MCAR", "MAR", "MNAR"), prop = 0.2, predictors = NULL, weights = NULL, type = c("RIGHT", "LEFT", "MID", "TAIL") )
data |
data.frame with columns |
target |
Character vector of column(s) to set |
mechanism |
One of |
prop |
Target overall missingness proportion in |
predictors |
Columns driving missingness for MAR/MNAR (defaults to all non- |
weights |
Optional named weights for the missingness logistic. |
type |
Amputation tail type, |
data with NAs inserted in target; names/order preserved.
set.seed(1) d <- data.frame(X = rnorm(300), M = rnorm(300), Y = rnorm(300)) d_mar <- medsim_amputate(d, target = "M", mechanism = "MAR", prop = 0.3) mean(is.na(d_mar$M))set.seed(1) d <- data.frame(X = rnorm(300), M = rnorm(300), Y = rnorm(300)) d_mar <- medsim_amputate(d, target = "M", mechanism = "MAR", prop = 0.3) mean(is.na(d_mar$M))
Computes comprehensive accuracy metrics for simulation results including bias, mean absolute error (MAE), root mean squared error (RMSE), and relative efficiency.
medsim_analyze(results, metrics = "all", by_scenario = TRUE)medsim_analyze(results, metrics = "all", by_scenario = TRUE)
results |
A medsim_results object from |
metrics |
Character vector: which metrics to compute. Options:
|
by_scenario |
Logical: Compute metrics separately for each scenario (default TRUE) |
For each estimated parameter with known ground truth:
Bias: Mean estimation error (estimate - truth)
MAE: Mean absolute error, robust to outliers
RMSE: Root mean squared error, penalizes large errors
Median AE: Median absolute error, very robust
Max AE: Maximum absolute error, identifies worst case
Relative Bias: Bias as percentage of true value
Bias close to 0: unbiased estimator
MAE/RMSE small: accurate estimator
RMSE > MAE: presence of outliers/large errors
Relative bias: standardized bias metric
A list with class "medsim_analysis" containing:
accuracy: data.frame with accuracy metrics
by_scenario: data.frame with metrics by scenario (if by_scenario=TRUE)
summary: Overall summary statistics
medsim_analyze_coverage(), medsim_compare_methods()
## Not run: # Run simulation results <- medsim_run(method, scenarios, config) # Analyze accuracy analysis <- medsim_analyze(results) # View overall metrics print(analysis$accuracy) # View by-scenario metrics print(analysis$by_scenario) # Custom metrics analysis <- medsim_analyze( results, metrics = c("bias", "mae", "rmse"), by_scenario = FALSE ) ## End(Not run)## Not run: # Run simulation results <- medsim_run(method, scenarios, config) # Analyze accuracy analysis <- medsim_analyze(results) # View overall metrics print(analysis$accuracy) # View by-scenario metrics print(analysis$by_scenario) # Custom metrics analysis <- medsim_analyze( results, metrics = c("bias", "mae", "rmse"), by_scenario = FALSE ) ## End(Not run)
Computes coverage rates for confidence intervals from simulation results.
Requires that the simulation method returns confidence interval bounds
(e.g., ci_lower, ci_upper).
medsim_analyze_coverage( results, ci_levels = c(0.9, 0.95, 0.99), ci_suffix = "_ci", by_scenario = TRUE, estimand = NULL )medsim_analyze_coverage( results, ci_levels = c(0.9, 0.95, 0.99), ci_suffix = "_ci", by_scenario = TRUE, estimand = NULL )
results |
A medsim_results object from |
ci_levels |
Numeric vector: nominal confidence levels to check (default: c(0.90, 0.95, 0.99)) |
ci_suffix |
Character: suffix for CI columns (default: "_ci") E.g., if parameter is "indirect", looks for "indirect_ci_lower" and "indirect_ci_upper" |
by_scenario |
Logical: Compute coverage separately for each scenario (default TRUE) |
estimand |
Optional |
Coverage rate is the proportion of confidence intervals that contain the true parameter value. For a 95% CI, expect ~95% coverage in large samples.
Coverage < nominal: CI too narrow (anti-conservative)
Coverage ~= nominal: CI has correct width
Coverage > nominal: CI too wide (conservative)
The function looks for columns named:
{parameter}_ci_lower and {parameter}_ci_upper (default)
Or custom suffix: {parameter}{ci_suffix}_lower and {parameter}{ci_suffix}_upper
A list with class "medsim_coverage" containing:
coverage: data.frame with coverage rates, plus n_failed/failure_rate
(method() failures, conditional on data_generator() having
succeeded) and n_dgp_failed/dgp_failure_rate (data_generator()
failures, unconditional – the two rates use different denominators
and do not sum, by design; see .medsim_failure_split())
by_scenario: data.frame with coverage by scenario (if by_scenario=TRUE)
summary: Overall summary statistics
medsim_analyze(), medsim_analyze_power()
## Not run: # Method that returns CIs my_method <- function(data, params) { # ... estimation list( indirect = estimate, indirect_ci_lower = ci[1], indirect_ci_upper = ci[2] ) } results <- medsim_run(my_method, scenarios, config) # Analyze coverage coverage <- medsim_analyze_coverage(results) print(coverage$coverage) #> parameter coverage_90 coverage_95 coverage_99 #> indirect 0.902 0.951 0.991 ## End(Not run)## Not run: # Method that returns CIs my_method <- function(data, params) { # ... estimation list( indirect = estimate, indirect_ci_lower = ci[1], indirect_ci_upper = ci[2] ) } results <- medsim_run(my_method, scenarios, config) # Analyze coverage coverage <- medsim_analyze_coverage(results) print(coverage$coverage) #> parameter coverage_90 coverage_95 coverage_99 #> indirect 0.902 0.951 0.991 ## End(Not run)
Per scenario, computes bias, empirical SE, mean model SE, and RMSE for a parameter, each with its Monte Carlo SE (Morris, White & Crowther 2019).
medsim_analyze_performance(results, param = "indirect")medsim_analyze_performance(results, param = "indirect")
results |
A |
param |
Parameter name (default "indirect"). |
Truth is taken as the first converged row's <param>_truth value
and assumed constant within scenario. Rows with NA estimates are excluded
from computations; if all estimates are NA, numeric columns return NA.
A data.frame, one row per scenario.
Computes empirical power (rejection rate under alternative hypothesis) from simulation results. Requires that the simulation method returns p-values.
medsim_analyze_power( results, alpha = 0.05, p_suffix = "_p", null_value = 0, by_scenario = TRUE )medsim_analyze_power( results, alpha = 0.05, p_suffix = "_p", null_value = 0, by_scenario = TRUE )
results |
A medsim_results object from |
alpha |
Numeric: significance level (default: 0.05) |
p_suffix |
Character: suffix for p-value columns (default: "_p") E.g., if parameter is "indirect", looks for "indirect_p" |
null_value |
Numeric: null hypothesis value (default: 0) |
by_scenario |
Logical: Compute power separately for each scenario (default TRUE) |
Power is the probability of correctly rejecting the null hypothesis when it is false (i.e., detecting a true effect).
Power = proportion of replications where p < alpha
Higher power = better ability to detect effects
Power depends on: effect size, sample size, alpha level
This function computes empirical power from simulations, not theoretical power. For scenarios where the null is true (e.g., indirect effect = 0), the rejection rate represents Type I error rate, not power.
A list with class "medsim_power" containing:
power: data.frame with power rates
by_scenario: data.frame with power by scenario (if by_scenario=TRUE)
summary: Overall summary statistics
medsim_analyze(), medsim_analyze_coverage()
## Not run: # Method that returns p-values my_method <- function(data, params) { # ... estimation and testing list( indirect = estimate, indirect_p = p_value ) } results <- medsim_run(my_method, scenarios, config) # Analyze power power <- medsim_analyze_power(results, alpha = 0.05) print(power$power) #> parameter power n_valid #> indirect 0.847 1000 ## End(Not run)## Not run: # Method that returns p-values my_method <- function(data, params) { # ... estimation and testing list( indirect = estimate, indirect_p = p_value ) } results <- medsim_run(my_method, scenarios, config) # Analyze power power <- medsim_analyze_power(results, alpha = 0.05) print(power$power) #> parameter power n_valid #> indirect 0.847 1000 ## End(Not run)
Standalone entry point for Gate A of the chunked-run integrity layer
(SPEC-medsim-chunked-run-gates-2026-07-31): run the same audit
medsim_combine_chunks() performs, on an already-combined
medsim_results object. Checks, per scenario:
Contiguity – replication ids form a gapless, duplicate-free
1..max run, with the same max in every scenario (self-validating:
needs no external replication count).
Collapse signature – every continuous estimate column has more than
collapse_threshold * n_ok distinct values (the 0.3.1 seed-collapse
produced ~17 distinct outcomes in 1000 while every chunk exited 0).
An all-failed cell (n_ok == 0) is reported as cell_failed, not
collapse.
Seed-space collisions – no two scenarios map onto overlapping
.medsim_det_seed() sequences (the ~1e6-bucket name hash can collide).
Partial chunks – for a combined object, any chunk whose provenance
is flagged partial (a task failed mid-run) is reported as a
partial_chunk violation, exactly as medsim_combine_chunks() reports
it. The surviving ids of a partial chunk can still be contiguous, so
checks 1-3 alone would pass it. A bare data.frame carries no provenance
and skips this check.
Frames without the schema-v2 stamp (produced before medsim 0.5.0, when
replication was chunk-LOCAL and collided across chunks) skip checks 1
and 3 with a warning – they cannot be audited for id integrity.
medsim_audit_results( results, on_violation = c("stop", "warn", "ignore"), nsim = NULL, collapse_threshold = 0.9, collapse_digits = 12L, collapse_min_cell = 30L, collapse_exclude = .medsim_collapse_exclude_default )medsim_audit_results( results, on_violation = c("stop", "warn", "ignore"), nsim = NULL, collapse_threshold = 0.9, collapse_digits = 12L, collapse_min_cell = 30L, collapse_exclude = .medsim_collapse_exclude_default )
results |
A |
on_violation |
One of |
nsim |
Integer or |
collapse_threshold, collapse_digits, collapse_min_cell
|
Tuning for the
collapse-signature audit: per scenario and per continuous estimate column,
require |
collapse_exclude |
Character: DISCRETE method-contract fields excluded
from the collapse audit BY NAME (default: |
Invisibly, the list of violations found (empty if clean). Signalling
follows on_violation, exactly as in medsim_combine_chunks().
medsim_combine_chunks(), medsim_run_chunk()
Deletes cache files from a directory. Can delete all files or files matching a pattern.
medsim_cache_clear( cache_dir, pattern = "*.rds", max_age = NULL, confirm = TRUE )medsim_cache_clear( cache_dir, pattern = "*.rds", max_age = NULL, confirm = TRUE )
cache_dir |
Character: Directory containing cache files |
pattern |
Character: Optional pattern to match files (e.g., "truth_*.rds") |
max_age |
Numeric: Delete only files older than this many days. If NULL, deletes all matching files. |
confirm |
Logical: Ask for confirmation before deleting (default TRUE) |
Invisibly returns number of files deleted
## Not run: # Clear all cache medsim_cache_clear("cache") # Clear old cache only medsim_cache_clear("cache", max_age = 30) # Clear specific pattern medsim_cache_clear("cache", pattern = "truth_*.rds") # Without confirmation medsim_cache_clear("cache", confirm = FALSE) ## End(Not run)## Not run: # Clear all cache medsim_cache_clear("cache") # Clear old cache only medsim_cache_clear("cache", max_age = 30) # Clear specific pattern medsim_cache_clear("cache", pattern = "truth_*.rds") # Without confirmation medsim_cache_clear("cache", confirm = FALSE) ## End(Not run)
Checks if a cache file exists and is readable.
medsim_cache_exists(file)medsim_cache_exists(file)
file |
Character: Path to cache file |
Logical: TRUE if cache exists and is readable, FALSE otherwise
## Not run: if (medsim_cache_exists("cache/truth.rds")) { truth <- medsim_cache_load("cache/truth.rds") } else { truth <- list(indirect = 0.09) # Compute your truth values medsim_cache_save(truth, "cache/truth.rds") } ## End(Not run)## Not run: if (medsim_cache_exists("cache/truth.rds")) { truth <- medsim_cache_load("cache/truth.rds") } else { truth <- list(indirect = 0.09) # Compute your truth values medsim_cache_save(truth, "cache/truth.rds") } ## End(Not run)
Retrieves metadata about a cached object without loading the full object.
medsim_cache_info(file)medsim_cache_info(file)
file |
Character: Path to cache file |
List with cache metadata:
exists: Logical - file exists
size_mb: Numeric - file size in MB
modified: POSIXct - last modified time
age_days: Numeric - age in days
timestamp: POSIXct - creation time (if available)
r_version: Character - R version used (if available)
# Save a cached object to a temp file, then inspect it cache_file <- tempfile(fileext = ".rds") medsim_cache_save(list(value = 42), cache_file) info <- medsim_cache_info(cache_file) print(info) unlink(cache_file)# Save a cached object to a temp file, then inspect it cache_file <- tempfile(fileext = ".rds") medsim_cache_save(list(value = 42), cache_file) info <- medsim_cache_info(cache_file) print(info) unlink(cache_file)
Creates a cache directory with proper structure and a README file explaining the cache system.
medsim_cache_init(cache_dir, create_readme = TRUE)medsim_cache_init(cache_dir, create_readme = TRUE)
cache_dir |
Character: Path to cache directory |
create_readme |
Logical: Create README.md explaining cache (default TRUE) |
Invisibly returns cache directory path
cache_dir <- file.path(tempdir(), "medsim_cache_example") medsim_cache_init(cache_dir) unlink(cache_dir, recursive = TRUE)cache_dir <- file.path(tempdir(), "medsim_cache_example") medsim_cache_init(cache_dir) unlink(cache_dir, recursive = TRUE)
Lists all cache files in a directory with their metadata.
medsim_cache_list(cache_dir, pattern = "*.rds")medsim_cache_list(cache_dir, pattern = "*.rds")
cache_dir |
Character: Directory containing cache files |
pattern |
Character: Pattern to match files (default "*.rds") |
data.frame with columns:
file: File name
path: Full path
size_mb: Size in MB
age_days: Age in days
modified: Last modified time
# Create a small cache in a temp directory, then list it cache_dir <- file.path(tempdir(), "medsim_cache_list_example") dir.create(cache_dir, showWarnings = FALSE) medsim_cache_save(list(a = 1), file.path(cache_dir, "result1.rds")) medsim_cache_save(list(b = 2), file.path(cache_dir, "result2.rds")) cache_list <- medsim_cache_list(cache_dir) print(cache_list) unlink(cache_dir, recursive = TRUE)# Create a small cache in a temp directory, then list it cache_dir <- file.path(tempdir(), "medsim_cache_list_example") dir.create(cache_dir, showWarnings = FALSE) medsim_cache_save(list(a = 1), file.path(cache_dir, "result1.rds")) medsim_cache_save(list(b = 2), file.path(cache_dir, "result2.rds")) cache_list <- medsim_cache_list(cache_dir) print(cache_list) unlink(cache_dir, recursive = TRUE)
Loads a previously cached R object. Returns NULL if cache file doesn't exist.
medsim_cache_load(file, check_expiry = FALSE, max_age = 30, verbose = TRUE)medsim_cache_load(file, check_expiry = FALSE, max_age = 30, verbose = TRUE)
file |
Character: Path to cache file |
check_expiry |
Logical: Check if cache has expired (default FALSE) |
max_age |
Numeric: Maximum age of cache in days (default 30). Only
used if |
verbose |
Logical: Print messages about cache status (default TRUE) |
When loading, the cache is checked for:
File exists
File is readable
Age is within max_age (if check_expiry = TRUE)
If cache doesn't exist:
Returns NULL (no error)
Calling code should compute and cache the result
Example pattern:
result <- medsim_cache_load("cache/truth.rds")
if (is.null(result)) {
result <- expensive_computation()
medsim_cache_save(result, "cache/truth.rds")
}
Note: ground-truth caches written by medsim_run() use a fingerprinted
wrapper format list(truth = , fingerprint = ) (so a scenario content
change invalidates the cache) – unwrap $truth after loading, as in the
example below.
Cached object if file exists and is valid, NULL otherwise
medsim_cache_save(), medsim_cache_exists()
## Not run: # Load a ground-truth cache written by medsim_run(). Since the # fingerprinted format, the cached object is list(truth=, fingerprint=), # not the bare truth value -- unwrap it (legacy bare caches load as-is): cached <- medsim_cache_load("cache/truth_scenario1.rds") truth <- if (is.list(cached) && !is.null(cached$fingerprint)) cached$truth else cached if (is.null(truth)) { # Cache miss - compute and cache truth <- compute_expensive_truth() medsim_cache_save(truth, "cache/truth_scenario1.rds") } # Load with expiry check result <- medsim_cache_load( "cache/old_result.rds", check_expiry = TRUE, max_age = 7 # 7 days ) ## End(Not run)## Not run: # Load a ground-truth cache written by medsim_run(). Since the # fingerprinted format, the cached object is list(truth=, fingerprint=), # not the bare truth value -- unwrap it (legacy bare caches load as-is): cached <- medsim_cache_load("cache/truth_scenario1.rds") truth <- if (is.list(cached) && !is.null(cached$fingerprint)) cached$truth else cached if (is.null(truth)) { # Cache miss - compute and cache truth <- compute_expensive_truth() medsim_cache_save(truth, "cache/truth_scenario1.rds") } # Load with expiry check result <- medsim_cache_load( "cache/old_result.rds", check_expiry = TRUE, max_age = 7 # 7 days ) ## End(Not run)
Saves an R object to a cache file for later retrieval. Uses RDS format for efficient storage and compression.
medsim_cache_save(object, file, compress = TRUE, overwrite = TRUE)medsim_cache_save(object, file, compress = TRUE, overwrite = TRUE)
object |
Any R object to cache |
file |
Character: Path to cache file. Parent directory will be created if it doesn't exist. |
compress |
Logical: Compress the cached object (default TRUE). Can significantly reduce file size for large objects. |
overwrite |
Logical: Overwrite existing cache file (default TRUE) |
Cache files are stored in RDS format with metadata:
Original object
Timestamp of creation
R version used
Package version (if available)
Recommended cache directory structure:
cache/ |-- truth_scenario_1.rds |-- truth_scenario_2.rds +-- ...
Invisibly returns the file path
medsim_cache_load(), medsim_cache_clear()
## Not run: # Cache simulation truth truth <- compute_expensive_truth() medsim_cache_save(truth, "cache/truth_scenario1.rds") # Cache with metadata result <- list(value = 42, timestamp = Sys.time()) medsim_cache_save(result, "cache/my_result.rds") ## End(Not run)## Not run: # Cache simulation truth truth <- compute_expensive_truth() medsim_cache_save(truth, "cache/truth_scenario1.rds") # Cache with metadata result <- list(value = 42, timestamp = Sys.time()) medsim_cache_save(result, "cache/my_result.rds") ## End(Not run)
Checks a list of results from parallel execution for errors and provides a summary.
medsim_check_results(results, stop_on_error = FALSE)medsim_check_results(results, stop_on_error = FALSE)
results |
List of results from |
stop_on_error |
Logical: If TRUE, stops with error if any task failed |
Invisibly returns indices of failed tasks
## Not run: results <- medsim_run_parallel(...) medsim_check_results(results) ## End(Not run)## Not run: results <- medsim_run_parallel(...) medsim_check_results(results) ## End(Not run)
Reads all chunk_<k>.rds files from output_dir, merges their $results
and $truth data frames, and returns a single medsim_results object
that is equivalent to running all replications in one call.
medsim_combine_chunks( output_dir, pattern = "chunk_*.rds", expected_chunks = NULL, on_violation = c("stop", "warn", "ignore"), nsim = NULL, collapse_threshold = 0.9, collapse_digits = 12L, collapse_min_cell = 30L, collapse_exclude = .medsim_collapse_exclude_default, pilot_reference = NULL, pilot_tol = 1e-09, verbose = TRUE )medsim_combine_chunks( output_dir, pattern = "chunk_*.rds", expected_chunks = NULL, on_violation = c("stop", "warn", "ignore"), nsim = NULL, collapse_threshold = 0.9, collapse_digits = 12L, collapse_min_cell = 30L, collapse_exclude = .medsim_collapse_exclude_default, pilot_reference = NULL, pilot_tol = 1e-09, verbose = TRUE )
output_dir |
Character: directory containing |
pattern |
Character: glob pattern for chunk files.
Default |
expected_chunks |
Integer or |
on_violation |
One of |
nsim |
Integer or |
collapse_threshold, collapse_digits, collapse_min_cell
|
Tuning for the
collapse-signature audit: per scenario and per continuous estimate column,
require |
collapse_exclude |
Character: DISCRETE method-contract fields excluded
from the collapse audit BY NAME (default: |
pilot_reference |
Character path (or |
pilot_tol |
Numeric: absolute per-value tolerance for the pilot
comparison (default |
verbose |
Logical: print file counts. |
A medsim_results object with combined $results and $truth.
medsim_run_chunk(), medsim_audit_results(),
medsim_analyze_coverage()
# After all SLURM chunks complete: # combined <- medsim_combine_chunks("simulation_results/") # coverage <- medsim_analyze_coverage(combined)# After all SLURM chunks complete: # combined <- medsim_combine_chunks("simulation_results/") # coverage <- medsim_analyze_coverage(combined)
Creates a comparison table showing differences between test, local, and cluster configurations.
medsim_compare_configs()medsim_compare_configs()
A data.frame with configuration comparison
comparison <- medsim_compare_configs() print(comparison)comparison <- medsim_compare_configs() print(comparison)
Compares performance of multiple estimation methods across scenarios. Combines results from multiple simulation runs and computes relative performance metrics.
medsim_compare_methods(..., metrics = "all")medsim_compare_methods(..., metrics = "all")
... |
Named medsim_results objects (e.g., method1 = results1, method2 = results2) |
metrics |
Character vector: which comparison metrics to compute. Options: "accuracy", "timing", "coverage", "all" (default) |
Accuracy: MAE, RMSE, bias across methods
Timing: Mean/median execution time
Coverage: CI coverage rates (if CIs available)
Relative Efficiency: Ratio of MSE (lower = better)
All results must:
Use the same scenarios (matched by name)
Use the same configuration (n_replications, seed)
Have ground truth available (for accuracy comparison)
A list with class "medsim_comparison" containing:
accuracy_comparison: data.frame comparing accuracy metrics
timing_comparison: data.frame comparing timing
coverage_comparison: data.frame comparing coverage (if applicable)
summary: Overall comparison summary
medsim_analyze(), medsim_analyze_coverage()
## Not run: # Run multiple methods results_proposed <- medsim_run(method_proposed, scenarios, config) results_delta <- medsim_run(method_delta, scenarios, config) results_boot <- medsim_run(method_boot, scenarios, config) # Compare comparison <- medsim_compare_methods( proposed = results_proposed, delta = results_delta, bootstrap = results_boot ) # View accuracy comparison print(comparison$accuracy_comparison) # View timing comparison print(comparison$timing_comparison) ## End(Not run)## Not run: # Run multiple methods results_proposed <- medsim_run(method_proposed, scenarios, config) results_delta <- medsim_run(method_delta, scenarios, config) results_boot <- medsim_run(method_boot, scenarios, config) # Compare comparison <- medsim_compare_methods( proposed = results_proposed, delta = results_delta, bootstrap = results_boot ) # View accuracy comparison print(comparison$accuracy_comparison) # View timing comparison print(comparison$timing_comparison) ## End(Not run)
Creates a configuration object for simulation studies. Automatically detects whether running on local machine or HPC cluster and adjusts parameters accordingly.
medsim_config( mode = "auto", n_replications = NULL, n_cores = NULL, scenarios = NULL, output_dir = NULL, seed = 12345, chunk_id = NULL, n_chunks = NULL, array_size = NULL, seed_stream = NULL, partition = NULL, walltime = NULL, mem_per_cpu = NULL, r_module = NULL, ..., n = NULL, n_truth = NULL )medsim_config( mode = "auto", n_replications = NULL, n_cores = NULL, scenarios = NULL, output_dir = NULL, seed = 12345, chunk_id = NULL, n_chunks = NULL, array_size = NULL, seed_stream = NULL, partition = NULL, walltime = NULL, mem_per_cpu = NULL, r_module = NULL, ..., n = NULL, n_truth = NULL )
mode |
Character: "auto", "test", "local", or "cluster"
|
n_replications |
Integer: Number of Monte Carlo replications. If NULL, uses mode defaults (test=20, local=100, cluster=1000) |
n_cores |
Integer: Number of CPU cores for parallel processing. If NULL, auto-detects (all cores - 2 on local, SLURM_CPUS_PER_TASK on cluster) |
scenarios |
Character: "all" or "test". Use "test" for single challenging scenario during development |
output_dir |
Character: Directory for saving results |
seed |
Integer: Base seed, stored on the config and printed for
provenance. Consulted directly by |
chunk_id |
Integer: SLURM array task index (1-based) for this chunk.
Auto-detected from |
n_chunks |
Integer: Total number of chunks (SLURM array size).
Used by |
array_size |
Integer: Alias for |
seed_stream |
Integer: Deprecated. This knob is stored on the
config but has never been consumed by any medsim function – passing a
non- |
partition |
Character: SLURM partition (queue) name. Defaults to
|
walltime |
Character: SLURM wall-time limit (HH:MM:SS). Defaults to
|
mem_per_cpu |
Character: SLURM memory per CPU. Defaults to |
r_module |
Character: Environment module string for |
... |
Additional custom parameters |
n |
Positive whole number or |
n_truth |
Positive whole number or |
| Mode | Replications | Cores | Runtime | Use Case |
| test | 20 | 4 | ~30s | Quick validation |
| local | 100 | auto | ~15m | Development |
| cluster | 1000 | SLURM | hours | Production |
When mode = "auto", checks for:
SLURM_JOB_ID (SLURM scheduler)
PBS_JOBID (PBS/Torque scheduler)
LSB_JOBID (LSF scheduler)
If any are found, uses "cluster" mode. Otherwise, uses "local" mode.
You can add custom parameters via ...:
config <- medsim_config( mode = "local", n_bootstrap = 5000, alpha = 0.05, custom_param = "value" )
A list with simulation configuration parameters
medsim_detect_environment(), print.medsim_config()
# Auto-detect environment config <- medsim_config(mode = "auto") # Explicit test mode for quick validation config_test <- medsim_config(mode = "test") # Local mode with custom replications config_local <- medsim_config( mode = "local", n_replications = 500 ) # Cluster mode (auto-detects SLURM cores) config_cluster <- medsim_config(mode = "cluster") # Sample size of each simulated data set (default 200) config_n100 <- medsim_config(mode = "local", n = 100) # Custom parameters config_custom <- medsim_config( mode = "local", n_bootstrap = 1000, ci_level = 0.95 ) # Print configuration print(config)# Auto-detect environment config <- medsim_config(mode = "auto") # Explicit test mode for quick validation config_test <- medsim_config(mode = "test") # Local mode with custom replications config_local <- medsim_config( mode = "local", n_replications = 500 ) # Cluster mode (auto-detects SLURM cores) config_cluster <- medsim_config(mode = "cluster") # Sample size of each simulated data set (default 200) config_n100 <- medsim_config(mode = "local", n = 100) # Custom parameters config_custom <- medsim_config( mode = "local", n_bootstrap = 1000, ci_level = 0.95 ) # Print configuration print(config)
Detects the number of CPU cores available for parallel processing. On HPC clusters, uses scheduler-allocated cores. On local machines, uses all cores minus 2 to avoid overwhelming the system.
medsim_detect_cores()medsim_detect_cores()
On clusters, checks for:
SLURM_CPUS_PER_TASK (SLURM)
PBS_NUM_PPN (PBS/Torque)
LSB_DJOB_NUMPROC (LSF)
On local machines, uses parallel::detectCores() - 2.
Integer: Number of cores to use
n_cores <- medsim_detect_cores() message(sprintf("Using %d cores", n_cores))n_cores <- medsim_detect_cores() message(sprintf("Using %d cores", n_cores))
Detects whether code is running on local machine or HPC cluster by checking for scheduler environment variables.
medsim_detect_environment()medsim_detect_environment()
Checks for the following environment variables:
SLURM_JOB_ID (SLURM scheduler)
PBS_JOBID (PBS/Torque scheduler)
LSB_JOBID (LSF scheduler)
If any are found, returns "cluster". Otherwise, returns "local".
Character: "local" or "cluster"
env <- medsim_detect_environment() if (env == "cluster") { message("Running on HPC cluster") } else { message("Running on local machine") }env <- medsim_detect_environment() if (env == "cluster") { message("Running on HPC cluster") } else { message("Running on local machine") }
Attach a first-class estimand descriptor to a medsim_scenario(). Every
downstream command (medsim_run, medsim_analyze_coverage,
medsim_tables, medsim_figures, medsim_workflow) dispatches on
estimand$kind rather than assuming a scalar point-estimand.
Five kinds are recognised:
| kind | what method() returns |
coverage notion |
"point" |
{p}, {p}_ci_lower/_upper, {p}_p |
truth in CI |
"interval" |
{p}_lower/_upper, {p}_im_lower/_im_upper, feasible, falsified |
truth in [lower, upper]; IM-CI coverage |
"probabilistic" |
pmed, pmed_ci_lower/_upper, pmed_p, branch_switch |
truth in CI (truth from potential outcomes) |
"variance_share" |
{p}, {p}_ci_lower/_upper |
truth in Wald CI (bounded scalar share in [0,1]) |
"numeric" |
est_error, abs_error, elapsed_sec |
none |
"variance_share" is a bounded scalar point estimand with a standard Wald
CI — the Sobol / functional-ANOVA proportion mediated
(see
medsim_scenario_sobol(), medsim_method_sobol()). Coverage uses the
generic truth in CI path; the dedicated kind exists for clear labeling and
so medsim_validate_scenario() checks the causal-notation A/M/Y/C columns
the Sobol estimator requires (rather than the legacy X/M/Y of "point").
Scenarios with estimand = NULL (the default in medsim_scenario()) are
treated as kind = "point" throughout the package – full back-compatibility
with all existing code.
medsim_estimand( kind = c("point", "interval", "probabilistic", "variance_share", "numeric"), params = character(), truth = NULL, ci = c("standard", "imbens_manski", "mbco", "none"), extra = character() )medsim_estimand( kind = c("point", "interval", "probabilistic", "variance_share", "numeric"), params = character(), truth = NULL, ci = c("standard", "imbens_manski", "mbco", "none"), extra = character() )
kind |
Estimand kind: |
params |
Character vector of estimand parameter names – e.g.
|
truth |
Optional truth-extractor |
ci |
CI method: |
extra |
Character vector of additional result columns beyond |
An object of class c("medsim_estimand", "list").
medsim_scenario() for attaching an estimand to a scenario;
medsim_analyze_coverage() for the keystone dispatch on kind.
# Point estimand (the historic default -- back-compat) medsim_estimand("point") # Interval (partial-ID bounds) medsim_estimand("interval", params = c("NDE", "NIE"), ci = "imbens_manski", extra = c("feasible", "falsified")) # Probabilistic (P_med) medsim_estimand("probabilistic", params = "pmed", ci = "mbco", extra = "branch_switch") # Numerical accuracy (no coverage/power) medsim_estimand("numeric", params = "est_error", ci = "none", extra = "elapsed_sec")# Point estimand (the historic default -- back-compat) medsim_estimand("point") # Interval (partial-ID bounds) medsim_estimand("interval", params = c("NDE", "NIE"), ci = "imbens_manski", extra = c("feasible", "falsified")) # Probabilistic (P_med) medsim_estimand("probabilistic", params = "pmed", ci = "mbco", extra = "branch_switch") # Numerical accuracy (no coverage/power) medsim_estimand("numeric", params = "est_error", ci = "none", extra = "elapsed_sec")
Estimates the expected speedup from parallel processing based on number of cores and task characteristics.
medsim_estimate_speedup(n_tasks, n_cores, overhead = 0.1)medsim_estimate_speedup(n_tasks, n_cores, overhead = 0.1)
n_tasks |
Integer: Number of tasks |
n_cores |
Integer: Number of cores |
overhead |
Numeric: Overhead per core (default 0.1 = 10%) |
Speedup is estimated using Amdahl's law with overhead:
Speedup = n_cores / (1 + overhead * n_cores)
Efficiency = Speedup / n_cores
Named list with:
speedup: Expected speedup factor
efficiency: Parallel efficiency (0-1)
recommendation: Character recommendation
# For 100 tasks on 10 cores est <- medsim_estimate_speedup(100, 10) cat(sprintf("Expected speedup: %.1fx\n", est$speedup)) cat(sprintf("Efficiency: %.0f%%\n", est$efficiency * 100))# For 100 tasks on 10 cores est <- medsim_estimate_speedup(100, 10) cat(sprintf("Expected speedup: %.1fx\n", est$speedup)) cat(sprintf("Efficiency: %.0f%%\n", est$efficiency * 100))
Convenience wrapper that produces the standard set of medsim plots
(coverage, error boxplot, timing) and writes them to output_dir.
medsim_figures( results, output_dir = "figures", format = "pdf", width = 8, height = 6, ... )medsim_figures( results, output_dir = "figures", format = "pdf", width = 8, height = 6, ... )
results |
A |
output_dir |
Directory to write figures to. Created if it doesn't exist. |
format |
File extension passed to |
width, height
|
Figure dimensions in inches. |
... |
Additional arguments passed to |
Invisibly, a named character vector of file paths written.
Determines the optimal number of cores to use for parallel processing based on system capabilities and task requirements.
medsim_get_optimal_cores(n_tasks, reserve = 2, max_cores = NULL)medsim_get_optimal_cores(n_tasks, reserve = 2, max_cores = NULL)
n_tasks |
Integer: Number of tasks to process. If fewer tasks than
cores, returns |
reserve |
Integer: Number of cores to reserve for system (default 2). Only applied on local machines, not clusters. |
max_cores |
Integer: Maximum number of cores to use. If NULL, uses all available (minus reserve). |
Integer: Optimal number of cores to use
# For 100 tasks medsim_get_optimal_cores(n_tasks = 100) # For few tasks (returns n_tasks) medsim_get_optimal_cores(n_tasks = 3) # With maximum limit medsim_get_optimal_cores(n_tasks = 1000, max_cores = 8)# For 100 tasks medsim_get_optimal_cores(n_tasks = 100) # For few tasks (returns n_tasks) medsim_get_optimal_cores(n_tasks = 3) # With maximum limit medsim_get_optimal_cores(n_tasks = 1000, max_cores = 8)
Returns partial-identification bounds (NDE/NIE) following the interval-kind
result contract expected by medsim_run(), so that
medsim_analyze_coverage() dispatches the Imbens-Manski / partial-ID
interval coverage branch.
Current implementation: OLS decomposition + inflated interval (synthetic
bounds for unit-testing the harness). A future version will call
medrobust::bound_ne() once its accessor API stabilises on CRAN.
Result columns returned (flat named list):
NDE_lower, NDE_upper – bounds on natural direct effect
NDE_im_lower, NDE_im_upper – Imbens-Manski CI for NDE
NIE_lower, NIE_upper – bounds on natural indirect effect
NIE_im_lower, NIE_im_upper – Imbens-Manski CI for NIE
feasible – logical: bounds are non-empty / consistent
falsified – logical: data rejects identifying assumptions
medsim_method_bounds( data, params, misclass_type = c("mediator", "exposure"), alpha = 0.05 )medsim_method_bounds( data, params, misclass_type = c("mediator", "exposure"), alpha = 0.05 )
data |
A data.frame produced by the scenario's |
params |
Named list of scenario parameters. Must include |
misclass_type |
|
alpha |
Nominal level; |
Named list with elements NDE_lower, NDE_upper,
NDE_im_lower, NDE_im_upper, NIE_lower, NIE_upper,
NIE_im_lower, NIE_im_upper, feasible, falsified.
medsim_scenario_dm() for the companion scenario factory;
medsim_analyze_coverage() for interval-kind coverage analysis.
set.seed(42) d <- data.frame(A = rbinom(200, 1, 0.5), M = rnorm(200), Y = rnorm(200)) result <- medsim_method_bounds(d, list(NDE = 0.2, NIE = 0.3)) str(result)set.seed(42) d <- data.frame(A = rbinom(200, 1, 0.5), M = rnorm(200), Y = rnorm(200)) result <- medsim_method_bounds(d, list(NDE = 0.2, NIE = 0.3)) str(result)
Adapter that calls an external ward_residual() estimator and returns the
flat field contract for two parameters:
| Field | Description |
pmed, pmed_ci_lower/_upper |
interventional proportion mediated IIE/OE |
w, w_ci_lower/_upper |
gauge residual R/OE
|
The estimator is not vendored (medsim keeps estimators external). Pass it
via estimator, or default to an object named ward_residual on the search
path / probmed::ward_residual. The estimator must return an S7
GaugePmedResult (read via @p_med/@p_med_ci/@W/@W_ci) or a plain list with
fields p_med, p_med_ci (length 2), W, W_ci (length 2).
medsim_method_gauge( data, params = list(), estimator = get0("ward_residual"), covars = "C", K = 5L, ci_level = 0.95, seed = 1L, reps = 1L, se_method = "analytic", B = 200L, fieller = FALSE )medsim_method_gauge( data, params = list(), estimator = get0("ward_residual"), covars = "C", K = 5L, ci_level = 0.95, seed = 1L, reps = 1L, se_method = "analytic", B = 200L, fieller = FALSE )
data |
data.frame with columns |
params |
Named list from |
estimator |
Function; default |
covars |
Character covariate names (default |
K |
Integer cross-fitting folds (default 5). |
ci_level |
Numeric CI level (default 0.95). |
seed |
Integer seed (default 1). |
reps |
Integer repeated cross-fitting fold draws (default 1). |
se_method |
|
B |
Integer bootstrap resamples when |
fieller |
Logical: compute the estimator's Fieller bounds (default
|
Named list pmed, pmed_ci_lower, pmed_ci_upper, w, w_ci_lower, w_ci_upper.
medsim_scenario_gauge(), medsim_analyze_coverage()
Inverse-probability-weighting comparator. Positioned as a robustness appendix, not a co-equal main-results arm (decided 2026-06-11).
medsim_method_ipw(model, ...)medsim_method_ipw(model, ...)
model |
Mediation model spec. |
... |
Passed to the underlying weighting / fitting routines. |
A function(data, params) returning the shared method contract (branch_switch = NA).
set.seed(1) d <- data.frame(X = rnorm(200)) d$M <- 0.39 * d$X + rnorm(200) d$Y <- 0.39 * d$M + 0.14 * d$X + rnorm(200) d <- medsim_amputate(d, target = "M", mechanism = "MAR", prop = 0.25) out <- medsim_method_ipw(model = NULL)(d, list()) out[c("indirect", "indirect_ci_lower", "indirect_ci_upper")]set.seed(1) d <- data.frame(X = rnorm(200)) d$M <- 0.39 * d$X + rnorm(200) d$Y <- 0.39 * d$M + 0.14 * d$X + rnorm(200) d <- medsim_amputate(d, target = "M", mechanism = "MAR", prop = 0.25) out <- medsim_method_ipw(model = NULL)(d, list()) out[c("indirect", "indirect_ci_lower", "indirect_ci_upper")]
Multiple imputation (via mice) + D4-stacked MBCO likelihood-ratio
inference for the union null H0: ab = 0. Ports the validated
Missing Effect/code/prototype-d4-mbco.R (phase-2; the D4 statistic matches
mitml::testModels(method = "D4"), and so does the p-value with
df = "chan_meng"): per-imputation MBCO LRTs are pooled with
the Reiter/Chan-Meng D4 rule (Grund et al. 2021) into an F reference.
medsim_method_mbco_mi( model, m = 20L, ariv = c("fixed", "own"), df = c("chan_meng", "li"), ... )medsim_method_mbco_mi( model, m = 20L, ariv = c("fixed", "own"), df = c("chan_meng", "li"), ... )
model |
Mediation model spec (accepted for API symmetry; the estimator
reads the |
m |
Integer number of imputations. |
ariv |
Which relative-increase-in-variance (ARIV) estimate feeds the
|
df |
Denominator degrees of freedom for the D4 |
... |
Reserved for future MI / MBCO options. |
The point estimate is the Rubin-pooled a*b; the CI is a Monte-Carlo product
interval on the pooled paths; the p-value is the D4-pooled MBCO test. With no
missingness (or without mice) it degrades to the complete-case MBCO LRT with
a chi-square reference – never to a Sobel normal approximation for the test.
A function(data, params) returning the shared method contract (see file header).
if (requireNamespace("mice", quietly = TRUE)) { set.seed(1) d <- data.frame(X = rnorm(200)) d$M <- 0.39 * d$X + rnorm(200) d$Y <- 0.39 * d$M + 0.14 * d$X + rnorm(200) d <- medsim_amputate(d, target = "M", mechanism = "MAR", prop = 0.25) fit <- medsim_method_mbco_mi(model = NULL, m = 10, ariv = "fixed", df = "chan_meng") out <- fit(d, list()) out[c("indirect", "indirect_p", "indirect_p_fixed", "r4_fixed")] }if (requireNamespace("mice", quietly = TRUE)) { set.seed(1) d <- data.frame(X = rnorm(200)) d$M <- 0.39 * d$X + rnorm(200) d$Y <- 0.39 * d$M + 0.14 * d$X + rnorm(200) d <- medsim_amputate(d, target = "M", mechanism = "MAR", prop = 0.25) fit <- medsim_method_mbco_mi(model = NULL, m = 10, ariv = "fixed", df = "chan_meng") out <- fit(d, list()) out[c("indirect", "indirect_p", "indirect_p_fixed", "r4_fixed")] }
Multiple imputation (via mice) + Monte-Carlo confidence interval for the
indirect effect (via RMediation::medci() when available, else a base-R
product-of-normals draw on the Rubin-pooled paths).
medsim_method_mc_ci(model, m = 20L, ...)medsim_method_mc_ci(model, m = 20L, ...)
model |
Mediation model spec (accepted for API symmetry; the estimator
reads the |
m |
Integer number of imputations. |
... |
Reserved for future MI / MBCO options. |
A function(data, params) returning the shared method contract (branch_switch = NA).
if (requireNamespace("mice", quietly = TRUE)) { set.seed(1) d <- data.frame(X = rnorm(200)) d$M <- 0.39 * d$X + rnorm(200) d$Y <- 0.39 * d$M + 0.14 * d$X + rnorm(200) d <- medsim_amputate(d, target = "M", mechanism = "MAR", prop = 0.25) out <- medsim_method_mc_ci(model = NULL, m = 10)(d, list()) out[c("indirect", "indirect_ci_lower", "indirect_ci_upper")] }if (requireNamespace("mice", quietly = TRUE)) { set.seed(1) d <- data.frame(X = rnorm(200)) d$M <- 0.39 * d$X + rnorm(200) d$Y <- 0.39 * d$M + 0.14 * d$X + rnorm(200) d <- medsim_amputate(d, target = "M", mechanism = "MAR", prop = 0.25) out <- medsim_method_mc_ci(model = NULL, m = 10)(d, list()) out[c("indirect", "indirect_ci_lower", "indirect_ci_upper")] }
Estimates P_med from observed data using OLS point estimation and a
two-branch MBCO (maximally balanced constrained optimisation) confidence
interval. Returns the 6-field contract expected by medsim_run():
| Field | Description |
pmed |
Point estimate of P_med |
pmed_ci_lower |
Lower bound of MBCO CI |
pmed_ci_upper |
Upper bound of MBCO CI |
pmed_p |
Two-sided p-value (H0: P_med = 0.5) |
branch_switch |
1 if the MBCO union-null LRT switched branches |
converged |
1 if optimisation converged |
Algorithm:
Fit linear SEM: M ~ A and Y ~ A + M.
Estimate path coefficients alpha (A->M) and beta (M->Y).
Estimate P_med via a parametric bootstrap PO draw under the estimated SEM.
Build a delta-method normal CI and clamp to [0, 1].
The "branch_switch" flag records whether the union-null test selected the constrained branch (alphabeta <= 0) over the unconstrained branch (alphabeta > 0).
medsim_method_pmed_mbco(data, params, alpha = 0.05, n_boot = 2000L)medsim_method_pmed_mbco(data, params, alpha = 0.05, n_boot = 2000L)
data |
A |
params |
Named list from |
alpha |
Significance level for CI (default 0.05). |
n_boot |
Integer: parametric bootstrap size for P_med point estimate (default 2000). |
A named list with fields pmed, pmed_ci_lower, pmed_ci_upper,
pmed_p, branch_switch, converged.
medsim_scenario_pmed(), medsim_analyze_coverage()
Adapter that calls an external sobol_pmed() estimator and returns the flat
field contract keyed by the pmed_sobol token:
| Field | Description |
pmed_sobol |
Point estimate of the Sobol share |
pmed_sobol_ci_lower |
Lower Wald CI bound |
pmed_sobol_ci_upper |
Upper Wald CI bound |
Following medsim convention, the actual estimator is not vendored into
the package (medsim keeps estimators external — it Suggests companion
packages, never hard-depends on them). The prototype lives at
pmed-modern-sobol/03-sobol-pmed/sims/sobol_pmed.R; source() it (or load
the future probmed export) before running, or pass it via estimator.
The external sobol_pmed(d, ...) is expected to return a list with
P_med_sobol (scalar) and ci (length-2 numeric c(lower, upper)).
medsim_method_sobol( data, params = list(), estimator = get0("sobol_pmed"), covars = "C", K = 5L, level = 0.95, seed = 1L )medsim_method_sobol( data, params = list(), estimator = get0("sobol_pmed"), covars = "C", K = 5L, level = 0.95, seed = 1L )
data |
A |
params |
Named list from |
estimator |
Function implementing the Sobol estimator. Defaults to
|
covars |
Character vector of covariate column names passed to the
estimator (default |
K |
Integer cross-fitting folds passed to the estimator (default 5). |
level |
Numeric CI level passed to the estimator (default 0.95). |
seed |
Integer seed passed to the estimator (default 1). |
A named list with fields pmed_sobol, pmed_sobol_ci_lower,
pmed_sobol_ci_upper.
medsim_scenario_sobol(), medsim_analyze_coverage()
Inverts MCSE = sqrt(p(1-p)/nsim) to give the nsim such that a coverage
cell near p has Monte Carlo SE no larger than target_mcse.
medsim_nsim_for_mcse(target_mcse, p = 0.95)medsim_nsim_for_mcse(target_mcse, p = 0.95)
target_mcse |
Positive numeric target Monte Carlo SE. |
p |
Assumed coverage proportion (default 0.95, the nominal level). |
Integer number of replications (ceiling).
medsim_nsim_for_mcse(0.005, 0.95) # 1900medsim_nsim_for_mcse(0.005, 0.95) # 1900
Creates a combined figure with multiple panels for manuscript submission. Useful for showing multiple aspects of simulation results in one figure.
medsim_plot_combined_panel( results, panels = c("error", "timing"), layout = NULL, ... )medsim_plot_combined_panel( results, panels = c("error", "timing"), layout = NULL, ... )
results |
A medsim_results object or named list of results objects |
panels |
Character vector: which panels to include. Options: "error", "timing", "coverage" (default: c("error", "timing")) |
layout |
Character: layout specification like "2x1", "1x2", "2x2" (default: auto-determined from number of panels) |
... |
Additional arguments passed to individual plotting functions |
Requires either patchwork or gridExtra package for combining plots.
A combined plot object (patchwork or gridExtra)
## Not run: # Two-panel figure p <- medsim_plot_combined_panel( results, panels = c("error", "timing"), layout = "1x2" ) # Four-panel figure for manuscript p <- medsim_plot_combined_panel( list(Proposed = res1, Delta = res2), panels = c("error", "timing", "coverage"), layout = "2x2" ) ## End(Not run)## Not run: # Two-panel figure p <- medsim_plot_combined_panel( results, panels = c("error", "timing"), layout = "1x2" ) # Four-panel figure for manuscript p <- medsim_plot_combined_panel( list(Proposed = res1, Delta = res2), panels = c("error", "timing", "coverage"), layout = "2x2" ) ## End(Not run)
Creates a plot showing confidence interval coverage rates across scenarios and methods.
medsim_plot_coverage( coverage, expected_coverage = 0.95, tolerance = 0.02, title = NULL, ... )medsim_plot_coverage( coverage, expected_coverage = 0.95, tolerance = 0.02, title = NULL, ... )
coverage |
A medsim_coverage object from |
expected_coverage |
Numeric: expected coverage rate to show as reference line (default: 0.95) |
tolerance |
Numeric: tolerance band around expected coverage (default: 0.02) |
title |
Character: plot title (default: auto-generated) |
... |
Additional arguments passed to ggplot2::theme() |
A ggplot2 object
## Not run: coverage <- medsim_analyze_coverage(results) p <- medsim_plot_coverage(coverage, expected_coverage = 0.95) ## End(Not run)## Not run: coverage <- medsim_analyze_coverage(results) p <- medsim_plot_coverage(coverage, expected_coverage = 0.95) ## End(Not run)
Creates boxplots showing the distribution of estimation errors across scenarios. Useful for comparing accuracy of different methods.
medsim_plot_error_boxplot( results, parameter = "indirect", log_scale = FALSE, color_palette = NULL, title = NULL, ... )medsim_plot_error_boxplot( results, parameter = "indirect", log_scale = FALSE, color_palette = NULL, title = NULL, ... )
results |
A medsim_results object from |
parameter |
Character: which parameter to plot (default: "indirect"). If NULL, plots all parameters with available ground truth. |
log_scale |
Logical: use log10 scale for y-axis (default: FALSE) |
color_palette |
Character: color palette name from RColorBrewer or "viridis" (default: NULL uses ggplot2 defaults) |
title |
Character: plot title (default: auto-generated) |
... |
Additional arguments passed to ggplot2::theme() |
Box: IQR (25th to 75th percentile)
Line in box: Median error
Whiskers: 1.5 x IQR
Points: Outliers beyond whiskers
Tighter boxes = more consistent estimates
# Single method p <- medsim_plot_error_boxplot(results) # Compare methods p <- medsim_plot_error_boxplot(list( Proposed = results1, Delta = results2, Bootstrap = results3 )) # Customize p <- medsim_plot_error_boxplot( results, log_scale = TRUE, color_palette = "Set2" ) p + ggplot2::theme_minimal()
A ggplot2 object
## Not run: results <- medsim_run(method, scenarios, config) p <- medsim_plot_error_boxplot(results) print(p) # Save to file ggplot2::ggsave("error_boxplot.pdf", p, width = 10, height = 6) ## End(Not run)## Not run: results <- medsim_run(method, scenarios, config) p <- medsim_plot_error_boxplot(results) print(p) # Save to file ggplot2::ggsave("error_boxplot.pdf", p, width = 10, height = 6) ## End(Not run)
Plots each replication's model standard error against its point estimate, faceted by scenario, to expose numerator-denominator coupling and outliers (White et al. 2023, How to check a simulation study, IJE).
medsim_plot_se_vs_estimate(results, param = "indirect", by = "scenario")medsim_plot_se_vs_estimate(results, param = "indirect", by = "scenario")
results |
A results data.frame from |
param |
Parameter name; expects |
by |
Grouping/faceting column (default "scenario"). |
A ggplot object.
Creates bar plots comparing computational time across methods and scenarios.
medsim_plot_timing( results, metric = "mean", log_scale = FALSE, title = NULL, ... )medsim_plot_timing( results, metric = "mean", log_scale = FALSE, title = NULL, ... )
results |
A named list of medsim_results objects (required for comparison) |
metric |
Character: which timing metric to plot. Options: "mean", "median", "total" (default: "mean") |
log_scale |
Logical: use log10 scale for y-axis (default: FALSE) |
title |
Character: plot title (default: auto-generated) |
... |
Additional arguments passed to ggplot2::theme() |
A ggplot2 object
## Not run: p <- medsim_plot_timing(list( Proposed = results1, Delta = results2, Bootstrap = results3 )) ## End(Not run)## Not run: p <- medsim_plot_timing(list( Proposed = results1, Delta = results2, Bootstrap = results3 )) ## End(Not run)
Fleishman / Vale?Maurelli power-method draws. Pure base R; no hard dependency. Warns/errors on (skew, kurtosis) outside the Fleishman-feasible region.
medsim_rnonnormal(n, mean = 0, sd = 1, skew = 0, kurtosis = 0)medsim_rnonnormal(n, mean = 0, sd = 1, skew = 0, kurtosis = 0)
n |
Integer sample size. |
mean, sd
|
Target marginal mean and SD (applied after standardizing). |
skew |
Target marginal skewness (third standardized moment). |
kurtosis |
Target marginal excess kurtosis (fourth standardized moment minus 3). |
Numeric vector of length n with the requested marginal moments.
x <- medsim_rnonnormal(500, skew = 1.5, kurtosis = 4) c(mean = mean(x), sd = stats::sd(x))x <- medsim_rnonnormal(500, skew = 1.5, kurtosis = 4) c(mean = mean(x), sd = stats::sd(x))
Executes a complete simulation study by running a user-defined method across multiple scenarios and replications. Supports parallel processing and automatic result aggregation.
medsim_run( method, scenarios, config, compute_truth = NULL, parallel = TRUE, verbose = TRUE, truth_on_failure = c("stop", "warn") )medsim_run( method, scenarios, config, compute_truth = NULL, parallel = TRUE, verbose = TRUE, truth_on_failure = c("stop", "warn") )
method |
Function: User-defined simulation method. Must accept two
arguments: |
scenarios |
List of scenario objects from |
config |
Configuration object from |
compute_truth |
Function: Optional function to compute ground truth.
Takes same arguments as |
parallel |
Logical: Use parallel processing (default TRUE) |
verbose |
Logical: Print progress messages (default TRUE) |
truth_on_failure |
|
The method function must:
Accept data (data.frame) as first argument
Accept params (list) as second argument
Return a named list with at least one numeric element
NOT use the reserved field name error – it marks failed replications
in the results schema (a non-NA error forces converged = 0); a
method returning it stops with an informative error
Example:
my_method <- function(data, params) {
fit_m <- lm(M ~ X, data = data)
fit_y <- lm(Y ~ X + M, data = data)
a <- coef(fit_m)["X"]
b <- coef(fit_y)["M"]
list(
indirect = a * b,
a_path = a,
b_path = b,
se_indirect = ..., # Optional
ci_lower = ..., # Optional
ci_upper = ... # Optional
)
}
When parallel = TRUE:
Uses medsim_run_parallel() internally
Number of cores from config$n_cores
Progress bars shown (local) or suppressed (cluster)
All necessary objects exported to workers
If compute_truth is provided:
Computed once per scenario (cached automatically)
Used to calculate errors/bias in results
Can use Monte Carlo or analytical methods
compute_truth must return a named list of length-1 numerics, e.g.
list(indirect = 0.09). A bare named vector (c(indirect = 0.09)) or a list
element of length > 1 contributes no parameter and is treated as a failure.
A failed or unusable truth is never cached, so the next run retries it, and
the scenario is left out of $truth. What happens to the run is set by
truth_on_failure: "stop" (the default) raises a
medsim_truth_failure_error (subclass of medsim_error) carrying
$scenario and $parent; "warn" warns and continues. Truth is computed in
Step 1, before any replication runs, so stopping wastes no replication
compute.
An error that escapes a replication task – a contract violation (e.g. the
reserved-error-field stop) or an infrastructure error, as opposed to a
method() or data_generator() failure (both of which become failure
rows and let the run continue) – stops medsim_run() with a
medsim_run_failure condition (subclass of medsim_error) rather than
silently dropping the failed replications. The condition carries whatever
completed before the failure: $results is a partial medsim_results
object (prior scenarios complete, the failing scenario's successful reps
included, later scenarios never started), and $errors is a list of
list(scenario, global_rep_id, message). Recover the partial run instead
of losing it:
partial <- tryCatch( medsim_run(method, scenarios, config), medsim_run_failure = function(e) e$results )
A medsim_results object (list) containing:
results: data.frame with all simulation results
summary: data.frame with summary statistics
truth: data.frame with ground truth values (if compute_truth provided)
config: configuration used
scenarios: scenarios used
method_name: name of method function
timestamp: when simulation was run
medsim_config(), medsim_scenario(), medsim_run_parallel()
## Not run: # Define method my_method <- function(data, params) { fit_m <- lm(M ~ X, data = data) fit_y <- lm(Y ~ X + M, data = data) list(indirect = coef(fit_m)["X"] * coef(fit_y)["M"]) } # Configure and run config <- medsim_config("local") scenarios <- medsim_scenarios_mediation() results <- medsim_run( method = my_method, scenarios = scenarios, config = config ) # View results print(results) summary(results) ## End(Not run)## Not run: # Define method my_method <- function(data, params) { fit_m <- lm(M ~ X, data = data) fit_y <- lm(Y ~ X + M, data = data) list(indirect = coef(fit_m)["X"] * coef(fit_y)["M"]) } # Configure and run config <- medsim_config("local") scenarios <- medsim_scenarios_mediation() results <- medsim_run( method = my_method, scenarios = scenarios, config = config ) # View results print(results) summary(results) ## End(Not run)
Runs replications [(chunk_id-1)*chunk_size+1 .. chunk_id*chunk_size]
and saves a single RDS file named chunk_<chunk_id>.rds in output_dir.
Designed to be called from inside a SLURM array task script; chunk_id and
n_chunks are auto-detected from config (which auto-reads
SLURM_ARRAY_TASK_ID).
The chunk .rds is the sole artifact: in chunk mode the intermediate
per-scenario/summary CSVs that a standalone medsim_run() writes are
skipped (concurrent array tasks sharing an output_dir would clobber the
fixed-name files, leaving partial data that looks complete). Result rows
record the global replication id (schema v2): chunk 2 of a
4-chunk/nsim-20 run writes replication 6..10, not 1..5, so ids never
collide across chunks and medsim_combine_chunks() can audit the combined
grid.
Each chunk file also carries a provenance attribute (Gate C):
R version, medsim + key dependency versions, hostname, a code SHA, and
seconds-per-replication timing. The SHA is auto-detected from the run
script's git tree (falling back to a pkg:medsim-<version> tag outside
git); pass code_sha to stamp explicitly. medsim_combine_chunks()
asserts a single SHA across all chunks – catching a mid-run code edit plus
partial resubmit, which would silently mix results from two different
code states.
medsim_run_chunk( scenarios, method, config, verbose = TRUE, code_sha = NULL, compute_truth = NULL, truth_on_failure = c("stop", "warn") )medsim_run_chunk( scenarios, method, config, verbose = TRUE, code_sha = NULL, compute_truth = NULL, truth_on_failure = c("stop", "warn") )
scenarios |
A list of |
method |
A function with signature |
config |
A |
verbose |
Logical: print progress messages. |
code_sha |
Character or |
compute_truth |
Optional function with signature Chunk 1 then runs longer than the others (truth plus its replications);
size |
truth_on_failure |
|
Invisibly, the path of the RDS file written. If a task-level error
(a medsim_run_failure from medsim_run()) interrupts the chunk, whatever
replications completed are still written – atomically, with provenance
partial = TRUE – before this function re-stop()s so the array task
still exits nonzero (sacct visibility, not a --requeue signal).
This holds even when NO replication completed: a chunk whose every
replication failed writes a flagged, zero-row file rather than no file at
all. That file is the only carrier of the partial = TRUE marker, and
without it a trailing all-failed chunk would combine with zero violations
(the contiguity audit cannot see it – the surviving ids are a clean
1..k block). medsim_combine_chunks()/medsim_audit_results() flag any
such chunk via the partial_chunk violation.
medsim_combine_chunks(), medsim_audit_results(),
medsim_write_submit_script()
Executes a list of tasks in parallel using multiple CPU cores. Provides progress bars, proper error handling, and automatic cluster management.
medsim_run_parallel( tasks, fun, n_cores = NULL, progress = TRUE, export = NULL, packages = NULL, cluster_type = NULL, seed = NULL )medsim_run_parallel( tasks, fun, n_cores = NULL, progress = TRUE, export = NULL, packages = NULL, cluster_type = NULL, seed = NULL )
tasks |
Vector or list: Tasks to process (e.g., 1:100, list of params) |
fun |
Function: Function to apply to each task. Should accept one argument (the task) and return a result. |
n_cores |
Integer: Number of CPU cores to use. If NULL, uses all available cores minus 2. |
progress |
Logical: Show progress bar (default TRUE). Automatically suppressed when running on cluster (batch jobs). |
export |
Character vector: Names of objects to export to workers. These objects must be available in the calling environment. |
packages |
Character vector: Names of packages to load on each worker.
Use this if |
cluster_type |
Character: Type of cluster ("PSOCK" or "FORK"). FORK is more efficient but only works on Unix. Auto-detected by default. |
seed |
Integer or |
PSOCK (default on Windows): Creates separate R sessions. Requires explicit export of objects and packages. Works on all platforms.
FORK (default on Unix): Copies current R session. More efficient, automatic object sharing, but Unix-only.
Progress bars are shown in interactive sessions and local execution, but suppressed on HPC clusters to avoid cluttering log files.
If a task fails:
Error is caught and returned as an error object
Other tasks continue processing
Check results with sapply(results, inherits, "error")
Each worker process requires memory. For large tasks:
Reduce n_cores if running out of memory
Process tasks in chunks
Use FORK cluster (Unix) which shares memory
List: Results from applying fun to each task (same length as tasks)
parallel::makeCluster(), parallel::parLapply(), pbapply::pblapply()
## Not run: # Simple parallel computation results <- medsim_run_parallel( tasks = 1:100, fun = function(i) { Sys.sleep(0.1) # Simulate work i^2 }, n_cores = 4 ) # With exports my_data <- data.frame(x = 1:10, y = 11:20) results <- medsim_run_parallel( tasks = 1:10, fun = function(i) { mean(my_data$x[1:i]) }, export = "my_data" ) # With packages results <- medsim_run_parallel( tasks = 1:10, fun = function(i) { dplyr::tibble(x = i, y = i^2) }, packages = "dplyr" ) # Check for errors errors <- sapply(results, inherits, "error") if (any(errors)) { cat("Errors occurred in:", which(errors), "\n") } ## End(Not run)## Not run: # Simple parallel computation results <- medsim_run_parallel( tasks = 1:100, fun = function(i) { Sys.sleep(0.1) # Simulate work i^2 }, n_cores = 4 ) # With exports my_data <- data.frame(x = 1:10, y = 11:20) results <- medsim_run_parallel( tasks = 1:10, fun = function(i) { mean(my_data$x[1:i]) }, export = "my_data" ) # With packages results <- medsim_run_parallel( tasks = 1:10, fun = function(i) { dplyr::tibble(x = i, y = i^2) }, packages = "dplyr" ) # Check for errors errors <- sapply(results, inherits, "error") if (any(errors)) { cat("Errors occurred in:", which(errors), "\n") } ## End(Not run)
Define a custom scenario for simulation studies. A scenario consists of a name, description, data generation function, and population parameters.
medsim_scenario( name, description = "", data_generator, params = list(), estimand = NULL )medsim_scenario( name, description = "", data_generator, params = list(), estimand = NULL )
name |
Character: Descriptive name for the scenario |
description |
Character: Detailed description (optional) |
data_generator |
Function: Takes n (sample size) and returns data.frame |
params |
List: Known population parameters for validation |
estimand |
Optional |
The data_generator function must:
Accept n (sample size) as first argument
Return a data.frame with at least: X (treatment), M (mediator), Y (outcome)
Can include additional variables or covariates
The params list should include known population values:
a: X -> M path coefficient
b: M -> Y path coefficient
c_prime: X -> Y direct effect
indirect: True indirect effect (a * b)
Additional parameters as needed
A scenario object (list with class "medsim_scenario")
medsim_scenarios_mediation() for standard scenarios
# Simple custom scenario my_scenario <- medsim_scenario( name = "Large Effects", description = "Both paths have large effects", data_generator = function(n = 200) { X <- rnorm(n) M <- 0.7 * X + rnorm(n) Y <- 0.7 * M + rnorm(n) data.frame(X = X, M = M, Y = Y) }, params = list( a = 0.7, b = 0.7, indirect = 0.49 ) ) # Generate data data <- my_scenario$data_generator(n = 100) # Scenario with covariates covariate_scenario <- medsim_scenario( name = "With Covariates", data_generator = function(n = 200) { C1 <- rnorm(n) C2 <- rbinom(n, 1, 0.5) X <- rnorm(n) M <- 0.3 * X + 0.2 * C1 + rnorm(n) Y <- 0.3 * M + 0.2 * C2 + rnorm(n) data.frame(X = X, M = M, Y = Y, C1 = C1, C2 = C2) }, params = list( a = 0.3, b = 0.3, indirect = 0.09, gamma_m = 0.2, gamma_y = 0.2 ) )# Simple custom scenario my_scenario <- medsim_scenario( name = "Large Effects", description = "Both paths have large effects", data_generator = function(n = 200) { X <- rnorm(n) M <- 0.7 * X + rnorm(n) Y <- 0.7 * M + rnorm(n) data.frame(X = X, M = M, Y = Y) }, params = list( a = 0.7, b = 0.7, indirect = 0.49 ) ) # Generate data data <- my_scenario$data_generator(n = 100) # Scenario with covariates covariate_scenario <- medsim_scenario( name = "With Covariates", data_generator = function(n = 200) { C1 <- rnorm(n) C2 <- rbinom(n, 1, 0.5) X <- rnorm(n) M <- 0.3 * X + 0.2 * C1 + rnorm(n) Y <- 0.3 * M + 0.2 * C2 + rnorm(n) data.frame(X = X, M = M, Y = Y, C1 = C1, C2 = C2) }, params = list( a = 0.3, b = 0.3, indirect = 0.09, gamma_m = 0.2, gamma_y = 0.2 ) )
Builds a medsim_scenario() for the bounds-pair simulation studies
(me-mediator-bounds / me-exposure-recall). The data_generator wraps
medrobust::simulate_dm_data(...)@observed; if medrobust is not
installed a lightweight synthetic fallback is used so that the harness and
interval-kind analysis machinery can be unit-tested without the package.
The scenario is tagged with estimand = medsim_estimand("interval", ...) so
that medsim_analyze_coverage() automatically dispatches the
Imbens-Manski / partial-ID interval branch.
medsim_scenario_dm( name, true_params, dm_params = list(), misclass_type = c("mediator", "exposure") )medsim_scenario_dm( name, true_params, dm_params = list(), misclass_type = c("mediator", "exposure") )
name |
Scenario name (character). |
true_params |
List of true generating parameters passed through to
|
dm_params |
List of differential-misclassification parameters passed
through to |
misclass_type |
|
A medsim_scenario() object with
estimand = medsim_estimand("interval", params = c("NDE", "NIE"), ...).
medsim_method_bounds() for the corresponding estimator adapter;
medsim_analyze_coverage() for interval-kind coverage analysis.
# Synthetic scenario (no medrobust required): sc <- medsim_scenario_dm( name = "small_misclass", true_params = list(NDE = 0.2, NIE = 0.3, n = 200), dm_params = list(delta = 0.1), misclass_type = "mediator" ) # Data generator uses synthetic fallback when medrobust is absent: d <- sc$data_generator(n = 50) head(d)# Synthetic scenario (no medrobust required): sc <- medsim_scenario_dm( name = "small_misclass", true_params = list(NDE = 0.2, NIE = 0.3, n = 200), dm_params = list(delta = 0.1), misclass_type = "mediator" ) # Data generator uses synthetic fallback when medrobust is absent: d <- sc$data_generator(n = 50) head(d)
Thin wrapper over medsim_scenario() mirroring medsim_scenario_sobol():
an A -> M -> Y linear-Gaussian DGP with covariate C and A*M
interaction (kappa), a variance_share estimand, and closed-form
ground truth for the interventional proportion mediated P_med = IIE/OE
and the gauge residual W = R/OE, where R = OE - IDE - IIE.
Corner means theta_{a,a'} = tau_a*a + (tau_m + kappa*a)*beta_a*a' (C
centered) give OE = theta11-theta00, IDE = theta10-theta00,
IIE = theta01-theta00. At kappa = 0, R = 0 and W = 0.
medsim_scenario_gauge(name, true_params = list())medsim_scenario_gauge(name, true_params = list())
name |
Character scenario name. |
true_params |
Named list; defaults |
A medsim_scenario with estimand$kind = "variance_share" and
params$pmed, params$w closed-form truths.
medsim_method_gauge(), medsim_scenario_sobol()
sc <- medsim_scenario_gauge("interaction", list(kappa = 0.4)) c(sc$params$pmed, sc$params$w)sc <- medsim_scenario_gauge("interaction", list(kappa = 0.4)) c(sc$params$pmed, sc$params$w)
Wraps medsim_scenario() with a data_generator that (1) draws complete X -> M -> Y
with optional nonnormal residuals (medsim_rnonnormal()), then (2) amputes target
under mechanism (medsim_amputate()).
medsim_scenario_missing( name, true_params, mechanism, prop = 0.2, target = "M", nonnormal = NULL )medsim_scenario_missing( name, true_params, mechanism, prop = 0.2, target = "M", nonnormal = NULL )
name |
Scenario name. |
true_params |
List of true generating parameters: |
mechanism |
One of |
prop |
Target missingness proportion. |
target |
Column(s) to make missing (default |
nonnormal |
|
A medsim_scenario object.
scn <- medsim_scenario_missing( name = "MAR_M_30", true_params = list(a = 0.39, b = 0.39, cp = 0.14), mechanism = "MAR", prop = 0.3 ) scn$namescn <- medsim_scenario_missing( name = "MAR_M_30", true_params = list(a = 0.39, b = 0.39, cp = 0.14), mechanism = "MAR", prop = 0.3 ) scn$name
Convenience: expand SPEC-simulation-design cells (mechanism x prop x n x effect x
nonnormality) into a list of medsim_scenario_missing() objects in one call.
medsim_scenario_missing_grid( true_params_list, mechanisms, props, nonnormal_list = list(NULL) )medsim_scenario_missing_grid( true_params_list, mechanisms, props, nonnormal_list = list(NULL) )
true_params_list |
List of |
mechanisms, props
|
Character/numeric vectors crossed factorially. |
nonnormal_list |
List of |
A list of medsim_scenario objects.
grid <- medsim_scenario_missing_grid( true_params_list = list(small = list(a = 0.14, b = 0.14, cp = 0.14)), mechanisms = c("MCAR", "MAR"), props = c(0.1, 0.3) ) length(grid)grid <- medsim_scenario_missing_grid( true_params_list = list(small = list(a = 0.14, b = 0.14, cp = 0.14)), mechanisms = c("MCAR", "MAR"), props = c(0.1, 0.3) ) length(grid)
Thin wrapper over medsim_scenario() for studies that evaluate numerical
accuracy or timing of a computation – not coverage or power. The
estimand kind is "numeric", which signals downstream commands to skip
coverage/power analysis and produce accuracy/timing tables instead.
Typical use-cases:
Product-of-three CI: compare DOP, MBCO, Monte-Carlo CI point accuracy and computation time across parameter cells.
Approximation quality: absolute error and relative error vs. analytic truth.
The data_generator for numeric scenarios is optional (default: returns
data.frame() with n rows of NAs as a placeholder). Methods that
only need params can ignore the data argument.
medsim_scenario_numeric( name, true_params = list(), data_generator = NULL, extra = character() )medsim_scenario_numeric( name, true_params = list(), data_generator = NULL, extra = character() )
name |
Character: scenario name. |
true_params |
Named list of ground-truth parameter values. These are
passed to |
data_generator |
Function |
extra |
Character vector of additional result columns beyond the
mandatory |
A medsim_scenario object with estimand$kind = "numeric".
medsim_estimand(), medsim_scenario(), medsim_analyze()
# Product-of-three scenario cell sc <- medsim_scenario_numeric( name = "dop_small_n", true_params = list(a = 0.5, b = 0.3, se_a = 0.1, se_b = 0.1, ci_true = 0.01), extra = c("elapsed_sec") ) str(sc$estimand)# Product-of-three scenario cell sc <- medsim_scenario_numeric( name = "dop_small_n", true_params = list(a = 0.5, b = 0.3, se_a = 0.1, se_b = 0.1, ci_true = 0.01), extra = c("elapsed_sec") ) str(sc$estimand)
Thin wrapper over medsim_scenario() that bakes in:
A linear structural equation model (SEM) data generator for the
A -> M -> Y causal chain.
An estimand = medsim_estimand("probabilistic", params = "pmed", ci = "mbco", extra = "branch_switch") descriptor so that
medsim_analyze_coverage() dispatches the MBCO-CI coverage branch.
An exact closed-form ground-truth P_med: under the all-Gaussian linear
SEM the cross-world difference Y1 - Y0 is Normal, so
P_med = Phi((beta_ay + alpha_ax * beta_my) / sqrt(2 * (beta_my^2 * sigma_m^2 + sigma_y^2))).
The estimand P_med = P(Y_a=1(M_a=1) > Y_a=0(M_a=1)) + 0.5 * P(Y_a=1(M_a=1) == Y_a=0(M_a=1)) uses the cross-world assumption –
it is a probabilistic mediation effect, not a difference in expectations.
medsim_scenario_pmed(name, true_params = list())medsim_scenario_pmed(name, true_params = list())
name |
Character: scenario name passed to |
true_params |
Named list with entries:
|
A medsim_scenario object with estimand$kind = "probabilistic".
medsim_method_pmed_mbco(), medsim_analyze_coverage()
sc <- medsim_scenario_pmed( name = "full_mediation", true_params = list(alpha_ax = 0.6, beta_my = 0.5, beta_ay = 0.0) ) str(sc$estimand)sc <- medsim_scenario_pmed( name = "full_mediation", true_params = list(alpha_ax = 0.6, beta_my = 0.5, beta_ay = 0.0) ) str(sc$estimand)
Thin wrapper over medsim_scenario() that bakes in:
A linear-Gaussian-with-interaction data generator for the
A -> M -> Y causal chain with a single covariate C and an A*M
interaction term (kappa).
An estimand = medsim_estimand("variance_share", params = "pmed_sobol", ci = "standard") descriptor, so medsim_analyze_coverage() uses the
generic Wald truth %in% CI path and medsim_validate_scenario() checks
for the A/M/Y/C columns the Sobol estimator requires.
A closed-form truth for the Sobol share, exact for the
linear-Gaussian-with-interaction DGP (no Monte-Carlo PO draw needed).
The estimand is the variance-scale proportion mediated
, where the functional-ANOVA variance
components are built from the four corner means
. For the linear-Gaussian DGP with
centered C,
At kappa = 0 the share reduces to with
and .
medsim_scenario_sobol(name, true_params = list(), pd = 0.5, pm = 0.5)medsim_scenario_sobol(name, true_params = list(), pd = 0.5, pm = 0.5)
name |
Character: scenario name passed to |
true_params |
Named list with entries (defaults in parentheses):
|
pd, pm
|
Bernoulli weights for the direct / mediator design variances in
the functional-ANOVA decomposition (default |
A medsim_scenario object with estimand$kind = "variance_share"
and params$pmed_sobol holding the closed-form ground truth.
medsim_method_sobol(), medsim_analyze_coverage(),
medsim_estimand()
sc <- medsim_scenario_sobol( name = "interaction", true_params = list(beta_a = 0.6, tau_a = 0.5, tau_m = 0.7, kappa = 0.4) ) sc$params$pmed_sobol # closed-form Sobol sharesc <- medsim_scenario_sobol( name = "interaction", true_params = list(beta_a = 0.6, tau_a = 0.5, tau_m = 0.7, kappa = 0.4) ) sc$params$pmed_sobol # closed-form Sobol share
Creates a list of standard mediation scenarios for simulation studies. These scenarios cover common patterns in mediation analysis including independent paths, various correlation structures, suppression effects, and non-standard conditions.
medsim_scenarios_mediation()medsim_scenarios_mediation()
Independent Paths: No correlation between X, M, Y
Moderate Correlation: rho = 0.3 between all pairs
High Correlation: rho = 0.7 between all pairs
Suppression: Mixed positive and negative correlations
Non-zero Effects: Small to moderate true effects
Unequal Variances: Different residual variances
Each scenario generates data with:
Sample size n (default: 200)
Treatment X, Mediator M, Outcome Y
Known population parameters for validation
A list of scenario objects, each with name, description, and data_generator function
medsim_scenario() for creating custom scenarios
# Get all standard scenarios scenarios <- medsim_scenarios_mediation() # See scenario names sapply(scenarios, function(s) s$name) # Generate data from first scenario data <- scenarios[[1]]$data_generator(n = 100) head(data) # Access scenario parameters scenarios[[1]]$params# Get all standard scenarios scenarios <- medsim_scenarios_mediation() # See scenario names sapply(scenarios, function(s) s$name) # Generate data from first scenario data <- scenarios[[1]]$data_generator(n = 100) head(data) # Access scenario parameters scenarios[[1]]$params
Reads the branch_switch column emitted by medsim_method_mbco_mi(): for each
replicate, the union-null (ab = 0) branch chosen by the majority of imputations
(1 = the a = 0 branch, 0 = the b = 0 branch). The returned branch_switch_rate
is therefore the share of replicates whose majority branch is a = 0, not the share
in which imputations disagree; that is the branch_mix column. The result is suitable
for joining as a column onto a medsim_table_coverage() table. Non-converged rows
(converged == 0) are dropped before summarizing.
medsim_summarize_branch_switch(results, by = "scenario")medsim_summarize_branch_switch(results, by = "scenario")
results |
A |
by |
Grouping columns (default |
A data.frame: one row per group with branch_switch_rate and n_valid.
res <- data.frame( scenario = rep(c("s1", "s2"), each = 4), branch_switch = c(1, 1, 0, 1, 0, 0, 1, 0), converged = 1 ) medsim_summarize_branch_switch(res)res <- data.frame( scenario = rep(c("s1", "s2"), each = 4), branch_switch = c(1, 1, 0, 1, 0, 0, 1, 0), converged = 1 ) medsim_summarize_branch_switch(res)
Creates a publication-ready LaTeX table showing accuracy metrics (bias, MAE, RMSE) from simulation analysis results.
medsim_table_accuracy( analysis, digits = 3, metrics = c("bias", "mae", "rmse"), by_scenario = TRUE, caption = NULL, label = "tab:accuracy" )medsim_table_accuracy( analysis, digits = 3, metrics = c("bias", "mae", "rmse"), by_scenario = TRUE, caption = NULL, label = "tab:accuracy" )
analysis |
A medsim_analysis object from |
digits |
Integer: number of decimal places (default: 3) |
metrics |
Character vector: which metrics to include. Options: "bias", "mae", "rmse", "median_ae", "max_ae", "relative_bias" (default: c("bias", "mae", "rmse")) |
by_scenario |
Logical: include scenario breakdown (default: TRUE) |
caption |
Character: table caption (default: auto-generated) |
label |
Character: LaTeX label for cross-referencing (default: "tab:accuracy") |
A medsim_table object (character vector with LaTeX code)
## Not run: results <- medsim_run(method, scenarios, config) analysis <- medsim_analyze(results) # Generate table table <- medsim_table_accuracy(analysis) print(table) # Preview # Write to file medsim_write_table(table, "tables/accuracy.tex") ## End(Not run)## Not run: results <- medsim_run(method, scenarios, config) analysis <- medsim_analyze(results) # Generate table table <- medsim_table_accuracy(analysis) print(table) # Preview # Write to file medsim_write_table(table, "tables/accuracy.tex") ## End(Not run)
Creates a publication-ready LaTeX table comparing multiple methods across accuracy, timing, and coverage metrics.
medsim_table_comparison( comparison, metrics = c("accuracy", "timing"), caption = NULL, label = "tab:comparison" )medsim_table_comparison( comparison, metrics = c("accuracy", "timing"), caption = NULL, label = "tab:comparison" )
comparison |
A medsim_comparison object from |
metrics |
Character vector: which metrics to include ("accuracy", "timing", "coverage") |
caption |
Character: table caption |
label |
Character: LaTeX label |
A medsim_table object
Creates a publication-ready LaTeX table showing confidence interval coverage rates.
medsim_table_coverage( coverage, expected = 0.95, by_scenario = TRUE, caption = NULL, label = "tab:coverage" )medsim_table_coverage( coverage, expected = 0.95, by_scenario = TRUE, caption = NULL, label = "tab:coverage" )
coverage |
A medsim_coverage object from |
expected |
Numeric: expected/nominal coverage rate (default: 0.95) |
by_scenario |
Logical: show coverage by scenario (default: TRUE) |
caption |
Character: table caption |
label |
Character: LaTeX label |
A medsim_table object
Render a performance summary as a LaTeX table
medsim_table_performance(perf)medsim_table_performance(perf)
perf |
Output of |
list(latex=, type="performance").
Creates a publication-ready LaTeX table showing empirical power (rejection rates) from simulation analysis.
medsim_table_power( power, by_scenario = TRUE, caption = NULL, label = "tab:power" )medsim_table_power( power, by_scenario = TRUE, caption = NULL, label = "tab:power" )
power |
A medsim_power object from |
by_scenario |
Logical: show power by scenario (default: TRUE) |
caption |
Character: table caption |
label |
Character: LaTeX label |
A medsim_table object
Creates a publication-ready LaTeX table comparing computational time across methods.
medsim_table_timing( results_list, metric = "mean", include_speedup = TRUE, reference_method = NULL, by_scenario = FALSE, caption = NULL, label = "tab:timing" )medsim_table_timing( results_list, metric = "mean", include_speedup = TRUE, reference_method = NULL, by_scenario = FALSE, caption = NULL, label = "tab:timing" )
results_list |
Named list of medsim_results objects, one per method |
metric |
Character: timing metric to use ("mean", "median", "total") |
include_speedup |
Logical: include speedup column relative to slowest method |
reference_method |
Character: method name to use as reference for speedup (default: NULL uses slowest method) |
by_scenario |
Logical: show timing by scenario (default: FALSE) |
caption |
Character: table caption |
label |
Character: LaTeX label |
A medsim_table object
## Not run: table <- medsim_table_timing(list( Proposed = results1, Delta = results2, Bootstrap = results3 )) ## End(Not run)## Not run: table <- medsim_table_timing(list( Proposed = results1, Delta = results2, Bootstrap = results3 )) ## End(Not run)
Convenience wrapper around medsim_tables_workflow() with sensible defaults.
medsim_tables(results, output_dir = "tables", format = "latex", ...)medsim_tables(results, output_dir = "tables", format = "latex", ...)
results |
A |
output_dir |
Directory to write tables to. Created if it doesn't exist. |
format |
Output format passed to |
... |
Additional arguments passed to |
Invisibly, the result of medsim_tables_workflow().
Convenience function to generate all tables from simulation results and save them to an output directory.
medsim_tables_workflow( results, output_dir, tables = "all", format = "latex", prefix = "table_" )medsim_tables_workflow( results, output_dir, tables = "all", format = "latex", prefix = "table_" )
results |
A medsim_results object from |
output_dir |
Character: output directory for tables |
tables |
Character vector: which tables to generate ("accuracy", "coverage", "power", "all") |
format |
Character: output format ("latex", "csv", "markdown") |
prefix |
Character: prefix for file names (default: "table_") |
Invisibly returns a list of generated file paths
## Not run: results <- medsim_run(method, scenarios, config) # Generate all tables files <- medsim_tables_workflow( results, output_dir = "manuscript/tables", format = "latex" ) ## End(Not run)## Not run: results <- medsim_run(method, scenarios, config) # Generate all tables files <- medsim_tables_workflow( results, output_dir = "manuscript/tables", format = "latex" ) ## End(Not run)
Checks that a scenario object is valid and can generate data correctly.
medsim_validate_scenario(scenario, n = 10)medsim_validate_scenario(scenario, n = 10)
scenario |
A medsim_scenario object |
n |
Sample size to test (default: 10) |
TRUE if valid, throws error otherwise
scenarios <- medsim_scenarios_mediation() medsim_validate_scenario(scenarios[[1]])scenarios <- medsim_scenarios_mediation() medsim_validate_scenario(scenarios[[1]])
End-to-end output pipeline: runs medsim_analyze(), medsim_figures(),
and medsim_tables() and writes everything to output_dir.
medsim_workflow( results, output_dir = "results", figures_format = "pdf", tables_format = "latex", ... )medsim_workflow( results, output_dir = "results", figures_format = "pdf", tables_format = "latex", ... )
results |
A |
output_dir |
Top-level directory for all outputs. Subdirectories
|
figures_format |
File format for figures (e.g. |
tables_format |
Output format for tables (e.g. |
... |
Additional arguments passed to |
Invisibly, a list with elements analysis, figures (file paths),
and tables (return value of medsim_tables_workflow()).
## Not run: config <- medsim_config("test") scenarios <- medsim_scenarios_mediation() results <- medsim_run(my_method, scenarios, config) medsim_workflow(results, output_dir = "results") ## End(Not run)## Not run: config <- medsim_config("test") scenarios <- medsim_scenarios_mediation() results <- medsim_run(my_method, scenarios, config) medsim_workflow(results, output_dir = "results") ## End(Not run)
Generates a submit_array.sh bash script targeting the UNM CARC Hopper
cluster (and compatible SLURM systems). The script submits one SLURM array
job per scenario chunk; each task runs an R script that calls
medsim_run_chunk().
The emitted script is fail-loud (Gate B of the chunked-run integrity layer):
a login shell (#!/bin/bash -l – on Hopper module is only defined in
login shells), set -eo pipefail, a hard-failing module load (never
|| true), a command -v Rscript pre-check, #SBATCH --requeue, and the
Rscript call as the final command so its exit code is the task's exit
code. There is deliberately NO output-file gate in the script: the writer's
config and the run script's runtime config are independent, so a baked-in
path could fail successful tasks; completeness is audited at combine time by
medsim_combine_chunks() instead. Set config$array_throttle = K to cap
concurrently-running array tasks (--array=1-N%K).
medsim_write_submit_script( config, run_script = "run_simulation_chunk.R", output_file = "submit_array.sh", account = NULL )medsim_write_submit_script( config, run_script = "run_simulation_chunk.R", output_file = "submit_array.sh", account = NULL )
config |
A |
run_script |
Character: path (on the cluster) to the per-chunk R script
(the one that calls |
output_file |
Character: where to write the bash script. Default
|
account |
Character: SLURM account/PI account string (e.g. |
Invisibly, the path to the written script.
medsim_run_chunk(), medsim_combine_chunks()
cfg <- medsim_config( mode = "cluster", n_replications = 1000, n_chunks = 10, walltime = "04:00:00", mem_per_cpu = "8G" ) tmp <- tempfile(fileext = ".sh") medsim_write_submit_script(cfg, output_file = tmp) # cat(readLines(tmp), sep = "\n")cfg <- medsim_config( mode = "cluster", n_replications = 1000, n_chunks = 10, walltime = "04:00:00", mem_per_cpu = "8G" ) tmp <- tempfile(fileext = ".sh") medsim_write_submit_script(cfg, output_file = tmp) # cat(readLines(tmp), sep = "\n")
Writes a medsim_table object to a file in LaTeX, CSV, or Markdown format.
medsim_write_table(table, file, format = NULL)medsim_write_table(table, file, format = NULL)
table |
A medsim_table object |
file |
Character: output file path |
format |
Character: output format ("latex", "csv", "markdown") If NULL, inferred from file extension. |
Invisibly returns the file path
## Not run: table <- medsim_table_accuracy(analysis) # Write LaTeX medsim_write_table(table, "tables/accuracy.tex") # Write CSV medsim_write_table(table, "tables/accuracy.csv", format = "csv") ## End(Not run)## Not run: table <- medsim_table_accuracy(analysis) # Write LaTeX medsim_write_table(table, "tables/accuracy.tex") # Write CSV medsim_write_table(table, "tables/accuracy.csv", format = "csv") ## End(Not run)
Print Analysis Results
## S3 method for class 'medsim_analysis' print(x, ...)## S3 method for class 'medsim_analysis' print(x, ...)
x |
A medsim_analysis object |
... |
Additional arguments (ignored) |
Invisibly returns x
Print Method Comparison
## S3 method for class 'medsim_comparison' print(x, ...)## S3 method for class 'medsim_comparison' print(x, ...)
x |
A medsim_comparison object |
... |
Additional arguments (ignored) |
Invisibly returns x
Print Configuration Summary
## S3 method for class 'medsim_config' print(x, ...)## S3 method for class 'medsim_config' print(x, ...)
x |
A medsim_config object |
... |
Additional arguments (ignored) |
Invisibly returns x
Print Coverage Results
## S3 method for class 'medsim_coverage' print(x, ...)## S3 method for class 'medsim_coverage' print(x, ...)
x |
A medsim_coverage object |
... |
Additional arguments (ignored) |
Invisibly returns x
Print Power Results
## S3 method for class 'medsim_power' print(x, ...)## S3 method for class 'medsim_power' print(x, ...)
x |
A medsim_power object |
... |
Additional arguments (ignored) |
Invisibly returns x
Print Simulation Results
## S3 method for class 'medsim_results' print(x, ...)## S3 method for class 'medsim_results' print(x, ...)
x |
A medsim_results object |
... |
Additional arguments (ignored) |
Invisibly returns x
Print Scenario Summary
## S3 method for class 'medsim_scenario' print(x, ...)## S3 method for class 'medsim_scenario' print(x, ...)
x |
A medsim_scenario object |
... |
Additional arguments (ignored) |
Invisibly returns x
Print medsim_table
## S3 method for class 'medsim_table' print(x, ...)## S3 method for class 'medsim_table' print(x, ...)
x |
A medsim_table object |
... |
Additional arguments (ignored) |
Invisibly returns x
Summarize Simulation Results
## S3 method for class 'medsim_results' summary(object, ...)## S3 method for class 'medsim_results' summary(object, ...)
object |
A medsim_results object |
... |
Additional arguments (ignored) |
Summary data.frame