--- name: simulation-architect description: Design and implementation of comprehensive simulation studies --- # Simulation Architect You are an expert in designing Monte Carlo simulation studies for statistical methodology research. ## Morris et al Guidelines ### The ADEMP Framework (Morris et al., 2019, Statistics in Medicine) The definitive guide for simulation study design requires five components: | Component | Question | Documentation Required | |-----------|----------|----------------------| | **A**ims | What are we trying to learn? | Clear research questions | | **D**ata-generating mechanisms | How do we create data? | Full DGP specification | | **E**stimands | What are we estimating? | Mathematical definition | | **M**ethods | What estimators do we compare? | Complete algorithm description | | **P**erformance measures | How do we evaluate? | Bias, variance, coverage | ### Morris et al. Reporting Checklist ```markdown □ Aims stated clearly □ DGP fully specified (all parameters, distributions) □ Estimand(s) defined mathematically □ All methods described with sufficient detail for replication □ Performance measures defined □ Number of replications justified □ Monte Carlo standard errors reported □ Random seed documented for reproducibility □ Software and version documented □ Computational time reported ``` --- ## Replication Counts ### How Many Replications Are Needed? **Monte Carlo Standard Error (MCSE)** formula: $$\text{MCSE}(\hat{\theta}) = \frac{\hat{\sigma}}{\sqrt{B}}$$ where $B$ is the number of replications and $\hat{\sigma}$ is the estimated standard deviation. ### Recommended Replications by Purpose | Purpose | Minimum B | Recommended B | MCSE for proportion | |---------|-----------|---------------|---------------------| | Exploratory | 500 | 1,000 | ~1.4% at 95% coverage | | Publication | 1,000 | 2,000 | ~1.0% at 95% coverage | | Definitive | 5,000 | 10,000 | ~0.4% at 95% coverage | | Precision | 10,000+ | 50,000 | ~0.2% at 95% coverage | ### MCSE Calculation ```r # Calculate Monte Carlo standard errors calculate_mcse <- function(estimates, coverage_indicators = NULL) { B <- length(estimates) list( # MCSE for mean (bias) mcse_mean = sd(estimates) / sqrt(B), # MCSE for standard deviation mcse_sd = sd(estimates) / sqrt(2 * (B - 1)), # MCSE for coverage (proportion) mcse_coverage = if (!is.null(coverage_indicators)) { p <- mean(coverage_indicators) sqrt(p * (1 - p) / B) } else NA ) } # Rule of thumb: B needed for desired MCSE replications_needed <- function(desired_mcse, estimated_sd) { ceiling((estimated_sd / desired_mcse)^2) } ``` --- ## R Code Templates ### Complete Simulation Template ```r # Full simulation study template following Morris et al. guidelines run_simulation_study <- function( n_sims = 2000, n_vec = c(200, 500, 1000), seed = 42, parallel = TRUE, n_cores = parallel::detectCores() - 1 ) { set.seed(seed) # Define parameter grid params <- expand.grid( n = n_vec, effect_size = c(0, 0.14, 0.39), model_spec = c("correct", "misspecified") ) # Setup parallel processing if (parallel) { cl <- parallel::makeCluster(n_cores) doParallel::registerDoParallel(cl) on.exit(parallel::stopCluster(cl)) } # Run simulations results <- foreach( i = 1:nrow(params), .combine = rbind, .packages = c("tidyverse", "mediation") ) %dopar% { scenario <- params[i, ] sim_results <- replicate(n_sims, { data <- generate_dgp(scenario) estimates <- apply_methods(data) evaluate_performance(estimates, truth = scenario$effect_size) }, simplify = FALSE) summarize_scenario(sim_results, scenario) } # Add MCSE results <- add_monte_carlo_errors(results, n_sims) results } # Summarize with MCSE add_monte_carlo_errors <- function(results, B) { results %>% mutate( mcse_bias = empirical_se / sqrt(B), mcse_coverage = sqrt(coverage * (1 - coverage) / B), mcse_rmse = rmse / sqrt(2 * B) ) } ``` ### Parallel Simulation Template ```r # Memory-efficient parallel simulation run_parallel_simulation <- function(scenario, n_sims, n_cores = 4) { library(future) library(future.apply) plan(multisession, workers = n_cores) results <- future_replicate(n_sims, { data <- generate_dgp(scenario$n, scenario$params) est <- estimate_effect(data) list( estimate = est$point, se = est$se, covered = abs(est$point - scenario$truth) < 1.96 * est$se ) }, simplify = FALSE) plan(sequential) # Reset # Aggregate estimates <- sapply(results, `[[`, "estimate") ses <- sapply(results, `[[`, "se") covered <- sapply(results, `[[`, "covered") list( bias = mean(estimates) - scenario$truth, empirical_se = sd(estimates), mean_se = mean(ses), coverage = mean(covered), mcse_bias = sd(estimates) / sqrt(n_sims), mcse_coverage = sqrt(mean(covered) * (1 - mean(covered)) / n_sims) ) } ``` --- ## Core Principles (Morris et al., 2019) ### ADEMP Framework 1. **Aims**: What question does the simulation answer? 2. **Data-generating mechanisms**: How is data simulated? 3. **Estimands**: What is being estimated? 4. **Methods**: What estimators are compared? 5. **Performance measures**: How is performance assessed? ## Data-Generating Process Design ### Standard Mediation DGP ```r generate_mediation_data <- function(n, params) { # Confounders X <- rnorm(n) # Treatment (binary) ps <- plogis(params$gamma0 + params$gamma1 * X) A <- rbinom(n, 1, ps) # Mediator M <- params$alpha0 + params$alpha1 * A + params$alpha2 * X + rnorm(n, sd = params$sigma_m) # Outcome Y <- params$beta0 + params$beta1 * A + params$beta2 * M + params$beta3 * X + params$beta4 * A * M + rnorm(n, sd = params$sigma_y) data.frame(Y = Y, A = A, M = M, X = X) } ``` ### DGP Variations to Consider - **Linearity**: Linear vs nonlinear relationships - **Model specification**: Correct vs misspecified - **Error structure**: Homoscedastic vs heteroscedastic - **Interaction**: No interaction vs A×M interaction - **Confounding**: Measured vs unmeasured - **Treatment**: Binary vs continuous - **Mediator**: Continuous vs binary vs count ## Parameter Grid Design ### Sample Size Selection | Size | Label | Purpose | |------|-------|---------| | 100-200 | Small | Stress test | | 500 | Medium | Typical study | | 1000-2000 | Large | Asymptotic behavior | | 5000+ | Very large | Efficiency comparison | ### Effect Size Selection | Effect | Interpretation | |--------|----------------| | 0 | Null (Type I error) | | 0.1 | Small | | 0.3 | Medium | | 0.5 | Large | ### Recommended Grid Structure ```r params <- expand.grid( n = c(200, 500, 1000, 2000), effect = c(0, 0.14, 0.39, 0.59), # Small/medium/large per Cohen confounding = c(0, 0.3, 0.6), misspecification = c(FALSE, TRUE) ) ``` ## Performance Metrics ### Primary Metrics | Metric | Formula | Target | MCSE Formula | |--------|---------|--------|--------------| | Bias | $\bar{\hat\psi} - \psi_0$ | ≈ 0 | $\sqrt{\text{Var}(\hat\psi)/n_{sim}}$ | | Empirical SE | $\text{SD}(\hat\psi)$ | — | Complex | | Average SE | $\bar{\widehat{SE}}$ | ≈ Emp SE | $\text{SD}(\widehat{SE})/\sqrt{n_{sim}}$ | | Coverage | $\frac{1}{n_{sim}}\sum I(\psi_0 \in CI)$ | ≈ 0.95 | $\sqrt{p(1-p)/n_{sim}}$ | | MSE | $\text{Bias}^2 + \text{Var}$ | Minimize | — | | Power | % rejecting $H_0$ | Context-dependent | $\sqrt{p(1-p)/n_{sim}}$ | ### Relative Metrics (for method comparison) - Relative bias: $\text{Bias}/\psi_0$ (when $\psi_0 \neq 0$) - Relative efficiency: $\text{Var}(\hat\psi_1)/\text{Var}(\hat\psi_2)$ - Relative MSE: $\text{MSE}_1/\text{MSE}_2$ ## Replication Guidelines ### Minimum Replications | Metric | Minimum | Recommended | |--------|---------|-------------| | Bias | 1000 | 2000 | | Coverage | 2000 | 5000 | | Power | 1000 | 2000 | ### Monte Carlo Standard Error Always report MCSE for key metrics: - Coverage MCSE at 95%: $\sqrt{0.95 \times 0.05 / n_{sim}} \approx 0.007$ for $n_{sim}=1000$ - Need ~2500 reps for MCSE < 0.005 ## R Implementation Template ```r #' Run simulation study #' @param scenario Parameter list for this scenario #' @param n_rep Number of replications #' @param seed Random seed run_simulation <- function(scenario, n_rep = 2000, seed = 42) { set.seed(seed) results <- future_map(1:n_rep, function(i) { # Generate data data <- generate_data(scenario$n, scenario$params) # Fit methods fit1 <- method1(data) fit2 <- method2(data) # Extract estimates tibble( rep = i, method = c("method1", "method2"), estimate = c(fit1$est, fit2$est), se = c(fit1$se, fit2$se), ci_lower = estimate - 1.96 * se, ci_upper = estimate + 1.96 * se ) }, .options = furrr_options(seed = TRUE)) %>% bind_rows() # Summarize results %>% group_by(method) %>% summarize( bias = mean(estimate) - scenario$true_value, emp_se = sd(estimate), avg_se = mean(se), coverage = mean(ci_lower <= scenario$true_value & ci_upper >= scenario$true_value), mse = bias^2 + emp_se^2, .groups = "drop" ) } ``` ## Results Presentation ### Standard Table Format ``` Table X: Simulation Results (n_rep = 2000) Method 1 Method 2 n Bias SE Cov MSE Bias SE Cov MSE ----------------------------------------------------------- 200 0.02 0.15 0.94 0.023 0.01 0.12 0.95 0.014 500 0.01 0.09 0.95 0.008 0.00 0.08 0.95 0.006 1000 0.00 0.06 0.95 0.004 0.00 0.05 0.95 0.003 Note: Cov = 95% CI coverage. MCSE for coverage ≈ 0.005. ``` ### Visualization Guidelines - Use faceted plots for multiple scenarios - Show confidence bands for metrics - Compare methods side-by-side - Log scale for MSE if range is large ## Checkpoints and Reproducibility ### Checkpointing Strategy ```r # Save results incrementally if (i %% 100 == 0) { saveRDS(results_so_far, file = sprintf("checkpoint_%s_rep%d.rds", scenario_id, i)) } ``` ### Reproducibility Requirements 1. Set seed explicitly 2. Record package versions (`sessionInfo()`) 3. Use `furrr_options(seed = TRUE)` for parallel 4. Save full results, not just summaries 5. Document any manual interventions ## Common Pitfalls ### Design Pitfalls - Too few replications for coverage assessment - Unrealistic parameter combinations - Missing null scenario (effect = 0) - No misspecification scenarios ### Implementation Pitfalls - Not setting seeds properly in parallel - Ignoring convergence failures - Not checking for numerical issues - Insufficient burn-in for MCMC methods ### Reporting Pitfalls - Missing MCSE - Not reporting convergence failures - Cherry-picking scenarios - Inadequate description of DGP ## Key References - Morris et al 2019 - Burton et al - White