# Dependence-Robust Inference in gdpar (Axis 2)

### Making the uncertainty honest under temporal and spatial dependence, without modelling it (Block 9, Axis 2)

#### José Mauricio Gómez Julián

#### 2026-07-06

# 1. What this vignette covers

gdpar’s Path 1 estimates — Empirical-Bayes (`gdpar_eb`) or full-Bayes (`gdpar`) — assume the observations are conditionally independent given the mean structure. When the data carry temporal (serial) or spatial autocorrelation, that assumption is violated: the point estimates remain consistent if the mean structure is correctly specified, but the model-based (posterior / Laplace) standard errors are too narrow and the intervals under-cover. Sections 2–3 use the Empirical-Bayes path; Section 4 shows the identical workflow on a full-Bayes fit.

Axis 2 of Block 9 is the inferential response to this. It is deliberately not a model of the dependence:

> Honest scope. gdpar does not model the correlated noise or the spatial random field (that is Axis 1 / a later modelling block, deferred and evidence-gated). Axis 2 only makes the inference robust to dependence. The stance is the working-independence + robust-variance estimator of Liang & Zeger (1986): the point estimates are unchanged; only the reported uncertainty is re-estimated to be dependence-robust.

Two symmetric pairs of functions, sharing one internal refit engine:

Dependence | Diagnostic | Robust SE / intervals
--- | --- | ---
Temporal | `gdpar_dependence_diagnostic()` | `gdpar_dependence_robust()`
Spatial | `gdpar_spatial_dependence_diagnostic()` | `gdpar_spatial_dependence_robust()`

The workflow is always the same: diagnose first, and only if dependence is flagged, re-estimate the uncertainty.

# 2. The temporal pair (recap)

    library(gdpar)

    set.seed(1)
    n <- 150
    x <- rnorm(n)
    # AR(1) errors: serial dependence the working-independence fit ignores.
    y <- 1 + 0.5 * x + as.numeric(stats::arima.sim(list(ar = 0.6), n))
    df <- data.frame(x = x, y = y, t = seq_len(n))

    fit <- gdpar_eb(y ~ x, amm = amm_spec(a = ~ x), data = df,
                    chains = 2, iter_warmup = 300, iter_sampling = 300)

    # Step 1 — diagnose. Lag-1 autocorrelation, Durbin-Watson, Ljung-Box.
    gdpar_dependence_diagnostic(fit, index = df$t)

    # Step 2 — if flagged, re-estimate the uncertainty by a temporal block bootstrap.
    gdpar_dependence_robust(fit, data = df, index = df$t, B = 199, seed = 1)

The robust table reports, per coefficient, the `estimate` (unchanged), the `model_se`, the bootstrap `robust_se`, their ratio `se_ratio` and percentile `ci_lower` / `ci_upper`. A `se_ratio > 1` means the model-based SE understates the dependence-robust SE.

The block length. By default `block_length = NULL` uses the rate-optimal `round(n^(1/3))` (it fixes only the rate). You may fix an integer manually, or let the data choose it:

    # Data-driven block length: the Politis-White (2004) automatic selector,
    # computed from the residuals (no extra refit), with the rate as fallback.
    gdpar_dependence_robust(fit, data = df, index = df$t,
                            block_length = "auto", B = 199, seed = 1)

`"auto"` runs the canonical Politis & White (2004) rule (Patton, Politis & White 2009 correction): it reads the residual autocorrelations, finds the lag beyond which they are negligible, and returns `b_opt = (2 ĝ² / D)^(1/3) n^(1/3)` with the overlapping-block constant `D = (4/3) spec²`. Stronger serial dependence gives a longer block; white-noise residuals give unit blocks; a degenerate series falls back to the rate. The chosen value and method are reported in `block_length` and `block_length_method`. (The selector assumes the fitted parameter count is small relative to `n`.)

# 3. The spatial pair

The spatial functions replace exactly two pieces of the temporal machinery — the statistic (lag-1 autocorrelation becomes Moran’s I) and the resampling (1-D contiguous blocks become 2-D spatial blocks) — and reuse everything else.

## 3.1. Diagnostic: Moran’s I

    set.seed(2)
    n <- 200
    gx <- runif(n); gy <- runif(n)            # spatial coordinates
    x <- rnorm(n)
    # An omitted smooth spatial trend lands in the residuals.
    y <- 1 + 0.5 * x + 3 * (gx + gy) + rnorm(n, sd = 0.3)
    df <- data.frame(x = x, y = y)

    fit_sp <- gdpar_eb(y ~ x, amm = amm_spec(a = ~ x), data = df,
                       chains = 2, iter_warmup = 300, iter_sampling = 300)

    gdpar_spatial_dependence_diagnostic(fit_sp, coords = cbind(gx, gy), seed = 1)

Moran’s I is

\[ I \;=\; \frac{n}{S_0}\, \frac{\sum_i \sum_j w_{ij}\,(r_i - \bar r)(r_j - \bar r)} {\sum_i (r_i - \bar r)^2}, \qquad S_0 = \sum_i \sum_j w_{ij}, \]

with \(r\) the residuals and \(W = (w_{ij})\) a spatial weight matrix. Under the null of spatial exchangeability \(\mathbb{E}[I] = -1/(n-1)\).

Choosing the weights. By default gdpar builds a row-standardized \(k\)-nearest-neighbour graph (so \(S_0 = n\)), with the declared heuristic \(k = \max(4, \min(\lfloor \log n \rceil, n-1))\) — robust to irregular spacing and guaranteeing no isolated point. Alternatives:

- `weights = "distance"` — a distance band whose threshold is the smallest that isolates no location.
- `W = <your matrix>` — supply your own neighbourhood (adjacency, flow network, administrative contiguity). This is the right choice whenever domain knowledge defines the structure; it is row-standardized internally.

Significance. The default is a two-sided permutation test (`n_perm = 999`): the residuals are relabelled across locations, \(I\) is recomputed, and \(p = (1 + \#\{|I_{\text{perm}} - \mathbb{E}[I]| \ge |I - \mathbb{E}[I]|\}) / (n_{\text{perm}} + 1)\). It assumes neither normal residuals nor a symmetric \(W\), so it is safe with Dunn-Smyth residuals and the asymmetric kNN graph. The analytic Cliff-Ord normal approximation is available via `test = "analytic"` (cheaper, but it warns under an asymmetric \(W\)).

## 3.2. Robust SE: spatial block bootstrap

    gdpar_spatial_dependence_robust(fit_sp, data = df, coords = cbind(gx, gy),
                                    B = 199, seed = 1)

The bounding box of `coords` is tiled into a `g × g` grid; non-empty cells are resampled with replacement and concatenated to length `n`. By default the grid origin is randomized per replicate (`random_origin = TRUE`, Politis-Romano-Lahiri), which breaks the deterministic cell-boundary artifact. An overlapping `scheme = "moving"` is also available. The output table has the same columns as the temporal one.

## 3.3. The default block size, and why `n^(1/4)`

The block side per axis defaults to `g = max(2, round(n^(1/4)))`. This is the \(d = 2\) case of the rate that minimises the mean-squared error of the block-bootstrap variance estimator. Writing \(M\) for the points per block (so the block has linear extent \(M^{1/d}\) per axis):

- bias from dependence broken at block edges: \(O(M^{-1/d})\) (Künsch 1989; Hall, Horowitz & Jing 1995),
- estimator variance: \(O(M/n)\),

so \(\mathrm{MSE}(M) \sim M^{-2/d} + M/n\) is minimised at \(M \sim n^{\,d/(d+2)}\). At \(d = 1\) this gives \(M \sim n^{1/3}\) points per block — exactly the temporal `block_length = round(n^(1/3))` default — and at \(d = 2\) it gives \(M \sim n^{1/2}\), i.e. \(g^2 = n/M \sim n^{1/2}\) cells, hence \(g \sim n^{1/4}\) cells per axis. The exponent is therefore the variance-optimal rate that reduces correctly to the canonical temporal rate.

> Registered dissent (D100). A decorrelating cross-lineage review argued for the \(n^{1/(d+4)}\) rate (\(n^{1/6}\) at \(d = 2\)). That rate governs a different estimand — the second-order bias / two-sided distribution-function coverage, which gives \(n^{1/5}\) at \(d = 1\) and so does not reduce to the variance default’s \(n^{1/3}\). It is recorded as a dissent, not adopted.

## 3.4. Data-driven block size (`block_size = "auto"`)

The default above fixes only the rate. Its constant can be chosen from the data — but, unlike the temporal case, Politis & White (2004) has no established spatial plug-in (its flat-top spectral-density-at-zero estimator does not extend cleanly to a field in the plane). gdpar therefore calibrates the cells-per-axis `g` over a grid of candidates:

    gdpar_spatial_dependence_robust(fit_sp, data = df, coords = cbind(gx, gy),
                                    block_size = "auto", B = 199, seed = 1)

For each candidate `g`, cheap (no-refit) spatial block resamples give the bootstrap variance \(V(g)\) of the design-weighted residual functionals \((1/n)\,[1, \tilde{gx}, \tilde{gy}]^\top z\) — these are the influence directions of the coefficient, so their optimal block size matches the coefficient’s, which the residual mean alone would not. The block size minimises an empirical mean-squared error,

\[ g^\ast \;=\; \arg\min_g\;\Big[\underbrace{(\tilde V(g) - \tilde V(g_{\min}))^2} _{\text{squared bias}} \;+\; \underbrace{c\,V(g)^2 / n_{\text{tiles}}(g)} _{\text{variance}}\Big], \]

with the bias anchored at the largest blocks \(g_{\min}\) (the least biased, since the dependence-breaking bias of the variance estimator grows like \(g/\sqrt n\)) and the variance scaling like the inverse number of blocks (Lahiri 2003). Stronger / longer-range spatial dependence is captured by larger blocks (smaller `g`); a short-range or near-independent field by smaller blocks (larger `g`); a degenerate calibration falls back to the \(n^{1/4}\) rate. The choice and method are reported in `block_size` and `block_size_method`.

> Provenance and honest scope (D101). A decorrelating cross-lineage review supplied the empirical-MSE skeleton; two of its concrete choices — the bias anchor and the variance term — were corrected after an audit and empirical validation, because the proposed forms would have made the selector anticonservative (anchoring at the smallest blocks) or non-adaptive (a Monte-Carlo jackknife variance that vanishes with the resample count). A single isotropic `g` is used; strongly anisotropic residual dependence is a documented limitation (the minimal fix, two coordinate-wise calibrations, is deferred). As everywhere in Axis 2, this only sizes the resampling blocks — it does not model the dependence.

# 4. The full-Bayes path

Everything above was shown on the Empirical-Bayes path (`gdpar_eb`). The exact same four functions also accept a scalar full-Bayes fit (`gdpar`), so the EB/FB asymmetry is closed: diagnose, then — if flagged — re-estimate, with no change of API.

    # A full-Bayes fit instead of an Empirical-Bayes one:
    fit_fb <- gdpar(y ~ x, amm = amm_spec(a = ~ x), data = df,
                    chains = 2, iter_warmup = 300, iter_sampling = 300)

    # Same diagnostic, same robust SE -- only the object class changes.
    gdpar_dependence_diagnostic(fit_fb, index = df$t)
    gdpar_dependence_robust(fit_fb, data = df, index = df$t, B = 199, seed = 1)

What differs is internal, and worth understanding before you read the table:

- Point estimate and model SE. On the full-Bayes path the per-coefficient point estimate is the posterior mean and `model_se` is the posterior standard deviation (vs. the Laplace mode / conditional posterior SD on the EB path). `robust_se` is, as always, the block-bootstrap SD of the per-refit point estimate, so `se_ratio = robust_se / model_se` stays a like-for-like SD-vs-SD ratio. The posterior mean / SD are chosen (over median / IQR) to keep that parity and to avoid an undeclared normal-scaling constant.
- Cost. Each refit re-runs the full HMC, so a full-Bayes `gdpar_dependence_robust()` is markedly more expensive than the Empirical-Bayes one. Keep `B` modest, and raise `iter_sampling` rather than `B` if the refit ESS (reported in the result’s `refit_diagnostics`) is low — a finite-iteration refit slightly and conservatively inflates `robust_se`.
- Reading `se_ratio` under an informative prior. Unlike the EB path, a full-Bayes `se_ratio < 1` is usually benign: an informative prior concentrates the posterior beyond what the data alone support, so the posterior SD is smaller than the bootstrap SD even under correct independent specification (the conjugate-Gaussian algebra gives \(\text{se\_ratio}^2 = n\tau/(n\tau+\tau_0) < 1\) with prior precision \(\tau_0\)). Only `se_ratio` clearly above 1 flags dependence / misspecification.
- EB vs FB are different estimands for `theta_ref`. The EB point estimate is the Laplace mode; the full-Bayes one is the posterior mean. They coincide asymptotically (Bernstein–von Mises) but may differ in finite samples, so the two paths’ tables need not agree to the last digit on the same data.
- Convergence accounting. The result carries `refit_diagnostics` (aggregate `max_rhat`, `min_ess_bulk`, divergent / high-R-hat refit counts). Under-converged refits are never excluded or down-weighted — that screen would be non-random and would bias the SE — but a single note is emitted when a refit’s R-hat clearly exceeds 1.05.

A bagged / widened posterior (BayesBag) would be a different object — a re-architected estimator, not a robust variance for the same one — and is a deliberately deferred lateral; Axis 2 keeps its honest scope, robust variance, not better estimates, on both paths.

# 5. Reading the result, and the caveats

- `se_ratio > 1` ⇒ the model-based uncertainty understated the dependence-robust uncertainty; trust the robust interval. `se_ratio ≈ 1` ⇒ little correction needed.
- Coordinates are Euclidean. Project lon/lat first (e.g. to UTM); raw lon/lat distorts the neighbour graph, severely at high latitudes. Great-circle distance is deliberately not supported, to avoid a heavy sf / geosphere dependency.
- Moran’s I tests residual exchangeability, not its cause. A significant \(I\) may reflect true spatial dependence or model misspecification (e.g. an omitted nonlinear effect). Interpret accordingly — a misspecified mean structure is better fixed than bootstrapped around.
- Honest scope, restated. The bootstrap delivers robust variance, not better point estimates, and is valid for weak / short-range dependence relative to the block size; it does not rescue long-memory or strong long-range dependence.

# References

- Cliff, A. D. & Ord, J. K. (1981). Spatial Processes: Models and Applications. Pion, London.
- Hall, P., Horowitz, J. L. & Jing, B.-Y. (1995). On blocking rules for the bootstrap with dependent data. Biometrika 82(3), 561-574.
- Künsch, H. R. (1989). The jackknife and the bootstrap for general stationary observations. Annals of Statistics 17(3), 1217-1241.
- Lahiri, S. N. (2003). Resampling Methods for Dependent Data. Springer.
- Liang, K.-Y. & Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika 73(1), 13-22.
- Patton, A., Politis, D. N. & White, H. (2009). Correction to “Automatic block-length selection for the dependent bootstrap”. Econometric Reviews 28(4), 372-375.
- Politis, D. N. & White, H. (2004). Automatic block-length selection for the dependent bootstrap. Econometric Reviews 23(1), 53-70.
- Moran, P. A. P. (1950). Notes on continuous stochastic phenomena. Biometrika 37(1/2), 17-23.
- Nordman, D. J. & Lahiri, S. N. (2004). On optimal spatial subsample size for variance estimation. Annals of Statistics 32(5), 1981-2027.
- Politis, D. N. & Romano, J. P. (1992). A circular block resampling procedure for stationary data. In Exploring the Limits of Bootstrap, 263-270. Wiley.