Class WildBootstrapMarkovCheck

java.lang.Object
edu.cmu.tetrad.search.WildBootstrapMarkovCheck
public final class WildBootstrapMarkovCheck extends Object
Joint, dependence-aware Markov check for a DAG or CPDAG.

For each target X we test the local constraints X _||_ Y | MB(X) for every Y outside MB(X) U {X}, using a GCM-style residual-product statistic 

   T_j = sqrt(n) * mean_i( rX_i * rY_i ) / sd_i( rX_i * rY_i ),
where rX, rY are residuals of X and Y regressed on MB(X). Under the null each T_j is approximately N(0,1), but the T_j are correlated (they share rX, and the non-blanket variables are correlated through the rest of the graph). We therefore do NOT convert to p-values and assume independence, and we do NOT flip the signs of the averaged T_j (a per-coordinate flip silently assumes the T_j are uncorrelated and badly over-rejexiects).

Instead we calibrate the whole vector jointly with a per-observation wild / multiplier bootstrap: a single sign (or Gaussian) multiplier vector e_1..e_n is drawn per replicate and applied to the centered per-observation products, shared across all constraints. This reproduces the joint covariance of the T_j, so the omnibus statistic is calibrated under dependence.

Row dependence. The dependence handled above is among the CONSTRAINTS; the rows are assumed exchangeable. Under serial or spatial dependence among observations the i.i.d. multiplier bootstrap understates the null variance of mean_i(rX_i*rY_i) and the check is anticonservative. Two dependence-aware variants are provided, both bootstrap-t: the observed and replicate statistics are studentized by a dependence-consistent variance (re-studentization is essential; with the plain statistic, global centering removes exactly the low-frequency variance that dependent multipliers are supposed to carry, and the level error does not vanish at any block size or bandwidth):

  • setBlocks(int[]) — wild CLUSTER bootstrap-t (Cameron, Gelbach & Miller 2008): one multiplier per block per replicate, cluster-robust variance in the denominator. Preferred whenever a grouping exists or can be constructed: towns for tract data, subjects for repeated measures, or spatial grid cells via gridBlocks(double[][], double). Choose blocks at least as large as the dependence range; too-small blocks leave an anticonservative residue. In simulation (GP rows, range ~9 in row units, n=120), i.i.d. multipliers rejected a true null 98% of the time at alpha = .05; block bootstrap-t with adequate blocks rejected ~6%.
  • setKernel(double[][], double) — dependent wild bootstrap-t (Shao 2010; self-normalized): Gaussian multipliers with covariance exp(-d^2/(2 h^2)) over row coordinates, HAC-type kernel variance in the denominator. For dependence without a clean block structure. More bandwidth-sensitive than blocks (same simulation: ~11% at the best h); run a small sensitivity grid over h. medianNearestNeighborDistance(double[][]) times 2–4 is a starting scale.
Defaults are unchanged: with neither option set, multipliers are i.i.d. per row and the statistic is the original one.

Two omnibus statistics are reported:

  • max |T_j| — power against a few badly-violated constraints (sparse alternatives);
  • sum T_j^2 — power against many small violations (diffuse alternatives).

The residual-product statistic is a conditional-covariance test, so like GCM it has no power against dependence that lives only in the conditional variance / higher moments. The bootstrap fixes the calibration under dependence, not the power direction. Swap in a richer WildBootstrapMarkovCheck.Residualizer (e.g. RFF / feature regression) to relax the functional form of the conditional mean.

Typical use: 

  DataSet data = ...;           // continuous (or use a mixed-data Residualizer)
  Graph cpdag  = ...;           // candidate model to check
  WildBootstrapMarkovCheck mc = new WildBootstrapMarkovCheck(data, cpdag)
          .setNumBootstraps(1000)
          .setSeed(13L)
          .setMultiplier(Multiplier.RADEMACHER);
  WildBootstrapMarkovCheck.Result r = mc.check();
  System.out.println(r);                       // p-values + worst constraints
  boolean passes = r.pMax > 0.05;              // thumbs-up / thumbs-down

Run main (within the Tetrad classpath) for a self-contained calibration demo that contrasts the per-observation bootstrap against the broken per-coordinate flip.

  • Constructor Details

    • WildBootstrapMarkovCheck

      public WildBootstrapMarkovCheck(DataSet data)
      Constructs a check over the given data. Configure it with the setters, then call checkMarkovBlanket(Graph) or checkFacts(List).
      Parameters:
      data - the dataset to check; columns are treated as continuous by the default residualizer.

  • Method Details

    • runEngine

      public static WildBootstrapMarkovCheck.Result runEngine(double[][] product, String[] labels, int numBootstraps, long seed, WildBootstrapMarkovCheck.Multiplier mult) throws InterruptedException
      Run the joint check directly on a precomputed product matrix.
      Parameters:
      product - n x K matrix; column j holds the per-observation products rX_i*rY_i for one constraint. Pass uncentered products; centering is handled here.
      labels - length-K constraint labels.
      numBootstraps - number of bootstrap replicates B.
      seed - RNG seed (reproducibility).
      mult - multiplier distribution.
      Returns:
      the omnibus p-values and per-constraint statistics. If K == 0 both p-values are 1.0.
      Throws:
      InterruptedException - if the calling thread is interrupted while the bootstrap runs.

    • runEngineBlocks

      public static WildBootstrapMarkovCheck.Result runEngineBlocks(double[][] product, String[] labels, int numBootstraps, long seed, WildBootstrapMarkovCheck.Multiplier mult, int[] blocks) throws InterruptedException
      Wild CLUSTER bootstrap-t engine for row dependence with a block (cluster) structure. The observed statistic per constraint is the cluster-studentized mean, 
         T_j = sqrt(n) * mean_i(c_ij) / s_j,   s_j^2 = (1/n) sum_b ( sum_{i in b} (c_ij - mean_j) )^2,
      and each replicate multiplies the centered products by one Rademacher/Gaussian draw per BLOCK (shared across constraints, preserving cross-constraint dependence) and is re-studentized by the same cluster-variance formula. Blocks should be at least as large as the row-dependence range.
      Parameters:
      product - n x K matrix of per-observation products (uncentered).
      labels - length-K constraint labels.
      numBootstraps - number of bootstrap replicates B.
      seed - RNG seed.
      mult - multiplier distribution (per block).
      blocks - length-n block assignment; rows with equal values share a multiplier.
      Returns:
      the omnibus p-values and per-constraint statistics.
      Throws:
      InterruptedException - if the calling thread is interrupted while the bootstrap runs.

    • runEngineKernel

      public static WildBootstrapMarkovCheck.Result runEngineKernel(double[][] product, String[] labels, int numBootstraps, long seed, double[][] coords, double bandwidth) throws InterruptedException
      Dependent wild bootstrap-t engine (Shao 2010, self-normalized) for row dependence indexed by coordinates (spatial locations; a single column of time indices for serial data). The observed statistic per constraint is the HAC-studentized mean, 
         T_j = sqrt(n) * mean_i(c_ij) / s_j,   s_j^2 = (1/n) sum_ik W_ik (c_ij - mean_j)(c_kj - mean_j),
      with W_ik = exp(-d_ik^2 / (2 h^2)), and each replicate multiplies the centered products by a Gaussian vector with covariance W (shared across constraints) and is re-studentized by the same kernel variance. O(n^2) work per constraint per replicate. Level is more bandwidth-sensitive than the block engine; prefer blocks when a grouping exists, and run a small sensitivity grid over h otherwise.
      Parameters:
      product - n x K matrix of per-observation products (uncentered).
      labels - length-K constraint labels.
      numBootstraps - number of bootstrap replicates B.
      seed - RNG seed.
      coords - n rows of coordinates aligned with the data rows.
      bandwidth - the kernel bandwidth h in coordinate units; must be positive.
      Returns:
      the omnibus p-values and per-constraint statistics.
      Throws:
      InterruptedException - if the calling thread is interrupted while the bootstrap runs.

    • gridBlocks

      public static int[] gridBlocks(double[][] coords, double cellSize)
      Bins coordinates into a square grid and returns one block id per row: rows falling in the same grid cell share a block. A convenient way to build spatial clusters for setBlocks(int[]) from locations (e.g. LON/LAT): choose the cell size at least as large as the suspected spatial dependence range.
      Parameters:
      coords - n rows of coordinates (any dimension).
      cellSize - the grid cell edge length, in coordinate units; must be positive.
      Returns:
      length-n block ids (arbitrary ints; equal id = same cell).

    • medianNearestNeighborDistance

      public static double medianNearestNeighborDistance(double[][] coords)
      Median Euclidean distance from each row to its nearest neighbor; a scale reference for setKernel(double[][], double) bandwidths and gridBlocks(double[][], double) cell sizes. Multiply by 2–4 to cover a dependence range extending a few neighbor-spacings.
      Parameters:
      coords - n rows of coordinates.
      Returns:
      the median nearest-neighbor distance (0.0 if fewer than 2 rows or if duplicated coordinates dominate; callers should then choose a scale by hand).

    • setResidualizer

      public WildBootstrapMarkovCheck setResidualizer(WildBootstrapMarkovCheck.Residualizer r)
      Sets the residualization used to regress X and Y on the conditioning set. Default is WildBootstrapMarkovCheck.OlsResidualizer.
      Parameters:
      r - the residualizer.
      Returns:
      this, for chaining.

    • setNumBootstraps

      public WildBootstrapMarkovCheck setNumBootstraps(int b)
      Sets the number of bootstrap replicates B. Default is 1000.
      Parameters:
      b - the number of replicates; the smallest attainable p-value is 1 / (b + 1).
      Returns:
      this, for chaining.

    • setSeed

      public WildBootstrapMarkovCheck setSeed(long s)
      Sets the RNG seed for the bootstrap, for reproducibility. Default is 0.
      Parameters:
      s - the seed.
      Returns:
      this, for chaining.

    • setMultiplier

      public WildBootstrapMarkovCheck setMultiplier(WildBootstrapMarkovCheck.Multiplier m)
      Sets the multiplier distribution for the wild bootstrap. Default is WildBootstrapMarkovCheck.Multiplier.RADEMACHER.
      Parameters:
      m - the multiplier distribution.
      Returns:
      this, for chaining.

    • setBlocks

      public WildBootstrapMarkovCheck setBlocks(int[] blocks)
      Dependent wild bootstrap, cluster form: switches the check to the wild cluster bootstrap-t of runEngineBlocks(double[][], String[], int, long, WildBootstrapMarkovCheck.Multiplier, int[]). Use when rows are grouped (tracts within towns, repeated measures within subjects) or a spatial grouping can be constructed via gridBlocks(double[][], double), and within-group row dependence is the concern. Blocks should be at least as large as the dependence range. Clears any kernel setting; pass null to return to i.i.d. multipliers.
      Parameters:
      blocks - length-n block assignment aligned with the data rows; rows with equal values share a multiplier.
      Returns:
      this, for chaining.

    • setKernel

      public WildBootstrapMarkovCheck setKernel(double[][] coords, double bandwidth)
      Dependent wild bootstrap, kernel form: switches the check to the dependent wild bootstrap-t of runEngineKernel(double[][], String[], int, long, double[][], double) (Shao 2010, self-normalized), with multiplier covariance exp(-d^2/(2 h^2)) over the given row coordinates. Use for spatial or serial dependence without a clean block structure; multipliers are necessarily Gaussian and the WildBootstrapMarkovCheck.Multiplier setting is ignored while a kernel is set. Level is more bandwidth-sensitive than the block form; run a small sensitivity grid over h. Clears any block setting; pass null coords to return to i.i.d. multipliers.
      Parameters:
      coords - n rows of coordinates aligned with the data rows (1 column for time order, 2 for LON/LAT, ...).
      bandwidth - the kernel bandwidth h in coordinate units; see medianNearestNeighborDistance(double[][]) for a scale reference.
      Returns:
      this, for chaining.

    • checkMarkovBlanket

      public WildBootstrapMarkovCheck.Result checkMarkovBlanket(Graph graph) throws InterruptedException
      Build the stacked product matrix from the data and graph, then run the joint check. One constraint X _||_ Y | MB(X) is contributed for each target X and each Y outside MB(X) U {X}. Graph nodes absent from the data are skipped.
      Parameters:
      graph - the DAG or CPDAG to check.
      Returns:
      the omnibus p-values and per-constraint statistics.
      Throws:
      InterruptedException - if the calling thread is interrupted while the bootstrap runs.

    • checkFacts

      public WildBootstrapMarkovCheck.Result checkFacts(List<IndependenceFact> facts) throws InterruptedException
      Run the joint wild-bootstrap check against an explicit list of independence facts.

      Each fact X _||_ Y | Z contributes one product column: residualize X and Y on Z (via the configured WildBootstrapMarkovCheck.Residualizer) and take the per-observation product rX_i*rY_i. The facts need not come from a Markov blanket; any constraint set works, e.g. the m-separation implications of a DAG / MAG / PAG. Calibration is the same joint wild bootstrap, so dependence among the T_j (including redundancy among implied facts) is handled.

      Parameters:
      facts - the independence facts to test; one product column each. Members of Z equal to X or Y are ignored.
      Returns:
      the omnibus p-values and per-constraint statistics.
      Throws:
      IllegalArgumentException - if a fact references an unknown variable or has X == Y.
      InterruptedException - if the calling thread is interrupted while the bootstrap runs.

    • main

      public static void main(String[] args)
      Self-contained calibration demo. Simulates null data in which every constraint holds, at several levels of cross-constraint correlation rho, and prints the rejection rate of the per-observation bootstrap against that of the broken per-coordinate sign flip. The former should sit near alpha at every rho; the latter inflates as rho grows.
      Parameters:
      args - ignored.