gdpar

installed. This is a packaging-only change with no effect on results when

cmdstanr is available.

This release consolidates the development cycles documented in detail below

(Sessions 8.3--8.6, the Block 8.4 audit, and the Block 9 re-validation):

• **Distributional regression, K >= 1 slots.** Gaussian, Student-t, Gamma, Beta, Tweedie, Poisson, negative-binomial, zero-inflated and hurdle families, plus heterogeneous per-slot families, with B-spline W bases and arbitrary covariate dimension p.
• **Empirical-Bayes and full-Bayes estimation**, with EB-vs-FB comparators (marginal total variation and joint KSD) and a unified Path C bit-exact comparator for the K > 1 and p > 1 regime.
• **Causal positioning.** A T-learner CATE/ITE bridge on the AMM side and an external meta-learner comparator (grf, EconML via reticulate).
• **Geometric-robustness sampling engine (Block RG).** An opt-in, default-bit-identical Riemannian / sub-Riemannian / Finsler metric stack with a certifying orchestrator and a learned expected-Fisher metric, for stiff and near-deterministic posteriors.
• **Dependence-robust inference (Block 9, Axis 2).** Diagnostics (lag-1 / Durbin-Watson; Moran's I) and block-bootstrap robust standard errors for temporally and spatially dependent data, with a data-driven (Politis-White) block length, on both the Empirical-Bayes and full-Bayes paths. gdpar does not model the dependence; only its inference is made robust to it.
• **Identifiability diagnostics and validity tests** for the population reference, residual diagnostics (Dunn-Smyth quantile residuals; G1/G2/G3), and a golden-regression test suite.

**Validation.** Block 9 closed adversarial internal re-validation, a synthetic

recovery benchmark against mgcv / brms / INLA / rstanarm with known ground truth

(3200 cells), and an organic eBird benchmark (80 cells). gdpar leads on the

distributional structure it models (heavy-tail quantiles, zero-inflated count

means, heteroscedastic scale), is statistically indistinguishable from the

strongest competitor (mgcv) on generic predictive accuracy, is robust where it

makes no modelling claim (e.g. autocorrelation), and is the most computationally

expensive method in the roster. Reports under inst/benchmarks/results/.

Detailed per-session development notes follow.

Axis 2 — full-Bayes dependence-robust inference (Session B9.39, 2026-06-26; D102)

• **The dependence diagnostics and the block-bootstrap robust standard errors now accept a scalar full-Bayes fit (gdpar_fit), not only a scalar Empirical-Bayes fit (gdpar_eb_fit).** gdpar_dependence_diagnostic(), gdpar_dependence_robust(), gdpar_spatial_dependence_diagnostic() and gdpar_spatial_dependence_robust() are unchanged in name and signature; they simply work on either path now, closing the EB/FB asymmetry left open at the start of Axis 2. The K > 1 / p > 1 paths remain deferred on both EB and FB.
• **One shared engine, three class-dispatched touchpoints.** The refit loop, the Politis-White / spatial block-length selectors, the temporal and spatial diagnostics and the honesty contract are path-agnostic; only the per-fit extraction of the point estimate, the model SE and the Dunn-Smyth residuals is dispatched by object class. The Empirical-Bayes code paths are byte-identical, so the existing bit-exact EB regression gate is preserved.
• **Full-Bayes point estimate / model SE = posterior mean / posterior SD** of each AMM coefficient (theta_ref, a_coef, b_coef, W_raw, the last on its raw scale, for parity with the Empirical-Bayes extractor; the theta_ref hyperparameters are excluded). robust_se is the block-bootstrap SD of the per-refit posterior means and se_ratio = robust_se / model_se remains a like-for-like (SD-vs-SD) ratio. Each refit re-runs the full HMC and is markedly more expensive than an Empirical-Bayes refit (the opt-in cost message says so).
• **Per-refit convergence accounting** (new refit_diagnostics field on the robust result: max_rhat, min_ess_bulk, n_divergent_refits, n_high_rhat_refits). Refits are never excluded or down-weighted (that would be a non-random, biasing screen); a single informational note is emitted only when a refit's R-hat clearly exceeds 1.05. This applies to both paths.
• **Documented full-Bayes caveats:** the finite-iteration Monte-Carlo error of each refit's posterior mean conservatively (slightly) inflates robust_se; under an informative prior se_ratio < 1 is benign prior regularization, not an overstated model SE; and the Empirical-Bayes theta_ref point estimate (the Laplace mode) and the full-Bayes one (the posterior mean) are different estimands that coincide only asymptotically. A bagged / widened posterior (BayesBag) is recorded as a deferred lateral, not adopted. The full-Bayes design was cross-lineage-reviewed; the review's concrete proposals (a hard minimum-ESS abort, an estimator-replacing bag) were audited and **not** adopted in favour of the more robust informational accounting, while its documentation findings (the se_ratio < 1 and mode-vs-mean subtleties) were adopted.

Axis 2 — data-driven block length (Session B9.38, 2026-06-26; D101)

• **gdpar_dependence_robust() and gdpar_spatial_dependence_robust() gain an opt-in data-driven block size, block_length = "auto" / block_size = "auto".** The default behaviour is unchanged (the rate-only n^(1/3) / n^(1/4)); "auto" replaces the rate's *constant* with a value selected from the fitted residuals, with no extra Stan refit, and the rate as the fallback. Both functions also gain residual_type / randomize_seed, used only to feed the "auto" selector. The chosen value and method are reported (block_length / block_size and the new block_length_method / block_size_method).
• **Temporal: the canonical Politis & White (2004) selector** (with the Patton, Politis & White 2009 correction), hand-rolled in base R (no **np** dependency): the adaptive flat-top spectral estimate of the residual autocovariances gives b_opt = (2 ghat^2 / D)^(1/3) n^(1/3) with the overlapping moving/circular-block constant D = (4/3) spec^2; longer dependence yields a longer block, white noise yields unit blocks, and a degenerate series falls back to the rate. A p << n caveat is documented.
• **Spatial: a data-driven calibration over the cell grid** (Politis & White has no established spatial plug-in). For each candidate g, cheap (no-refit) spatial block resamples give the bootstrap variance of the design-weighted residual functionals (the coefficient's influence directions); g minimises an empirical mean-squared error — squared bias anchored at the *largest* blocks (least biased) plus a 1/n_tiles sampling-variance term (Lahiri 2003). A decorrelating cross-lineage review supplied the MSE skeleton; two of its choices were corrected after an audit and empirical validation (the bias anchor, and the variance term), because the proposed forms would have made the selector either anticonservative or non-adaptive. A single isotropic g is used (anisotropy is a documented, deferred limitation).

Axis 2 — spatial dependence-robust inference (Session B9.37, 2026-06-26; D100)

• **gdpar_spatial_dependence_diagnostic() (new, exported) — Moran's I for the residuals of a scalar Empirical-Bayes fit.** The spatial sibling of gdpar_dependence_diagnostic(): it builds a row-standardized spatial weight matrix (default k-nearest-neighbour, the data-driven k = max(4, min(round( log n), n - 1)); or a distance band; or a user-supplied W carrying domain knowledge), hand-rolls Moran's I in base R (no **spdep** / **sf** dependency), and reports a significance test — a two-sided permutation test by default (robust to non-normal Dunn-Smyth residuals and asymmetric weights), or the cheaper analytic Cliff-Ord normal approximation (which warns under an asymmetric W). Guards cover isolated locations (zero-weight rows warn and return NA), duplicate coordinates (permitted), small n (hard/soft warnings), and the lon/lat (project-first) and model-misspecification caveats are documented. gdpar does **not** model the spatial dependence; this only makes the violation visible.
• **gdpar_spatial_dependence_robust() (new, exported) — spatial block-bootstrap-by-refit standard errors.** Refits the model on B spatial block-bootstrap resamples and reports bootstrap SEs and percentile intervals alongside the model-based ones, in the working-independence + robust-variance spirit of Liang & Zeger (1986) (point estimates unchanged). Default scheme is non-overlapping **tiled** blocks with a **randomized grid origin** per replicate (Politis-Romano-Lahiri), plus an opt-in overlapping **moving** scheme. The default block side is g = max(2, round(n^(1/4))) cells per axis, the d = 2 case of the variance-MSE-optimal rate M ~ n^(d/(d+2)) points per block, which reduces **exactly** to the temporal n^(1/3) block length at d = 1 (derivation and the registered cross-lineage dissent toward n^(1/6) are in ?gdpar_spatial_dependence_robust). Collinear coordinates abort; a single-cell collapse warns.
• **Internal refactor (regression-gated).** The temporal gdpar_dependence_ robust() and the new spatial function share one extracted engine (.gdpar_dependence_robust_engine) and one residual extractor (.gdpar_dependence_residuals); the temporal path is verified **bit-identical** pre/post refactor by a fixed-seed regression gate (a frozen pre-refactor copy vs. the live engine-backed function, same robust-SE table) plus a Stan-free golden that locks the engine's RNG-consumption order. No change to the temporal API or output.

Bloque RG — RG.7 (step 4): the automated Laplace fallback (Session B9.35, 2026-06-25; D99)

• **gdpar_geom_laplace() (new, exported) — the Laplace approximation as a first-class capability.** Given a gdpar_geom_target, it climbs to the posterior mode (L-BFGS-B warm start + a modified-Newton polish when an exact Hessian is exposed, reading the **same** target and gradient the sampler uses), forms the precision M = -Hessian (the observed information, exact cmdstan Hessian or finite differences), reports its covariance M^{-1}, optional iid draws mode + M^{-1/2} z, and — crucially — a **fidelity diagnostic** of the Gaussian against the true posterior over the same draws: the importance-sampling ESS, the PSIS Pareto-k (when **loo** is installed), and the mean/max log-density drop against its Gaussian expectation d/2. These are distilled into a single scalar label "good" / "poor" / "very_poor" so the approximation is never mistaken for exact MCMC. The **un-floored** condition number is reported and a saddle (non-positive-definite curvature) is flagged loudly. This promotes the validated extractor core of RG.7 (inst/benchmarks/scripts/rg7_laplace_elpd.R) into R/.
• **gdpar_geom_orchestrate() and gdpar_geom_fit() gain laplace_fallback = FALSE (opt-in) and laplace_draws = 0L.** When opted in and a run ends in a **certified limit** (the genuinely non-Gaussian canyon the geometry ladder cannot sample — RG.7's eBird Tweedie count), the orchestrator attaches the Laplace approximation ($laplace) and relabels the status "certified_limit_laplace": the sampling limit still stands **and** a labelled, fidelity-diagnosed Laplace posterior is provided — exactly the mgcv/REML and INLA/Laplace competitor-parity regime. On the out-of-scope path the fallback is attached only when the curvature at the mode is genuinely positive-definite (status "out_of_scope_laplace"); otherwise the certificate is left untouched (no overreach). **The default FALSE leaves the output bit-identical** (no $laplace, the original status set), so the goldens stay bit-identical.

Bloque RG — RG.6 (part ii): integration of the geometry-adaptive orchestrator + the pedagogical vignette (Session B9.31, 2026-06-03)

• **gdpar_geom_bridge() (new, exported) — the durable, path-agnostic core.** Turns an already-fitted gdpar object into the inputs the geometry-adaptive controller gdpar_geom_orchestrate() consumes, **without touching gdpar()**. It reads the fit's compiled cmdstan model and Stan data, exposes the standalone log_prob / grad_log_prob / hessian methods, derives the unconstrained dimension and a posterior-mean warm-start, and returns a (target, geom_target, fisher, reference) tuple. Because the diagnostic needs a *re-samplable* model and a CmdStanMCMC object cannot be re-sampled, the bridge recompiles one from the fit's own Stan source ($code(), with the standalone methods; cmdstanr's content-hash cache makes it a cache hit), while the engine target reuses the fitted object directly. The Hessian compilation is **best-effort**: if higher-order autodiff will not compile (models with custom densities, the Tweedie lpdf being the case RG.7 targets), it falls back to gradient-only methods, which still serve the Euclidean, dense and sub-Riemannian (expected-Fisher) levels. This is the tool RG.7 points at the real Tweedie count of benchmark 9.2.O.
• **gdpar_geom_fit() (new, exported) — the one-call ergonomic entry.** A *sister* of gdpar() (not an internal branch) for the K-individual path: it builds and compiles the model through the **shared seam .gdpar_K_build()**, then runs the orchestrator instead of the default NUTS fit, returning a rich gdpar_geom_fit object. Print methods for both new classes are provided.
• **.gdpar_K_build() (internal) extracted from .gdpar_K().** A behaviour-preserving refactor of the model-building phase (everything up to but not including cs_model$sample): .gdpar_K() and gdpar_geom_fit() now share a **single source** (no duplication, no throwaway computation). The default branch (compile_model_methods = FALSE) issues the byte-identical cmdstanr::cmdstan_model(stan_path) call used before, so the K-path goldens stay bit-identical (verified: the K=2 Gaussian golden re-fits with max|delta| = 0); gdpar_geom_fit() passes compile_model_methods = TRUE to expose the standalone methods.
• **Pedagogical vignette vop08_geometric_robustness.Rmd (new).** An abundant, two-level narrative of the whole Block RG (RG.0 → RG.6): the conceptual bridge to the user's geometry manuscript (with an organic critical reading of its §12 hierarchy and its §14 "future statistical line", and of the ORPHEUS §16.3 demonstrated-vs-conjectured honesty), the diagnostic with size-invariant signals, the five sampler levels, the orchestrator with its armour and its certified limit, Option A, and the integration. The **evaluated** chunks are cheap (suite metadata, the budget/criteria/thresholds, and the geometry engine on pure-R closure targets — no Stan); the heavy cmdstan-backed runs are shown with eval = FALSE and reproduced end to end by the gated script inst/scripts/geometry_pilots_deep.R (GDPAR_RUN_GEOMETRY_PILOTS=1). Opt-in throughout: the default gdpar() fit path is untouched, so the goldens remain bit-identical (canonized as **D97**, session B9.31).

Bloque RG — RG.6 (part 1): Option A folded in — sigma_a_k