param_specs as the Per-Slot Resolution
Layerlognormal_loc_scalestan_id = 8)stan_id = 9)
predict.gdpar_fit for \(K
> 1\)
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):
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:
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).
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.
The framework’s slot convention is positional, mirroring the parameter order of the underlying density:
mu (identity link), slot 2 =
sigma (log link).mu (log link); for NB, slot 2 = phi (log
link).mu (logit link, support \((0,1)\)), slot 2 = phi (log
link, precision).mu (log link), slot 2 = phi (log
link, shape under mean-shape parametrization).mu
(identity), slot 2 = sigma (log), slot 3 = nu
(log).mu (log), slot 2 =
phi (log), slot 3 = p (identity, bounded \((1.01, 1.99)\)).mu (log), slot 2 =
pi (logit).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.
param_specs as the
Per-Slot Resolution LayerThe 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\):
"mu",
"sigma", "phi", "nu",
"p", "pi"),link_id used by the Stan helper
apply_inv_link_by_id (Sub-phase 8.3.7; see §9),THETA_REF_PRIOR_BLOCK (per-family slot expansion; see
§7),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.
Sub-phase 8.3.2 introduced the per-slot identifiability check (D-ID) as a three-layer guard:
did_layer1 field (a logical scalar; default
TRUE).R/check_identifiability.R::.check_per_slot_rank and
reported by
gdpar_check_identifiability(rigor = "full").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.
lognormal_loc_scaleSub-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.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.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.).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.
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.
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.
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.
THETA_REF_PRIOR_BLOCK per-family
expansionDecision 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:
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.phi, log link, \(\mathbb{R}_+\)):
weakly-informative half-normal prior on the log-scale with scale 1.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.
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.
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) |
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.
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.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.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.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.
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.
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:
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).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.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.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.
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).
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.
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.
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.
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.
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.
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.
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.
Sub-phase 8.3.10 is the release-gate that closes Session 8.3 as a coherent extension cycle. Six decisions canonized:
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.).vop04_amm_intermediate.Rmd for B-spline + heterogeneous
families, vop05_distributional_K_dharma.Rmd for the \(K > 1\) distributional path + DHARMa
diagnostics).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.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.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.0 ERRORs / 0 WARNINGs / 0 NOTEs,
Status: OK. The release-gate verifies clean R CMD
check from /tmp/ at the time of closure.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:
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).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.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.
The release-gate of Sub-phase 8.3.10 recorded the following debts for Sub-phase 8.4 (the audit phase):
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+.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.gdpar_compare_eb_fb reports total-variation (TV) on
marginal posteriors; the joint KSD (kernel Stein discrepancy) deferred
to Block 9.x.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:
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.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.param_specs. Adding a new family adds a new slot’s residual
contribution without touching the diagnostic code.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.
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:
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.dharma = TRUE argument.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):
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.
End of v08d.