This report documents the adversarial code review and empirical
validation of bqmm. It is intended as a standing record of
what was checked, how, and with what result — the evidence
behind the claim that the package’s estimates and uncertainty can be
trusted. All scripts referenced live under tools/ and their
raw outputs under tools/validation/.
1. Method
The review combined two tracks:
-
Static review. Three independent adversarial passes
— over (a) the asymmetric-Laplace density, random-number generator and
Stan model, (b) the
lme4→ Stan random-effect mapping, and (c) the Yang–Wang–He (YWH) variance correction. Each finding had to be demonstrated (derived or reproduced in R), not asserted, and was re-verified against the source before being accepted. -
Empirical validation. Simulation studies run
against the compiled package: parameter recovery, simulation-based
calibration (SBC), a variance-estimator coverage “bake-off”, a
correlated-random-effect recovery study, and head-to-head comparisons
with
quantregandlqmm.
2. Verified correct
| Component | Check | Result |
|---|---|---|
| ALD density / quantile | integrates to 1; τ-quantile at mu; R vs Stan agree |
exact |
| Kozumi–Kobayashi RNG | σ-scaling derivation; empirical variance & quantile | correct |
lme4 → Stan sd_map
|
structural check across slope / 3-coef / unbalanced / crossed / interaction designs | 100% correct |
| Parameter recovery (ALD DGP, 100 reps) | bias of β, σ, σ_u | β bias ≈ −0.04 / +0.02; σ̂ = 0.998; σ_u = 0.739 |
| Koenker sandwich | 1/n scaling; vs quantreg nid
|
scaling exact; matches quantreg |
| Simulation-based calibration (300 sims) | rank uniformity for β, σ, sd_re | χ² p = 0.39 / 0.52 / 0.84 / 0.34 — calibrated |
3. Bugs found and fixed
| Severity | Issue | Fix |
|---|---|---|
| Severe | Singular-bread fallback in the variance code was O(1/n²), able to collapse intervals to ~0% coverage | Removed; bandwidth-grow + ridge keep the bread full rank at the correct O(1/n) scale; regression tests added |
| Major (method) | The pure Koenker sandwich under-covers mixed-model fixed effects (intercept coverage 0.72–0.92) because it discards the between-cluster variance | Default replaced with the YWH multiplicative form
Σ_post · G · Σ_post (cluster meat); coverage 0.95–1.00 |
| Minor |
match.call() captured the wrong call |
capture in bqmm()
|
| Minor (hardening) | predictive draws combined parameters from separate
extract() calls |
single joint extract |
4. Coverage bake-off (the headline)
Two-level location-scale designs with known τ-varying coefficients, 60 reps/cell. Frequentist coverage of nominal-95% intervals (intercept / slope):
| DGP / τ | koenker | naive | ywh_cluster | ywh_indep |
|---|---|---|---|---|
| homo 0.25 | 0.80/0.92 | 0.93/0.95 | 1.00/0.98 | 1.00/0.98 |
| homo 0.50 | 0.87/0.92 | 0.93/0.97 | 1.00/0.98 | 1.00/0.98 |
| homo 0.75 | 0.72/0.83 | 0.92/0.85 | 0.97/0.97 | 1.00/0.97 |
| hetero 0.25 | 0.87/0.85 | 0.92/0.90 | 0.98/0.95 | 1.00/0.98 |
| hetero 0.50 | 0.92/0.90 | 0.97/0.93 | 0.98/0.95 | 1.00/0.97 |
| hetero 0.75 | 0.87/0.88 | 0.93/0.87 | 0.95/0.98 | 0.97/0.97 |
ywh_cluster is the only estimator never below nominal.
Confirmed end-to-end through the package’s
confint(adjusted = TRUE) (80 reps): adjusted coverage
0.94–0.99 vs naive 0.90–0.91.
Comparison to lqmm (random-intercept
QR, τ = 0.5): mean |bqmm − lqmm| = 0.10 (intercept), 0.04 (slope) — the
fixed-effect estimates agree with the established frequentist
package.
5. Correlated random effects
A second Stan model adds LKJ-correlated random effects
(cov = "unstructured"). A well-identified, well-converged
diagnostic fit (0 divergences, all R̂ ≤ 1.01) recovers the covariance: β
= (0.93, 2.13) vs (1, 2); SDs = (1.45, 1.04) vs (1.5, 1.0);
correlation ρ = 0.60, 95% CrI [0.36, 0.76] vs true 0.5;
σ = 0.503 vs 0.5.
6. Verification status
- Unit tests: 74 passed / 0 failed, including regression tests for every fix above plus the correlated-model path.
- SBC, recovery, bake-off, and cross-package comparisons all pass.
7. Remaining limitations (documented, not defects)
- Correlated random effects are limited to a single grouping factor; multiple or crossed correlated terms use the diagonal covariance.
- The YWH correction is a large-sample / many-clusters argument and is mildly conservative under weak misspecification.
- Recommended future work: a larger correlation-coverage study and SBC for the correlated model.