Code
mcse_prop <- function(p, R) sqrt(p * (1 - p) / R)
R_needed <- function(p, half_width) ceiling(p * (1 - p) / (half_width / 2)^2)
mcse_prop(0.05, R = 1000)[1] 0.006892024Code
R_needed(0.05, half_width = 0.01)[1] 1900Simulation studies for methods research
ADEMP, Monte Carlo error, and the protocol you write before running anything.
Before reading: Bias, standard errors, tests, and enough R to write functions. Evaluating a statistical method defines the performance summaries. The executable network example uses Association targets, pooling, and imperfect adjustment and Jackknife pseudo-values; the protocol and Monte Carlo precision sections can be read before that example.
1 Why simulation matters
A simulation study can do more than illustrate a theory after the fact. It checks the theory, shows what matters in practice, and often reveals the problem in the first place. Writing down a data-generating process forces you to say exactly what you believe about the data, and the moment you cannot write it down is the moment you have found the gap in your understanding.
A simulation study is also an experiment, and it should be designed like one. The framework below is ADEMP, from Morris, White, and Crowther (2019). Write the protocol before you run anything, because a protocol written afterwards only describes the study and cannot protect against choices made after seeing the results.
One replicate is one complete synthetic study: generate the people and all their features, apply each analysis, and record its results. Repeating an analysis 1,000 times on one fixed dataset is a different experiment. So is generating 1,000 genes for one cohort. The replicate must reproduce the source of variation relevant to the performance claim.
2 Connect the evaluation to the estimand
For every scenario, record the scientific target \theta(F) and, when derivable, the procedure’s probability limit \theta^*(F). Bias that persists as sample size grows can indicate \theta^*(F)\ne\theta(F); poor small-sample coverage can also arise from a variance approximation, nuisance estimation, or an implementation error. Vary sample size and the specific misspecification separately to distinguish these mechanisms.
For a cross-fitted estimator, vary nuisance complexity and the independent-unit sample size, and keep all repeated observations of a subject in the same fold. Compare estimated standard errors with empirical standard deviations, not just point-estimate bias. Record fitting failures, Monte Carlo uncertainty, and computation as separate outcomes.
Keep the protocol and replicate outputs in the research record described in Research practice and reproducibility. Saving only a final summary table makes later checks and paired method comparisons harder.
3 ADEMP
A simulation protocol has five parts, abbreviated ADEMP.
The aims say what question the study answers. A good aim is a question with a yes-or-no or numeric answer, such as “Does method A hold its level when groups are unbalanced?” An aim like “compare methods” is too vague to answer.
The data-generating mechanisms are the processes that produce the simulated datasets, with every parameter named. The protocol says which factors vary (sample size, imbalance, effect size, correlation structure) and which are fixed, and why. A factorial design can separate factors and their interactions; choose the grid for the scientific aim and computational budget, and state any deliberate omissions. At least one mechanism should have a known answer (a null, or a case where the method is known to work), so that a bug shows up as a wrong answer there.
The estimands are what each method is trying to estimate or test. This is where the What is your estimator actually estimating? page applies. The estimand and the limit of a method often differ, and the simulation can measure that difference, so the estimand has to be defined independently of any method.
The methods are the procedures being compared, specified precisely enough to reimplement: the test statistic, the reference distribution, the number of permutations, the multiple-testing adjustment, and the software and version.
The performance measures are what you compute from each replicate and how you summarize over replicates. For tests, they are the rejection rate under the null (level) and under alternatives (power). For estimators, they are bias, empirical SD, mean squared error, and coverage of a confidence interval. Every performance measure is reported with its Monte Carlo standard error.
3.1 Make the protocol concrete before coding
The running example asks whether methods falsely report a between-group correlation difference when only the group means change. Its scientific target is
\begin{aligned} \Delta_\rho&=\rho_{\rm patient}-\rho_{\rm control},\\ H_0&:\Delta_\rho=0. \end{aligned}
Each group has unit feature variances and correlation 0.5; the patient means shift together. Group sizes are 20/80 or 50/50, always totaling 100 people. Shifts of 0, 1, and 2 give six scenarios. Every scenario is a scientific null for the correlation contrast, including those with very different expression means. Rejection rates in those scenarios assess false-positive behavior for this stated scientific claim, rather than power to detect mean differences.
We compare two candidate tests based on cohort-derived pseudo-values with a direct Fisher-z correlation comparison. The candidate scores are computed with and without within-group centering. A pseudo-value compares the pooled correlation with the correlation after deleting one person; it is not a raw independent measurement of that person’s correlation. This is the target-mismatch example developed in Jackknife pseudo-values.
The direct test compares the two within-group sample correlations after the Fisher transformation. Its approximate normal reference is justified asymptotically under the generator’s independent bivariate-normal groups. It is a useful benchmark for this mechanism, with its own finite-sample approximation. The two score t-tests require a separate argument: centering alone does not validate t inference for shared-cohort scores. The evaluation is designed to investigate their behavior, rather than to assume all three tests are calibrated.
Before running, predict the mechanisms. A mean shift can create a difference in raw pooled scores despite equal within-group correlations; imbalance can amplify the discrepancy. Centering removes that particular mean-mixture mechanism, but need not solve the score variance problem. A direct within-group correlation comparison should be unaffected by adding constants to a group’s features. These predictions give the numerical results something specific to confirm or challenge.
4 Monte Carlo error, and how many replicates you need
A simulation estimates a performance measure from R replicates, so it has sampling error of its own. Reporting a level of 0.052 without saying how many replicates produced it is reporting a number without a standard error.
Distinguish n, the independent people in each dataset; R, the independent simulated datasets; and B, any bootstrap or permutation draws used inside one analysis. Increasing n changes the method’s statistical problem. Increasing R improves precision of the evaluation. Increasing B refines an internal resampling calculation, subject to its validity. None is a substitute for the others.
Let J_r=1 if the method rejects on dataset r and zero otherwise. Its rejection probability is p=E[J_r]. The simulated rate is \hat p=R^{-1}\sum_rJ_r. Independence between complete datasets gives \operatorname{Var}(\hat p)=p(1-p)/R, motivating the estimated MCSE
\operatorname{MCSE}(\hat p) = \sqrt{\frac{\hat p(1 - \hat p)}{R}}.
To estimate a nominal 5% level to within about \pm 0.01 (two MCSEs), you need 0.05 \times 0.95 / 0.005^2 \approx 1900 replicates. To within \pm 0.005, about 7600. Power near 0.5 has the largest MCSE, 0.5/\sqrt{R}. Choose the number of replicates from the precision needed for the claim, not a universal minimum; smaller pilot studies can diagnose problems before a final evaluation.
More generally, for an approximate two-MCSE half-width h, plan R\approx4p(1-p)/h^2. This uses a planning value of p, not a promise that the true rejection rate equals the nominal level. If p is poorly known, p(1-p)\le1/4 gives a conservative precision calculation. A 0.052 rate from 1,000 datasets has MCSE about 0.007: a difference of 0.002 from 0.05 is small relative to that simulation uncertainty.
The same calculation starts from a dataset-level contribution for other summaries, with independent replicates and finite variance of that contribution. To estimate bias, average \hat\theta_r-\theta and use its sample SD divided by \sqrt R. For MSE use (\hat\theta_r-\theta)^2. For FDR use the entire dataset’s FDP. FDP is a fraction, so its MCSE is generally not the binomial formula above. To compare paired methods, average their within-dataset contribution differences and use the SD of those differences divided by \sqrt R, as worked through in Evaluating a statistical method. An empirical SD is itself a nonlinear summary; it needs a variance formula or a justified resampling calculation, rather than treating it as an average of estimates.
At a boundary the plug-in binomial MCSE can be misleading. Zero rejections in 200 datasets gives an estimated MCSE of zero, but does not establish a zero rejection probability. For independent Bernoulli replicates, an exact one-sided 95% upper bound after zero events is 1-0.05^{1/R}, about 0.015 when R=200. This follows from \Pr(\text{zero events})=(1-p)^R. It concerns binary rejection events; do not automatically reuse it for a general FDP.
A binomial band is the cleanest way to show level results. If a method truly has level 0.05, the observed rejection rate over R replicates lies between qbinom(0.025, R, 0.05) / R and qbinom(0.975, R, 0.05) / R with at least about 95% probability under that point null, allowing for binomial discreteness. Draw a reference band on level plots. Falling inside it is not proof of calibration, and interpretation across many cells needs to account for their joint examination. Across many calibrated cells, some points can fall outside by chance; discreteness and dependence affect the pattern. A collection of marginal bands is not a simultaneous calibration test.
Code
R <- 2000
c(lower = qbinom(0.025, R, 0.05) / R, upper = qbinom(0.975, R, 0.05) / R) lower upper
0.0405 0.0600 5 A minimal, correct skeleton
The code separates generation, analysis, and summaries. One function generates a dataset. A named list defines the methods. A runner applies every method to the same dataset and retains one row per method and replicate, including failures. A separate step summarizes those rows. The data generator is the two-group mechanism used in Association targets, pooling, and imperfect adjustment. Here it is used to evaluate tests of a within-group contrast.
Code
# Data-generating mechanism: two groups, same within-group correlation,
# patients shifted by `shift` SDs on both genes. Bivariate normal within group.
SimulateCohort <- function(n_ctrl, n_pat, shift, rho = 0.5) {
group <- rep(c("control", "patient"), c(n_ctrl, n_pat))
n <- length(group)
z1 <- rnorm(n)
z2 <- rho * z1 + sqrt(1 - rho^2) * rnorm(n)
data.frame(group = group,
gene_a = z1 + shift * (group == "patient"),
gene_b = z2 + shift * (group == "patient"))
}
# Per-person LIONESS-style scores for the correlation of x and y.
PersonScores <- function(x, y) {
n <- length(x)
r_all <- cor(x, y)
sapply(seq_len(n), function(i) n * r_all - (n - 1) * cor(x[-i], y[-i]))
}
# Each function returns one p-value; truth is not passed to the methods.
method_functions <- list(
"Raw scores" = function(d) {
score <- PersonScores(d$gene_a, d$gene_b)
t.test(score ~ group, data = d)$p.value
},
"Centered scores" = function(d) {
a <- d$gene_a - ave(d$gene_a, d$group)
b <- d$gene_b - ave(d$gene_b, d$group)
score <- PersonScores(a, b)
t.test(score ~ group, data = d)$p.value
},
"Fisher z" = function(d) {
by_group <- split(d, d$group)
n <- vapply(by_group, nrow, integer(1))
r <- vapply(by_group, function(a) cor(a$gene_a, a$gene_b), numeric(1))
stopifnot(all(n > 3), all(abs(r) < 1))
statistic <- diff(atanh(r)) / sqrt(sum(1 / (n - 3)))
2 * pnorm(-abs(statistic))
}
)
ApplyMethods <- function(d, method_set) {
do.call(rbind, lapply(names(method_set), function(name) {
error_message <- ""; warnings <- character()
elapsed <- system.time({
p <- tryCatch(withCallingHandlers({
p <- method_set[[name]](d)
stopifnot(length(p) == 1, is.finite(p), p >= 0, p <= 1)
unname(p)
}, warning = function(w) {
warnings <<- c(warnings, conditionMessage(w))
invokeRestart("muffleWarning")
}), error = function(e) {
error_message <<- conditionMessage(e)
NA_real_
})
}, gcFirst = FALSE)[["elapsed"]]
data.frame(method = name, p_value = p, success = !nzchar(error_message),
stage = "analysis", error = error_message,
warning = paste(warnings, collapse = " | "), elapsed)
}))
}
# One scenario: preserve every method attempt, even when another method fails.
RunCell <- function(cell, n_ctrl, n_pat, shift, R, seed,
method_set = method_functions) {
set.seed(seed)
rows <- lapply(seq_len(R), function(b) {
d <- tryCatch(SimulateCohort(n_ctrl, n_pat, shift), error = identity)
if (inherits(d, "error")) {
value <- data.frame(method = names(method_set), p_value = NA_real_,
success = FALSE, stage = "generation", error = conditionMessage(d),
warning = "", elapsed = NA_real_)
} else {
value <- ApplyMethods(d, method_set)
}
data.frame(cell, n_ctrl, n_pat, shift, replicate = b, value)
})
do.call(rbind, rows)
}The method names label actual functions, not truth labels. A failing method produces a failed row while the others still run; a failed generator marks every method for that replicate as unavailable. Warnings are retained even when a finite p-value is returned, so success means computational completion under these checks, not verified statistical validity. This is a serial teaching runner. Production work should also retain the generator state or datasets needed to reproduce a specific attempt.
The design is then a table of cells, and the study is a loop over its rows. Give each cell a stable identifier and recorded seed so it can be rerun. Distinct integer seeds alone do not prove nonoverlapping or independent streams; parallel execution needs a justified stream scheme. Deliberate common random numbers can also be useful, provided the pairing is recorded.
Code
design <- expand.grid(n_ctrl = c(20, 50), n_pat = c(80, 50), shift = c(0, 1, 2))
design <- design[design$n_ctrl + design$n_pat == 100, ] # keep total N fixed
design$cell <- seq_len(nrow(design))
design$seed <- 1000 + seq_len(nrow(design))
design n_ctrl n_pat shift cell seed
1 20 80 0 1 1001
4 50 50 0 2 1002
5 20 80 1 3 1003
8 50 50 1 4 1004
9 20 80 2 5 1005
12 50 50 2 6 1006Code
# Production run: R = 2000 per cell. Save per-replicate results, never summaries only.
results <- do.call(rbind, lapply(seq_len(nrow(design)), function(i) {
with(design[i, ], RunCell(cell, n_ctrl, n_pat, shift, R = 2000, seed = seed))
}))
dir.create("results", showWarnings = FALSE)
saveRDS(results, "results/level_study.rds")Code
# Teaching run at R = 200. Report the MCSE with it.
results <- do.call(rbind, lapply(seq_len(nrow(design)), function(i) {
with(design[i, ], RunCell(cell, n_ctrl, n_pat, shift, R = 200, seed = seed))
}))
cells <- split(results, interaction(results$cell, results$method, drop = TRUE))
summ <- do.call(rbind, lapply(cells, function(d) {
valid <- d[d$success, ]; M <- nrow(valid)
rate <- if (M) mean(valid$p_value <= 0.05) else NA_real_
data.frame(d[1, c("cell", "n_ctrl", "n_pat", "shift", "method")],
attempted = nrow(d), successful = M, failures = sum(!d$success),
rejection = rate, MCSE = if (M) sqrt(rate*(1-rate)/M) else NA_real_)
}))
knitr::kable(summ[c("n_ctrl", "n_pat", "shift", "method", "successful",
"failures", "rejection", "MCSE")], digits = 3, row.names = FALSE)| n_ctrl | n_pat | shift | method | successful | failures | rejection | MCSE |
|---|---|---|---|---|---|---|---|
| 20 | 80 | 0 | Centered scores | 200 | 0 | 0.080 | 0.019 |
| 50 | 50 | 0 | Centered scores | 200 | 0 | 0.050 | 0.015 |
| 20 | 80 | 1 | Centered scores | 200 | 0 | 0.085 | 0.020 |
| 50 | 50 | 1 | Centered scores | 200 | 0 | 0.065 | 0.017 |
| 20 | 80 | 2 | Centered scores | 200 | 0 | 0.060 | 0.017 |
| 50 | 50 | 2 | Centered scores | 200 | 0 | 0.055 | 0.016 |
| 20 | 80 | 0 | Fisher z | 200 | 0 | 0.060 | 0.017 |
| 50 | 50 | 0 | Fisher z | 200 | 0 | 0.030 | 0.012 |
| 20 | 80 | 1 | Fisher z | 200 | 0 | 0.060 | 0.017 |
| 50 | 50 | 1 | Fisher z | 200 | 0 | 0.055 | 0.016 |
| 20 | 80 | 2 | Fisher z | 200 | 0 | 0.040 | 0.014 |
| 50 | 50 | 2 | Fisher z | 200 | 0 | 0.055 | 0.016 |
| 20 | 80 | 0 | Raw scores | 200 | 0 | 0.070 | 0.018 |
| 50 | 50 | 0 | Raw scores | 200 | 0 | 0.050 | 0.015 |
| 20 | 80 | 1 | Raw scores | 200 | 0 | 0.200 | 0.028 |
| 50 | 50 | 1 | Raw scores | 200 | 0 | 0.035 | 0.013 |
| 20 | 80 | 2 | Raw scores | 200 | 0 | 0.835 | 0.026 |
| 50 | 50 | 2 | Raw scores | 200 | 0 | 0.015 | 0.009 |
The displayed run is a 200-replicate teaching evaluation; the 2,000-replicate code above is a proposed larger run and is not executed when rendering. Read the actual rates with their MCSEs, rather than treating an expected pattern as a numerical result. Mean shifts can make the raw scores answer a different question, especially with unbalanced groups. Centered scores can behave differently in this generator, but a rate inside a reference band does not prove their calibration elsewhere.
The summary reports method-specific rates conditional on computational success, with their actual denominators, and retains every attempt in results. If failures occur, investigate them and label comparisons on jointly successful replicates explicitly. Near level 0.05, a 200-replicate rate has MCSE about 0.015. Use a prespecified larger run when the scientific claim needs greater precision.
Read the plot against the scientific null: shifting means never changes either within-group correlation in this generator. A rising rejection curve therefore cannot be called power for a correlation change. Balance and shift are crossed to reveal their interaction, rather than changing sample size at the same time. The direct comparison helps distinguish a problem with the score construction from a general inability to estimate correlations in these samples. Agreement with a reference band in six settings is still a finite evaluation, not a theorem covering every distribution.
For the 20/80 split with shift 2, the raw-score test rejects 167 of 200 datasets: 0.835 with MCSE about 0.026. The direct test rejects 8 of 200: 0.040 with MCSE about 0.014. The scientific null is true for both. In the balanced setting the raw-score artifact can cancel in the between-group comparison, even though the pooled scores still do not become individual correlation measurements. This is why a single balanced null scenario would be an incomplete evaluation.
The functions here use no analysis randomness. If a new method uses random folds, bootstrap draws, or stochastic fitting, give those operations separate recorded random-number streams. Otherwise adding that method can consume random draws and change later generated datasets even though the cell seed is unchanged. For parallel work, use a stream scheme such as the one provided by R’s parallel RNG tools, rather than assuming nearby integer seeds establish independent streams.
6 Rules for a simulation study
- Write the protocol first, as a file in the project repository (for example
protocol.md) in ADEMP form, before the first production run. Pilot runs to debug the code are fine. If the design changes after seeing production results, record the change in the protocol with a date and a reason. - Give every cell a seed derived from the design table, so that the production run is reproducible from the repository alone.
- Save the per-replicate results to a results folder. Summaries are computed from them when the manuscript renders.
- Report MCSEs for bias and other numerical performance summaries. For rejection and coverage rates across independent datasets, use a binomial MCSE or an equivalent uncertainty band.
- Include a cell you can check: a null with no heterogeneity, or a case with a known closed form. If that cell is wrong, stop.
- Keep the total N fixed when varying balance. Otherwise a balance effect is confounded with a sample size effect.
- Record the environment by saving
sessionInfo()with the production run. - For factorial designs, draw one plot panel per combination of the two most important factors, with methods distinguished by color and the appropriate Monte Carlo uncertainty shown. Use tables for the numbers and plots for the pattern.
7 Common failures
- Making precise calibration claims from a small pilot. At level 0.05, 100–200 replicates have MCSE about 0.015–0.022; use the precision needed for the claim rather than banning pilot runs.
- Varying several things at once between two cells so no single factor can be blamed.
- Simulating only from the working model the method assumes. Many studies exist to see what happens when that assumption fails, so the data-generating mechanisms must include a richer model with the departure (see the What is your estimator actually estimating? page).
- Giving only one procedure favorable generating assumptions, without justifying those assumptions from the research question. A scenario can deliberately match a method’s model, but include relevant departures and competing-model settings. If a generator is fitted to motivating data, state that evaluation is conditional on that fitted law and what features it omits.
- Confusing a permutation p-value’s grid spacing 1/(B+1) with precision of a rejection-rate estimate across independent datasets. The finite grid can affect the test’s attainable level; replicate count controls MCSE of that actual rejection rate.
- Summaries only. When a result looks odd, the per-replicate results are the only way to find out why.
8 Exercises
Checkpoints are completion targets for this page. Additional practice below is optional. Each checkpoint gives a completion standard, a hint, and a solution outline.
Checkpoint · SIM-01 — Write a target-based simulation protocol
Task type: Design and interpretation.
Write an ADEMP protocol for the study in this page. Specify the within-group edge target and its null, the balance/mean-shift factors, both score-based methods and the direct correlation benchmark, rejection rates and failure reporting, and the Monte Carlo precision that determines the replicate count. Explain why an inflated rejection rate can reveal a target mismatch rather than evidence of a true within-group edge difference. Distinguish increasing n, R, and an internal resampling count B.
Ready to continue when: another reader could implement your protocol and distinguish the scientific null, the procedure’s target, and simulation uncertainty.
Transfer: What new mechanism would test power against an actual within-group correlation difference while keeping the target definition unchanged?
Hint
The data generator holds each within-group correlation at 0.5, even when pooled means shift. Near a rejection rate of 0.05, MCSE is about \sqrt{0.05(0.95)/R}.
Solution outline — open after trying
One protocol asks whether the methods reject the null of equal within-group correlations too often. Use the stated bivariate-normal mechanisms with both correlations 0.5, total sample size 100, group sizes 20/80 or 50/50, and shifts 0, 1, or 2. Compare the raw and within-group-centered score tests and the direct Fisher-z test on the same simulated datasets. Report rejection rates, their MCSEs, and method-specific failures. To target a Monte Carlo half-width of roughly 0.01 at two MCSEs near 0.05, use about 1,900 or more replicates; 2,000 is a reasonable prespecified choice. Neither score test’s calibration should be presumed from its label. A mean-shift artifact can generate apparent score differences while the within-group correlation null remains true. The count n sets each study’s independent-unit sample size; R sets precision of the performance evaluation; B sets an internal resampling calculation if a method uses one.
Additional practice (optional)
These exercises extend the ideas through further derivations, simulations, or applications. Choose those useful for your work.
- Add a method to
method_functionswithout changingSimulateCohort()orRunCell(). If you try a Wilcoxon test on the scores, state its own null and explain why a rank test does not automatically fix a mismatch with within-group correlation. Check that one method’s deliberate computational failure is retained while the others still run. - Run one cell at R = 200, 1000, and 5000, and plot the estimated level with its MCSE at each. Watch the interval shrink and decide what R you would use in a paper.
- Suppose no rejections occur in 200 replicates. Derive the one-sided upper bound from (1-p)^R=0.05, then compare it with the plug-in MCSE. How many zero-event replicates would be needed to put that upper bound below 0.005?
9 Reading
- Morris, T. P., White, I. R., and Crowther, M. J. (2019). Using simulation studies to evaluate statistical methods. Statistics in Medicine, 38(11), 2074 to 2102. Practical guidance on protocol design, performance measures, and Monte Carlo uncertainty.
- Burton, A., Altman, D. G., Royston, P., and Holder, R. L. (2006). The design of simulation studies in medical statistics. Statistics in Medicine, 25(24), 4279 to 4292.
- Siepe, B. S., Bartoš, F., Morris, T. P., Boulesteix, A.-L., Heck, D. W., and Pawel, S. (2024). Simulation studies for methodological research in psychology: A standardized template for planning, preregistration, and reporting. Psychological Methods. doi:10.1037/met0000695. A template for writing a simulation protocol, useful as an example of the form in use.