2.1. Three equivalent input forms

gdpar() accepts three syntactically equivalent ways of declaring a K > 1 distributional regression. All three canonicalise to the same internal gdpar_formula_set object (sub-phase 8.3.3, decision E):

library(gdpar)

# (a) brms-style `bf()` sugar
fit <- gdpar(
  gdpar_bf(y ~ a(x1), sigma ~ a(x2)),
  data = d, family = gdpar_family("gaussian")
)

# (b) Named list of formulas
fit <- gdpar(
  list(mu = y ~ a(x1), sigma = ~ a(x2)),
  data = d, family = gdpar_family("gaussian")
)

# (c) Named list of amm_spec (low-level, bypasses formula parsing)
fit <- gdpar(
  list(
    mu    = amm_spec(a = ~ x1),
    sigma = amm_spec(a = ~ x2)
  ),
  data = d, family = gdpar_family("gaussian")
)

Three contract notes:

• The slot names of the gdpar_formula_set must match the eligible param_specs of the family. For gaussian, the eligibles are {mu, log_sigma}; the alias sigma is canonicalised to log_sigma at construction time. For other families see gdpar_family(name)$param_specs.
• The left-hand-side response y appears only in the first formula (the observation slot, slot 1). Subsequent slots use one-sided formulas (~ a(x2)).
• Suppressing the anchor via - 1 or + 0 is an error: gdpar canonicalises theta_ref as a structural anchor, not as an optional parameter (sub-phase 8.3.3, decision 5).

2.2. Choosing K

The number of slots K is determined by the input. K = 1 retains the legacy path; K = 2 adds dispersion / scale modelling; K = 3 adds shape / weight modelling. The minimum K per family is enforced by .gdpar_guard_K_below_family_min:

A K = 1 fit on a beta family aborts with gdpar_input_error pointing to elevation to K = 2.

The pattern name lognormal_loc_scale is not part of the enum of gdpar_family(name): the package registers it as a K = 2 custom-family pattern (canonised in Sub-phase 8.3.4), accessed via gdpar_family_custom_K(stan_lpdf_id = "lognormal_loc_scale", ...). See §2.4 below for the literal recipe.

2.3. End-to-end example: Gaussian K = 2

set.seed(2026L)
n <- 100L
x1 <- rnorm(n); x2 <- rnorm(n)
mu_true       <- 0.4 + 0.6 * (x1 - mean(x1))
log_sigma_eta <- -0.2 + 0.4 * (x2 - mean(x2))
y <- rnorm(n, mu_true, exp(log_sigma_eta))
d <- data.frame(y = y, x1 = x1, x2 = x2)

library(gdpar)
fit_K2 <- gdpar(
  gdpar_bf(y ~ a(x1), sigma ~ a(x2)),
  data   = d,
  family = gdpar_family("gaussian"),
  chains = 2L, iter_warmup = 400L, iter_sampling = 400L,
  refresh = 0L
)

co <- coef(fit_K2)
co$mu
co$sigma

coef.gdpar_fit for K > 1 returns a named list of gdpar_coef objects (decision E4.A, sub-phase 8.3.10). Each entry follows the scalar gdpar_coef contract: posterior summaries of theta_ref, the additive a, the multiplicative b/c_b, and the modulating W. The modulating block is globally shared across slots (replicated identically in every slot’s W component).

2.4. Custom K > 1 families via gdpar_family_custom_K()

The constructor gdpar_family_custom_K() exposes the K-family custom-pattern registry opened in Sub-phase 8.3.4 of Block 8. Each registered pattern is identified by a stan_lpdf_id (the name of a pre-validated Stan _lpdf function shipped with the package) and carries its own minimum K. The first pattern registered is lognormal_loc_scale (min_K = 2); subsequent sub-phases extend the whitelist.

Signature:

gdpar_family_custom_K(
  name,                 # character scalar; must not collide with a built-in
  stan_lpdf_id,         # character scalar; key in the registry
  did_holds     = TRUE, # logical; user declaration of D-ID
  did_condition = NULL, # character scalar describing any conditional D-ID
  did_reference = NULL  # citation supporting did_holds
)

Literal recipe for lognormal_loc_scale (a K = 2 location-scale family on the log scale; slot 1 carries the location and slot 2 carries the log-scale):

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

fit_lognorm <- gdpar(
  gdpar_bf(y ~ a(x1), sigma ~ a(x2)),
  data   = d,
  family = my_lognorm,
  chains = 2L, iter_warmup = 400L, iter_sampling = 400L,
  refresh = 0L
)

The user is responsible for asserting that the chosen pattern is identifiable in its parameter (did_holds = TRUE); the package does not test identifiability from data, only registers the declaration. Attempting to use an unregistered stan_lpdf_id aborts with a gdpar_input_error that enumerates the allowed patterns. The general (K = 1) custom-family constructor gdpar_family_custom() is documented in §6 below.

2.5. Prediction

# In-sample prediction (theta_i_k draws)
pred_in <- predict(fit_K2, summary = "mean_se")
str(pred_in, max.level = 1L)

# Out-of-sample prediction on new covariates
new_d <- data.frame(x1 = c(-1, 0, 1), x2 = c(-1, 0, 1))
pred_new <- predict(fit_K2, newdata = new_d, summary = "mean_se")
str(pred_new, max.level = 1L)

Three points of contract:

• summary = "draws" returns the posterior array of shape \((S, n, K)\), with the third dimension named after the slots.
• summary = "mean_se" and summary = "quantiles" return a named list of length K with one summary data frame per slot.
• type = "response" applies the slot’s canonical inverse link from object$family$param_specs[[k]]$inv_link. Slot 1 uses the location’s primary link; slot \(k > 1\) uses its own canonical link (log for sigma / phi, logit for pi, etc.). Applying the location link to all slots indiscriminately would be incorrect, so the implementation iterates the per-slot link.
• newdata on K > 1 polynomial W fits is supported (sub-phase 8.3.10, decision Ruta B + newdata extension). For B-spline W on K > 1 and for grouping (J_groups > 1), the function raises gdpar_unsupported_feature_error and points to Session 8.4 deudas.

3.1. API

# G1: deviance and Pearson (frequentist canonical)
r_dev  <- residuals(fit_K2, type = "deviance")
r_pear <- residuals(fit_K2, type = "pearson")

# G2: Bayesian quantile residuals (Dunn-Smyth)
r_q <- residuals(fit_K2, type = "quantile", randomize_seed = 1L)

# Response residuals (y_obs - mean of y_pred draws)
r_resp <- residuals(fit_K2, type = "response")

head(data.frame(deviance = r_dev, pearson = r_pear,
                quantile = r_q, response = r_resp))

The signature residuals.gdpar_fit(object, type, coord = NULL, randomize_seed = NULL, ...) lets the user pin the randomisation seed for reproducible G2 residuals across runs. For multi-coordinate fits (p > 1), coord selects which coordinate is summarised.

3.2. Posterior predictive draws and PPC

# Posterior-predictive draws (S x n matrix for K=1 or K>1 with p=1)
pp <- gdpar_posterior_predict(fit_K2)
dim(pp)

# Visual PPCs via bayesplot::pp_check generic
if (requireNamespace("bayesplot", quietly = TRUE)) {
  pp_check(fit_K2, type = "dens_overlay", ndraws = 30L)
}

gdpar_posterior_predict is the exported posterior-predictive draws extractor; pp_check.gdpar_fit is an S3 method off the bayesplot::pp_check generic and supports five PPC types: dens_overlay, hist, ecdf_overlay, stat, intervals. Loading bayesplot makes pp_check(fit_K2) work directly; without it the user can still call pp_check.gdpar_fit(fit_K2) if bayesplot is installed.

4.1. API

dh <- gdpar_dharma_object(fit_K2)
class(dh)
#> [1] "DHARMa"
DHARMa::testResiduals(dh)

#> $uniformity
#> 
#>  Asymptotic one-sample Kolmogorov-Smirnov test
#> 
#> data:  simulationOutput$scaledResiduals
#> D = 0.085, p-value = 0.7789
#> alternative hypothesis: two-sided
#> 
#> 
#> $dispersion
#> 
#>  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
#>  simulated
#> 
#> data:  simulationOutput
#> dispersion = 1.0076, p-value = 0.86
#> alternative hypothesis: two.sided
#> 
#> 
#> $outliers
#> 
#>  DHARMa outlier test based on exact binomial test with approximate
#>  expectations
#> 
#> data:  simulationOutput
#> outliers at both margin(s) = 0, observations = 60, p-value = 1
#> alternative hypothesis: true probability of success is not equal to 0.004987531
#> 95 percent confidence interval:
#>  0.00000000 0.05962949
#> sample estimates:
#> frequency of outliers (expected: 0.00498753117206983 ) 
#>                                                      0

The returned object is a standard DHARMa::createDHARMa() simulationOutput with:

• nSim: number of posterior-predictive draws used for the empirical CDF (default: all available).
• simulatedResponse: matrix of posterior-predictive draws, shape \((n, \text{nSim})\).
• observedResponse: the original y_obs.
• fittedPredictedResponse: posterior mean of y_pred per observation.

All DHARMa post-processing functions (testUniformity, testDispersion, testOutliers, testZeroInflation, plotResiduals, plotQQunif) work off this object.

4.2. When to use DHARMa vs the built-in G2

The two paths agree on methodology (Bayesian randomized quantile residuals à la Dunn-Smyth 1996). They differ in scope:

• Built-in residuals(fit, type = "quantile") returns the residuals as a numeric vector mapped through qnorm() so they are approximately standard normal under the correct model. Useful as input to qqnorm(), shapiro.test(), or any downstream tool that expects standard-normal residuals.
• gdpar_dharma_object(fit) returns the residuals as a DHARMa simulationOutput in the \([0, 1]\) scale. Useful when the user wants DHARMa’s plotting machinery and the formal tests for over- / under-dispersion, outliers, and zero-inflation.

Both paths are reproducible: pass randomize_seed to residuals() or set set.seed() before gdpar_dharma_object().