# Theoretical Addendum – Block 8 (Sub-phases 8.3.1–8.3.10):

### Per-Slot Parametric Resolution, Identifiability under Heterogeneity, Custom Families, \(K \geq 1\) Distributional Slots, Mixtures and Hurdles, B-spline Basis, Residual Diagnostics, and S3 Methods

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

#### 2026-07-06

# 1. Purpose and Relation to the v08* Family

The Block 8 vignette family canonizes the AMM extensions that live around but outside the Block 1 core (Asymmetric Multiplicative Model with \(\theta_i = \theta_{\mathrm{ref}} + \Delta(x_i, \theta_{\mathrm{ref}})\) and scalar-slot Gaussian response):

- v08 positions AMM relative to the CATE / ITE literature (Block 8 conceptual frame).
- v08b implements the causal bridge T-learner standalone (Sub-phase 8.5.A).
- v08c introduces the meta-learner comparator as an external bench (Sub-phase 8.5.B).
- This vignette (v08d) canonizes the internal per-slot parametric extension of AMM, summarized as Sub-phases 8.3.1 to 8.3.10 of the development log.

The decisions canonized here previously lived only in development handoffs and project memories. The release-gate of Sub-phase 8.3.10 (development version `0.0.0.9001`) closed the corresponding R-side and Stan-side implementation cycle, but the theoretical canonization — what each sub-phase decided, why, and how it composes with the others — was deferred to this addendum.

The structure follows the three-layer rigor standard used throughout the gdpar documentation:

- (L1) Algebraic / structural. What the per-slot architecture is: \(K\) distributional slots, \(K \geq 1\), each carrying its own family, link, parameter set, and identifiability constraint.
- (L2) Statistical. What it means to fit a model with \(K \geq 1\) heterogeneous slots (priors, sampling, posterior identifiability).
- (L3) Numerical / operational. How the codegen, the per-slot inverse-link dispatch, the B-spline basis dispatch, and the residual diagnostics actually compose in `R/` and `inst/stan/`.

This vignette does not claim the per-slot architecture is unique or novel. The design choices are documented for auditability (what was decided, with which alternatives present, under which constraint).

# 2. Setting and Notation

## 2.1. AMM with \(K \geq 1\) Distributional Slots

Let the response \(Y_i \in \mathcal{Y}\) for observation \(i = 1, \dots, n\) follow a conditional law indexed by a \(K\)-tuple of slot parameters:

\[Y_i \mid X_i \;\sim\; \mathcal{D}\bigl(\eta_{i,1}, \dots, \eta_{i,K}\bigr),\]

where each slot index \(k \in \{1, \dots, K\}\) carries its own natural (linear-predictor) scale \(\eta_{i,k} \in \mathbb{R}\), its own response (parameter) scale \(\theta_{i,k} = g_k^{-1}(\eta_{i,k})\) under an inverse-link \(g_k^{-1}\), and its own AMM decomposition:

\[\eta_{i,k} \;=\; \eta_{\mathrm{ref}, k} \;+\; \Delta_k(x_i, \eta_{\mathrm{ref}, k}),\]

with \(\eta_{\mathrm{ref}, k}\) the per-slot population reference and \(\Delta_k\) the per-slot AMM deviation. The scalar single-slot case \(K = 1\) recovers the canonical Block 1 setup; this is enforced as a corollary in §3.

The multivariate-reference dimension \(p\) (Sub-phase 8.6.A) is orthogonal to \(K\): \(p\) enumerates the components of \(\theta_{\mathrm{ref}}\) that an anchored hierarchical prior pools across grouping levels, whereas \(K\) enumerates the slots of the response distribution. The full \((p \geq 1, K \geq 1)\) regime — including the \(p > 1 \wedge K > 1\) Path C — is treated in v07b and v08; the present vignette focuses on the \(K\)-axis.

## 2.2. Slot Conventions

The framework’s slot convention is positional, mirroring the parameter order of the underlying density:

- Gaussian (\(K = 2\)): slot 1 = `mu` (identity link), slot 2 = `sigma` (log link).
- Poisson / NB (\(K = 1\) or \(K = 2\)): slot 1 = `mu` (log link); for NB, slot 2 = `phi` (log link).
- Beta (\(K = 2\)): slot 1 = `mu` (logit link, support \((0,1)\)), slot 2 = `phi` (log link, precision).
- Gamma (\(K = 2\)): slot 1 = `mu` (log link), slot 2 = `phi` (log link, shape under mean-shape parametrization).
- Student-\(t\) (\(K = 3\)): slot 1 = `mu` (identity), slot 2 = `sigma` (log), slot 3 = `nu` (log).
- Tweedie (\(K = 3\)): slot 1 = `mu` (log), slot 2 = `phi` (log), slot 3 = `p` (identity, bounded \((1.01, 1.99)\)).
- ZIP / Hurdle-Poisson (\(K = 2\)): slot 1 = `mu` (log), slot 2 = `pi` (logit).
- ZINB / Hurdle-NB (\(K = 3\)): slot 1 = `mu` (log), slot 2 = `phi` (log), slot 3 = `pi` (logit).

The per-slot conventions are not arbitrary: each link is the canonical link for the corresponding scale (location-on-real, scale-on-positive, probability-on-unit-interval), and the framework’s residual diagnostics (G1 / G2 / G3 of Sub-phase 8.3.9) rely on this convention for unit-coherent contributions.

# 3. Sub-phase 8.3.1 – `param_specs` as the Per-Slot Resolution Layer

The release-gate of Sub-phase 8.3 (Decision 1C, canonized in Session 1) introduces an internal data structure named `param_specs`: a list of length \(K\) whose \(k\)-th element is a `gdpar_param_spec` object recording, for slot \(k\):

- the name of the parameter (e.g., `"mu"`, `"sigma"`, `"phi"`, `"nu"`, `"p"`, `"pi"`),
- the link function and its inverse (canonical, dispatched by name),
- the integer `link_id` used by the Stan helper `apply_inv_link_by_id` (Sub-phase 8.3.7; see §9),
- the support of the parameter on the response scale (e.g., \(\mathbb{R}\), \(\mathbb{R}_+\), \((0, 1)\), \((1.01, 1.99)\)),
- a default weakly-informative prior block written in Stan as `THETA_REF_PRIOR_BLOCK` (per-family slot expansion; see §7),
- the minimum admissible \(K\) for the family (e.g., \(K_{\min} = 2\) for Gaussian, \(K_{\min} = 3\) for Student-\(t\) or Tweedie).

`param_specs` replaces what previous iterations carried as scattered per-family `switch`-statements in the codegen and as ad-hoc per-slot hyperparameter blocks in the Stan template. The resolution layer is single-source-of-truth: any addition of a new family (whether built-in or custom; see §5) extends `param_specs` and the codegen consumes the extension uniformly.

The structural alternative considered (Session 1) was to keep per-family branches in the codegen and let each branch own its slot layout. That alternative was rejected on the ground that it scales as \(O(\text{families})\) in maintenance burden and as \(O(\text{families} \times K_{\max})\) in code-path duplication; `param_specs` scales as \(O(\text{families})\) in data and as \(O(1)\) in code paths.

The `param_specs` object is exposed in the public S3 print methods (`print.gdpar_param_spec`) for diagnostic use and is returned as part of `print.gdpar_fit` summaries when `level = "spec"` is requested.

# 4. Sub-phase 8.3.2 – The (D-ID) Per-Slot Identifiability Check (4C)

Sub-phase 8.3.2 introduced the per-slot identifiability check (D-ID) as a three-layer guard:

- (D-ID layer 1, algebraic.) For each slot \(k\), the parameter \(\theta_{i,k}\) must be uniquely determined by the per-slot natural parameter \(\eta_{i,k}\) via the canonical link \(g_k\). This is a property of the family, not of the dataset: it is verified at design time by the family’s `did_layer1` field (a logical scalar; default `TRUE`).
- (D-ID layer 2, statistical.) For each slot \(k\), the empirical design \(Z_k\) (the per-slot covariate matrix) must satisfy a per-slot rank condition: \(\mathrm{rank}(Z_k) = \mathrm{ncol}(Z_k)\) over the active grouping level. The check is performed pre-fit in `R/check_identifiability.R::.check_per_slot_rank` and reported by `gdpar_check_identifiability(rigor = "full")`.
- (D-ID layer 3, numerical / posterior.) For each slot \(k\), the marginal posterior of the slot coefficients must concentrate (verified post-fit by `gdpar_bvm_check` and the contraction diagnostic of `gdpar_contraction_diagnostic`).

The release-gate decision is that (D-ID) is run per slot, not pooled: a Gaussian fit with \(K = 2\) where slot 1 (`mu`) is identifiable but slot 2 (`sigma`) is collinear under the active grouping is flagged at slot 2, not at the fit as a whole. This per-slot granularity is structurally necessary under Sub-phase 8.3.7 heterogeneity (different families per slot impose different (D-ID) checks per slot) and is preserved across all subsequent sub-phases.

The `did_override` parameter of `gdpar_family()` admits family-side plasticity. When a custom family violates a default `(D-ID layer 1)` assumption — for example, when introducing a new link with non-standard monotonicity — the user can override the default via `did_override = function(slot_idx, family_spec) {...}` and supply a slot-specific check. This was introduced in Sub-phase 8.3.5a to avoid duplicating the entire family registry for the Student-\(t\) extension.

# 5. Sub-phase 8.3.4 – Custom Families and `lognormal_loc_scale`

Sub-phase 8.3.4 introduced the (D-A3.B) custom-family descriptor pattern: the public function `gdpar_family_custom_K(stan_lpdf_id, links, supports, min_K, ...)` constructs a `gdpar_family` object from a descriptor whose Stan-side log-density is identified by `stan_lpdf_id` and whose R-side metadata mirrors the built-in registry. The descriptor pattern subsumes the previous `gdpar_family_custom()` (single slot) by reading \(K\) from the user-supplied `min_K` and dispatching the per-slot links via the same mechanism as built-in families.

Three concrete consequences:

- `lognormal_loc_scale` is the canonical example of a family that is only accessible via `gdpar_family_custom_K(stan_lpdf_id = "lognormal_loc_scale", ...)`: it is not exposed through `gdpar_family(name)` because the framework’s enumerated registry treats it as a custom-family example (the Stan-side `lognormal_loc_scale` lpdf lives in `inst/stan/amm_distrib_K.stan`, and the R-side metadata is in `R/families.R::.gdpar_family_lognormal_loc_scale_K`). The custom-family path documents how new families can be added without touching the built-in registry; `lognormal_loc_scale` itself was added during Sub-phase 8.3.4 to validate the pattern end-to-end.
- Beta (\(K = 2\), `stan_id = 5`) and Gamma (\(K = 2\), `stan_id = 6`) were added as built-in \(K = 2\) families using the same descriptor mechanism internally, with mean-parametrized lpdfs and the canonical logit / log link decomposition.
- D-B1 / D-B2 / D-B3 normalization (Sub-phase 8.3.4): the Stan-side wires the slot parameters via `apply_inv_link_by_id` (see §9 for the helper) and the R-side wires the `link_id` integer (1 = identity, 2 = log, 3 = logit, 4 = bounded-identity-Tweedie, etc.).
- D-A2 cleanup (Sub-phase 8.3.4): the previously inconsistent `family_name` slot of `gdpar_fit` was harmonized to dispatch through `param_specs[[k]]$family_name` per slot, eliminating the ambiguity that arose under future heterogeneity (Sub-phase 8.3.7).

The custom-family descriptor pattern is reused in Sub-phase 8.3.5a and 8.3.5b without further structural changes.

# 6. Sub-phase 8.3.5a – Student-\(t\) (\(K = 3\), `stan_id = 8`)

The Student-\(t\) family (\(K = 3\)) was added in Sub-phase 8.3.5a with slots `(mu, sigma, nu)` under links `(identity, log, log)`. The lpdf is the standard location-scale Student-\(t\) density (`student_t_lpdf` in Stan) with the slot-\(3\) `nu` interpreted as the degrees-of-freedom parameter.

The sub-phase emerged the `did_override` mechanism described in §4: the (D-ID layer 1) default for `nu` (positive real, log link) was correct for the family, but the registry needed a per-family table of minimum \(K\) values (`min_K = 3` for Student-\(t\)) to gate the public API. The helper `.gdpar_guard_K_below_family_min` was renamed from a session-local name to its current canonical form to make the gate explicit.

The release-gate of 8.3.5a (Session 16 of Block 8.3) was a non-trivial test of the per-slot architecture under \(K > 2\): the codegen had to expand `THETA_REF_PRIOR_BLOCK` over three slots instead of two, and the residual diagnostics (introduced later in Sub-phase 8.3.9) had to compose with a three-slot inverse-link dispatch. The test passed (suite delta `+86 PASS / 0 FAIL / +1 SKIP` for the Stan-gated smoke); the registry entry for `stan_id = 8` is canonical from that session forward.

# 7. Sub-phase 8.3.5b – Tweedie (\(K = 3\), `stan_id = 9`)

The Tweedie family (\(K = 3\)) was added in Sub-phase 8.3.5b with slots `(mu, phi, p)` under links `(log, log, identity-bounded-on-(1.01, 1.99))`. The bounded-identity link on slot 3 is necessary because the Tweedie exponent \(p\) has a known degenerate behaviour at \(p = 1\) (Poisson limit) and \(p = 2\) (Gamma limit); the bounded interval excludes both degeneracies while preserving the continuous compound-Poisson regime.

## 7.1. Hybrid lpdf: Dunn-Smyth series + saddlepoint

The Tweedie density for \(1 < p < 2\) has no closed-form expression. The framework’s Stan-side implementation uses a hybrid of two complementary approximations:

- Dunn-Smyth series. For \(1 < p < 2\), the density at \(y > 0\) admits the series representation \[f(y \mid \mu, \phi, p) \;=\; \frac{1}{y} \sum_{j=1}^{\infty} W_j(y, \phi, p), \qquad W_j(y, \phi, p) = \frac{y^{-j\alpha} (p - 1)^{j\alpha}}{\phi^{j(1-\alpha)} (2 - p)^{j} j! \Gamma(-j\alpha)},\] where \(\alpha = (2 - p) / (1 - p) \in (-\infty, 0)\). The series converges for all \((y, \mu, \phi, p)\) in the admissible region but converges slowly when \(\phi\) is small or \(y / \mu\) is large.

- Saddlepoint approximation. When the Dunn-Smyth series converges slowly, the saddlepoint approximation provides high accuracy: \[f(y \mid \mu, \phi, p) \;\approx\; \frac{1}{\sqrt{2 \pi \phi y^{p}}} \exp\!\left(-\frac{(y - \mu)^{2}}{2 \phi y \mu^{p-1}}\right).\] The saddlepoint is essentially exact in the regime where the series struggles, which is precisely the high-\(y / \mu\) tail.

The cut-over between the two approximations is implicit (Decision D5 of 8.3.5b): the function selects the series at moderate \(y / \mu\) and falls back to the saddlepoint otherwise, with the threshold derived from the Stan-side numerical tolerance rather than from a user-tunable parameter. The decision rationale is that exposing the threshold would invite mis-specification (users would tune it on a per-fit basis without principled grounds); the implicit hybrid is conservative in both regimes.

## 7.2. `THETA_REF_PRIOR_BLOCK` per-family expansion

Decision E7 of Sub-phase 8.3.5b introduced `THETA_REF_PRIOR_BLOCK` as a per-family placeholder in the Stan template, expanded at codegen time according to the slot’s link and support:

- Tweedie slot 1 (`mu`, log link, \(\mathbb{R}_+\)): weakly-informative log-normal prior centered at \(\log(\bar{y})\) (the empirical log-mean of the response), with scale 1.
- Tweedie slot 2 (`phi`, log link, \(\mathbb{R}_+\)): weakly-informative half-normal prior on the log-scale with scale 1.
- Tweedie slot 3 (`p`, bounded-identity, \((1.01, 1.99)\)): uniform prior on the bounded interval.

The placeholder mechanism is uniform across all built-in families: each family’s `param_specs[[k]]$theta_ref_prior_block` field carries the slot-specific prior expansion as a literal Stan code chunk, and the codegen splices it into the model block at compile time. The mechanism is reused unchanged in Sub-phases 8.3.6 and 8.3.7.

## 7.3. Initial-value helper

The release-gate of 8.3.5b also introduced a per-family initial-value helper for the Stan optimizer / sampler. Tweedie’s slot 3 requires an initial value strictly inside the bounded interval (a naive zero-initialization fails because \(p = 0 \notin (1.01, 1.99)\)). The helper `gdpar:::.tweedie_init_p` returns a stratified initial value (default: the midpoint \(1.5\) jittered by \(\mathrm{Uniform}(-0.1, 0.1)\) per chain). The mechanism is exposed via the `init` argument of `gdpar()` for any family that needs it; users who want to control initialization explicitly can pass a list of per-chain initial values.

# 8. Sub-phase 8.3.6 – Mixtures and Hurdle Models

Sub-phase 8.3.6 added four mixture / hurdle families:

Family | \(K\) | `stan_id` | Slot 1 | Slot 2 | Slot 3
--- | --- | --- | --- | --- | ---
ZIP | 2 | 10 | `mu` (log) | `pi` (logit) | –
ZINB | 3 | 11 | `mu` (log) | `phi` (log) | `pi` (logit)
Hurdle-Poisson | 2 | 12 | `mu` (log) | `pi` (logit) | –
Hurdle-NB | 3 | 13 | `mu` (log) | `phi` (log) | `pi` (logit)

## 8.1. Algebraic decompositions

The four families’ log-densities decompose as follows. Let \(\mu_i = \exp(\eta_{i,1})\) denote the count-process mean, \(\phi_i = \exp(\eta_{i,2})\) (for NB-based families) the dispersion, and \(\pi_i\) the slot-\(K\) inflation / hurdle probability on the unit-interval scale.

- ZIP. \(\Pr(Y_i = y \mid \mu_i, \pi_i) = \pi_i \cdot \mathbf{1}\{y = 0\} + (1 - \pi_i) \cdot \mathrm{Poisson}(y \mid \mu_i)\). The log-density at \(y_i\) is \[\log p(y_i) = \begin{cases} \log\bigl(\pi_i + (1 - \pi_i) \exp(-\mu_i)\bigr), & y_i = 0, \\ \log(1 - \pi_i) + y_i \log \mu_i - \mu_i - \log y_i!, & y_i > 0. \end{cases}\] The Stan-side implementation uses `log_mix(pi, 0, poisson_log_lpmf(0 | log_mu))` for the zero case and `log1m(pi) + poisson_log_lpmf(y | log_mu)` for the positive case.

- ZINB. Replace \(\mathrm{Poisson}(y \mid \mu_i)\) by \(\mathrm{NB}(y \mid \mu_i, \phi_i)\) in the ZIP decomposition. The Stan-side uses `neg_binomial_2_log_lpmf` instead of `poisson_log_lpmf`; the rest of the structure is identical.

- Hurdle-Poisson. The hurdle parametrization separates the zero process from the positive-count process: \[\Pr(Y_i = y \mid \mu_i, \pi_i) = \begin{cases} \pi_i, & y_i = 0, \\ (1 - \pi_i) \cdot \dfrac{\mathrm{Poisson}(y \mid \mu_i)}{1 - \exp(-\mu_i)}, & y_i > 0. \end{cases}\] The truncation factor \(1 - \exp(-\mu_i)\) in the denominator normalizes the positive-count distribution to integrate to one on the positive integers. In log-space, \(\log\bigl(1 - \exp(-\mu_i)\bigr) = \mathtt{log1m\_exp}(-\mu_i) = \mathtt{log1m\_exp}\bigl(\mathtt{poisson\_log\_lpmf}(0 \mid \log \mu_i)\bigr)\).

- Hurdle-NB. Replace the Poisson zero-probability \(\exp(-\mu_i)\) by the NB zero-probability \(\bigl(\phi_i / (\phi_i + \mu_i)\bigr)^{\phi_i}\) (or equivalently `neg_binomial_2_log_lpmf(0 | log_mu, phi)` in log-space). The hurdle truncation factor and the per-observation log-density compose identically to Hurdle-Poisson otherwise.

## 8.2. Structural decisions

- (8.3.6 D1) Stan-side log-density via `target +=`. The mixture / hurdle lpdfs are not built-in to Stan as single function calls; they require the explicit decomposition above. The decomposition is written as `target += ...` in the model block rather than as a custom `_lpdf` function, because Stan’s built-in vectorized lpdfs for `poisson_log` and `neg_binomial_2_log` are used as building blocks and the mixture structure is composed at the model-block level.
- (8.3.6 D2) Hurdle truncated-at-one via `log1m_exp`. The hurdle-component conditional (positive count given \(> 0\)) requires the truncation factor \(1 - \Pr(Y = 0)\), computed in log-space as `log1m_exp(...)` to avoid underflow when \(\mu_i\) is large. The Hurdle parametrization treats \(\pi\) as the probability of the zero outcome (not as a mixture weight on the zero component), preserving the orthogonality between the zero-process and the count-process that distinguishes hurdle from zero-inflation.
- (8.3.6 D6 = a, non-generalization) The placeholder `THETA_REF_PRIOR_BLOCK` of §7 does not generalize automatically to the mixture / hurdle case: the zero-inflation / hurdle slot `pi` (logit link, support \((0, 1)\)) receives a Beta-prior block written as a literal expansion in the per-family template rather than via the helper. This was an explicit decision (D6 = a) to avoid coupling the placeholder mechanism to the mixture-component algebra; the placeholder remains the canonical mechanism for continuous-link slots.

## 8.3. Deferred cases

The Hurdle-continuous families (Hurdle-Gamma, Hurdle-Lognormal) and the \(K < K_{\min}\) degenerate case (e.g., a user supplying `K = 1` for ZIP, in which case the inflation slot is absent) were deferred to a future sub-phase as out-of-scope for 8.3.6. The guard `.gdpar_guard_K_below_family_min` rejects \(K < K_{\min}\) configurations at the public API level.

# 9. Sub-phase 8.3.7 – Heterogeneous Families per Slot

Sub-phase 8.3.7 introduced heterogeneous families per slot: a single fit in which slot 1 uses one family and slot 2 uses another (e.g., Beta in slot 1 + Gamma in slot 2, both \(K = 2\)). The public API takes a named list in the `family` argument: `gdpar(y ~ x, family = list(slot1 = gdpar_family("beta", K = 2), slot2 = gdpar_family("gamma", K = 2)))`. The flat-family API (`family = gdpar_family("gaussian", K = 2)`) is preserved as a corollary.

The architectural alpha-decision (L1, refined in Sub-phase 8.3.7) is that heterogeneity is per slot, not per coordinate: under multivariate \(\theta_{\mathrm{ref}} \in \mathbb{R}^p\) (Sub-phase 8.6.A), the slot-\(k\) family applies uniformly across the \(p\) coordinates of slot \(k\). The per-coordinate heterogeneity (different families across the \(p\) axis within a slot) was considered and rejected on identifiability grounds: the per-coordinate decomposition \(\theta_{\mathrm{ref}, k} \in \mathbb{R}^p\) requires a coherent prior structure across coordinates, and mixing families across coordinates within a slot fragments that coherence.

The implementation introduces two Stan-side conveniences:

- `apply_inv_link_by_id(link_id, eta)` (data-block integer dispatch). The helper, defined at the top of `inst/stan/amm_distrib_K.stan`, takes a `link_id` integer and a real `eta`, returning the inverse-link-applied response-scale parameter. The dispatch is `if (link_id == 1) identity; else if (link_id == 2) exp; else if (link_id == 3) inv_logit; else if (link_id == 4) bounded_identity_tweedie; else reject`. The Stan-side `inv_link_id_per_slot` data field is a length-\(K\) integer array that codegen populates at compile time.
- `apply_W_basis_diff(W_type_id, ...)` (Sub-phase 8.3.8 sibling). The same integer-dispatch pattern is reused for the per-slot W-basis dispatch; see §10.

The release-gate guard rejects \(K \geq 3\) heterogeneous configurations (e.g., heterogeneous over three slots) as out-of-scope for 8.3.7. Five \(K = 2\) branches were refactored to the named-list API: `stan_id` 1 (Gaussian), 3 (NB), 5 (Beta), 6 (Gamma), 7 (custom-via-lognormal-loc-scale).

The bit-exactness preservation at the release-gate is structurally non-trivial: the goldens for \(K = 2\) Beta + Gamma (from Sub-phase 8.3.4) had to be preserved bit-exact under the refactor (no re-bootstrap). The audit-compare path (`gdpar_compare_eb_fb`) registered zero FAIL on the bit-exact compare for the refactored families.

# 10. Sub-phase 8.3.8 – B-spline W Basis (D-D3)

Sub-phase 8.3.8 introduced the B-spline W basis as a canonical alternative to the polynomial W basis previously used in the AMM linear predictor. Six decisions canonized:

- (D1) Alcance A: B-spline canonical. The B-spline basis is added as a first-class W-basis option (not a custom-only extension). The polynomial basis remains canonical for the scalar single-slot case; the B-spline becomes canonical for the multi-knot non-polynomial case.
- (D2) Stan-side Cox-de Boor recurrence. The Stan-side evaluation of the B-spline basis uses the Cox-de Boor recurrence (rather than a precomputed design matrix passed as data), because the basis must be evaluated at the per-slot inverse-link-transformed natural parameter at every iteration, and pre-computing the design at every iteration is wasteful when the recurrence is cheap.
- (D3) Signature `W_basis(x, type, degree, knots, boundary_knots)`. The R-side public constructor takes the explicit boundary knots and a default `degree = 3` (cubic B-spline). The default cubic was chosen because it is the lowest degree under which the second derivative is continuous (a property needed by the Sub-phase 8.3.9 residuals; see §11).
- (D4) R-side strict, Stan-side liberal. The R-side `W_basis` constructor validates the inputs strictly (knot monotonicity, boundary inclusion, degree \(\geq 0\)); the Stan-side accepts whatever the R-side passes without re-validating, on the assumption that the R-side is the single source of truth.
- (D5) Helper `apply_W_basis_diff` (Sub-phase 8.3.7 pattern). The Stan-side helper `apply_W_basis_diff(W_type_id, ...)` replicates the integer-dispatch pattern of `apply_inv_link_by_id` (see §9). The dispatch is `if (W_type_id == 1) polynomial; else if (W_type_id == 2) bspline; else reject`. The dispatch is per slot, mirroring the per-slot family heterogeneity of Sub-phase 8.3.7.
- (D6) Three smokes gated + goldens deferred. Three integration smokes were added under environment-variable gates (`GDPAR_RUN_STAN_SMOKE_BSPLINE=1`); the corresponding goldens were deferred to Sub-phase 8.3.9 (bootstrapped roster) to avoid duplicating the bootstrap effort across two sub-phases.

A side-effect of the refactor: the Stan-side `W_is_polynomial` boolean data field was eliminated and replaced by the `W_type_id` integer. The elimination is canonical and the field does not return.

## 10.1. Cox-de Boor recurrence (D2 detailed)

Given a non-decreasing knot vector \(\boldsymbol{\tau} = (\tau_0 \leq \tau_1 \leq \dots \leq \tau_{m})\) and degree \(d \geq 0\), the B-spline basis functions \(B_{j, d}(x)\) for \(j = 0, \dots, m - d - 1\) are defined by the Cox-de Boor recurrence:

\[B_{j, 0}(x) \;=\; \begin{cases} 1, & \tau_j \leq x < \tau_{j+1}, \\ 0, & \text{otherwise}, \end{cases}\]

\[B_{j, d}(x) \;=\; \frac{x - \tau_j}{\tau_{j+d} - \tau_j} B_{j, d-1}(x) \;+\; \frac{\tau_{j+d+1} - x}{\tau_{j+d+1} - \tau_{j+1}} B_{j+1, d-1}(x),\]

with the convention \(0/0 = 0\) for the boundary fractions when consecutive knots coincide. For the default cubic case (\(d = 3\)), the recurrence evaluates 4 basis function values at each \(x\) (the non-zero ones overlapping the active knot interval), independent of \(m\).

The Stan-side implementation in `inst/stan/functions/W_basis.stan` evaluates the recurrence in a single forward pass over \(d\) levels, with a local array of length \(d + 1\) holding the partial values. The numerical cost per evaluation is \(O(d^2)\), which is \(O(1)\) at the default cubic. The R-side `W_basis()` constructor returns the basis evaluation matrix (for sampling-free contexts like coefficient interpretation) using the equivalent `splines::splineDesign` call from base R; the two evaluations agree to machine precision (verified by the regression test `test-bspline_W.R`).

## 10.2. Boundary knot semantics

The `boundary_knots` argument of `W_basis(x, type, degree, knots, boundary_knots)` specifies the two outermost knots of the knot vector. The interior knots are passed via `knots`. The full knot vector is then constructed as

\[\boldsymbol{\tau} \;=\; \bigl(\underbrace{\tau_{\mathrm{lo}}, \dots, \tau_{\mathrm{lo}}}_{d + 1}, \mathtt{knots}, \underbrace{\tau_{\mathrm{hi}}, \dots, \tau_{\mathrm{hi}}}_{d + 1}\bigr),\]

with the boundary knots repeated \(d + 1\) times to enforce zero-derivative-of-order \(d\) at the boundaries (the “clamped” B-spline convention). This convention ensures that the basis evaluation at the boundary is well-defined and equal to one for the boundary basis function (a useful property for interpretation: the boundary coefficient is the response at the boundary).

The R-side `W_basis()` constructor rejects configurations where the boundary knots are not strictly outside the range of `knots` (D4 strict validation), preventing the degenerate case where the boundary repetition overlaps the interior knots.

# 11. Sub-phase 8.3.9 – Residuals G1 / G2 / G3 + S3 Methods Re-validation + Route B `predict.gdpar_fit` for \(K > 1\)

Sub-phase 8.3.9 introduced the three residual diagnostics that compose under \(K \geq 1\). Each diagnostic is defined per slot and is unit-coherent with the slot’s response scale.

## 11.1. (G1) Pearson-type residual per slot

For slot \(k\) and observation \(i\), denote by \(\mu_{i,k}(\theta_i)\) and \(s_{i,k}(\theta_i)\) the slot-implied response mean and standard deviation (both functions of the full \(\theta_i = (\theta_{i,1}, \dots, \theta_{i,K})\) in the general case, because slot statistics may depend on cross-slot quantities — e.g., for ZINB the slot-1 implied mean is the unconditional mean \(\mu_i (1 - \pi_i)\) rather than just \(\mu_i\)). The Pearson-type residual is \[r_{i,k}^{(G1)} \;=\; \frac{y_i - \mu_{i,k}(\theta_i)}{s_{i,k}(\theta_i)}.\] For \(K = 1\) Gaussian this reduces to the standard Pearson residual. For mixture / hurdle slots the residual uses the unconditional moments of the slot’s marginal contribution.

The (G1) residual is fast to compute (no auxiliary sampling) and is the canonical first-pass diagnostic. Its limitation is that it is symmetric in the location-scale family and does not detect distributional mis-specification beyond the second moment.

## 11.2. (G2) Randomized quantile residual (Dunn-Smyth) per slot

For slot \(k\) and observation \(i\), the (G2) residual is \[r_{i,k}^{(G2)} \;=\; \Phi^{-1}\bigl(F_{i,k}(y_i^*)\bigr),\] where \(F_{i,k}\) is the slot-conditional CDF and \(y_i^*\) is the realized response, possibly jittered for discrete cases following Dunn & Smyth (1996). For continuous slots, \(y_i^* = y_i\) and \(F_{i,k}(y_i)\) is uniform on \([0, 1]\) under correct specification, so \(r_{i,k}^{(G2)}\) is standard normal under the null. For discrete slots (Poisson, NB), \(F_{i,k}\) is jittered as \(u_i F_{i,k}(y_i) + (1 - u_i) F_{i,k}(y_i - 1)\) with \(u_i \sim \mathrm{Uniform}(0, 1)\), restoring the continuous-uniform property under the null.

For mixture / hurdle slots, the jitter respects the zero-component and the count-component separately: the zero-component is treated as the lump-mass at 0 (the jitter range is \([0, F_{i,k}(0)] = [0, \pi_i]\) for ZIP / Hurdle), and the positive-count component is jittered as in the count-only case.

The (G2) residual detects distributional mis-specification beyond the second moment (skewness, kurtosis, tail behaviour). Its computational cost is moderate; it requires evaluating the slot-conditional CDF at every observation, which is closed-form for the built-in families but may require numerical integration for custom families.

## 11.3. (G3) DHARMa-compatible scaled residual per slot

For slot \(k\) and observation \(i\), the (G3) residual is the rank of the realized response among \(M\) posterior predictive draws, scaled to \([0, 1]\): \[r_{i,k}^{(G3)} \;=\; \frac{1}{M} \sum_{m=1}^{M} \mathbf{1}\bigl\{y_{i,k}^{(m)} \leq y_i\bigr\} \;+\; \frac{u_i}{M}, \qquad u_i \sim \mathrm{Uniform}(0, 1),\] where \(y_{i,k}^{(m)}\) is the \(m\)-th posterior predictive draw for observation \(i\) in slot \(k\) and \(u_i\) is a continuity-correction jitter. The diagnostic is compatible with the DHARMa R package’s `simulationOutput` interface and is opt-in through the `dharma = TRUE` argument of `gdpar_residuals()`.

The (G3) residual is the most expensive of the three (it requires \(M\) posterior predictive draws per observation) but is the most general: it works for any family with a valid posterior predictive distribution, including custom families and complex mixture / hurdle compositions where the slot-conditional CDF (needed by G2) is unavailable in closed form.

## 11.4. S3 method re-validation under \(K > 1\)

The S3 methods `coef.gdpar_fit`, `print.gdpar_fit`, `summary.gdpar_fit`, `predict.gdpar_fit`, and `gdpar_residuals` were re-validated under \(K > 1\) in Sub-phase 8.3.9. Each method had been written under the assumption of \(K = 1\) and required extension to the per-slot architecture.

The release-gate decision is that the methods return named lists per slot under \(K > 1\) (preserving access to per-slot diagnostics) and flat vectors under \(K = 1\) (preserving backward compatibility with the Block 1 default). The dispatch is via a `is_K_gt_one(fit)` helper at the top of each method; the per-slot iteration is then uniform across families.

The Route B `predict.gdpar_fit` for \(K > 1\) is the most non-trivial of the S3 re-validations: the predicted response under \(K > 1\) is an array of dimension \(S \times n \times K\) (posterior draws \(\times\) observations \(\times\) slots), and the per-slot inverse-link must be applied dimension-wise. The implementation in `R/methods.R::predict.gdpar_fit` dispatches to `predict_from_newdata_K` (the \(K > 1\) helper) when `K > 1` and to `predict_from_newdata` (the \(K = 1\) helper) otherwise. The slot-wise inverse-link dispatch uses `apply_inv_link_by_id` (see §9) on the R side, mirroring the Stan-side implementation.

## 11.5. Bootstrapped goldens and the 16-column manifest

The release-gate of 8.3.9 introduced 14 bootstrapped goldens with a 16-column manifest CSV recording the goldens’ metadata: `stan_id`, \(K\), family name, seed, bootstrap iterations, posterior draws, observations, slot link IDs, golden RDS path, golden hash, generation timestamp, framework version, R version, Stan version, gdpar git commit, and per-golden notes.

The goldens are stored as compressed RDS files under `tests/testthat/_snaps/` and verified bit-exact in the regression tests `test-golden_regression_K2.R` and `test-golden_regression_8_3_9.R`. The bit-exactness check is at the per-slot coefficient level (machine precision); the test fails if any coefficient diverges by more than \(10^{-12}\) relative to the golden, which catches both implementation drift and Stan-version-induced numerical changes.

The 14 goldens cover the 7 mixture / hurdle and heterogeneous configurations plus 7 distributional / B-spline configurations selected to exercise the codegen paths that emerged in 8.3.4 to 8.3.8. The 16-column manifest is read by the regression tests via `manifest <- read.csv(system.file("extdata", "golden_manifest_8_3_9.csv", package = "gdpar"))` to dispatch each test case to its corresponding RDS file.

# 12. Sub-phase 8.3.10 – Release-Gate for Session 8.3

Sub-phase 8.3.10 is the release-gate that closes Session 8.3 as a coherent extension cycle. Six decisions canonized:

- (E4.A) `coef.gdpar_fit` under \(K > 1\) returns a named list of `gdpar_coef` per slot. Each per-slot `gdpar_coef` is a `gdpar_coef` object with the slot’s coefficient table, posterior summaries, and metadata. The list is named by the slot parameter names (`mu`, `sigma`, `phi`, etc.).
- (E5.C) Two vignettes added (`vop04_amm_intermediate.Rmd` for B-spline + heterogeneous families, `vop05_distributional_K_dharma.Rmd` for the \(K > 1\) distributional path + DHARMa diagnostics).
- (E6.C) Compare-path tiered Tier 1 / Tier 2. The audit-compare path between two fits (`gdpar_compare_fits`) is tiered: Tier 1 is the structural compare (same `stan_id`, same \(K\), same `inv_link_id_per_slot`) verifiable bit-exact; Tier 2 is the semantic compare across paths and bridges (Path A vs Path B; EB vs FB), verifiable on summary statistics. The tiering is a release-gate decision to prevent the bit-exact gate from over-fitting to numerical noise across paths.
- (V.B) Development version `0.0.0.9001`. The version bump from `0.0.0.9000` to `0.0.0.9001` marks the closure of Session 8.3 as a coherent development cycle. The version stays at `0.0.0.9001` through Sub-phases 8.4, 8.5.A, 8.5.B, 8.6 (entire), and the present audit (Block 8.4 Stage 2); a release bump (to `0.1.0` or similar) is deferred to post-Block-9 validation.
- (Fuzz S3 PASS + Tier 1 PASS + Tier 2 PASS bit-exact.) The release-gate verified `coef` / `predict` / `summary` / `print` under a randomized fuzz over 108 fits (varying \(K\), `stan_id`, family, grouping, paths); 36 Tier 1 compares; 6 Tier 2 compares. All passed.
- R CMD check `0 ERRORs / 0 WARNINGs / 0 NOTEs`, `Status: OK`. The release-gate verifies clean R CMD check from `/tmp/` at the time of closure.

## 12.1. Fuzz testing protocol

The release-gate fuzz exercises the S3 methods (`coef`, `predict`, `summary`, `print`, `gdpar_residuals`) under randomized configurations to catch silent regressions that fixed-configuration smokes might miss. The protocol is:

1. Configuration sampler. For each fuzz run, sample \((K, \texttt{stan\_id}, \texttt{family\_args}, \texttt{group\_size}, \texttt{path})\) from the admissible product space: \(K \in \{1, 2, 3\}\), \(\texttt{stan\_id} \in \{1, \dots, 13\}\) subject to \(K \geq K_{\min}(\texttt{stan\_id})\), `family_args` drawn from a per-family valid range, `group_size` \(\in \{1, 5, 25\}\), `path` \(\in \{\)Path A, Path B\(\}\). Sampling rejects degenerate configurations (e.g., \(K = 1\) for ZIP).
2. Synthetic data generator. For the sampled configuration, generate \(n = 200\) observations from the corresponding distribution under known true parameters. The true parameters are stored as ground truth for the post-fit comparison.
3. Fit and call all S3 methods. Run `gdpar()` with the sampled configuration; call `coef(fit)`, `predict(fit)`, `summary(fit)`, `print(fit)`, `gdpar_residuals(fit)`. Each call must return an object of the expected class, satisfy the K-dispatch contract (named list per slot for \(K > 1\); flat for \(K = 1\)), and pass schema-validation expectations.
4. Tier 1 bit-exact compare. For 36 of the 108 fits, re-fit under identical seed and assert bit-exact equality of the coefficient table, the predicted means, and the residuals (machine precision: \(10^{-12}\) relative tolerance).
5. Tier 2 semantic compare. For 6 of the 108 fits, fit twice under different paths (Path A and Path B) and assert that the posterior summary statistics agree within statistical sampling tolerance (Monte Carlo standard error \(\times 3\)).

The fuzz protocol is implemented in `tests/testthat/test-fuzz_s3.R`, gated by the environment variable `GDPAR_RUN_STAN_S3_FUZZ=1` because each fuzz fit takes 10-60 seconds and the full sweep is several CPU-hours. The release-gate ran the full fuzz once at closure; subsequent CI runs sample a subset.

## 12.2. Debts recorded at release-gate

The release-gate of Sub-phase 8.3.10 recorded the following debts for Sub-phase 8.4 (the audit phase):

- Vignette additions for the Sub-phase 8.3.x decisions. The present vignette canonizes that debt.
- Uniform Stan-template consolidation. Originally registered as `project_gdpar_deuda_8_4_unificacion_stan`. Canonized in Sub-sub-fase 9.3.a (Block 9, Session B9.3, 2026-05-27) under the lateral decision B.iv: the unification is realized at the R-side codegen layer rather than at the `.stan` file layer. The R-side dispatcher `.gdpar_emit_canonical_stan(spec)` is the single source of truth for the K = 1 FB Stan emissions; the two legacy templates `amm_main.stan` (p = 1) and `amm_distrib_multi.stan` (p \(\geq\) 1) are renamed to `amm_canonical_p1.stan` and `amm_canonical_pmulti.stan` and relocated to `inst/stan/_canonical_pieces/`. The public-internal wrappers `generate_stan_code()` and `generate_stan_code_multi()` are reduced to thin delegators (~10 lines each) that construct a codegen spec and call the dispatcher. Bit-exactness of the substituted Stan source is preserved by construction (the dispatcher emits byte-identical strings to the legacy paths for the same inputs); no golden re-bootstrap was required. The dedicated sub-debt of factoring the duplicated `bspline_basis_eval` and `apply_W_basis_diff` helpers into a single shared piece (the explicit deuda f originally noted inline in the multivariate template’s header) was closed in Sub-sub-fase 9.3.a colateral (Block 9, Session B9.4, 2026-05-27) under the lateral decision G.iv: the helpers live in a dedicated `inst/stan/_canonical_pieces/amm_canonical_helpers.stan` Stan source FRAGMENT (no surrounding `functions { }` block; not standalone-parseable) and the dispatcher inserts them via R-side textual substitution at the `// {{CANONICAL_HELPERS}}` placeholder inside the body pieces’ `functions { }` block. No cmdstan `#include` semantics are used (the dispatcher is the sole assembly site). Stan semantics are preserved bit-exact (helper-block comments unified across the body pieces, taking the detailed algorithmic version from `canonical_p1`); goldens remain valid by-construction; `cmdstanr` cache hash recompiles once on the first call post-B9.4. The same helpers piece is consumed by the EB-side cascade pending in B9.5+.
- Path C bit-exact comparator template (originally projected as `amm_distrib_KxP.stan`): codegen + dispatcher canonized in Sub-bloque 9.3.d (Block 9, Session B9.5, 2026-05-27) under the lateral decision I.iv + the piece arquitectonica J.iv.A (DESIGN_9_3_D_PATH_C.md sub-decision 3.2). The realization avoids creating a new top-level `.stan` file: a dedicated canonical piece `inst/stan/_canonical_pieces/amm_canonical_pmulti_KxP.stan` is dispatched via `.gdpar_emit_canonical_stan(spec)` with `spec$p_class = "pmulti_KxP"`, reusing the canonical helpers piece infrastructure of G.iv. The new piece mirrors `amm_eb_marginal_KxP.stan` (the EB-side architecture canonized in Sub-fase 8.6.D) extended with the globally shared W modulating component and the full Path B family set \(\{1, 3, 5, 6, 7, 8, 9, 10, 11, 12, 13\}\). The fit harness was the deuda D69 deferred from B9.5; operationally closed in Session B9.6 (2026-05-27) under the canonized decisions M.iv unified + N.a: the unified assembler `.assemble_stan_data_KxP(spec)` now dispatches on `path = c("EB", "FB")` (the EB branch preserves byte-identical the 8.6.D first-iteration restrictions; the FB branch lifts W disabled + restricted stan_id to the full Path B set with the W block metadata fields appended) and the new internal driver `.gdpar_fb_KxP_fit()` composes design + assembler + cmdstanr sampling, reusing the canonical EB-side init helper `.gdpar_eb_make_random_init_KxP()` (Sub-fase 8.6.D Session 13c). The 4 reserved canonical seeds \(91001 \ldots 91004\) for the K.c roster (gaussian / beta / gamma / student_t) are now produced as Tier 1 fitable goldens in `tests/testthat/data/golden_fb_KxP_<family>_polynomial_K<K>_p<p>.rds` with the canonical 16-column manifest row.
- KSD-joint for EB-vs-FB comparison. The current `gdpar_compare_eb_fb` reports total-variation (TV) on marginal posteriors; the joint KSD (kernel Stein discrepancy) deferred to Block 9.x.

# 13. Synthesis: The Per-Slot Architecture as a Coherent Design

Read sequentially, Sub-phases 8.3.1 to 8.3.10 implement a single architectural decision applied consistently: the per-slot resolution layer is the single source of truth for everything that depends on the slot index \(k\). Five corollaries follow:

- (Corollary 1, codegen.) Stan-side per-slot expansions (`THETA_REF_PRIOR_BLOCK`, `apply_inv_link_by_id`, `apply_W_basis_diff`) read their per-slot configuration from `param_specs` rather than from per-family branches. The codegen is therefore family-agnostic at the slot-iteration level.
- (Corollary 2, identifiability.) The (D-ID) check (§4) is naturally per slot because the architecture is per slot. Heterogeneous families per slot (§9) inherit (D-ID) without further work.
- (Corollary 3, custom families.) A custom family added via `gdpar_family_custom_K` (§5) participates in the per-slot architecture immediately, without modifications to the codegen, the (D-ID) check, the residual diagnostics, or the S3 methods.
- (Corollary 4, residuals.) The G1 / G2 / G3 diagnostics (§11) compose per slot under any \(K \geq 1\) family registered in `param_specs`. Adding a new family adds a new slot’s residual contribution without touching the diagnostic code.
- (Corollary 5, S3 methods.) The S3 methods (§11, §12) dispatch on `K > 1` once at the top level and then iterate over slots. The dispatch logic is uniform across families.

The architecture’s structural cost is the up-front investment in `param_specs` and the integer-dispatch helpers (§5, §9, §10). The structural benefit is that the cost is paid once and amortized across all families and all future extensions.

# 14. Implementation Status (Release-Gate Context)

The present vignette canonizes decisions implemented in development version `0.0.0.9001` and verified at the release-gate of Sub-phase 8.3.10. Concretely:

- The 13 built-in families with `stan_id` in \(\{1, \dots, 13\}\) are implemented and tested. The custom-family path via `gdpar_family_custom_K` is operational; `lognormal_loc_scale` is the canonical example.
- The B-spline W basis is implemented and gated by the Sub-phase 8.3.8 smokes; the bootstrapped goldens of Sub-phase 8.3.9 cover the regression cases.
- The G1 / G2 / G3 residuals are implemented under any registered family; the DHARMa-compatible (G3) is opt-in through the `dharma = TRUE` argument.
- The S3 methods (`coef`, `predict`, `summary`, `print`, `gdpar_residuals`) are implemented under any \(K \geq 1\) and verified bit-exact for the registered families.

The following items are deferred at the time of this vignette (Block 8.4 Stage 2 audit closure):

- Hurdle-continuous families (Hurdle-Gamma, Hurdle-Lognormal): deferred to a future sub-phase as out-of-scope for Sub-phase 8.3.6.
- \(K \geq 3\) heterogeneous configurations: rejected at Sub-phase 8.3.7 release-gate as out-of-scope; opt-in via the public named-list API rejects these configurations at the guard.
- Multi-parametric extension (\(p \geq 1\) Block 8 + \(K \geq 1\) Block 8.3 coupled to anchored hierarchical pooling): implemented at Sub-phase 8.6 (Blocks 8.6.A / B / C / D); the present vignette is restricted to the \(K\)-axis.
- External validation against TOP-3 competitors (mgcv, brms, INLA on the eBird benchmark): scheduled for Block 9 (re-validation; the original Block 7 internal benchmark validated \(K = 1\) Gaussian / Poisson / NB).
- Path C (\(K > 1 \wedge p > 1\)) bit-exact comparator template: codegen + dispatcher + piece architecture closed in Sub-bloque 9.3.d (Block 9, Session B9.5) per the I.iv + J.iv.A canonization; the FB fit harness (assembler + driver + 4 fitable goldens with bootstrapped manifest rows) closed in Session B9.6 under the M.iv unified + N.a canonizations. Cluster K4 of Sub-bloque 9.3 is now fully closed (FB-side + EB-side cascade L heredada + fit harness D69).
- Internal synthetic adversarial re-validation (Sub-bloque 9.1): opened in Session B9.7 under Decision C.iv (layered defence with drift-gated reseed). Layers 1-3 (structural adversarial fuzz, deterministic property-based invariances, NA / outlier robustness) plus a bit-exact drift smoke executed substantively; the Stan-bound layers (statistical fit fuzz, 12-family outcome-outlier stress, full bit-exact compare-path) are gated to the nocturnal B9.8 sub-unit. Report: `inst/benchmarks/results/block9_internal.md`.

The release version bump from `0.0.0.9001` to a stable `0.1.0` or similar is scheduled for post-Block-9 closure.

# 15. References (selective)

- Dunn, P. K. & Smyth, G. K. (1996). Randomized quantile residuals. Journal of Computational and Graphical Statistics 5(3), 236–244.
- Dunn, P. K. & Smyth, G. K. (2005). Series evaluation of Tweedie exponential dispersion model densities. Statistics and Computing 15(4), 267–280.
- Hartig, F. (2024). DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models. R package, CRAN.
- Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., Brubaker, M., Guo, J., Li, P., & Riddell, A. (2017). Stan: A probabilistic programming language. Journal of Statistical Software 76(1), 1–32.
- Cox, M. G. (1972). The numerical evaluation of B-splines. IMA Journal of Applied Mathematics 10(2), 134–149.
- de Boor, C. (1972). On calculating with B-splines. Journal of Approximation Theory 6(1), 50–62.
- Mullahy, J. (1986). Specification and testing of some modified count data models. Journal of Econometrics 33(3), 341–365 (hurdle models).
- Lambert, D. (1992). Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics 34(1), 1–14.

End of v08d.