Simulation Diagnostics for High-Dimensional Test Size

Why Diagnose Test Size by Simulation?

A high-dimensional test can have a theoretically specified significance level and still reject too often or too rarely in the sample sizes and dependence regimes used in practice. Simulation provides a direct diagnostic: generate data under a null hypothesis whose properties are controlled, apply the entire testing procedure, and estimate how often it rejects. The crucial qualification is that an estimated rejection rate is itself random. A rate slightly above the target may be ordinary Monte Carlo variation rather than evidence of a defect.

A useful diagnostic therefore has two parts: a design that covers scientifically relevant null conditions, and an uncertainty calculation for every estimated rejection probability. The goal is not to declare a procedure valid from one favorable simulation cell. It is to locate where size departs from its nominal level and to judge whether that departure is larger than simulation error.

Definition: For a test with rejection event $A$, its type-I error probability under a particular null distribution $P_0$ is $P_0(A)$. If the test is designed for nominal level $\alpha$, its size under that distribution is compared with $\alpha$. For a composite null hypothesis $\mathcal{P}_0$, the size is the supremum $\sup_{P\in\mathcal{P}_0}P(A)$; a simulation at one chosen null distribution estimates only that distribution's rejection probability, not the supremum.

In a simulation study, each replicate is one independent data set generated under a specified null condition. Write $R$ for the number of replicates and $I_r=1$ if the procedure rejects in replicate $r$, with $I_r=0$ otherwise. The estimated rejection rate is $\hat{\pi}=R^{-1}\sum_{r=1}^R I_r$. When replicates are independent and identically generated under a fixed null distribution, the $I_r$ are Bernoulli with success probability $\pi=P_0(A)$. This simple structure is what makes binomial uncertainty calculations appropriate.

Designing a Diagnostic Grid

One simulation cell should specify the null distribution, the sample size, the dimension, the test implementation and its nominal level. Change design features systematically rather than changing several unrecorded details at once. For high-dimensional mean tests, a practical grid can vary:

Dimension and sample size: include dimensions below, near and above relevant sample sizes, as well as values representative of the intended application.
Covariance dependence: compare independence with structured dependence, such as block dependence, autoregressive correlation or an admissible equicorrelation model.
Tail behavior: include a Gaussian baseline and suitable non-Gaussian null distributions. State whether they have finite variance and how their covariance is matched to the baseline.
Other implementation conditions: vary nuisance-parameter estimation or missingness only when those features are part of the intended procedure, and label the resulting cells clearly.

Set the null mean to zero for a mean test, but do not assume that this alone specifies the null: covariance and tail behavior can affect calibration. A test evaluated under a Gaussian null with independent coordinates has not thereby been checked under dependent or heavy-tailed observations. Conversely, a single difficult null does not describe behavior over every distribution allowed by a composite null.

Each replicate must execute the procedure that would actually be used. If the test selects a threshold, chooses a tuning parameter, splits the sample, or uses random permutation or other randomized calibration, repeat those steps inside every replicate. Omitting a data-adaptive step evaluates a different procedure. If a calibration object is deliberately held fixed, describe the result as conditional on that object rather than as the performance of a fully rerun pipeline.

Specify a null cell.
Record $n$, dimension, mean, covariance structure, distributional family and nominal $\alpha$.

Generate independent replicates.
Use a reproducible random-number procedure and generate a fresh data set for each replicate.

Run the complete test.
Apply the same estimation, selection, splitting and calibration steps as in the intended analysis.

Count rejections and quantify uncertainty.
Report the rejection count, rate, replicate count and a binomial uncertainty measure for every cell.

Compare patterns across cells.
Look for stable departures from $\alpha$, and account for the uncertainty and number of cells before interpreting the largest deviation.

Monte Carlo Uncertainty for a Rejection Rate

The rejection count $K=\sum_{r=1}^R I_r$ has a binomial distribution with parameters $R$ and $\pi$. Therefore, $\hat{\pi}=K/R$ is unbiased and its variance is known exactly.

Result: If the $R$ simulation replicates are independent and have the same rejection probability $\pi$, then $$ E[\hat{\pi}]=\pi,\qquad \operatorname{Var}(\hat{\pi})=\frac{\pi(1-\pi)}{R}. $$ The Monte Carlo standard error is therefore $\sqrt{\pi(1-\pi)/R}$. Since $\pi$ is unknown, it is commonly estimated by $\sqrt{\hat{\pi}(1-\hat{\pi})/R}$.

Proof. Since each $I_r$ is Bernoulli with probability $\pi$, $E[I_r]=\pi$ and $\operatorname{Var}(I_r)=\pi(1-\pi)$. Linearity of expectation gives $E[\hat{\pi}]=R^{-1}\sum_{r=1}^R E[I_r]=\pi$. Independence makes the covariance terms zero, so $$ \operatorname{Var}(\hat{\pi}) =\frac{1}{R^2}\sum_{r=1}^R\operatorname{Var}(I_r) =\frac{R\pi(1-\pi)}{R^2} =\frac{\pi(1-\pi)}{R}. $$ Replacing $\pi$ by its estimate gives the stated estimated standard error. $\square$

A standard error is a useful scale, not a proof that a particular deviation is harmless. For a confidence interval with better behavior than a simple normal interval near zero or one, the Wilson interval can be obtained by inverting the binomial score inequality. If $z$ is the desired standard-normal critical value, solve $$ \frac{(\hat{\pi}-\pi)^2}{\pi(1-\pi)/R}\le z^2 $$ for $\pi$. The resulting interval has center and half-width $$ \frac{\hat{\pi}+z^2/(2R)}{1+z^2/R} \quad\text{and}\quad \frac{z\sqrt{\hat{\pi}(1-\hat{\pi})/R+z^2/(4R^2)}}{1+z^2/R}, $$ respectively. It is a binomial interval for the rejection probability in that cell, not a guarantee about unexamined null distributions.

Result: Under the same independent-replicate assumptions, for every $\varepsilon>0$, $$ P\bigl(|\hat{\pi}-\pi|\ge\varepsilon\bigr)\le\frac{1}{4R\varepsilon^2}. $$ This conservative bound guarantees that Monte Carlo error exceeds $\varepsilon$ with probability at most the displayed quantity.

Proof. For any random variable $W$ with finite variance, Markov's inequality applied to the nonnegative variable $(W-E[W])^2$ gives $P(|W-E[W]|\ge\varepsilon)\le\operatorname{Var}(W)/\varepsilon^2$. Apply this inequality to $W=\hat{\pi}$, use $E[\hat{\pi}]=\pi$, and note that $\pi(1-\pi)\le 1/4$, since $(\pi-\tfrac12)^2\ge0$. Thus $$ P(|\hat{\pi}-\pi|\ge\varepsilon) \le\frac{\pi(1-\pi)}{R\varepsilon^2} \le\frac{1}{4R\varepsilon^2}. $$ $\square$

Worked Examples

Worked Example: Is a Rejection Rate of 6.8% Persuasive?

A test has nominal level $\alpha=0.05$. In $R=1000$ independent null simulations, it rejects $K=68$ times. Estimate the rejection probability and quantify its uncertainty.

The estimated rate is $\hat{\pi}=68/1000=0.068$, or $6.8\%$. Its estimated standard error is $$ \sqrt{\frac{0.068(1-0.068)}{1000}} =\sqrt{0.000063376} \approx 0.00796. $$ Thus the observed excess over $0.05$, which is $0.018$, is about $0.018/0.00796\approx2.26$ estimated standard errors.

For a 95% Wilson interval, take $z=1.96$. The denominator is $1+3.8416/1000=1.0038416$; the center is $(0.068+3.8416/2000)/1.0038416\approx0.06965$; and the half-width is approximately $0.01566$. The interval is therefore about $(0.0540,0.0853)$. Since $0.05$ lies just below this interval, the result is evidence of excess rejection for this particular null cell, but it is not a diagnosis of every possible null condition. More replicates or a repeated study would help assess how stable the departure is.

Worked Example: Choosing Replicates for One-Percentage-Point Precision

Plan a simulation intended to estimate a rejection probability near $0.05$ with an approximate 95% margin of error no larger than $0.01$. Use the binomial standard-error calculation to choose $R$.

Near $\pi=0.05$, the standard error is approximately $\sqrt{0.05(0.95)/R}$. Setting $1.96$ standard errors equal to $0.01$ gives $$ 1.96\sqrt{\frac{0.0475}{R}}\le0.01 \quad\Longrightarrow\quad R\ge\frac{1.96^2(0.0475)}{0.01^2} =\frac{0.182476}{0.0001} =1824.76. $$ Rounding up gives $R=1825$ as the approximate planning minimum for a rate near $0.05$. If the true rejection rate is unknown and one wants a conservative normal-approximation calculation valid over all rates, use $\pi(1-\pi)\le0.25$: then $R\ge1.96^2(0.25)/0.01^2=9604$. The larger choice protects against the greatest binomial variance, which occurs at $\pi=0.5$.

Worked Example: A Matched Gaussian and Heavy-Tailed Null Grid

Design two null cells for a mean test with $n=40$ observations and dimension $d=100$. Use an equicorrelation covariance matrix with diagonal entries 1 and off-diagonal entries $\rho=0.2$. Compare Gaussian observations with a multivariate $t$-based construction having five degrees of freedom, while matching their covariance matrices.

The covariance matrix is $\Sigma=(1-\rho)I_d+\rho\mathbf{1}\mathbf{1}^{\mathsf T}$, where $\mathbf{1}$ is the vector of $d$ ones. Its eigenvalues are $1-\rho=0.8$, with multiplicity $99$, and $1+(d-1)\rho=1+99(0.2)=20.8$, in the direction of $\mathbf{1}$. Both are positive, so $\Sigma$ is positive definite and is a valid covariance matrix.

For the Gaussian cell, generate each observation as $X\sim N(0,\Sigma)$. For the heavy-tailed cell, independently generate $Z\sim N(0,\Sigma)$ and $U\sim\chi^2_5$, then set $X=\sqrt{3/U}\,Z$. The mean is zero because $Z$ has mean zero and is independent of $U$. Since $E[1/U]=1/(5-2)=1/3$, the covariance is $$ \operatorname{Cov}(X)=E[3/U]\Sigma=3(1/3)\Sigma=\Sigma. $$ Thus the two cells have the same mean and covariance but different tail behavior. Generate 40 independent observations per replicate, apply the complete test, and repeat each cell, for example, $R=5000$ times. Report the two rejection counts, rates, and binomial intervals separately. This design isolates a tail comparison at fixed $n$, $d$, mean and covariance; it does not presume what either simulated rejection rate will be.

Interpreting the Pattern, Not Just the Largest Number

A grid is useful because size failures can appear only in particular regimes. For example, rates might remain near the nominal level at smaller dimensions but depart as $d$ grows relative to $n$, or behave differently under dependence than under independence. Treat these as questions for the simulation, not conclusions to assume in advance. A table or plot should show the rate and its uncertainty for every cell, alongside $R$, $\alpha$, and the generating conditions.

When many cells are examined, the largest estimated excess is especially likely to be affected by chance: even if every cell has correct size, some rates will fluctuate upward. Use intervals or error bars to show cellwise uncertainty, and avoid treating a single interval crossing or missing $\alpha$ as a definitive global verdict. If the scientific question is whether any cell in a prespecified collection has substantial distortion, plan a suitable multiple-comparison or simultaneous assessment rather than selecting the most striking cell after inspection.

Also distinguish Monte Carlo uncertainty from variation in the data-generation design. The binomial interval quantifies uncertainty from using a finite number of independent simulation replicates in one fixed cell. It does not account for having chosen an unrepresentative covariance matrix, an unrealistic tail model, or too narrow a range of dimensions. Nor does it repair dependent simulation replicates: if replicates share random data or otherwise are not independent, the binomial standard error need not apply.

Common Pitfall: Reading a Small Excess as Proof of Invalid Size

Wrong reasoning: “The estimated rejection rate is 5.4% rather than 5%, so the test is oversized.” This ignores that the estimated rate varies from run to run even when the true rejection probability is exactly $0.05$. With a modest number of replicates, a difference of a few tenths of a percentage point can be ordinary simulation noise.

Correction: Report the number of rejections, $R$, and a binomial uncertainty measure, and judge the difference in that context. Then repeat the assessment over a prespecified grid of dimensions, sample sizes, dependence structures and tail conditions. A consistent, practically important departure that exceeds plausible Monte Carlo fluctuation is more informative than one isolated rate. Even then, the conclusion is specific to the simulated null conditions unless the study supports a broader claim.

Key takeaway: A simulated rejection rate is an estimate, not a fixed property of the run. Specify the null cell, rerun the complete testing pipeline in independent replicates, report binomial uncertainty, and interpret patterns across the grid rather than over-reading one percentage.

Check Your Understanding

Use the definitions and calculations above to assess each design choice.

In one null cell, a procedure rejects 42 times in 800 independent replicates. What is the estimated rejection rate, and what formula gives its estimated Monte Carlo standard error?
Why does simulating a zero mean with independent Gaussian coordinates not establish size control under dependent heavy-tailed observations?
For a composite null, what does the rejection rate from one specified covariance matrix estimate, and what does it not establish?
Why should sample splitting, tuning selection or randomized calibration be repeated inside each replicate when those steps are part of the procedure?
Explain why the largest rate among many simulated cells can be high by chance even if all cell-specific tests have correct size.