---
title: "Evaluating a statistical method"
subtitle: "Targets, limits, uncertainty, and computational reliability"
guide-status: "Available"
---
**Before reading:** Read [What is your estimator actually estimating?](../questions-and-design/estimands-and-misspecification.qmd); basic means, variance, and confidence intervals are assumed.
## Evaluate a claim, not a method name
A method can have a small standard error and still answer the wrong question. Before comparing procedures, write the population target $\psi(P)$, the limit of each procedure, and the claim being evaluated: consistency, confidence-interval coverage, efficiency, robustness, prediction, or computational reliability. Those claims require different evidence.
For an omics example, a pooled coefficient, an adjusted group contrast, and a standardized population contrast may differ even with infinitely many participants. First use [What is your estimator actually estimating?](../questions-and-design/estimands-and-misspecification.qmd) to determine whether a comparison has a common target. Then assess finite-sample behavior.
There are three sources of discrepancy to keep apart. **Target mismatch** means the procedure approaches a different population quantity. **Sampling error** is variation caused by which people enter a finite dataset. **Numerical error** arises when the implementation cannot compute its intended answer accurately. Increasing the sample size, increasing simulation replicates, and tightening an optimizer's tolerance address different problems.
## Choose measures that answer the claim
Suppose a scenario has fixed population truth $\theta$ and produces R independent datasets. For replicate r, retain the estimate $\hat\theta_r$, reported SE $\hat s_r$, interval $I_r$, and fitting status. Then the main summaries are
$$
\begin{aligned}
\widehat{\operatorname{Bias}}&=\overline{\hat\theta}-\theta,\\
s_{\rm emp}^2&=\frac1{R-1}\sum_r
(\hat\theta_r-\overline{\hat\theta})^2,\\
\overline{\operatorname{SE}}&=\frac1R\sum_r\hat s_r,\\
\widehat{\operatorname{MSE}}&=\frac1R\sum_r(\hat\theta_r-\theta)^2,\\
\widehat{\operatorname{Coverage}}&=\frac1R\sum_r1\{\theta\in I_r\}.
\end{aligned}
$$
Bias measures the average signed error; positive and negative errors can cancel. MSE averages squared errors, penalizing both bias and variation. Empirical SD measures how estimates vary **between datasets**. Mean SE summarizes what the method reports **within datasets**. Coverage asks whether the intervals actually include the common target at the advertised frequency. Comparing interval length is meaningful alongside coverage: very short intervals are easy to produce if they miss the target.
For testing, rejection under a true scientific null estimates Type I error; rejection under a specified alternative estimates power. For a family of omics hypotheses, estimate FDR by averaging one false discovery proportion per dataset. Genes and edges from the same dataset are not independent simulation replicates. See [Multiple testing when the tests are not independent](../analysis-problems/multiple-testing.qmd) for that distinction.
For the formulas above assume all fits return usable results; the failure section below explains the denominators when they do not. If truth changes randomly between replicates, also state whether the summaries average over that random-truth distribution or concern one fixed parameter. Those are different evaluation claims.
## A precise example: shrinkage toward zero
Let independent $Y_i\sim N(\mu,1)$ and suppose the target is $\mu$. Compare the sample mean with $\tilde\mu=0.9\bar Y$. Their limits are $\mu$ and $0.9\mu$. For the second procedure,
$$
\begin{aligned}
\operatorname{Bias}(\tilde\mu)&=-0.1\mu,\\
\operatorname{Var}(\tilde\mu)&=0.81/n.
\end{aligned}
$$
Its sampling SD is smaller: multiplying the sample mean by 0.9 multiplies its variance by $0.9^2$. But
$$
\tilde\mu-\mu=0.9(\bar Y-\mu)-0.1\mu.
$$
The first term is sampling noise that shrinks with n; the second is a fixed displacement from the target. A smaller SE can describe the first term correctly while omitting the second completely. The bias measured in units of sampling SD is $-\mu\sqrt n/9$, whose magnitude grows when $\mu\ne0$.
Let $c_n$ be the 97.5th percentile of $t_{n-1}$. Compare the usual mean interval with the naively shrunken interval,
$$
\begin{aligned}
&\bar Y\pm c_n\frac{s}{\sqrt n},\\
&\tilde\mu\pm c_n\frac{0.9s}{\sqrt n}.
\end{aligned}
$$
Under the normal generator, the first interval has exact 95% coverage for $\mu$. The second is the first interval multiplied by 0.9, so it covers $0.9\mu$ at 95%, but generally does not cover $\mu$ at that rate. At $\mu=0$ those targets coincide. At $\mu=1$, the interval shrinks around 0.9 while the target stays at 1. Dividing the shrunken estimate and both endpoints by 0.9 recovers ordinary mean inference; simply shrinking the center and SE does not provide that correction.
The MSE decomposition makes a different tradeoff:
$$
\operatorname{MSE}(\tilde\mu)=0.01\mu^2+0.81/n.
$$
At $\mu=1,n=10$ this is 0.091, less than the sample mean's 0.1. That point-estimation advantage does not establish consistency or interval coverage. For a nonzero fixed $\mu$, the MSE advantage disappears once $n\ge19/\mu^2$. Always say which loss or inferential property makes one procedure preferable.
### Consistency alone is not enough for inference
Now let the multiplier approach one: $\tilde\mu_n=(1-n^{-1/4})\bar Y$. Both the multiplier and the sample mean converge to their desired limits, so this estimator is consistent. Its bias is nevertheless $-\mu n^{-1/4}$, larger than its root-n sampling SD. At a nonzero target, ordinary centered intervals can still fail increasingly badly. A multiplier $1-n^{-1}$ instead has bias negligible at the root-n scale. This is the same reason a semiparametric inference argument controls the **rate** of a nuisance remainder, rather than asking only whether it converges to zero.
## Compare procedures on the same datasets
The following experiment uses the same simulated sample for both procedures in each replicate. This keeps the comparison tied to the same scenarios and permits a paired analysis. Pairing reduces Monte Carlo uncertainty when the performance contributions have positive covariance; it is not a universal variance-reduction theorem for any two algorithms. All values are synthetic, with truth supplied by the generator.
```{r}
set.seed(61)
n_sim <- 1000
records <- list(); i <- 0
for (mu in c(0, 1)) for (n in c(20, 100, 500)) {
y <- matrix(rnorm(n * n_sim, mean = mu), nrow = n)
average <- colMeans(y)
sample_sd <- sqrt(colSums(sweep(y, 2, average)^2) / (n - 1))
for (shrink in c(1, 0.9)) {
estimate <- shrink * average
se <- shrink * sample_sd / sqrt(n)
covered <- abs(estimate - mu) <= qt(0.975, df = n - 1) * se
i <- i + 1
records[[i]] <- data.frame(mu = mu, n = n, replicate = seq_len(n_sim),
method = if (shrink == 1) "Mean" else "0.9 x mean",
estimate, se, covered, squared_error = (estimate - mu)^2)
}
}
records <- do.call(rbind, records)
cells <- split(records, interaction(records$mu, records$n, records$method))
performance <- do.call(rbind, lapply(cells, function(d) {
R <- nrow(d); coverage <- mean(d$covered)
data.frame(d[1, c("mu", "n", "method")],
bias = mean(d$estimate) - d$mu[1], empirical_sd = sd(d$estimate),
mean_se = mean(d$se), MSE = mean(d$squared_error), coverage,
bias_mcse = sd(d$estimate) / sqrt(R),
coverage_mcse = sqrt(coverage * (1 - coverage) / R))
}))
print(performance, digits = 3, row.names = FALSE)
```
Each row in `records` belongs to one method on one replicate; the replicate ID is shared across methods within the same scenario. Summaries are computed afterwards. **MCSE** describes uncertainty in such a summary because the simulation uses finitely many datasets. Here bias MCSE is empirical SD divided by $\sqrt R$, whereas the SE reported by a method scales with $\sqrt n$. Increasing R improves knowledge of the method's behavior at that n; it does not improve the method's interval for an individual dataset.
For example, an empirical SD of 0.10 across 1,000 datasets gives a bias MCSE of about 0.0032. It does not give each dataset an estimation SE of 0.0032. This distinction matters when reporting small bias or coverage differences.
```{r}
#| label: fig-limit-coverage
#| fig-cap: "A smaller sampling SD does not remove target bias. At mu = 1 the shrunken procedure's coverage deteriorates with n; the null scenario alone would conceal that problem. Bars show plus/minus two Monte Carlo standard errors."
#| fig-alt: "Two coverage panels, for population means zero and one, comparing ordinary and shrunken means as sample size increases."
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
for (target in c(0, 1)) {
plot(NA, xlim = c(15, 650), ylim = c(0.2, 1), log = "x",
xlab = "Independent observations n", ylab = "Coverage of mu",
main = paste("mu =", target))
abline(h = 0.95, lty = 2, col = "gray50")
for (j in 1:2) {
method_name <- c("Mean", "0.9 x mean")[j]
a <- performance[performance$mu == target & performance$method == method_name, ]
a <- a[order(a$n), ]
color <- c("#2F80C8", "#D55E00")[j]
lines(a$n, a$coverage, type = "b", pch = j + 14, col = color)
arrows(a$n, a$coverage - 2 * a$coverage_mcse,
a$n, a$coverage + 2 * a$coverage_mcse,
angle = 90, code = 3, length = 0.04, col = color)
}
legend("bottomleft", c("Mean", "0.9 x mean"), col = c("#2F80C8", "#D55E00"),
pch = 15:16, lty = 1, bty = "n", cex = 0.8)
}
```
At $\mu=0$, both methods have exactly the same coverage indicators in this construction: multiplying the entire interval by a positive constant does not change whether it includes zero. That reassuring null scenario conceals the failure at nonzero targets. In the right panel the shrunken estimator's SE tracks its sampling SD, yet its intervals increasingly miss $\mu=1$.
### Calculate uncertainty in a paired difference
Let $C_{A,r}$ and $C_{B,r}$ indicate coverage by methods A and B on dataset r. Form $d_r=C_{B,r}-C_{A,r}$, which is −1, 0, or 1. Then
$$
\begin{aligned}
\widehat\Delta_{\rm cov}&=\bar d,\\
\widehat{\operatorname{MCSE}}(\bar d)&=s_d/\sqrt R.
\end{aligned}
$$
A value −1 means B misses where A covers; zero means they agree; 1 means B alone covers. The variance is based on these **within-dataset differences**. It includes the methods' covariance, which would be lost by treating the marginal coverage estimates as independent.
```{r}
a <- subset(records, mu == 1 & n == 500 & method == "Mean")
b <- subset(records, mu == 1 & n == 500 & method == "0.9 x mean")
stopifnot(identical(a$replicate, b$replicate))
d <- as.integer(b$covered) - as.integer(a$covered)
c(coverage_difference = mean(d), paired_MCSE = sd(d) / sqrt(length(d)))
```
The same recipe works for differences in squared error, interval length, runtime, or dataset-level FDP: form the relevant pair of contributions, subtract, then summarize across independent datasets. Different performance measures require their own contributions, not a single universal binomial SE.
## Failures are part of performance
Suppose an algorithm returns valid intervals in 90 of 100 replicates, and 86 of those cover the target. Report success rate $90/100$, coverage conditional on success $86/90$, and joint success-and-coverage $86/100$. The last quantity treats failure as no usable covering interval and must be labeled that way; it is not a conventional coverage calculation over 100 returned intervals. Do not quietly delete failed replicates or report only the successful methods' fastest runs.
Record warnings, convergence, nonfinite results, nuisance clipping, runtime, and memory when relevant. Prespecify retry rules and record each retry. Separate statistical failure (for example, an unidentifiable target) from numerical failure to compute a defined estimate.
If A succeeds on every dataset but B fails on the difficult ones, comparing A on all datasets with B only on its successes favors B by changing the evaluation population. Report method-specific success rates and conditional summaries, then label any paired comparison restricted to **jointly successful** datasets as such. Keep an attempt-level record so an absent row cannot silently become a successful fit. Runtime should include the prespecified retries and failures when assessing the complete workflow.
## Build a useful evaluation grid
Vary the assumptions that support the claim: independent-unit count, cluster sizes, nuisance-model correctness, smoothness, overlap, sparsity, batch/design confounding, and measurement noise. Include settings where theory predicts success, settings near its boundary, and explicit counterexamples. Change one factor at a time to learn a mechanism, then investigate scientifically plausible combinations.
For the estimator in [From an influence function to an estimator](../semiparametric-methods/influence-function-to-estimator.qmd), compare bias and coverage while changing each nuisance model separately and both together. Record fitted probability ranges and clipping. Compare empirical variance with the known observed-data bound only in scenarios where that bound and the model are specified.
Real omics data usually do not supply the true effect or network. Agreement with another method, pathway enrichment, or reproducibility across cohorts answers a different question from simulated coverage. State what each piece of evidence can support. An empirical analysis can demonstrate feasibility and interpretation without establishing ground-truth accuracy.
## Exercises
Try each checkpoint before opening its hint or solution.
::: {#checkpoint-eval-01 .checkpoint}
### Checkpoint · EVAL-01 — Separate target bias from uncertainty calibration {.unnumbered}
**Task type:** Derivation and evaluation design.
For the shrunken mean with target $\mu=1$, derive its probability limit, bias, sampling variance, and MSE. Explain why a mean SE close to the empirical SD can coexist with poor coverage. Propose a comparison that would expose the failure, and specify how you would report 10 failed fits among 100 replicates.
For a coverage comparison on the same R datasets, define the dataset-level paired contribution and derive the MCSE of the estimated difference.
**Ready to continue when:** you can distinguish the target, the procedure’s limit, sampling uncertainty, and Monte Carlo uncertainty, and compute a paired performance difference.
**Transfer:** Would a method with lower MSE necessarily give a better confidence interval for the same target?
<details class="checkpoint-help">
<summary>Hint</summary>
Compare the fixed bias with the n-dependent interval width; keep the denominators explicit.
</details>
<details class="checkpoint-help">
<summary>Solution outline — open after trying</summary>
The limit is 0.9, bias is -0.1, variance is 0.81/n, and MSE is 0.01+0.81/n. A correct SE for sampling variability does not center the estimate at 1. At increasing n the bias dominates its SD. Include nonzero targets and increasing n, report bias, empirical SD, mean SE, coverage, and MCSE on paired datasets. Report success as 90/100; coverage among successful fits has denominator 90; joint success-and-coverage has denominator 100. The number covering must also be reported.
For coverage differences let $d_r=1\{\theta\in I_{B,r}\}-1\{\theta\in I_{A,r}\}$. Estimate the difference by $\bar d$ and its MCSE by $s_d/\sqrt R$. This retains the within-dataset pairing; a formula that assumes independent marginal coverages omits their covariance.
</details>
:::
### Additional practice (optional) {.unnumbered}
1. Extend the paired calculation to MSE and interval-length differences. State the sign that favors each method and explain why shorter intervals should be judged alongside coverage.
2. Evaluate the missing-outcome estimator under one misspecified nuisance and under two; write down the predicted limits first.
3. Compare multipliers $1-n^{-1/4}$ and $1-n^{-1}$ at a nonzero mean. Derive their biases relative to the root-n sampling SD before simulating coverage. How can both estimators be consistent but support different conclusions about centered intervals?
## Reading
- Morris, White, and Crowther, [Using simulation studies to evaluate statistical methods](https://pmc.ncbi.nlm.nih.gov/articles/PMC6492164/).
<span id="continue-through-the-collection"></span>
