Bootstrap Randomization Test Inference
Source:R/inference_all_abstract_rand_bootstrap.R
InferenceRandBootstrap.RdEach of the B null draws is generated by (1) resampling n subject rows
with replacement from the observed data and (2) drawing one fresh assignment
vector w from the actual experimental design run on the resampled covariates.
The test statistic computed on each such draw forms a null distribution against which
the observed statistic (actual data, actual w) is compared.
Details
Abstract class implementing the bootstrap randomization test (BRT), a hybrid of the
nonparametric bootstrap and the randomization test of Fisher's sharp null. This
construction is inspired by Kallus, N. (2018), "Optimal a priori balance in the design
of controlled experiments," Journal of the Royal Statistical Society: Series B,
80(1), 85-112, Section 5 ("Algorithms for inference"), Algorithm 4 – see
@references below for the full attribution and what EDI's implementation adds
beyond it.
Motivating heuristic: the sharp null makes the science table fully known. Under Fisher's sharp null \(H_0: y_i(0) = y_i(1) = y_i\) for all \(i\), the observed outcome is the outcome regardless of assignment, so each observed row \((x_i, y_i)\) is a complete description of that subject: it can be paired with any \(w\) and the outcome is still correct. Resampling rows and drawing a fresh \(w\) therefore simulates an entire new experiment — new subjects from \(\hat{F}_n\) (the empirical distribution of subjects), new assignment from the true, known assignment mechanism — under the assumption that outcomes do not respond to treatment. Both simulation ingredients are faithful to the real data-generating process under \(H_0\): the assignment mechanism is exact (drawn from the actual design, including sequential covariate-dependent designs, since the design depends only on the covariates), and the subject distribution is approximately right by the usual bootstrap argument \(\hat{F}_n \to F\). This is a motivating heuristic, not a proof — see "Known limitations and open theoretical question" below for exactly where it stops short of establishing asymptotic validity, and why. If \(H_0\) is false, the construction still generates the null distribution (it forcibly treats \(y\) as invariant to \(w\)) while the observed statistic drifts into the tail — that asymmetry is the intended source of power, independent of the open validity question below.
Why this outgrew its original motivation. This construction was originally
implemented for a narrow reason: it works even when \(w\) is deterministic or
near-deterministic, exactly the degenerate case Kallus (2018) needed it for (see
@references). It turns out to be useful well beyond that: (1) it is
design-agnostic for free, since it calls draw_ws_according_to_design() on the
resampled data rather than enumerating a permutation space – for the matching-on-the-fly
designs (DesignSeqOneByOneKK14/KK21/KK21stepwise), that permutation
space depends on the whole sequential arrival/matching history and has no closed form
worth enumerating, so a classic randomization test would need bespoke combinatorics this
construction avoids entirely; (2) it generalizes the inferential target from
finite-population (conditional on exactly these \(n\) subjects) to superpopulation,
independent of the degenerate-design motivation; (3) it inherits the rest of the
inference machinery (CI inversion, the estimand/testing-type axes) automatically by
living in the same class hierarchy as everything else, rather than being a bolted-on
special case usable only for the one scenario that motivated it.
How it differs from the pure randomization test. The classic randomization test conditions on the realized sample (the "science table") and is exact for any \(n\). The BRT is unconditional: it targets the distribution of the statistic over both subject sampling and randomization, and the null tested is the compound "sharp null AND subjects are i.i.d. draws from \(F\)". Exactness is traded for a population-level (superpopulation) interpretation. Equivalently, the BRT is a parametric bootstrap test whose "parameter" is the pair (\(F\), design): the design part is known exactly, \(F\) is plugged in via \(\hat{F}_n\), and the sharp null is precisely what makes \(\hat{F}_n\) estimable from one arm-agnostic dataset. The pure randomization test is the special case where \(F\) is conditioned away.
Known limitations and open theoretical question: asymptotic validity is not proven, for two distinct, identifiable reasons. Kallus (2018) himself states that asymptotic validity of this bootstrap test is "believed" rather than established, and explicitly calls it an open question. EDI's implementation does not resolve that question; the two specific gaps are:
Glivenko-Cantelli is not bootstrap consistency. The motivating heuristic above leans on \(\hat{F}_n \to F\) to justify substituting a bootstrap resample for a genuinely fresh i.i.d. draw from \(F\). But \(\hat{F}_n \to F\) (uniform convergence of the empirical CDF, a Glivenko-Cantelli statement) is a claim about the empirical distribution itself; it is not the same as consistency of the bootstrap, i.e. convergence of the sampling distribution of the test statistic computed under resampling from \(\hat{F}_n\) to the statistic's true sampling distribution under \(F\). That second, load-bearing claim is well known to require statistic-specific regularity conditions (smoothness/Hadamard differentiability of the statistic as a functional of the empirical process, adequate moment conditions, etc. — the classical Bickel-Freedman-type distinction) and is known to fail for some statistics even when \(\hat{F}_n \to F\) holds trivially (e.g. non-smooth statistics, extreme-value/max-type statistics, heavy-tailed distributions without enough moments, non-regular estimators). No such regularity condition is verified, case by case, for the estimators this class is composed into.
The naive version of this concern is already mitigated by design; a narrower residual question remains open. A naive i.i.d. row bootstrap would create exact duplicate rows that never occur in real data, and for a design whose assignment mechanism depends on the joint covariate configuration of the whole sample — exactly what EDI's matching-on-the-fly designs (
DesignSeqOneByOneKK14/KK21/KK21stepwise) and fixed matched-pair designs (DesignFixedBinaryMatch) do — those artificial exact ties would fabricate zero-distance ("perfect") matches that a genuinely fresh sample from \(F\) would essentially never produce. This package does not resample rows i.i.d. for matching-capable designs.DesignMatchingAbstract$draw_bootstrap_indices()(design_matching_abstract.R), inherited by every matching-capable design, dispatches to a pair-aware resampler (private$draw_matching_bootstrap_indices()) that resamples reservoir subjects i.i.d. and matched pairs as intact units — preserving each pair's true, historically realized within-pair covariate distance rather than fabricating an artificial zero-distance one. This is the correct fix for the naive-duplicate-row concern, and it is already active for every BRT draw on a matching-capable design (InferenceRandBootstrapinheritsbootstrap_sample_indices()from the same chain, so no separate wiring is needed). What remains open, narrower than the naive concern above: reservoir subjects are still resampled i.i.d. independently of the intact pairs, so a single reservoir subject can still appear more than once in one bootstrap draw; whether two duplicate copies of the same original reservoir subject can subsequently be matched to each other by the re-run sequential matching algorithm (fabricating a same-subject zero-distance pair as a second-order effect, distinct from the naive first-order concern this mitigation addresses) has not been analyzed in this package. This is a narrower, unquantified residual question, not a demonstrated bias.
Status: point 1 above is a documented open theoretical question, not a settled result; point 2's naive form is mitigated by the pair-aware resampler described above, with only the narrower residual question left open. Neither this package nor Kallus (2018) supplies a general proof of asymptotic validity covering point 1; both explicitly flag it as unresolved rather than claiming a proof. Treat the p-values and confidence intervals from this class as resting on a well-motivated but unproven asymptotic argument. A rigorous asymptotic proof (or a demonstrated counterexample) remains future work.
Other implementation notes. (1) For sequential designs the bootstrap sample needs an arrival order; the i.i.d. resampling order plays that role, matching the i.i.d.-arrivals assumption of the KK designs. (2) Studentization is not needed for the motivating heuristic above but, as with any bootstrap test, an asymptotically pivotal statistic improves the level's rate of convergence where the construction is valid.
Users do not instantiate this class directly: every concrete inference class in the
package inherits from it, so its methods (compute_rand_bootstrap_two_sided_pval,
approximate_rand_bootstrap_distribution_beta_hat_T) are available on any
inference object. See InferenceRandBootstrapCI for the companion
confidence interval.
References
Kallus, N. (2018), "Optimal a priori balance in the design of controlled
experiments," Journal of the Royal Statistical Society: Series B, 80(1), 85-112,
Section 5 ("Algorithms for inference"), Algorithm 4 – the nearest prior appearance of
this exact construction (resample subjects, then redraw one fresh assignment from the
design on the resampled covariates), introduced there for the degenerate case of his a
priori balancing designs (PSODs), where a design can admit only one or very few
distinct treatment permutations, so the classic Fisher randomization test has no power
(it always returns p-value 1). Kallus himself cites Efron, B. and Tibshirani, R.
(1993), An Introduction to the Bootstrap, Chapman & Hall, for the bootstrap
ingredient, and Good, P. (2005), Permutation, Parametric and Bootstrap Tests of
Hypotheses, Springer, for the test/confidence-interval duality used to invert his
Algorithm 4 into intervals – and states asymptotic validity of the bootstrap test as
an open question, not a proven result (see "Known limitations and open theoretical
question" above, which this package inherits and extends with the matching-design-
specific failure mode). EDI's contribution beyond Algorithm 4 is generalizing the
construction to arbitrary designs (fixed and sequential, including the matching-on-
the-fly family) and response types, and pairing it with a confidence-interval
inversion (InferenceRandBootstrapCI); it does not resolve the open
validity question.
Super classes
Inference -> InferenceRand -> InferenceRandCI -> InferenceNonParamBootstrap -> InferenceRandBootstrap
Methods
+ inherited public methods from InferenceNonParamBootstrap
InferenceNonParamBootstrap$approximate_bootstrap_distribution_beta_hat_T()InferenceNonParamBootstrap$approximate_m_out_of_n_bootstrap_distribution_beta_hat_T()InferenceNonParamBootstrap$approximate_subsampling_distribution_beta_hat_T()InferenceNonParamBootstrap$compute_bootstrap_confidence_interval()InferenceNonParamBootstrap$compute_bootstrap_two_sided_pval()InferenceNonParamBootstrap$compute_m_out_of_n_bootstrap_confidence_interval()InferenceNonParamBootstrap$compute_m_out_of_n_bootstrap_two_sided_pval()InferenceNonParamBootstrap$compute_subsampling_confidence_interval()InferenceNonParamBootstrap$compute_subsampling_sensitivity()InferenceNonParamBootstrap$compute_subsampling_two_sided_pval()InferenceNonParamBootstrap$get_supported_bootstrap_ci_types()InferenceNonParamBootstrap$get_supported_bootstrap_pval_types()InferenceNonParamBootstrap$select_optimal_b_subsampling()InferenceNonParamBootstrap$select_optimal_m_out_of_n_bootstrap()+ inherited public methods from InferenceRandCI
InferenceRandCI$compute_rand_confidence_interval()InferenceRandCI$compute_rand_two_sided_pval()+ inherited public methods from InferenceRand
+ inherited public methods from Inference
Inference$capabilities()Inference$compute_asymp_confidence_interval()Inference$compute_asymp_two_sided_pval()Inference$compute_estimate()Inference$compute_exact_confidence_interval()Inference$compute_exact_two_sided_pval_for_treatment_effect()Inference$duplicate()Inference$get_analysis_data()Inference$get_covariates()Inference$get_design_object()Inference$get_model_formula()Inference$get_nonestimable_reason()Inference$get_nonestimable_stage()Inference$get_optimization_alg()Inference$get_response()Inference$get_response_type()Inference$get_treatment()Inference$initialize()Inference$is_nonestimable()Inference$set_optimization_alg()Inference$set_seed()Inference$supports()
InferenceRandBootstrap$get_supported_rand_bootstrap_pval_types()
Returns the type values
compute_rand_bootstrap_two_sided_pval() accepts.
InferenceRandBootstrap$approximate_rand_bootstrap_distribution_beta_hat_T()
Computes the bootstrap randomization null distribution of the test
statistic under Fisher's sharp null (shifted by delta): each draw resamples
subject rows with replacement and draws one fresh assignment vector from the design.
Usage
InferenceRandBootstrap$approximate_rand_bootstrap_distribution_beta_hat_T(
B = 501,
delta = 0,
transform_responses = "none",
show_progress = TRUE,
debug = FALSE,
bootstrap_type = NULL,
rand_bootstrap_draws = NULL,
zero_one_logit_clamp = .Machine$double.eps
)Arguments
BNumber of bootstrap randomization draws. Default 501.
deltaThe null treatment effect (on the
transform_responsesscale). Default 0.transform_responsesType of response transformation used to impose the sharp null shift for nonzero
delta. Default "none".show_progressA flag indicating whether a progress bar should be displayed.
debugIf
TRUE, return a list with the distribution values and per-iteration diagnostics including error messages, warning messages, counts of each, and summary proportions for iterations with errors, warnings, and illegal (non-finite) values. Runs serially. DefaultFALSE.bootstrap_typeOptional bootstrap-resampling scheme; see
approximate_bootstrap_distribution_beta_hat_Tfor legal values. DefaultNULL.rand_bootstrap_drawsOptional pre-generated draws (as returned by the private method
generate_rand_bootstrap_draws) enabling common random numbers across calls with differentdeltavalues (used by the CI inversion). DefaultNULL.zero_one_logit_clampThe clamping amount for exact 0 and 1 values when logging.
InferenceRandBootstrap$compute_rand_bootstrap_two_sided_pval()
Computes a bootstrap randomization two-sided p-value for Fisher's sharp
null (shifted by delta): the observed statistic (actual data, actual w)
is compared against the null distribution generated by resampling rows and drawing
fresh assignments from the design.
Usage
InferenceRandBootstrap$compute_rand_bootstrap_two_sided_pval(
B = 501,
delta = 0,
transform_responses = "none",
na.rm = TRUE,
show_progress = TRUE,
bootstrap_type = NULL,
rand_bootstrap_draws = NULL,
zero_one_logit_clamp = .Machine$double.eps,
type = "percentile"
)Arguments
BNumber of bootstrap randomization draws. Default 501.
deltaThe null treatment effect. Default 0.
transform_responsesType of response transformation for the sharp null shift. The default "none" resolves by response type (logit for proportion, log for count and survival, identity otherwise), matching
compute_rand_two_sided_pval.na.rmRemove non-finite null draws. Default
TRUE.show_progressA flag indicating whether a progress bar should be displayed.
bootstrap_typeOptional bootstrap-resampling scheme; see
approximate_bootstrap_distribution_beta_hat_Tfor legal values. DefaultNULL.rand_bootstrap_drawsOptional pre-generated draws for common random numbers across
deltavalues. DefaultNULL.zero_one_logit_clampThe clamping amount for exact 0 and 1 values when logging.
typeTest statistic type.
"percentile"(default) uses the raw estimator as the BRT test statistic, giving an asymptotic p-value that inherits the unconditional superpopulation validity of the BRT."studentized"divides each null draw's signed deviation from the null by its per-draw standard error: \(p = 2\min(P(z^0_b \ge z), P(z^0_b \le z))\) where \(z^0_b = (t^0_b - \delta)/\hat{s}^0_b\); yields asymmetric CI under inversion."symmetric-percentile-t"uses the absolute pivot \(p = P(|t^0_b - \delta|/\hat{s}^0_b \ge |t - \delta|/\hat{s})\); CI is symmetric. Both SE-based types require the class to exposes_beta_hat_T; fall back to"percentile"if the SE is unavailable."smoothed"adds kernel noise \(\varepsilon_b \sim N(0, \hat{\sigma}/\sqrt{n})\) to each resampled draw before imposing the null shift, reducing discreteness in the null distribution. Only meaningful for continuous responses.Theoretical justification. Order-statistic/rank-based estimators (e.g. the Hodges-Lehmann pseudo-median) take only finitely many values, so their bootstrap/ randomization null distribution is a step function; this both coarsens p-values and destabilizes the
uniroot-baseddeltasearch inInferenceRandBootstrapCI'scompute_rand_bootstrap_confidence_interval, which assumes an approximately continuous, monotone p-value curve. Restoring continuity by convolving the resampling distribution with a shrinking-bandwidth kernel is the classical "smoothed bootstrap" device: Silverman (1981, "Density ratios, empirical likelihood and cot death", Applied Statistics 30(2):142-145) for kernel smoothing of a resampled empirical distribution; Silverman & Young (1987, "The bootstrap: To smooth or not to smooth?", Biometrika 74(3):469-479) for applying that smoothing directly to the bootstrap resampling scheme; and Hall, DiCiccio & Romano (1989, "On smoothing and the bootstrap", Annals of Statistics 17(2):692-704) for the bandwidth conditions under which smoothing improves, rather than degrades, coverage accuracy.Implementation caveat (ad hoc, not a certified instance of the above). The bandwidth used here, \(\hat{\sigma}/\sqrt{n}\) (the raw SE-of-the-mean scale of the response), is a pragmatic engineering choice, not one derived from or validated against the bandwidth-selection results in the sources above. Those references smooth the resampling distribution itself with a bandwidth chosen to trade off bias against variance (often shrinking slower than \(n^{-1/2}\), e.g. \(n^{-1/5}\)-type KDE rates); here, noise is instead added directly to each already-resampled response at a fixed \(n^{-1/2}\) rate. The bandwidth is not exposed as a parameter, has no zero-noise escape hatch (the only way to disable smoothing is to pick a different
type), and its coverage behavior has not been validated by simulation in this package. Treat it as a discreteness patch that is qualitatively motivated by the literature above, not a certified implementation of it.Performance. Every class with a C++
compute_fast_rand_bootstrap_distrkernel that operates on real-valued responses (Wilcox HL, simple mean difference, OLS, robust regression, CoxPH, Weibull marginal, log-rank, RMST, KM-diff) accepts the smoothing noise directly in its kernel. This only speeds upInferenceRandBootstrapCI'scompute_rand_bootstrap_confidence_interval(type = "smoothed"), since CI inversion pre-materializes one set of fresh assignments up front (common random numbers reused across everydeltaevaluated during root-finding) and the fast kernel can engage on each evaluation; measured onInferenceAllSimpleWilcox,n = 30,B = 99: the forced-slow-fallback CI took 25.5s, the fast-kernel CI took 0.5s (about 50x). A standalonecompute_rand_bootstrap_two_sided_pval(type = "smoothed")call (no CI inversion) is not accelerated by this: it deliberately draws the fresh assignment lazily per replicate (materialize_w = FALSE), sorand_bootstrap_draw_matrices()cannot build the matrices the fast kernels need and the R-level fallback still runs regardless of this fix. The two ordinal classes (InferenceOrdinalRidit,InferenceOrdinalJonckheereTerpstraTest) still use the slower R-level fallback in every case, because adding continuous Gaussian noise to integer category codes is not statistically meaningful — see the response-type caveat above.
Examples
if (FALSE) { # \dontrun{
seq_des = DesignSeqOneByOneKK14$new(n = 100, response_type = "continuous")
# ... run the experiment: add subjects and responses ...
seq_des_inf = InferenceAllSimpleAverageDiff$new(seq_des)
seq_des_inf$compute_rand_bootstrap_two_sided_pval(B = 501)
} # }