RATest / Khmaladze Permutation Tests: Errata and Bug Analysis
Status: complete. All quantitative claims are backed by full-scale simulations (lognormal design, 2000 replications, 499 permutations, N ∈ {100, 400}).
Scope. Two procedures are affected, both implemented in the RATest package and used in the associated papers:
- JoE procedure —
PT.Khmaladze.fit()(package).Tests equality of distributions via a martingale-transformed two-sample Kolmogorov–Smirnov statistic. - JAE procedure —
PTQTE.Khmaladze.fit()(package). Quantile-treatment-effect version.
This note documents a formal erratum-style statement of each error and its correction. The single most consequential error is related to how we set the evaluation grid, a bug that silently disables the score-direction projection of the martingale transform; it is discussed last because understanding it requires the others.
1 Erratum
Each item states the error, the mechanism, the observable consequence, the correction, and the current validation status.
1.1 E1. Two-sample empirical process was constructed incorrectly (JoE)
Error. In Khm.trans, the process was built from F1 <- ecdf(Y1)(Y1) and Vn.hat1 <- ecdf(F1), i.e. the probability-integral transform of the treatment sample by its own ECDF, together with g1 <- ecdf(quantile(Y0, t))(t), an ECDF of the control quantiles evaluated on the probability grid.
Mechanism. ecdf(Y1)(Y1) is, for distinct values, always the fixed set {1/m, 2/m, …, 1} regardless of the data, so Vn.hat1 carries no information about Y1. The g1 term mixes the outcome and probability scales. The maximand collapses to F̂₀(t) − t − compensator, a comparison of the control ECDF to the identity on the probability grid.
Consequence. The test statistic does not depend on the treatment sample. Holding the control sample fixed and replacing the treatment outcomes with substantially different values leaves the statistic unchanged to full numerical precision (verified: identical to 10 decimal places across a location shift, a different distribution, and a rescaling).
Correction. Evaluate both ECDFs at the control quantiles, on the outcome scale, per equation (7) of the paper: v(t) = ecdf(Y1)(q0) − ecdf(Y0)(q0) with q0 = quantile(Y0, t, type = 1). The generalized-inverse convention F⁻¹(τ) = inf{y : F(y) ≥ τ} corresponds to type = 1, not the R default type = 7.
Status. Established; confirmed numerically.
1.2 E2. Score evaluated on the wrong scale (JoE)
Error. The extended score was obtained via akj(Y0, t), evaluating the density estimate at the probability grid t.
Mechanism. akj’s evaluation points must be on the outcome scale; t lies in [0,1] while the control outcomes span the data range.
Correction. Evaluate at the control quantiles: akj(sort(Y0), z = q0) with q0 = quantile(Y0, t, type = 1).
Status. Established.
1.3 E3. Score sign/normalization (JoE)
Error. The score column used -score$psi * score$dens.
Mechanism. quantreg::akj returns psi = -f'/f. The paper’s extended score is ġ₂(s) = f₀'(Q₀(s)) / f₀(Q₀(s)) = -psi. The implemented expression -psi * dens = (f'/f) · f = f', i.e. the density derivative, not the score. The density normalization is missing.
Correction. Use gdot2 = -akj(...)$psi.
Status. Established.
1.4 E4. Observed statistic never included among the permutations (JAE)
Error. In PTQTE.Khmaladze.fit, the assignment matrix was generated from random permutations, and stat[1] was treated as the observed statistic; the original assignment vector Z was not included as the observed assignment.
Mechanism. With stat[1] itself a random permutation, it is exchangeable with stat[-1], so P(stat[1] ≥ quantile) = α by construction, independent of the data-generating process.
Consequence. The test rejects at exactly the nominal level under both the null and the alternative — i.e. no power. This matches the reported symptom of rejection rates near 0.05 under both.
Correction. Prepend the observed assignment before computing statistics: Sn <- cbind(Z, Z.pi), so stat[1] is the observed statistic and stat[-1] are the permutation statistics. This mirrors the (correct) structure in the research code main_sim_qp.R, which computes the observed statistic separately from the permutation statistics.
Status. Established.
1.5 E5. Grid resolution disables the score-direction projection (JoE and JAE) — dominant error
Error. The evaluation/integration grid was set to p = 3N.
Mechanism. The two-sample empirical process v(t) = F̂₁(Q̂₀(t)) − F̂₀(Q̂₀(t)) is evaluated at the control quantiles Q̂₀(t). Its resolution is therefore governed by the control sample: v can change only at the (at most n₀) distinct control order statistics, where n₀ = length(y0). This holds regardless of the sample size of the treatment group; that is, the treatment ECDF F̂₁, however many observations it has, is read off at the same n₀ control-quantile points, so the number of distinct increments of v is capped by n₀. On a grid with more than n₀ points, the additional points fall on flat stretches where diff(v) = 0; over those runs the score column of the recursive-least-squares design becomes collinear with the constant column, and the fitted score coefficient is pinned to exactly zero at every step. The compensator then projects out only the constant, discarding the score-direction (location-drift) component the transform is designed to remove.
The transition is sharp and reproducible: the score coefficient is nonzero for p ≤ n₀ and identically zero for p substantially above n₀. Since p = 3N = 6·(N/2) = 6n₀ in the balanced case (and p = 3N > n₀ in every case), the package default always lands in the degenerate region.
This was confirmed directly under both directions of group imbalance. With n₀ = 200 control and 20 treatment observations, p = 200 keeps the score column live (198 of 199 coefficients nonzero) while p = 660 = 3N kills it (0 nonzero). With the imbalance reversed — n₀ = 20 control, 200 treatment — the process still has only 19 distinct increments (set by the 20 control points, not the 200 treatment points), so p = 20 is live while p = 200 and p = 660 are both dead. The bottleneck is the control sample in every configuration; neither the treatment size nor min(n₀, n₁) is the governing quantity.
Consequence. The permutation test is rendered essentially inert. Under a known nuisance parameter (no estimation effect to remove) the transform still appears pivotal in isolation, which is why the defect survived simple size checks on the raw transform. But in the actual permutation procedure, relabeling induces a location drift on each permuted sample that the compensator cannot remove when the score is essentially disabled: the permutation statistics become strongly associated with the residual group-mean gap (measured correlation of approx. 0.71) and are inflated, while the observed (exactly recentered) statistic carries no such drift. The observed statistic is therefore almost never extreme relative to the permutation distribution, and the test rejects at essentially zero rate under both the null and the alternative.
A controlled contrast (identical data, recentering, and permutation scheme, differing only in the grid; lognormal design, 2000 replications, 499 permutations) quantifies this directly. With p = 3N the rejection rate is 0.000 under the null at both N = 100 and N = 400, and never exceeds 0.006 under the alternatives except at the largest configuration (N = 400, σ = 0.5), where it reaches only 0.104. The bug’s effect attenuates with sample size but is not close to resolved at any N the package supports, because the grid p = 3N remains well above n₀ (the control-sample resolution) for all N.
Note this corrects an earlier interim characterization: the degenerate transform does not merely lose power while preserving size, it actually collapses both. Size is 0.000, not the 0.05 one would expect.
Provenance. The 3N choice traces to Parker (2012, A Comparison of Alternative Approaches to Supremum-Norm Goodness of Fit Tests with Estimated Parameters, §6.2), which uses a grid of 3n points “to increase the precision of the compensator.” In that setting the score functions are analytic and the process is a single EDF, so a fine grid is beneficial. The choice is harmful only in a setting with an estimated score and a two-sample difference of ECDFs, where the fine grid produces the flat stretches that zero the score column, as is the papers behind the proposed methds.
Correction. Set the grid to p = length(y0) — the control-group size — equally spaced on (0,1). Because the process is evaluated at control quantiles, p ≤ n₀ guarantees the score projection never degenerates, in balanced and unbalanced designs (confirmed above in both imbalance directions). Note this is the control size specifically: neither the treatment size nor min(n₀, n₁) is correct. Actually, the former re-triggers the bug when control is the smaller group, the latter needlessly discards resolution when control is the larger group.
Two grid constructions were tested on balanced designs: 1. p = n₀, equally spaced on (0,1); 2. a data-adaptive grid at the distinct control order statistics (points only where the process changes).
A full comparison (lognormal design, 2000 replications, 499 permutations, N ∈ {100, 400}) found the two statistically indistinguishable at N = 400 and the order-statistic grid slightly worse at N = 100 (e.g. estimated-score power 0.190 vs 0.208 at σ = 0.5). The order-statistic grid places points densely where data is dense and sparsely in the tails, starving the tail resolution the recursive compensator needs. The equally spaced p = n₀ grid is therefore adopted: equal or better power, and a one-line rule rather than a data-driven grid with tie-handling and permutation-varying size.
Status. Established and validated. On the corrected p = n₀ grid the test has correct size (0.038–0.055 across N and arms) and substantial power. The controlled contrast gives the power attributable to the grid correction (fixed − broken, estimated score throughout):
N |
σ |
p = 3N (broken) |
p = n₀ (fixed) |
recovery |
|---|---|---|---|---|
| 100 | 0.2 | 0.001 | 0.099 | +0.098 |
| 400 | 0.2 | 0.001 | 0.288 | +0.287 |
| 100 | 0.5 | 0.006 | 0.198 | +0.192 |
| 400 | 0.5 | 0.104 | 0.823 | +0.719 |
A residual finite-sample gap remains between the test that estimates the score, and the analytic-score (“oracle”) reference on the corrected grid (e.g. 0.198 vs 0.314 at N = 100, σ = 0.5), but it closes as n grows (0.823 vs 0.935 at N = 400, σ = 0.5), consistent with a consistent score estimator. Attempts to reduce this gap by retuning akj (larger bandwidth, alpha = 1) lowered power on both grids. This is so, because minimizing score sup-norm error is not the same objective as maximizing test power. Therefore, the akj defaults are retained.
Caveat for citation. In this contrast the broken arm has E1–E3 corrected and only the grid left broken; it therefore isolates the grid error’s effect and is not a bit-for-bit reproduction of the originally published tables, which also contained E1–E3. It should be described as the controlled effect of the grid error, not as a reproduction of the published values.
1.6 Finite-sample size under heavier tails
The corrected test has nominal size under the design used in the papers (Lognormal(0, 1/4)), but its finite-sample size degrades under heavier-tailed designs, though things improve as n grows. This is a property of the martingale-transform asymptotics under heavy tails, not a residual error, and is documented here so that a reader re-running the corrected code under a heavier tail does not mistake it for an incomplete fix.
Holding the corrected code fixed and varying only the lognormal scale sdlog, the null rejection rate (σ = 0, homogeneous shift, so H0 holds) is:
N |
sdlog = 0.25 |
0.50 |
0.75 |
1.00 |
|---|---|---|---|---|
| JoE 100 | 0.045 | 0.050 | 0.076 | 0.186 |
| JoE 400 | 0.048 | 0.057 | 0.072 | 0.114 |
| JAE 100 | 0.041 | 0.051 | 0.071 | 0.100 |
| JAE 400 | 0.048 | 0.052 | 0.051 | 0.060 |
Two features confirm this is a finite-sample tail effect rather than a bug. First, at sdlog = 0.25 (the papers’ design) size is essentially exact for both methods at both N. Second, at any fixed heavier tail the over-rejection decreases toward nominal as N grows — e.g. JoE at sdlog = 1 falls from 0.186 (N = 100) to 0.114 (N = 400), and JAE from 0.100 to 0.060 — which is the signature of a consistent test whose asymptotic approximation improves with sample size, not of a coding defect (a defect would not systematically improve with N). JAE recovers faster than JoE, consistent with its default tail trimming at (0.05, 0.95) reducing exposure to the extreme upper tail that the full-range JoE KS statistic retains.