Package 'gdpar'

Description

Subset a gdpar_formula_set by slot name or position

Usage

## S3 method for class 'gdpar_formula_set'
x[i]

Arguments

Value

A named list of formulas (not a gdpar_formula_set; downstream re-validation requires reconstruction via gdpar_formula_set).

Description

Extract a single formula from a gdpar_formula_set

Usage

## S3 method for class 'gdpar_formula_set'
x[[i]]

Arguments

Value

The formula stored at slot i.

Description

Initialise an empty amm_builder object that can be incrementally configured through a sequence of setter calls and finalised into an amm_spec via as_amm_spec. The builder pattern is intended for low-verbosity programmatic construction of specifications in user scripts and serialisation routines (see upcoming amm_save_spec and amm_load_spec). Functionally equivalent to a direct call to amm_spec; the builder adds no new modelling capability.

Usage

amm_build(p = 1L)

Arguments

Details

The builder is a thin wrapper around the existing dims_spec semantics: every amm_set_a_uniform / amm_set_b_uniform call mutates the base of the embedded dims_spec, and every amm_set_a / amm_set_b call layers a per-dimension override on top. Overrides survive subsequent uniform changes, in direct correspondence to dimwise composed with override.

Construction is bifurcated by p at the finalisation step, not at builder time: when as_amm_spec is called with p = 1L, the embedded dims_spec is resolved to a length-one list and the resulting a and b formulas are passed to amm_spec on the scalar path; with p > 1L the embedded dims_spec is passed directly to the multivariate path. This keeps a single source of truth for per-dimension semantics.

Value

An object of class amm_builder with components p, dims (a dims_spec object initialised with NULL base and no overrides), W (NULL), and x_vars (NULL). The object is intended to be passed through a chain of amm_set_* calls and terminated by as_amm_spec.

Methodological notes

The builder exposes the same structural asymmetry between per-dimension components (a, b) and the cross-dimension component (W) as amm_spec: dedicated per-dimension setters are provided for a and b, while W is set globally via amm_set_W. Declaring W per-dimension would silently restrict the model class to the separable sub-class, which is rejected by construction.

See Also

amm_set_a_uniform, amm_set_b_uniform, amm_set_a, amm_set_b, amm_set_W, amm_set_x_vars, as_amm_spec, amm_spec, dimwise, override

Examples

spec <- amm_build(p = 1L) |>
  amm_set_a_uniform(~ x1 + x2) |>
  amm_set_b_uniform(~ x1) |>
  amm_set_W(W_basis(type = "polynomial", degree = 2)) |>
  as_amm_spec()
print(spec)

spec_mv <- amm_build(p = 3L) |>
  amm_set_a_uniform(~ x1 + x2) |>
  amm_set_a(k = 2L, ~ x1) |>
  amm_set_b_uniform(~ x1) |>
  as_amm_spec()
print(spec_mv)

Description

Parse a file written by amm_save_spec and reconstruct the corresponding amm_spec object. Parsing is purely lexical (no source or eval of the file contents) and validates the version header strictly against the running package version.

Usage

amm_load_spec(path)

Arguments

Value

An object of class amm_spec.

Errors

Aborts with a gdpar_input_error for: missing or unreadable file; absent or malformed version header; version mismatch; malformed record line (missing colon); unknown or missing required keys; invalid value grammar for a recognised key; W.type = "user" (not serialisable in the canonical form); dims.K.* records outside 1:p.

See Also

amm_save_spec, amm_spec

Examples

spec <- amm_spec(a = ~ x1 + x2,
                 W = W_basis(type = "polynomial", degree = 2))
tmp  <- tempfile(fileext = ".gdpar")
amm_save_spec(spec, tmp)
spec2 <- amm_load_spec(tmp)
print(spec2)

Description

Write an amm_spec object to a text file in a canonical, human-readable format suitable for version control, archival, and bit-exact reproducibility. The file uses a small key: value grammar with a mandatory version header and a hierarchical naming convention for the per-dimension entries on the multivariate path. The file is parsed by amm_load_spec via a dedicated parser (no source or eval of the file contents), so the serialized form is safe to load from untrusted locations.

Usage

amm_save_spec(spec, path)

Arguments

Details

The serialized format records exactly the constructor inputs of the specification, not derived state. In particular, when the modulating basis W has been materialised at a specific $\theta_{\mathrm{ref}}$ dimension (via the internal materialize_W_basis), the materialised fields are not written; the reconstructed object after amm_load_spec is the unmaterialised W_basis corresponding to the same constructor arguments, which is the form normally produced by amm_spec.

User-defined W_basis objects (type = "user") cannot be serialised because the evaluator is an arbitrary R function whose definition cannot be canonised into the file format. An attempt to serialise such a spec aborts with an informative error.

Value

Invisibly returns path.

File format

The first non-empty line must be the version header # gdpar_spec_version: <version>. Subsequent lines are either comments starting with #, empty lines (ignored), or records of the form key: value. Recognised keys and value grammars:

• p: positive integer.

• a, b: NULL or a one-sided formula literal such as ~ x1 + x2. Used on the scalar path; set to NULL on the multivariate path.

• x_vars: NULL or a literal of the form c("x1", "x2", ...).

• W.type: NULL, polynomial or bspline. user is not serialisable.

• W.degree: positive integer (required for polynomial and bspline).

• W.knots: c(...) of numerics (bspline with interior knots only).

• W.df: positive integer (bspline with df only).

• dims.K.a, dims.K.b: same grammar as a, b, for K in 1:p. Required when p > 1.

Version policy

The version field is checked strictly against the running package version. Until the package reaches its first stable release a mismatch is treated as an error; a subsequent release will introduce a forward-compatible upgrade path.

See Also

amm_load_spec, amm_spec, W_basis

Examples

spec <- amm_spec(a = ~ x1 + x2, b = ~ x1,
                 W = W_basis(type = "polynomial", degree = 2))
tmp  <- tempfile(fileext = ".gdpar")
amm_save_spec(spec, tmp)
spec2 <- amm_load_spec(tmp)
identical(spec$level, spec2$level)

Description

Register a per-dimension override of the additive component for index k. The override replaces the uniform base (set via amm_set_a_uniform, or NULL if never set) for dimension k only. Calling amm_set_a twice with the same k replaces the previous override.

Usage

amm_set_a(builder, k, a)

Arguments

Value

The modified amm_builder.

See Also

amm_build, amm_set_a_uniform, override

Examples

b <- amm_build(p = 3L) |>
  amm_set_a_uniform(~ x1 + x2) |>
  amm_set_a(k = 2L, ~ x1) |>
  amm_set_a(k = 3L, NULL)
print(b)

Description

Replace the base additive formula of the embedded dims_spec with the supplied one-sided formula (or NULL to disable the additive component on the base). Per-dimension overrides previously registered via amm_set_a are preserved and continue to take precedence at their respective indices.

Usage

amm_set_a_uniform(builder, a)

Arguments

Value

The modified amm_builder, returned invisibly for the pipe convention but suitable for direct inspection.

See Also

amm_build, amm_set_a

Examples

b <- amm_build(p = 2L) |>
  amm_set_a(k = 2L, ~ x1) |>
  amm_set_a_uniform(~ x1 + x2)
print(b)

Description

Register a per-dimension override of the multiplicative component for index k. The override replaces the uniform base (set via amm_set_b_uniform, or NULL if never set) for dimension k only. Calling amm_set_b twice with the same k replaces the previous override.

Usage

amm_set_b(builder, k, b)

Arguments

Value

The modified amm_builder.

See Also

amm_build, amm_set_b_uniform, override

Examples

b <- amm_build(p = 2L) |>
  amm_set_b_uniform(~ x1) |>
  amm_set_b(k = 2L, NULL)
print(b)

Description

Replace the base multiplicative formula of the embedded dims_spec with the supplied one-sided formula (or NULL to disable the multiplicative component on the base). Per-dimension overrides previously registered via amm_set_b are preserved.

Usage

amm_set_b_uniform(builder, b)

Arguments

Value

The modified amm_builder.

See Also

amm_build, amm_set_b

Examples

b <- amm_build(p = 2L) |>
  amm_set_b_uniform(~ x1)
print(b)

Description

Store a W_basis object as the modulating component of the specification under construction, or clear it by passing NULL. The modulating component is global to all dimensions of $\theta_i$ and is therefore stored as a single top-level slot of the builder.

Usage

amm_set_W(builder, W)

Arguments

Value

The modified amm_builder.

See Also

amm_build, W_basis

Examples

wb <- W_basis(type = "polynomial", degree = 2)
b <- amm_build(p = 1L) |>
  amm_set_a_uniform(~ x1) |>
  amm_set_W(wb)
print(b)

Description

Record the character vector identifying the covariates that enter the modulating component as the linear factor x in $W(\theta) x$, or clear it by passing NULL. The value is forwarded to amm_spec at finalisation time. When NULL, the package uses the covariates derived from the right-hand side of the model formula passed to gdpar.

Usage

amm_set_x_vars(builder, x_vars)

Arguments

Value

The modified amm_builder.

See Also

amm_build, amm_spec

Examples

b <- amm_build(p = 1L) |>
  amm_set_a_uniform(~ x1 + x2) |>
  amm_set_x_vars(c("x1", "x2"))
print(b)

Description

Declare which components of the Additive-Multiplicative-Modulated canonical form $\Delta(x, \theta) = a(x) + b(x) \odot \theta + W(\theta) x$ are active and how they are parametrized. The result is consumed by gdpar to assemble the design matrices and the Stan model.

Usage

amm_spec(a = NULL, b = NULL, W = NULL, x_vars = NULL, p = 1L, dims = NULL)

Arguments

Details

Two construction paths exist, selected by the dimension p of the per-individual parameter vector $\theta_i$: the scalar path (p = 1L, the default) uses a and b as one-sided formulas directly; the multivariate path (p > 1L) uses the dims argument together with the dimwise and override helpers to declare per-dimension specifications. The two paths are deliberately bifurcated to keep contract semantics clean and backward compatibility absolute.

The additive and multiplicative bases are declared as one-sided R formulas; the package uses model.matrix to build the design matrices from the data frame supplied to gdpar. Each formula is evaluated without an intercept column; the centering conditions (C2) and (C3) of Block 1 are enforced empirically by column-wise centering of the design matrices Z_a and Z_b inside the constructor (each column has its sample mean subtracted, so that colMeans(Z_a) = 0 and colMeans(Z_b) = 0 exactly). Because E_mu[a(X)] = colMeans(Z_a) * a_coef under the empirical mu used by Path 1, the centering of Z_a alone satisfies (C2) for any choice of a_coef in the full J_a-dimensional Euclidean space; analogously for (C3). No additional restriction on the basis coefficients is imposed.

On the multivariate path, the additive and multiplicative components factor per dimension of $\theta_i$ because they depend only on the covariates $x_i$ (and therefore on coordinates of $\mathcal{X}$, not of $\Theta$). The modulating component $W$ depends on the full vector $\theta_{\mathrm{ref}}$ and therefore couples the dimensions; declaring $W$ per-dimension would silently restrict the model class to the separable sub-class. The package therefore enforces the structural asymmetry between (a, b) and W at the API level: dims for the former, a single top-level W argument for the latter.

Value

An object of class amm_spec with components a, b, W, x_vars, level, p, and dims. On the scalar path (p = 1L), dims is NULL and a, b hold the formulas. On the multivariate path (p > 1L), a and b are NULL and dims holds the resolved per-dimension list of length p, each entry being a list with a and b.

Methodological notes

The AMM level implied by the specification is recorded on the returned object: Level 0 if every component is NULL across all dimensions, Level 1 if only the additive component is active in any dimension and no multiplicative or modulating term is active anywhere, and Level 2 otherwise. The identifiability theorems of Block 1 apply at each level under the linearity assumption (LIN), which holds automatically for formula-based bases (linear subspaces of L^2_0(\mu)) and for the polynomial and B-spline cases of W_basis.

Cross-component non-identifiability is not detected by this constructor; it is detected at the Gram-matrix diagnostic of gdpar_check_identifiability, which is called automatically before fitting. For p > 1L, an additional cross-dimension identifiability condition applies (see the planned vignette vop02_arbitrary_p).

Dependencies

Uses terms for parsing the formulas and model.matrix downstream for building the design matrices.

References

See vignette("v01_amm_identifiability", package = "gdpar"), Section 3 (the AMM hierarchy) and Section 5 (standing assumptions (C1)-(C6)).

See Also

W_basis, dimwise, override, gdpar_check_identifiability, gdpar

Examples

spec <- amm_spec(
  a = ~ x1 + x2,
  b = ~ x1,
  W = W_basis(type = "polynomial", degree = 2)
)
print(spec)

spec_mv <- amm_spec(
  p    = 2L,
  dims = dimwise(a = ~ x1 + x2, b = ~ x1),
  W    = W_basis(type = "polynomial", degree = 2)
)
print(spec_mv)

Description

Validate the accumulated state of an amm_builder and convert it to an amm_spec object. Dispatches on p: the scalar path (p = 1L) resolves the embedded dims_spec to the single per-dimension entry and forwards a, b, W, x_vars to amm_spec; the multivariate path (p > 1L) forwards the embedded dims_spec directly. All structural validation (range of override indices, consistency between p and the supplied components) is delegated to amm_spec.

Usage

as_amm_spec(builder)

Arguments

Value

An object of class amm_spec; see amm_spec for the slot layout.

See Also

amm_build, amm_spec

Examples

spec <- amm_build(p = 2L) |>
  amm_set_a_uniform(~ x1 + x2) |>
  amm_set_b_uniform(~ x1) |>
  as_amm_spec()
print(spec)

Description

For separable multivariate W bases (polynomial and bspline with $p > 1$), return a list of length p whose $k$-th entry is a univariate W_basis describing the contribution of coordinate $k$ of theta_ref. Useful for introspection, per-coordinate diagnostics (e.g., the per-$k$ CP/NCP toggle of the pre-flight pipeline), and downstream methods that operate on one coordinate at a time.

Usage

as_per_k(wb)

Arguments

Details

Polynomial and bspline bases are separable by construction: the total basis is the concatenation, in coordinate order, of the univariate bases applied independently to each coordinate of theta_ref. as_per_k() simply reconstructs that decomposition as explicit W_basis objects, each materialized with p = 1L.

This decomposition is the natural input to the per-coordinate CP/NCP toggle of the pre-flight pipeline (Path B'), which operates on each coordinate of theta_ref independently when p > 1.

For type = "user" the separability is unverifiable and the function returns NULL. Callers that need per-coordinate decompositions for user bases must construct them explicitly.

Value

A list of length wb$p. Each entry is itself a W_basis of the same type as wb, materialized with p = 1L. For type = "user" the function returns NULL: separability of a user-supplied evaluator cannot be inferred automatically. A warning is emitted in that case.

See Also

W_basis

Examples

wb <- W_basis(type = "polynomial", degree = 2, p = 3L)
subs <- as_per_k(wb)
length(subs)
subs[[1]]$dim

Description

Flattens the hierarchical per-component, per-coordinate structure into a single data.frame with columns (component, k, identifier, x_name, mean, q05, q50, q95). The identifier column carries the term for a and b slots, the basis_idx (formatted as a string) for W, and NA for theta_ref. The x_name column is NA except for W rows. Useful as input to dplyr / ggplot2 pipelines.

Usage

## S3 method for class 'gdpar_coef'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Value

A data.frame with one row per scalar coefficient summarized.

Description

Returns the tidy per-dim data frame. Suitable for use with subset, aggregate, dplyr::filter, etc.

Usage

## S3 method for class 'gdpar_preflight_report'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Value

Data frame with one row per (component, dim).

Description

Returns an object analogous to coef.gdpar_fit: a gdpar_coef list with components theta_ref, a, b, W. The EB version reports theta_ref$method == "EB" and includes a flag theta_ref$eb_correction_applied so downstream consumers can tell that the credible intervals are EB-corrected (or that the correction was disabled).

Usage

## S3 method for class 'gdpar_eb_fit'
coef(object, ...)

Arguments

Value

A list with class c("gdpar_coef_eb", "gdpar_coef", "list").

Description

Returns posterior summaries (mean and 5%/50%/95% quantiles) of theta_ref and the basis coefficients (additive a_coef, multiplicative b_coef/c_b, modulating W_raw) as an object of class gdpar_coef. The format is consistent across p = 1 (scalar) and p > 1 (multi) fits.

Usage

## S3 method for class 'gdpar_fit'
coef(object, ...)

Arguments

Details

For multi fits, the modulating slot uses the effective per-sample product W_raw * sigma_W (when cp_W is FALSE) so that the coefficients are reported on the natural modulating scale. The per-coordinate b slot draws from the multi parameter c_b (already on the theta_ref-multiplied scale; see Recovery 2, handoff 4).

For K-individual fits (K > 1, distributional regression on the per-slot template amm_distrib_K.stan) the function returns a named list of length K whose entries are gdpar_coef objects (each with p = 1L), one per slot (decision E4.A of sub-phase 8.3.10). The modulating block W_raw is globally shared across slots in the K-individual template and its contribution enters eta_k for every slot (see amm_distrib_K.stan, lines 500-523, and the canonical decision "Scope of W: global" in the likelihood K > 1 decision of handoff 28); accordingly, when any slot declared a non-NULL W the resulting W component is attached to every slot's gdpar_coef (and the per-sample product W_raw * sigma_W is reported when cp_W is FALSE, identically to the scalar path). K-individual fits with grouping (J_groups > 1) are queued for Session 8.4 and raise gdpar_unsupported_feature_error.

Value

For scalar (p = 1) and multi (p > 1) fits an object of class gdpar_coef; for K-individual fits (K > 1) a named list of length K of gdpar_coef objects. Use as.data.frame() on a single gdpar_coef to obtain a long-tidy table; for K-individual fits use lapply(coef(fit), as.data.frame).

Description

Convenience accessor for the convergence diagnostics computed at fit time and stored in a gdpar_fit object.

Usage

diagnostics(fit)

Arguments

Details

The diagnostics are computed once at fit time and stored. Calling this function does not re-run any computation.

Value

An object of class gdpar_diagnostics with R-hat, effective sample size, divergent-transition and tree-depth summaries, and a logical converged flag.

Dependencies

Diagnostics are computed with posterior at fit time.

See Also

gdpar

Examples

if (requireNamespace("cmdstanr", quietly = TRUE)) {
  df <- data.frame(x1 = rnorm(100), y = rnorm(100))
  fit <- gdpar(y ~ x1, amm = amm_spec(a = ~ x1), data = df,
               iter_warmup = 200, iter_sampling = 200, chains = 2)
  diagnostics(fit)
}

Description

Construct a dims_spec object that records the same additive and multiplicative formulas for every dimension $k = 1, \ldots, p$ of a multivariate theta_i vector. Used as the value of the dims argument of amm_spec when all dimensions share the same per-component specification. Per-dimension overrides can be layered on top via override.

Usage

dimwise(a = NULL, b = NULL)

Arguments

Details

This helper enables a low-verbosity declaration of p > 1 specs for the common case where all dimensions of the multivariate theta_i share the same a_k and b_k structures. Broadcasting is explicit: the user must wrap the per-component formulas in dimwise() to opt in. Bare formulas passed directly to the dims argument of amm_spec when p > 1 are rejected, to avoid silent recycling masking bugs.

The dimension p is not stored in the dims_spec; it is resolved at the point of consumption by amm_spec, which takes p as an explicit argument and validates coherence across the spec, the multivariate W basis, and any overrides.

The modulating component W is not part of dims_spec because W couples all dimensions of theta_ref in the canonical AMM form and therefore cannot be declared per-dimension. W is supplied to amm_spec as a separate top-level argument via a multivariate basis.

Value

An object of class dims_spec with components base (the uniform template, a list with a and b) and overrides (an initially empty named list of per-dimension overrides keyed by character form of the integer index k).

Methodological notes

The asymmetry between per-dimension components (a, b) and the cross-dimension component (W) reflects the canonical AMM form

$\theta_i[k] = \theta_{\mathrm{ref}}[k] + a_k(x_i) + b_k(x_i) \theta_{\mathrm{ref}}[k] + (W_k(\theta_{\mathrm{ref}}) - W_k(\theta_{\mathrm{anchor}})) x_i, \quad k = 1, \ldots, p,$

in which $a_k$ and $b_k$ depend only on the covariates $x_i$ (and therefore factor per $k$) while each $W_k$ depends on the full vector $\theta_{\mathrm{ref}}$ (and therefore couples the dimensions). Representing W as a list of independent per-dimension bases would silently restrict the model class to the separable sub-class. dims_spec therefore excludes W by construction.

See Also

override, amm_spec

Examples

base <- dimwise(a = ~ x1 + x2, b = ~ x1)
print(base)

with_override <- override(base, k = 2L, a = ~ x1)
print(with_override)

Description

Useful for logs and condensed printouts.

Usage

## S3 method for class 'gdpar_coef'
format(x, ...)

Arguments

Value

A length-1 character vector.

Description

One-line representation of the report. Used by format in contexts where a print is not desired.

Usage

## S3 method for class 'gdpar_preflight_report'
format(x, ...)

Arguments

Value

Character scalar.

Description

Main entry point of the package. Fits the AMM canonical decomposition of the individual parameter via one of three estimation paths. In this version of the package, only Path 1 (hierarchical Bayesian inference via Stan) is fully implemented; Path 2 and Path 3 wrappers are placeholders that signal the development status.

Usage

gdpar(
  formula,
  family = gdpar_family("gaussian"),
  amm = amm_spec(),
  W = NULL,
  data,
  prior = NULL,
  path = c("bayes", "vcm", "hyper"),
  anchor = "prior_mean",
  skip_id_check = FALSE,
  chains = 4L,
  iter_warmup = 1000L,
  iter_sampling = 1000L,
  adapt_delta = 0.95,
  max_treedepth = 12L,
  refresh = 100L,
  verbose = TRUE,
  seed = NULL,
  group = NULL,
  parametrization = c("auto", "ncp", "cp"),
  parametrization_a = NULL,
  parametrization_W = NULL,
  parametrization_aggregation = NULL,
  id_check_rigor = c("full", "fast"),
  ...
)

Arguments

Details

The function orchestrates the Path 1 fit in five steps: (1) input validation and standardization of the covariates entering the design matrices; (2) the basis-restricted identifiability diagnostic of gdpar_check_identifiability; (3) generation of the Stan model source by substituting the prior placeholders into the static template inst/stan/amm_main.stan; (4) compilation and sampling via cmdstanr; (5) collection of convergence diagnostics from posterior and assembly of the returned object.

Covariates entering the additive and multiplicative bases are centered before fitting, enforcing assumption (C1) of Block 1 at the empirical level. The covariates entering the modulating component are additionally scaled to unit standard deviation; this standardization is recorded so that predictions on new data can apply the same transformation consistently.

The anchor value theta_0 enters the parametrization W(theta) - W(theta_anchor) and is therefore a parametrization device, not an inferential statement about the posterior of theta_ref. Choosing a different anchor changes the parametrization but not the data-generating model.

Value

An object of class gdpar_fit with components fit (the underlying cmdstanr fit object), amm, family, prior, design, anchor, identifiability_report, diagnostics, parametrization (a list with the resolved CP/NCP flags per component plus the pre-flight diagnostic statistics when applicable), call and path.

Methodological notes

Path 1 supports finite-dimensional parametric AMM specifications in this version. Non-parametric extensions (Gaussian-process priors on a, b, W; adaptive splines with growing basis dimension) are deferred to a future version; users who need such flexibility are referred to Path 2 (varying-coefficient models via mgcv).

The polynomial restriction on W_basis in Path 1 is also a v0 limitation; B-spline and user-defined W bases require either a Stan-side basis evaluator that the current template does not implement, or a precomputed basis values approach that complicates the data block. Both extensions are planned.

Dependencies

Uses cmdstanr (Suggests) for compilation and sampling, and posterior (Imports) for convergence diagnostics on the resulting posterior draws. cmdstanr is loaded conditionally via requireNamespace; the function aborts with an informative error if the package is not installed.

References

See vignette("v01_amm_identifiability", package = "gdpar") for the canonical form and identifiability conditions; vignette("v04_asymptotics_path1_bayesian", package = "gdpar") for the asymptotic theory of Path 1; and vignette("vop01_parametrization_toggle", package = "gdpar") for the operational guide to the CP/NCP toggle, including the three-filter pre-flight diagnostic, introspection of fit$parametrization$meta, and known limitations under confounding.

Stan Development Team (2024). Stan User's Guide, version 2.35.

Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., and Riddell, A. (2017). Stan: A probabilistic programming language. Journal of Statistical Software, 76(1).

See Also

amm_spec, W_basis, gdpar_family, gdpar_prior, gdpar_check_identifiability, gdpar_bvm_check, gdpar_contraction_diagnostic

Examples

if (requireNamespace("cmdstanr", quietly = TRUE)) {
  set.seed(NULL)
  n <- 200
  df <- data.frame(
    x1 = rnorm(n),
    x2 = rnorm(n)
  )
  df$y <- with(df, 1 + 0.5 * x1 - 0.3 * x2 + rnorm(n, sd = 0.5))
  spec <- amm_spec(a = ~ x1 + x2)
  fit <- gdpar(
    formula       = y ~ x1 + x2,
    family        = gdpar_family("gaussian"),
    amm           = spec,
    data          = df,
    iter_warmup   = 200,
    iter_sampling = 200,
    chains        = 2
  )
  print(fit)
}

Description

Build a gdpar_meta_learner_adapter that wraps a Python-side EconML estimator (Chernozhukov et al., 2018) via reticulate for use with gdpar_compare_meta_learners. The default estimator is econml.dml.CausalForestDML, the orthogonal double-machine-learning causal forest of Athey, Tibshirani, and Wager (2019). The adapter exposes both the mandatory fit_predict_fun and the optional predict_fun; the latter reuses the fitted Python estimator on a fresh evaluation grid without a refit. Native CIs are produced by EconML's effect_interval(X, alpha = 1 - level) method.

Usage

gdpar_adapter_econml(
  estimator = "CausalForestDML",
  n_estimators = 1000L,
  model_y = NULL,
  model_t = NULL,
  seed = NULL
)

Arguments

Details

The Python module econml is in Suggests and must be installed by the user outside of the package (e.g. reticulate::py_install("econml") or manual installation in the active Python environment). The adapter aborts cleanly when reticulate or the econml module is unavailable.

Cached state caveat: the returned state carries a reference to a Python object. The reference is valid for the duration of the R session in which the bridge was built; serializing the comparison via saveRDS and reloading in a fresh R session invalidates the Python reference and the predict_fun aborts with gdpar_unsupported_feature_error when invoked on a restored state. Rebuild the comparison in such cases.

Value

A gdpar_meta_learner_adapter object with requires_r = "reticulate", requires_py = "econml", native_ci = TRUE, and both fit_predict_fun and predict_fun populated.

References

Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1), C1-C68.

Athey, S., Tibshirani, J., and Wager, S. (2019). Generalized random forests. The Annals of Statistics, 47(2), 1148-1178.

See Also

gdpar_compare_meta_learners, gdpar_meta_learner_adapter, gdpar_adapter_grf.

Examples

if (requireNamespace("reticulate", quietly = TRUE)) {
  adapter <- gdpar_adapter_econml(n_estimators = 200L)
  print(adapter)
}

Description

Build a gdpar_meta_learner_adapter that wraps the R-side causal forest of grf (Athey, Tibshirani, and Wager, 2019) for use with gdpar_compare_meta_learners. The adapter exposes both the mandatory fit_predict_fun and the optional predict_fun so the comparator's predict method can reuse the fitted causal forest on a fresh evaluation grid without a refit. Native CIs are obtained by the normal approximation $\widehat\tau(x) \pm z_{1-\alpha/2} \cdot \sqrt{\widehat{\mathrm{Var}}(\widehat\tau(x))}$ using grf's built-in variance estimator (predict(..., estimate.variance = TRUE)).

Usage

gdpar_adapter_grf(
  num_trees = 2000L,
  sample_fraction = 0.5,
  mtry = NULL,
  honesty = TRUE,
  seed = NULL
)

Arguments

Details

Categorical covariates are not handled by grf directly; the adapter coerces character columns to factors and then applies stats::model.matrix(~ . - 1, ...) to obtain a fully numeric design matrix. Numeric or factor inputs pass through unchanged.

Value

A gdpar_meta_learner_adapter object with requires_r = "grf", native_ci = TRUE, and both fit_predict_fun and predict_fun populated.

References

Athey, S., Tibshirani, J., and Wager, S. (2019). Generalized random forests. The Annals of Statistics, 47(2), 1148-1178.

Wager, S., and Athey, S. (2018). Estimation and inference of heterogeneous treatment effects using random forests. Journal of the American Statistical Association, 113(523), 1228-1242.

See Also

gdpar_compare_meta_learners, gdpar_meta_learner_adapter, gdpar_adapter_econml.

Examples

if (requireNamespace("grf", quietly = TRUE)) {
  adapter <- gdpar_adapter_grf(num_trees = 500L)
  print(adapter)
}

Description

Build a gdpar_formula_set from a sequence of two-sided formulas in the style of brms::bf. The first formula carries the outcome on its LHS and defaults the first slot name to "mu" (the canonical name of the location parameter across built-in families); each subsequent formula must carry the canonical parameter name as its LHS (e.g., sigma ~ a(x)). The result is an object of class gdpar_formula_set identical to what the explicit constructor produces.

Usage

gdpar_bf(...)

Arguments

Details

This is sugar: gdpar_bf(y ~ a(x1), sigma ~ a(x2)) is equivalent to gdpar_formula_set(mu = y ~ a(x1), sigma = ~ a(x2)). Naming the first argument overrides the default "mu": gdpar_bf(theta = y ~ a(x)) produces a single-slot set with name "theta", intended for custom families whose location parameter is not canonically called mu.

The brms package is not a runtime dependency of gdpar; the constructor returns a native gdpar_formula_set regardless of whether brms is installed.

Value

An object of class gdpar_formula_set (see gdpar_formula_set for components).

Dependencies

Internally calls gdpar_formula_set.

See Also

gdpar_formula_set, gdpar

Examples

fs <- gdpar_bf(y ~ a(x1) + b(z1), sigma ~ a(x2))
print(fs)

Description

Verify numerically the conclusion of the Bernstein-von Mises theorem (Theorem 4C of Block 4) for a fitted Path 1 model: that the posterior is asymptotically Gaussian around the maximum likelihood estimator with covariance equal to the inverse Fisher information matrix divided by n. The function refits the model by maximum likelihood (MLE) using a derived Stan model in which the prior block is stripped (see generate_stan_code(mle = TRUE) and the // BEGIN PRIORS / // END PRIORS markers in inst/stan/amm_main.stan); it computes a Hessian-based covariance estimate via a Laplace approximation around the MLE, both on the unconstrained scale (the Stan optimizer is invoked with jacobian = FALSE); and it compares the resulting interval coverage with the Bayesian posterior intervals reported by the fitted model.

Usage

gdpar_bvm_check(fit, parameters = NULL, level = 0.95, verbose = TRUE)

Arguments

Details

This function is opt-in and computationally expensive: it refits the model in MLE mode and inverts a Hessian matrix. It is intended as a methodological audit, not as part of the standard inference flow. The conclusions of gdpar are not affected by calling this function.

Theorem 4C of Block 4 establishes that, for finite-dimensional parametric AMM specifications under the (LAN) condition with non-singular Fisher information at the true parameter, the posterior distribution converges in total variation to the asymptotic-Gaussian distribution of the maximum likelihood estimator. Empirically, this entails that posterior credible intervals at any nominal level should agree with the Hessian-based asymptotic confidence intervals as the sample size grows.

This function performs the empirical comparison at the observed sample size. Substantial discrepancy between the two interval families at large n signals either (i) the limit has not yet been approached, requiring more data; (ii) the (LAN) condition fails (e.g., singular Fisher information at the true parameter); or (iii) the model is misspecified (Block 7).

Applies only to fits with finite-dimensional parametric AMM specifications. Non-parametric components (planned for a future version) are outside the scope of Theorem 4C; the function will abort if invoked on a non-parametric fit. The hierarchical regime activated by the group argument of gdpar() (Block 6.5) is likewise out of scope, because the classical asymptotic theory that underwrites Theorem 4C assumes a fixed-dimension parameter vector while the hierarchical case introduces an increasing number of random anchors; the function aborts with gdpar_unsupported_feature_error when invoked on a grouped fit.

Value

A list of class gdpar_bvm_report with components table (data frame comparing posterior intervals with Hessian-based intervals per parameter), discrepancy (numeric vector of relative interval-width differences), level and warnings. A print method provides a human-readable summary.

Methodological notes

This function is part of the methodological audit toolkit, not of the standard inference flow. It does not modify the fit object in any way; its output is informational. Users running large simulations are advised to call this function selectively rather than after every fit.

The MLE estimate uses Stan's optimize method with the LBFGS algorithm on the prior-stripped variant of the model produced by generate_stan_code(mle = TRUE); the optimizer is invoked with jacobian = FALSE so the maximum is taken on the constrained (natural) scale of the parameters rather than on the unconstrained scale that would carry the Jacobian of the transformation. The Hessian-based interval estimate comes from cmdstanr::laplace applied around the MLE on the same prior-stripped model and with the same jacobian = FALSE convention.

The comparison is restricted to the user-facing parameters that the likelihood identifies in MLE mode without the prior anchoring: the global reference parameter theta_ref, the response-level scale sigma_y (Gaussian families), and the dispersion phi (Negative Binomial). The hierarchical scales sigma_a, sigma_b, sigma_W are not identified by the likelihood alone (they enter only through the random-effect priors that have been stripped) and are excluded from the table.

Dependencies

Uses cmdstanr for the optimization run and posterior to extract the Bayesian intervals.

References

See vignette("v04_asymptotics_path1_bayesian", package = "gdpar"), Section 7 (Theorem 4C).

See Also

gdpar, gdpar_contraction_diagnostic

Examples

if (requireNamespace("cmdstanr", quietly = TRUE)) {
  df <- data.frame(x1 = rnorm(200), y = rnorm(200))
  fit <- gdpar(y ~ x1, amm = amm_spec(a = ~ x1), data = df,
               iter_warmup = 200, iter_sampling = 200, chains = 2)
  gdpar_bvm_check(fit)
}

Description

Estimate the conditional average treatment effect (CATE) from a pair of independent gdpar_fit objects fitted to disjoint arms of a treatment / control design. Implements the T-learner meta-learner of Kuenzel et al. (2019) on the AMM-side: each arm is fitted independently and the CATE is the per-observation difference of the two predictive distributions, evaluated on a common evaluation set.

Usage

gdpar_causal_bridge(
  fit_treat,
  fit_ctrl,
  newdata = NULL,
  type = c("response", "theta_i", "linear_predictor"),
  level = 0.95,
  ...
)

Arguments

Details

This function does not modify either fit. It assumes the two fits are independent posterior samples from disjoint subsets of the population (treatment arm and control arm). It does not perform any causal adjustment beyond what is encoded in the two AMM specifications: the assumption of no-unmeasured-confounding within each arm is the responsibility of the user (see Section 4 of the bridge vignette).

Structural compatibility. The function aborts with gdpar_unsupported_feature_error when the two fits differ in any of: path (both must be the Path 1 Bayesian fit), family identifier (family$name, or per-slot family identifiers when K > 1), K, p, AMM level (amm$level or the equivalent level inferred from each slot's spec), modulating basis type (amm$W$type), anchor value, or covariate column structure. The function also aborts when either fit was sampled in the hierarchical regime (stan_data$use_groups == 1L); the T-learner bridge for grouped fits is queued for a future sub-phase and would require careful treatment of the per-group anchors which is outside the scope of Sub-phase 8.5.A.

Identifiability per arm. The constructor records the identifiability report of each fit in the id_check slot. (C7) anti-aliasing of Block 6.5 is not invoked because the hierarchical guard above rules out the regime in which (C7) applies; this is documented for the eventual extension to grouped fits.

CATE estimator. For each posterior draw indexed by $s = 1, \dots, S$ and observation $i = 1, \dots, n_{\text{new}}$, the bridge computes

$\widehat{\tau}^{(s)}_i = \hat{\mu}^{(s)}_{\text{treat}}(x_i) - \hat{\mu}^{(s)}_{\text{ctrl}}(x_i),$

where $\hat{\mu}^{(s)}_{\text{arm}}(x)$ is the posterior prediction of the chosen type at $x$, drawn from the fit's predictive distribution. The marginal posterior of the CATE at each $x_i$ is summarized by the empirical mean and the $(\alpha/2, 1-\alpha/2)$ quantiles with $\alpha = 1 -$ level.

Independence of draws. The two fits are independent (they were sampled from disjoint data subsets), so the joint posterior of $(\theta_{\text{treat}}, \theta_{\text{ctrl}})$ factorizes and any pairing of marginal draws is a valid sample from the joint. The function trims to $S = \min(S_{\text{treat}}, S_{\text{ctrl}})$ when the two fits differ in number of draws and emits a gdpar_diagnostic_warning.

Multi-dimensional and K-individual fits. For p > 1 (multivariate) and K > 1 (distributional regression), predict.gdpar_fit returns a 3-array of shape [S, n, dim]; the CATE is computed elementwise and the per-coordinate or per-slot CATEs are returned as the last dimension of cate_draws. For type = "response", the canonical inverse link of each coordinate or slot is applied by predict.gdpar_fit before the difference is taken; the resulting CATE is therefore on the natural response scale of each slot, not a uniform link-transformed scale.

Value

An object of class gdpar_causal_bridge with components cate_draws (matrix [S, n] when both fits are scalar, or array [S, n, dim] when both fits are multivariate or K-individual), cate_mean, cate_ci, newdata, id_check, fits, type, level, n_draws, n_obs, call, warnings (character vector recording fallback notifications such as posterior-draw trimming; empty in the happy path), and meta. The companion S3 methods print and summary are documented in print.gdpar_causal_bridge and summary.gdpar_causal_bridge.

Methodological notes

The T-learner is the most direct meta-learner to map onto the gdpar pipeline: each arm is one gdpar_fit and the CATE reuses the posterior machinery of predict.gdpar_fit. S-learner and X-learner are queued for Block 9. The T-learner is known to suffer from regularization-induced bias in unbalanced samples (see Kuenzel et al. 2019, Section 3.4); the bridge vignette discusses the trade-off and the alternatives.

Dependencies

Inherits the posterior dependency of predict.gdpar_fit.

References

Kuenzel, S. R., Sekhon, J. S., Bickel, P. J., and Yu, B. (2019). Metalearners for estimating heterogeneous treatment effects using machine learning. Proceedings of the National Academy of Sciences, 116(10), 4156-4165.

See Also

gdpar, predict.gdpar_fit

Examples

if (requireNamespace("cmdstanr", quietly = TRUE)) {
  n <- 300L
  df <- data.frame(x1 = rnorm(2L * n))
  df$arm <- rep(c("treat", "ctrl"), each = n)
  df$y <- with(df, ifelse(arm == "treat", 0.5, 0) + 0.8 * x1 + rnorm(2L * n, sd = 0.5))
  df_treat <- subset(df, arm == "treat")
  df_ctrl  <- subset(df, arm == "ctrl")
  fit_t <- gdpar(y ~ x1, amm = amm_spec(a = ~ x1), data = df_treat,
                 iter_warmup = 200, iter_sampling = 200, chains = 2)
  fit_c <- gdpar(y ~ x1, amm = amm_spec(a = ~ x1), data = df_ctrl,
                 iter_warmup = 200, iter_sampling = 200, chains = 2)
  bridge <- gdpar_causal_bridge(fit_t, fit_c,
                                 newdata = data.frame(x1 = seq(-2, 2, length.out = 21L)))
  print(bridge)
  summary(bridge)
}

Description

Diagnose, before model fitting, whether the chosen finite parametric representation of the AMM canonical form satisfies the basis-restricted Functional Independence Condition at a candidate value of the population reference. Failure indicates that the model parameters are not identifiable in the chosen basis and that fitting should not proceed without revising the specification.

Usage

gdpar_check_identifiability(
  amm,
  data,
  theta_ref_init = NULL,
  formula_rhs = NULL,
  family = NULL,
  tol = 1e-08,
  rigor = c("full", "fast")
)

Arguments

Details

Proposition 1C of Block 1 establishes that, in a chosen finite basis B for the AMM components, the basis-restricted Functional Independence Condition at a value theta_ref holds if and only if the population Gram matrix of the extended design matrix Z_n(theta_ref) is non-singular. This function computes the empirical Gram matrix from the observed data and checks its condition number.

The diagnostic is by design local to the supplied theta_ref_init: it reports the condition at one point in the parameter space. The basis-restricted Functional Independence Condition is a property of the basis, not of the point; in practice the diagnostic at the prior mean is informative because non-identifiability that is structural to the basis manifests at typical reference points.

The function does not test the abstract Functional Independence Condition, which is a property of the full function classes F_a, F_b, F_W. When these classes are infinite-dimensional and the basis B is a finite truncation, abstract failure may occur in directions outside B and would not be detected here.

Value

An object of class gdpar_identifiability_report with components passed (logical), lambda_min, lambda_max, condition_number, collinear_directions (a list describing basis-function combinations corresponding to near-zero eigenvectors when the diagnostic fails, and NULL otherwise), theta_ref_used, tol_used and column_labels. A print method provides a human-readable summary.

Methodological notes

The threshold criterion is relative (smallest eigenvalue divided by the largest, compared against tol), not absolute. The relative criterion coincides with the inverse of the condition number, which is invariant to rescaling of the basis columns. An absolute eigenvalue threshold would depend on the scale of the covariates and would produce false positives when the basis matrix has columns of very different magnitudes. To make the diagnostic robust to scale heterogeneity across basis terms, the columns of the extended design matrix are normalized to unit norm before the Gram matrix is computed. Normalization affects only the diagnostic; model fitting uses the matrices in their natural scale.

When the diagnostic fails, collinear_directions is populated by projecting the eigenvectors associated with eigenvalues below tol times the largest eigenvalue onto the named basis columns, producing a human-readable description of which combinations of basis functions are linearly dependent.

Dependencies

This function calls eigen with symmetric = TRUE for the eigendecomposition of the Gram matrix, and model.matrix for evaluating the formula-based bases of the additive and multiplicative components.

References

Block 1 of the package theoretical addendum, Section 6.6, Proposition 1C. See vignette("01_amm_identifiability", package = "gdpar").

See Also

amm_spec, W_basis, gdpar

Examples

set.seed(1)
df <- data.frame(x1 = rnorm(50), x2 = rnorm(50))
spec <- amm_spec(
  a = ~ x1 + x2,
  b = ~ x1,
  W = W_basis(type = "polynomial", degree = 1),
  x_vars = "x2"
)
report <- gdpar_check_identifiability(spec, df, theta_ref_init = 0.5)
print(report)

Description

Sub-phase 8.6.E (Charter Section 3.5, decision 2.5 Trio of vignettes). Reports the operational comparison promised in v07 Section 11 between a gdpar_eb_fit (Empirical Bayes via gdpar_eb) and a gdpar_fit (Fully Bayes via gdpar) fitted on the same dataset: per-component differences in the population anchor $\theta_{\text{ref}}$, the empirical total variation distance between the lower-level posteriors of $\xi = (a, b, W, \text{dispersion})$ marginally per parameter, and the operational verification of the higher-order coverage discrepancy of v07 Section 6 (Proposition 7B scalar / 7B* matricial / 7B* tensorial) on the nominal EB and FB credible intervals.

Usage

gdpar_compare_eb_fb(eb_fit, fb_fit, level = 0.95, tv_bins = 30L, ...)

Arguments

Details

The comparator is descriptive: it does not assert algorithmic equivalence, nor does it test hypotheses across the EB and FB inferential frames. The TV distance is computed marginally parameter by parameter via histogram-based plug-in (relative bin counts at a common breakpoint grid); a finite-sample correction is not applied. Joint TV across the high-dimensional $\xi$ typically requires kernel Stein discrepancy or similar density-free metrics that are out of scope of the initial 8.6.E iteration; the marginal TV reported here is the operational proxy recommended in v07 Section 11.1.

Value

An object of class gdpar_eb_fb_comparison with components theta_diff_table, tv_table, coverage_table, level, tv_bins, n_common_params, path_eb, path_fb, call, warnings (character vector recording per-helper fallback notifications: silent extraction failures of the FB theta_ref draws, missing EB or FB $\xi$ draws, or zero-common-parameter TV inputs; empty in the happy path), and meta. See print.gdpar_eb_fb_comparison and summary.gdpar_eb_fb_comparison for the companion S3 methods.

Path coverage

The comparator handles all four EB regimes uniformly by extracting the per-element anchor estimate (vector / matrix / 3D array) and the corresponding lower-level posterior draws via the canonical posterior::as_draws_matrix interface. For Path C (K > 1 + p > 1) the theta_ref_kp_hat tensor is flattened to a length-K*p vector keyed by (slot, coord) for per-element comparison; the joint K x p inflation tensor is reported in the coverage_table per cell via its diagonal block entries.

Diagnostic value

Under the standing hypotheses of v07 Section 4 (EB-MARG-ID + PRIOR-FB-WEAK + HIER-COMPLEX), Theorem 7A predicts marginal TV -> 0 in probability as n -> Inf (specializing to Theorems 7A* / 7C* / 7C* compound multi-slot under Path A / Path B / Path C of v07b Sections 4-6). The empirical TV reported here is the operational diagnostic of this theoretical prediction. Persistent large marginal TV across $\xi$ suggests one of the discrepancy conditions of Proposition 7D (multi-modality of the marginal likelihood, near-singular Fisher information, informative prior on theta_ref, deep hierarchy). The coverage_table operationally verifies the $O(n^{-1})$ under-cover claim of Proposition 7B by comparing EB-nominal vs FB-nominal IC widths per anchor cell.

References

Petrone, S., Rousseau, J., and Scricciolo, C. (2014). Bayes and empirical Bayes: do they merge? Biometrika, 101(2), 285–302.

Rousseau, J., and Szabo, B. (2017). Asymptotic behaviour of the empirical Bayes posteriors associated to maximum marginal likelihood estimator. Annals of Statistics, 45(2), 833–865.

Carlin, B. P., and Gelfand, A. E. (1990). Approaches for empirical Bayes confidence intervals. JASA, 85(409), 105–114.

See Also

gdpar_eb, gdpar, vignette("v07_eb_vs_fb", package = "gdpar"), vignette("v07b_eb_multivariate", package = "gdpar"), vignette("vop07_eb_workflow", package = "gdpar").

Examples

if (requireNamespace("cmdstanr", quietly = TRUE) &&
    requireNamespace("posterior", quietly = TRUE)) {
  set.seed(20260526L)
  n <- 120L
  df <- data.frame(x = stats::rnorm(n))
  df$y <- 0.5 + 0.4 * df$x + stats::rnorm(n, sd = 0.3)
  spec <- amm_spec(a = ~ x)
  fit_eb <- gdpar_eb(
    formula = y ~ x, family = gdpar_family("gaussian"),
    amm = spec, data = df,
    iter_warmup = 200L, iter_sampling = 200L, chains = 2L,
    refresh = 0L, verbose = FALSE, seed = 1L
  )
  fit_fb <- gdpar(
    formula = y ~ x, family = gdpar_family("gaussian"),
    amm = spec, data = df,
    iter_warmup = 200L, iter_sampling = 200L, chains = 2L,
    refresh = 0L, verbose = FALSE, seed = 1L
  )
  cmp <- gdpar_compare_eb_fb(fit_eb, fit_fb)
  print(cmp)
}

Description

Evaluate a fitted gdpar_causal_bridge object against a user-supplied set of external meta-learner adapters (e.g. gdpar_adapter_grf for grf R-side, gdpar_adapter_econml for EconML Python-side) on a common evaluation grid. The function does not refit either of the two gdpar fits embedded in bridge; it only consumes the bridge's CATE estimates and reconstructs the (X, T, Y) dataset needed by the external adapters from the captured calls of the two fits, or from an explicit data argument when the captured calls cannot be resolved.

Usage

gdpar_compare_meta_learners(
  bridge,
  methods,
  newdata = NULL,
  data = NULL,
  seed = NULL,
  ...
)

Arguments

Details

The comparator is descriptive: it reports per-method posterior / point CATE estimates together with their native CIs (when the adapter exposes one) and three concordance metrics (RMSE, Pearson correlation, mean absolute discrepancy) between every ordered pair of methods over cate_mean. Tests of hypothesis and claims of algorithmic equivalence are deliberately out of scope (the inferential origin of each method differs); the interpretation of the discrepancy is left to the user.

Value

An object of class gdpar_meta_learner_comparison with components bridge_cate, external, comparison, newdata, level, n_obs, n_methods, call, meta. See print.gdpar_meta_learner_comparison and summary.gdpar_meta_learner_comparison for the companion S3 methods.

Scalar-outcome restriction

The current scope of Sub-phase 8.5.B supports scalar outcomes only. Bridges constructed from fits with K > 1 (distributional regression) or p > 1 (multivariate response) are rejected with gdpar_unsupported_feature_error; multi-output external adapters are queued for Block 9 (see vignette v08c_meta_learner_comparison, section "Limits").

Dataset reconstruction

When data is NULL, the helper .assemble_bridge_dataset recovers the training data of each arm via eval(fit$call$data, eval_env) (same mechanism used by the bridge constructor in gdpar_causal_bridge), identifies the outcome via the LHS of fit$call$formula, and assembles a single (X, T, Y) dataset with T = 1L for the treatment arm and T = 0L for the control arm. When the evaluations fail, the helper aborts with gdpar_input_error and instructs the user to pass data explicitly.

Concordance metrics

For every ordered pair of methods $(i, j)$ (including the bridge as a method indexed by "bridge"), the comparator reports

$\mathrm{RMSE}_{ij} = \sqrt{\mathrm{mean}((\widehat\tau_i - \widehat\tau_j)^2)}$

$\mathrm{Pearson}_{ij} = \mathrm{cor}(\widehat\tau_i, \widehat\tau_j)$

$\mathrm{MAD}_{ij} = \mathrm{mean}(\mid \widehat\tau_i - \widehat\tau_j \mid)$

computed on cate_mean only; CIs are not pooled across methods because the inferential origin of each CI is heterogeneous (posterior vs. asymptotic vs. bootstrap; see Appendix B of the bridge vignette).

References

Kuenzel, S. R., Sekhon, J. S., Bickel, P. J., and Yu, B. (2019). Metalearners for estimating heterogeneous treatment effects using machine learning. Proceedings of the National Academy of Sciences, 116(10), 4156-4165.

Athey, S., and Wager, S. (2019). Estimating treatment effects with causal forests: An application. Observational Studies, 5, 37-51.

Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018). Double/debiased machine learning for treatment and structural parameters. The Econometrics Journal, 21(1), C1-C68.

See Also

gdpar_causal_bridge, gdpar_meta_learner_adapter, gdpar_adapter_grf, gdpar_adapter_econml.

Examples

if (requireNamespace("cmdstanr", quietly = TRUE) &&
    requireNamespace("grf", quietly = TRUE)) {
  n <- 300L
  df <- data.frame(x1 = rnorm(2L * n))
  df$arm <- rep(c("treat", "ctrl"), each = n)
  df$y <- with(df, ifelse(arm == "treat", 0.5, 0) + 0.8 * x1 +
                   rnorm(2L * n, sd = 0.5))
  df_t <- subset(df, arm == "treat", select = -arm)
  df_c <- subset(df, arm == "ctrl",  select = -arm)
  fit_t <- gdpar(y ~ x1, amm = amm_spec(a = ~ x1), data = df_t,
                 iter_warmup = 200, iter_sampling = 200, chains = 2)
  fit_c <- gdpar(y ~ x1, amm = amm_spec(a = ~ x1), data = df_c,
                 iter_warmup = 200, iter_sampling = 200, chains = 2)
  bridge <- gdpar_causal_bridge(fit_t, fit_c,
                                newdata = data.frame(x1 = seq(-2, 2, 0.2)))
  cmp <- gdpar_compare_meta_learners(bridge,
                                     methods = list(gdpar_adapter_grf()))
  print(cmp); summary(cmp)
}

Description

Verify numerically the predicted posterior contraction rate (Theorem 4B of Block 4) for a fitted Path 1 model. Refits the model at multiple subsample sizes, records the median posterior credible interval width across the user-facing parameters, and fits the regression $\log(\text{width}) = \alpha + \beta \log(n)$ across the subsample sizes. The slope estimate is consistent with a parametric $n^{-1/2}$ contraction rate when its value is in the interval (-0.6, -0.4).

Usage

gdpar_contraction_diagnostic(
  fit,
  data,
  sizes = NULL,
  replicates = 1L,
  parameters = NULL,
  level = 0.95,
  iter_warmup = 500L,
  iter_sampling = 500L,
  chains = 2L,
  verbose = TRUE,
  ...
)

Arguments

Details

This function is opt-in and computationally expensive: it refits the model length(sizes) * replicates times. A cost message is printed at the start. The conclusions of gdpar are not affected by calling this function.

Theorem 4B of Block 4 establishes that, under the conditions of Theorem 4A plus the prior thickness condition (PRIOR-THICK) and the sieve condition (SIEVE), the posterior contracts at rate $\varepsilon_n$ with $n \varepsilon_n^2 \to \infty$. For finite-dimensional parametric AMM specifications, the rate is $n^{-1/2}$. This diagnostic checks the empirical slope of the log-width against log-n.

Deviations from the predicted slope can indicate (i) prior misspecification (the prior fails (PRIOR-THICK) at the true parameter), (ii) failure of the homogeneity (HOM) or regularity (REG) conditions of Block 2, or (iii) non-parametric components whose smoothness assumption does not match the truth. The diagnostic flags the discrepancy without diagnosing the cause.

Value

A list of class gdpar_contraction_report with components table (data frame with columns n, replicate, median_width), slope_estimate, slope_se, slope_ci_lower, slope_ci_upper, verdict (a character indicating whether the empirical slope is consistent with the parametric $n^{-1/2}$ rate), and warnings (character vector recording per-refit fallback notifications; empty when every refit succeeded). A print method provides a human-readable summary.

Methodological notes

This function is part of the methodological audit toolkit, not of the standard inference flow. It is computationally expensive: each subsample size requires a full refit. The default settings keep the refits short (500 + 500 iterations, 2 chains) to make the diagnostic affordable on moderate datasets; users with more computational budget should pass higher values via the relevant arguments.

Subsamples are drawn without replacement and stratified by row order; users with structured data (time series, clustered observations) should pass an explicit sizes vector that respects the structure if random subsampling is inappropriate.

Dependencies

Uses cmdstanr for the refits and posterior to extract credible-interval widths.

References

See vignette("v04_asymptotics_path1_bayesian", package = "gdpar"), Section 6.2 (numerical verification of contraction).

See Also

gdpar, gdpar_bvm_check

Examples

if (requireNamespace("cmdstanr", quietly = TRUE)) {
  df <- data.frame(x1 = rnorm(400), y = rnorm(400))
  fit <- gdpar(y ~ x1, amm = amm_spec(a = ~ x1), data = df,
               iter_warmup = 200, iter_sampling = 200, chains = 2)
  gdpar_contraction_diagnostic(fit, data = df, replicates = 1)
}

Description

Quantifies serial (temporal) dependence in the residuals of a fitted scalar Path 1 Empirical-Bayes model. gdpar assumes conditional independence; under temporal (or spatial) autocorrelation that assumption is violated and the model-based (posterior / Laplace) uncertainty is too narrow. This diagnostic makes the violation *visible and measurable* before any remedy is applied, and is the natural gate for gdpar_dependence_robust.

Usage

gdpar_dependence_diagnostic(
  object,
  index = NULL,
  residual_type = c("quantile", "response", "pearson", "deviance"),
  max_lag = NULL,
  level = 0.95,
  randomize_seed = NULL,
  ...
)

Arguments

Details

The Durbin-Watson statistic is reported descriptively as $DW = \sum_{t=2}^{n}(r_t - r_{t-1})^2 / \sum_{t=1}^{n} r_t^2$ ($DW \approx 2(1 - \hat{\rho}_1)$); values near 2 indicate no first-order autocorrelation. The Ljung-Box test (stats::Box.test) provides the omnibus p-value across lags and drives the verdict. The Ljung-Box degrees of freedom are not reduced by the number of estimated AMM coefficients (fitdf = 0); for residuals of a fitted model this makes the test mildly optimistic, a caveat stated honestly rather than masked.

Spatial dependence (Moran's I and a spatial weight structure) is handled by the sibling gdpar_spatial_dependence_diagnostic. Both the scalar Empirical-Bayes (gdpar_eb_fit) and the scalar full-Bayes (gdpar_fit) paths are supported (decision D102); the residuals are the Bayesian Dunn-Smyth residuals of whichever fit is supplied. The K > 1 / p > 1 paths remain deferred.

Value

A list of class gdpar_dependence_diagnostic with components residual_type, n, max_lag, lag1_autocorr, lag1_p_value (normal approximation $\sqrt{n}\,\hat{\rho}_1 \sim N(0, 1)$), durbin_watson, ljung_box_statistic, ljung_box_df, ljung_box_p_value, level, index_supplied and verdict. A print method provides a human-readable summary.

Methodological note

A flagged verdict says the model-based uncertainty is not trustworthy under the detected dependence, not that the point estimates are wrong. The companion remedy gdpar_dependence_robust re-estimates the uncertainty by a temporal block bootstrap.

References

Ljung, G. M. & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika 65(2), 297-303.

Durbin, J. & Watson, G. S. (1950). Testing for serial correlation in least squares regression. I. Biometrika 37(3/4), 409-428.

See Also

gdpar_dependence_robust, gdpar_eb

Examples

if (requireNamespace("cmdstanr", quietly = TRUE) &&
    requireNamespace("posterior", quietly = TRUE)) {
  n <- 100
  x <- rnorm(n)
  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 = 100, iter_sampling = 100)
  gdpar_dependence_diagnostic(fit, index = df$t)
}

Description

Re-estimates the uncertainty of a scalar Path 1 Empirical-Bayes fit so that it is robust to temporal (serial) dependence in the data, without modelling that dependence. It refits the model on B moving (or circular) block bootstrap resamples of the data ordered by index, and reports the bootstrap standard deviation and percentile intervals of each AMM coefficient alongside the model-based (Laplace / posterior) standard errors. This is the working-independence + robust-variance stance of Liang & Zeger (1986): the point estimates are unchanged (consistent when the mean structure is correct, not efficient), only the reported uncertainty is made dependence-robust.

Usage

gdpar_dependence_robust(
  object,
  data,
  index = NULL,
  block_length = NULL,
  residual_type = c("quantile", "response", "pearson", "deviance"),
  randomize_seed = NULL,
  type = c("moving", "circular"),
  B = 199L,
  level = 0.95,
  seed = NULL,
  iter_warmup = 500L,
  iter_sampling = 500L,
  chains = 2L,
  verbose = TRUE,
  ...
)

Arguments

Value

A list of class gdpar_dependence_robust with components table (data frame with one row per coefficient and columns estimate, model_se, robust_se, se_ratio, ci_lower, ci_upper), block_length, block_length_method ("rate", "fixed" or "auto"; "rate" also flags an "auto" request that fell back), type, B, B_ok (successful refits), level, index_supplied, seed, warnings and refit_diagnostics (aggregate per-refit convergence: max_rhat, min_ess_bulk, n_divergent_refits, n_high_rhat_refits). A print method provides a human-readable summary.

Honest scope

The bootstrap delivers robust variance, not better point estimates, and is valid for weak / short-range dependence relative to block_length; it does not rescue long-memory or unit-root processes. gdpar does not model the dependence (that is deferred to a future block); here it only makes its inference robust to it.

Empirical-Bayes and full-Bayes paths

Both the scalar Empirical-Bayes (gdpar_eb_fit) and the scalar full-Bayes (gdpar_fit) paths are supported (decision D102), through a single shared engine; only the per-fit extraction of the point estimate, the model SE and the residuals is class-dispatched. On the EB path the point estimate / model SE are the Laplace / conditional-posterior mean and SD; on the full-Bayes path they are the posterior mean and posterior SD of each AMM coefficient (theta_ref, a_coef, b_coef, W_raw, the latter on its raw scale, for parity with the EB extractor). In both cases robust_se is the block-bootstrap SD of the per-refit point estimate and se_ratio = robust_se / model_se is a like-for-like SD-vs-SD ratio: it contrasts the dependence-robust sampling variability of the point estimate against the within-model (posterior / Laplace) SD, and exceeds 1 when the latter understates the former. The posterior mean / SD choice (rather than median / IQR) preserves this parity and avoids an undeclared normal-scaling constant. (For a strongly skewed coefficient posterior – W_raw under sparse data or strong shrinkage is the usual culprit – the posterior mean can sit off the posterior mode; inspect object$fit$draws() directly in that case. The reported point estimate is always the posterior mean.)

Three honest full-Bayes caveats. (1) Each refit re-runs the full HMC (markedly more costly than an EB refit). (2) A finite-iteration refit carries Monte-Carlo error in its posterior mean that slightly and conservatively inflates robust_se (reducible by a larger iter_sampling; the aggregate refit ESS is reported in refit_diagnostics). (3) The se_ratio has a subtly different reading across paths: under an informative prior the full-Bayes posterior SD can be smaller than the bootstrap SD even under correct independent specification, because the prior concentrates the posterior beyond what the data alone support, giving se_ratio < 1. That is benign prior regularization, not the model SE overstating uncertainty; only se_ratio clearly above 1 signals dependence / misspecification (the analytic conjugate-Gaussian check gives $\mathrm{se\_ratio}^2 = n\tau / (n\tau + \tau_0) < 1$ with prior precision $\tau_0$). Note too that the EB and full-Bayes theta_ref point estimates are different estimands – the Laplace mode (EB) versus the posterior mean (full Bayes) – which coincide asymptotically (Bernstein-von Mises) but may differ in finite samples; their se_ratio values are therefore not expected to match to the last digit when the same data are run through both paths. A widened / bagged posterior (BayesBag; Huggins & Miller 2019) is a different object – a re-architected estimator rather than a robust variance for the same one – and is a documented deferred lateral, not adopted here.

Dependencies

Uses cmdstanr for the refits and posterior to extract the coefficient estimates (Empirical-Bayes or full-Bayes).

References

Liang, K.-Y. & Zeger, S. L. (1986). Longitudinal data analysis using generalized linear models. Biometrika 73(1), 13-22.

Kuensch, H. R. (1989). The jackknife and the bootstrap for general stationary observations. Annals of Statistics 17(3), 1217-1241.

Politis, D. N. & White, H. (2004). Automatic block-length selection for the dependent bootstrap. Econometric Reviews 23(1), 53-70.

Patton, A., Politis, D. N. & White, H. (2009). Correction to "Automatic block-length selection for the dependent bootstrap". Econometric Reviews 28(4), 372-375.

See Also

gdpar_dependence_diagnostic, gdpar_eb

Examples

if (requireNamespace("cmdstanr", quietly = TRUE) &&
    requireNamespace("posterior", quietly = TRUE)) {
  n <- 100
  x <- rnorm(n)
  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 = 100, iter_sampling = 100)
  # B is kept small here for a fast example; use B >= 199 in practice.
  gdpar_dependence_robust(fit, data = df, index = df$t, B = 10,
                          seed = 1, iter_warmup = 100,
                          iter_sampling = 100, chains = 2)
  # Data-driven block length (Politis-White), opt-in:
  gdpar_dependence_robust(fit, data = df, index = df$t,
                          block_length = "auto", B = 10, seed = 1,
                          iter_warmup = 100, iter_sampling = 100, chains = 2)
}

Description

If DHARMa is available (Suggests), returns a DHARMa::DHARMa object built from the posterior predictive draws of y_pred and the observed response y. The user can then call DHARMa::testResiduals(), DHARMa::testZeroInflation(), DHARMa::plotResiduals() and other DHARMa tests on the returned object.

Usage

gdpar_dharma_object(object, coord = NULL)

Arguments

Value

A DHARMa simulation object.

Dependencies

Requires DHARMa (Suggests). If the package is not installed the function raises gdpar_input_error pointing to residuals(., type = "quantile") as the built-in fallback.

Description

Path 1 Empirical-Bayes counterpart of gdpar: estimates the population reference $\theta_{ref}$ by maximizing the marginal likelihood (Type II ML) and samples the lower-level parameters $\xi = (a, b, W, \sigma_*, \phi)$ from the conditional posterior given $\widehat\theta_{ref}^{EB}$. See vignette("v07_eb_vs_fb", package = "gdpar") for the Fully-Bayes versus Empirical-Bayes asymptotic comparison ($p = 1$, $K = 1$) and vignette("v07b_eb_multivariate", package = "gdpar") for the multivariate extension that motivates the family of dedicated templates introduced in this sub-phase.

Usage

gdpar_eb(
  formula,
  family = gdpar_family("gaussian"),
  amm = amm_spec(),
  W = NULL,
  data,
  prior = NULL,
  anchor = "prior_mean",
  skip_id_check = FALSE,
  chains = 4L,
  iter_warmup = 1000L,
  iter_sampling = 1000L,
  adapt_delta = 0.95,
  max_treedepth = 12L,
  refresh = 100L,
  verbose = TRUE,
  seed = NULL,
  group = NULL,
  parametrization = c("auto", "ncp", "cp"),
  id_check_rigor = c("full", "fast"),
  eb_correction = TRUE,
  laplace_control = list(),
  ...
)

Arguments

Value

An object of class gdpar_eb_fit with components:

• theta_ref_hat: numeric vector of length J_groups with the EB point estimates.

• theta_ref_se: numeric vector of length J_groups with the marginal standard errors derived from the Laplace covariance.

• conditional_fit: the underlying cmdstanr fit object of Step (iii).

• correction_applied: logical scalar.

• eb_correction_constant: numeric scalar; the scalar Proposition 7B inflation constant when eb_correction = TRUE, NA_real_ otherwise.

• diagnostics_numerical: list with kappa, lm_perturbation, lm_n_iter, lm_status (one of "not_needed", "converged", "exhausted"), kappa_post_ridge, multi_start_dispersion, marginal_log_lik_history. For Path C (K > 1 and p > 1) the slot-vectorized counterparts replace the scalars: kappa_per_slot, lm_lambda_per_slot, lm_n_iter_per_slot, lm_status_per_slot.

• diagnostics: gdpar_diagnostics object from the conditional HMC fit (same shape as in gdpar_fit).

• amm, family, prior, design, anchor, stan_data, group_info, identifiability_report, parametrization, call, path ("eb").

Scope of Sub-phase 8.6.B

This version implements the base regime $K = 1$, $p = 1$ (single distributional slot per group, scalar $\theta_{ref}$). Multivariate $p > 1$ and multi-slot $K > 1$ are explicitly rejected by the input guard with gdpar_unsupported_feature_error; their habilitation is canonized in Charter Sub-phases 8.6.C and 8.6.D and exercises the K_slots / p_dim fields declared in the two dedicated templates inst/stan/amm_eb_marginal.stan and inst/stan/amm_eb_conditional.stan.

Dependencies

Uses cmdstanr (Suggests) for Laplace approximation (cmdstanr::laplace()) in Step (i) and HMC sampling in Step (iii); posterior (Imports) for diagnostics. The Laplace approximation requires the laplace() method of cmdstanr, available since cmdstanr 0.7.0.

References

Carlin, B. P., and Gelfand, A. E. (1990). Approaches for empirical Bayes confidence intervals. JASA 85(409), 105–114.

Petrone, S., Rousseau, J., and Scricciolo, C. (2014). Bayes and empirical Bayes: do they merge? Biometrika 101(2), 285–302.

Rousseau, J., and Szabo, B. (2017). Asymptotic behaviour of the empirical Bayes posteriors associated to maximum marginal likelihood estimator. Annals of Statistics 45(2), 833–865.

See Also

gdpar, amm_spec, gdpar_family, gdpar_prior.

Examples

if (requireNamespace("cmdstanr", quietly = TRUE)) {
  set.seed(NULL)
  n <- 200
  df <- data.frame(
    x1 = stats::rnorm(n),
    x2 = stats::rnorm(n)
  )
  df$y <- with(df, 1 + 0.5 * x1 - 0.3 * x2 + stats::rnorm(n, sd = 0.5))
  spec <- amm_spec(a = ~ x1 + x2)
  fit_eb <- gdpar_eb(
    formula       = y ~ x1 + x2,
    family        = gdpar_family("gaussian"),
    amm           = spec,
    data          = df,
    iter_warmup   = 200,
    iter_sampling = 200,
    chains        = 2
  )
  print(fit_eb)
}

Description

Define the response distribution that links the individual parameter theta_i to the observed outcome y_i. The family object carries the link function, the inverse link, the metadata for the parameter identifiability condition (D-ID) of Lemma 1B in Block 1, and the family identifier consumed by the Stan code generator.

Usage

gdpar_family(
  name = c("gaussian", "poisson", "neg_binomial_2", "bernoulli", "beta", "gamma",
    "student_t", "tweedie", "zip", "zinb", "hurdle_poisson", "hurdle_neg_binomial_2"),
  link = NULL,
  did_override = NULL
)

Arguments

Details

The Path 1 implementation in this version of the package fits the AMM canonical form on the linear-predictor scale of the family. The inverse link is applied at the likelihood block in Stan; the centering, anchoring and prior specifications all live on the linear-predictor scale, consistent with assumptions (C1)-(C6) of Block 1 and the prior conditions (PRIOR-KL), (PRIOR-THICK) of Block 4.

Built-in families and their D-ID status:

"gaussian"

Identity link by default. D-ID holds unconditionally for the location parameter when the variance is identifiable from the data.

"poisson"

Log link by default. D-ID holds unconditionally on the rate parameter.

"neg_binomial_2"

Log link by default; Stan's neg_binomial_2 parametrization with mean mu and dispersion phi such that variance equals mu + mu^2 / phi. D-ID holds unconditionally for mu when phi is identifiable from the data.

"bernoulli"

Logit link by default. D-ID holds unconditionally on the success probability.

"beta"

Logit link by default for the mean mu on (0, 1); precision phi on log link. Stan parametrization beta_proportion(mu, phi) with variance mu*(1-mu)/(1+phi). D-ID holds for mu unconditionally; for phi when the empirical dispersion identifies it (sub-phase 8.3.4 of Block 8).

"gamma"

Log link by default for the mean mu on positive reals; shape on log link. Stan parametrization gamma(shape, shape/mu) with variance mu^2/shape. D-ID holds for mu unconditionally; for shape when the empirical dispersion identifies it (sub-phase 8.3.4 of Block 8).

"student_t"

Identity link by default for the location mu on the real line; sigma on log link; nu (degrees of freedom) on log link. Stan parametrization student_t(nu, mu, sigma) with mean mu (for nu > 1) and variance sigma^2 * nu / (nu - 2) (for nu > 2). D-ID holds for mu unconditionally; for sigma when the empirical residual dispersion identifies it; for nu when the empirical tail thickness identifies it. Sub-phase 8.3.5a of Block 8 wires the family exclusively under the K = 3 distributional regression path of amm_distrib_K.stan; routings with K < 3 (population-scoped sigma and / or nu) are deferred.

"tweedie"

Log link by default for the mean mu on positive reals; dispersion phi on log link; power p on identity link with support on the open interval (1.01, 1.99) (canonical compound Poisson-gamma regime, with point mass at zero and continuous positive body). Stan parametrization tweedie(mu, phi, p) with mean mu and variance phi * mu^p. The log-pdf is not native in Stan and is implemented in the model's functions block as a hybrid: the Dunn–Smyth (2005) infinite-series expansion in the central region |p - 1.5| < tau (tau = 0.4 by default) and the saddlepoint approximation elsewhere, with the same dispatch applied to the random-number generator. D-ID holds for mu unconditionally; for phi when the empirical dispersion identifies it; for p when the empirical balance between zero mass and continuous body identifies it. Sub-phase 8.3.5b of Block 8 wires the family exclusively under the K = 3 distributional regression path of amm_distrib_K.stan.

"zip"

Zero-inflated Poisson. Log link by default for the count mean mu on positive reals; mixture probability pi on logit link with support on (0, 1). Likelihood: with probability pi the outcome is a structural zero, otherwise it is drawn from Poisson(mu). Stan implementation via log_sum_exp(bernoulli_logit_lpmf(1 | eta_pi), bernoulli_logit_lpmf(0 | eta_pi) + poisson_log_lpmf(0 | eta_mu)) for y = 0 and bernoulli_logit_lpmf(0 | eta_pi) + poisson_log_lpmf(y | eta_mu) for y > 0. D-ID holds for mu unconditionally; for pi when the empirical proportion of structural zeros, distinct from the sampling zeros of the count component, is sufficient to identify the zero-inflation probability. Sub-phase 8.3.6 of Block 8 wires the family exclusively under the K = 2 distributional regression path of amm_distrib_K.stan; routings with K < 2 are deferred.

"zinb"

Zero-inflated negative binomial. Log link by default for the count mean mu on positive reals; dispersion phi on log link; mixture probability pi on logit link. Likelihood: with probability pi the outcome is a structural zero, otherwise it is drawn from NegBinomial2(mu, phi) with variance mu + mu^2 / phi. D-ID holds for mu unconditionally; for phi when the empirical overdispersion identifies it; for pi as in the zero-inflated Poisson case. Sub-phase 8.3.6 wires the family exclusively under the K = 3 distributional regression path; routings with K < 3 are deferred.

"hurdle_poisson"

Hurdle Poisson. Log link by default for the truncated-count mean mu on positive reals; hurdle probability pi on logit link. Likelihood: a Bernoulli draw with probability pi decides whether the outcome equals zero or is strictly positive; the positive branch is drawn from a Poisson truncated at one. Stan implementation via bernoulli_logit_lpmf(1 | eta_pi) for y = 0 and bernoulli_logit_lpmf(0 | eta_pi) + poisson_log_lpmf(y | eta_mu) - log1m_exp(-exp(eta_mu)) for y > 0. Distinct from the zero-inflated Poisson in that the zero mass is a structural decision, not a mixture between a structural zero and a Poisson sampling zero. D-ID holds for mu unconditionally; for pi when the empirical proportion of zeros is sufficient to identify the hurdle probability. Sub-phase 8.3.6 wires the family exclusively under the K = 2 distributional regression path.

"hurdle_neg_binomial_2"

Hurdle negative binomial. Log link by default for the truncated-count mean mu; dispersion phi on log link; hurdle probability pi on logit link. Likelihood: Bernoulli draw decides zero vs. positive; the positive branch is drawn from a NegBinomial2(mu, phi) truncated at one. D-ID holds for mu unconditionally; for phi as in negative binomial; for pi as in hurdle Poisson. Sub-phase 8.3.6 wires the family exclusively under the K = 3 distributional regression path.

Value

An object of class gdpar_family with components name, link, inv_link, linkfun, stan_id, has_dispersion, did_status, did_condition and did_reference. A print method provides a human-readable summary.

Methodological notes

Identifiability of the response family in its parameter is a structural property of the model, not a property that can be tested from finite data. The package therefore documents this property through the did_status field but does not attempt to verify it at fitting time. When the family carries a conditional D-ID status, the fitting routine emits an informative message reminding the user of the relevant condition. See Lemma 1B and the commentary in Block 1, Section 6.4.

The returned family declares every structural parameter the distribution admits as eligible for an individual specification (e.g., mu and sigma for Gaussian; mu and phi for negative binomial; mu for Poisson and Bernoulli). Each eligible parameter carries a canonical scope in its gdpar_param_spec: the location parameter is per_observation (it is modeled through the AMM canonical form) and auxiliary parameters default to population. K-individual membership (which auxiliaries are promoted to per_observation) is declared exclusively at the entry of gdpar via a gdpar_formula_set (for the high-level formula path) or via a named list of amm_spec objects (for the low-level path). Parameters that are not named in that entry keep their canonical scope and are estimated as population-level constants. The family is the registry of eligibles, and the entry to gdpar() is the single source of truth for K.

Dependencies

This function uses stats family objects (gaussian, poisson, binomial) for the link metadata.

References

See vignette("01_amm_identifiability", package = "gdpar"), Section 6.4 (Lemma 1B) for D-ID; Lambert (1992), Greene (1994) and Mullahy (1986) for identifiability of zero-inflated and hurdle counts (sub-phase 8.3.6 of Block 8).

See Also

gdpar_family_custom, gdpar

Examples

fam <- gdpar_family("poisson")
print(fam)

Description

Build a user-defined family for use with gdpar when the built-in families of gdpar_family do not cover the application. The user is responsible for declaring whether the identifiability condition (D-ID) of Lemma 1B in Block 1 holds for the family.

Usage

gdpar_family_custom(
  name,
  link,
  did_holds,
  did_condition,
  stan_loglik_block,
  stan_log_lik_block,
  stan_y_pred_block,
  y_type,
  did_reference
)

Arguments

Details

Building a custom family is an advanced use of the package. The user assumes the responsibility of ensuring (i) that the Stan likelihood block is mathematically correct and (ii) that the family is identifiable in its parameter. The package emits an informative message when the custom family is created, restating these responsibilities.

Value

An object of class gdpar_family.

Methodological notes

The package never attempts to test identifiability from data; it only registers the user's declaration. See Lemma 1B in Block 1, Section 6.4.

Dependencies

None beyond the base R installation.

See Also

gdpar_family

Examples

my_family <- gdpar_family_custom(
  name               = "my_log_normal",
  link               = "log",
  did_holds          = TRUE,
  did_condition      = NA_character_,
  stan_loglik_block  =
    "target += normal_lpdf(log(y_real[i]) | eta[i], sigma_y[1]);",
  stan_log_lik_block =
    "log_lik[i] = normal_lpdf(log(y_real[i]) | eta[i], sigma_y[1]);",
  stan_y_pred_block  =
    "y_pred[i] = exp(normal_rng(eta[i], sigma_y[1]));",
  y_type             = "real",
  did_reference      = "User declaration"
)
print(my_family)

Description

Build a custom distributional regression family by selecting a canonical bi-parametric likelihood pattern from the registry returned by .gdpar_K_custom_patterns. The constructor is the descriptor-based wiring agreed in sub-phase 8.3.4 of Block 8 (D-A3.B; option (b) of the scoping): the user does not supply free-form Stan code; they choose a stan_lpdf_id from the whitelist and the family routes through the same amm_distrib_K.stan branch that the built-in distributional regression families use.

Usage

gdpar_family_custom_K(
  name,
  stan_lpdf_id,
  did_holds = TRUE,
  did_condition = NULL,
  did_reference = NULL
)

Arguments

Details

This is the K = 2 sibling of gdpar_family_custom (the K = 1 free-form custom path). The two coexist with distinct contracts: gdpar_family_custom() accepts user-authored Stan code blocks at K = 1 (user-validated identifiability and correctness); gdpar_family_custom_K() accepts only registered canonical patterns at K = 2 (descriptor-validated, no Stan injection surface).

The descriptor path enforces structural validation: name must be unique, stan_lpdf_id must be in the registry, and the slot configuration (links, supports, prior kinds) is fixed by the registry entry. Users who need to deviate from the canonical parametrization must request a new entry in the registry; this keeps the Stan-side surface auditable and bit-exact across versions of the package.

Sub-phase 8.3.4 (D-A3.B option (b)) decided this descriptor approach over (i) plug-in of free-form Stan code per family slot and (ii) a hybrid descriptor + escape-hatch. The decision applies [[feedback-max-robustness-priority]] (structural validation over surface flexibility) at the cost of forcing the user to contribute upstream when a new K = 2 family is needed.

Value

An object of class gdpar_family with two param_specs (slot 1 = location, slot 2 = scale), the stan_id of the pattern (e.g. 7 for lognormal_loc_scale), and is_custom = TRUE. The family routes through amm_distrib_K.stan via the family_id_k dispatcher.

Methodological notes

The K = 2 custom family inherits the canonical likelihood, links, and priors of the chosen registry pattern. The user only chooses the family name and (optionally) overrides the D-ID metadata. The AMM canonical form, the per-slot scope, and the dispatch in amm_distrib_K.stan are the same as for the built-in gdpar_family families wired in 8.3.4.

See Also

gdpar_family, gdpar_family_custom

Examples

my_lognorm <- gdpar_family_custom_K(
  name          = "my_lognormal_K2",
  stan_lpdf_id  = "lognormal_loc_scale",
  did_holds     = TRUE,
  did_reference = "User declaration"
)
print(my_lognorm)

Description

Build a per-coordinate family object for use with gdpar when the AMM specification has dimension p > 1L. The resulting object declares one univariate family per coordinate of the individual parameter vector $\theta_i \in \mathbb{R}^p$; the likelihood factorizes across coordinates as

$p(y_i \mid \theta_i) = \prod_{k=1}^{p} D_k(y_{ik} \mid \theta_i[k]),$

with cross-dimensional coupling carried exclusively by the modulating component $W(\theta_{\mathrm{ref}})$ of the AMM canonical form.

Usage

gdpar_family_multi(family, p, link = NULL)

Arguments

Details

The factorization is the canonization of architectural Option B of the Phase F decision (handoff 10 continuation, 2026-05-11): $y_i \in \mathbb{R}^p$ is multivariate with marginals independent conditional on $\theta_i$; cross-dimensional dependence enters the model through the coupling of $\theta_i[k]$ via $W(\theta_{\mathrm{ref}})$. Multi-parametric families (a single univariate outcome parametrized by the whole vector $\theta_i \in \mathbb{R}^p$, e.g., gaussian with $\theta_i = (\mu_i, \log \sigma_i)$ in the distributional regression sense) are deferred to a dedicated post-validation block; see the project memory entry project_gdpar_multiparametric_extension_postvalidation.

Value

An object of class gdpar_family_multi with components families (a list of p gdpar_family objects), p, homogeneous, stan_id (the common Stan family identifier when homogeneous), has_dispersion (the common dispersion flag), name (the common family name) and did_status (the common identifiability status). A print method provides a human-readable summary.

Methodological notes

The identifiability condition (D-ID) of Lemma 1B in Block 1 applies coordinate-wise under this factorization: each univariate marginal $D_k$ identifies $\theta_i[k]$ from $y_{ik}$ independently. The cross-dimensional identifiability condition (C4-bis), which guards against aliasing between coordinates of $\theta_{\mathrm{ref}}$ that share basis structure, is checked by gdpar_check_identifiability (Phase H pending).

Dependencies

Calls gdpar_family when family is supplied as a name.

See Also

gdpar_family, amm_spec

Examples

fam_mv <- gdpar_family_multi("gaussian", p = 2L)
print(fam_mv)

fam_mv2 <- gdpar_family_multi(gdpar_family("poisson"), p = 3L)
print(fam_mv2)

Description

Build the canonical internal representation that gdpar consumes when more than one structural parameter of the family is modeled with an AMM design (multi-parametric distributional regression in the sense of decision D-F1 of Block 8 Session 1, materialized in sub-phase 8.3.3 of the package). One slot per individual parameter, named with the canonical parameter name as declared by the family's param_specs; the first slot carries the outcome variable on its left-hand side, while subsequent slots are one-sided formulas that describe the AMM design of an auxiliary parameter (e.g., the Gaussian sigma or the negative-binomial phi).

Usage

gdpar_formula_set(...)

Arguments

Details

This constructor is the canonical low-cost entry point; the brms-style sugar gdpar_bf produces an equivalent object from a sequence of two-sided formulas.

Value

An object of class gdpar_formula_set with components outcome (character scalar with the outcome variable name), formulas (named list of formula objects; the first is two-sided, subsequent ones are one-sided), param_names (character vector of slot names, identical to names(formulas)) and env (the environment of the first formula, used downstream for evaluation).

Validation contract

• All formulas must be named with the canonical parameter name (no positional formulas).

• The first slot must be a two-sided formula whose left-hand side is a single symbol naming the outcome variable (e.g., mu = y ~ a(x)). Function calls on the LHS such as log(y) are rejected to keep the contract explicit; pre-transform the outcome in the data frame instead.

• Subsequent slots must be one-sided formulas (e.g., sigma = ~ a(x)).

• Slot names must be unique.

• No formula may suppress the intercept via -1 or +0: the AMM population anchor $\theta_{ref}$ is a structural term that cannot be removed. The constructor aborts with an informative gdpar_input_error when this pattern is detected.

Set membership against the family's eligible param_specs (slot names must be a subset of the eligible parameters of the family passed to gdpar) is enforced downstream by gdpar once both the formula set and the family are known.

Methodological notes

Block 8 Session 1 decision 1C established that a family is a list of parameter specifications. Block 8.3.3 closes the high-level API: the formula set is the single source of truth for the K-individual parameter set; any parameter named in this object is promoted to scope = "per_observation" in the family copy that gdpar uses internally. Parameters of the family that are not named in the formula set remain with their canonical scope (typically population). See memory entries project_gdpar_block_8_3_extended_plan and project_gdpar_block_8_session_1_decisions (2026-05-20).

Dependencies

Uses terms to detect intercept suppression.

See Also

gdpar_bf, gdpar, gdpar_family

Examples

fs <- gdpar_formula_set(
  mu    = y ~ a(x1) + b(z1),
  sigma = ~ a(x2)
)
print(fs)
fs[["mu"]]
names(fs)

Description

Turns an already-fitted gdpar object into the inputs the geometry-adaptive controller gdpar_geom_orchestrate consumes, without touching the fit path. This is the durable, path-agnostic core of the Block RG integration (RG.6 part ii): it reads the compiled cmdstan model and the Stan data carried by the fit, exposes the standalone log_prob / grad_log_prob / hessian methods, derives the unconstrained dimension and a posterior-mean warm-start, and packages a target (a re-samplable cmdstan model for the size-invariant diagnostic) together with a geom_target (the engine sampling target on the unconstrained scale).

Usage

gdpar_geom_bridge(
  object,
  fisher = NULL,
  reference = NULL,
  hessian = TRUE,
  methods_seed = 1L,
  ...
)

Arguments

Details

The orchestrator needs two things: a target it can re-sample for the diagnostic pilots, and an engine target exposing the unconstrained log-density and its gradient (and Hessian, for the Riemannian / SoftAbs level). A fitted gdpar object carries the posterior draws but its CmdStanMCMC object cannot be re-sampled. The bridge therefore recompiles a fresh CmdStanModel from the fit's own Stan source ($code()) with compile_model_methods = TRUE (cmdstanr's content-hash cache makes this a cache hit when the methods variant already exists), which serves both consumers, and reads the unconstrained dimension and posterior mean from the fit's draws.

The bridge is path-agnostic: it works for any gdpar_fit that carries $fit (a CmdStanMCMC) and $stan_data – the K-individual, multi-coordinate and single-coordinate paths alike. It is the tool RG.7 points at the real Tweedie count of benchmark 9.2.O.

Nothing here modifies gdpar; the returned object is plain data plus closures, so the default fit branch stays bit-identical and the goldens are untouched.

Value

An object of class gdpar_geom_bridge: a list with target (the cmdstan-model diagnostic target), geom_target (the engine target), fisher, reference, dim (the unconstrained dimension), model (the methods-enabled CmdStanModel) and stan_data. Feed these to gdpar_geom_orchestrate.

See Also

gdpar_geom_orchestrate, gdpar_geom_fit, gdpar_geom_target, gdpar_geom_fisher_simulator.

Examples

if (requireNamespace("cmdstanr", quietly = TRUE)) {
  set.seed(1)
  n <- 80
  x <- rnorm(n); z <- rnorm(n)
  y <- rnorm(n, 0.5 + 0.8 * (x - mean(x)), exp(-0.2 + 0.4 * (z - mean(z))))
  d <- data.frame(y = y, x = x, z = z)
  fit <- gdpar(gdpar_bf(y ~ a(x), sigma ~ a(z)), data = d,
               family = gdpar_family("gaussian"), chains = 1,
               iter_warmup = 200, iter_sampling = 200, refresh = 0,
               seed = 1, skip_id_check = TRUE, verbose = FALSE)
  bridge <- gdpar_geom_bridge(fit)
  bridge$dim
  # res <- gdpar_geom_orchestrate(bridge$target, bridge$geom_target,
  #                               reference = bridge$reference)
}

Description

Build the fisher(theta) function that feeds the learned metric gdpar_geom_metric_gp_fisher where the expected Fisher information has no closed form – the general realisation completing the Riemannian level of the Block RG hierarchy. The expected Fisher is $I(\theta) = \mathbb{E}_{y\sim p(\cdot\mid\theta)}\!\left[s(\theta, y)\, s(\theta, y)^\top\right]$ with $s$ the score (the gradient of the log-likelihood). The estimator is the average outer product of the scores over n_sim data sets simulated from the model at $\theta$: it is positive semi-definite by construction (a sum of rank-one outer products) and unbiased for the expected Fisher.

Usage

gdpar_geom_fisher_simulator(target, n_sim = 64L, seed = 1L, floor = 1e-08)

Arguments

Details

The estimate is a deterministic function of $\theta$: each evaluation reseeds the RNG from a key derived from $\theta$ and the base seed (and the RNG state is saved and restored around the call), so the surrogate trained on it is reproducible and independent of call order. The simulation noise is then absorbed by the Gaussian-process surrogate: its SoftAbs mean carries the bulk of the curvature deterministically and the process learns only the smooth residual, so the SoftAbs acts as a structural control variate and few replicates per site suffice (the deferred half of decision D91, now closed). Antithetic simulation is deliberately not used: for location families the score is odd in the centred data, so the antithetic partner has score $-s$ and $(-s)(-s)^\top = s s^\top$ leaves the outer product unchanged – no variance reduction.

For a well-conditioned estimate take n_sim comfortably larger than the dimension (the average of fewer than dim rank-one terms is singular and is only made strictly positive-definite by the eigenvalue floor).

Value

A function fisher(theta) returning the dim x dim estimated expected Fisher matrix (symmetric positive-definite), with the number of simulations recorded in its "n_sim" attribute. Pass it as the fisher argument of gdpar_geom_metric_gp_fisher.

See Also

gdpar_geom_metric_gp_fisher, gdpar_geom_target, gdpar_geom_rmhmc_adaptive.

Examples

# A bivariate normal location model y ~ N(theta, Sigma0): the expected Fisher
# is the constant precision Sigma0^{-1}, which the estimator recovers.
Sigma0 <- matrix(c(1, 0.3, 0.3, 2), 2, 2)
P0 <- solve(Sigma0)
tgt <- gdpar_geom_target(
  log_prob = function(theta) -0.5 * sum(theta^2),
  grad_log_prob = function(theta) -theta, dim = 2,
  simulate = function(theta)
    as.numeric(theta + t(chol(Sigma0)) %*% stats::rnorm(2)),
  score = function(theta, y) as.numeric(P0 %*% (y - theta)))
fisher <- gdpar_geom_fisher_simulator(tgt, n_sim = 4000, seed = 1)
round(fisher(c(0, 0)), 2)   # close to solve(Sigma0)