1. Who this vignette is for

This is the first vignette to read if you have never used gdpar before. It assumes only that you have a working R installation (R >= 4.2.0); it does not assume familiarity with Stan, Bayesian inference, or the AMM canonical form. The goal is a single end-to-end fit on a synthetic dataset, with one paragraph per command explaining what is happening and what to look at in the output.

After this quickstart you can move on, in order of operational depth, to:

• vop01_parametrization_toggle — how the parametrization auto-decision works and how to override it.
• vop02_arbitrary_p — fits with covariate dimension \(p > 1\).
• vop03_grouped_anchors — fits with grouped references via the group argument.
• vop04_amm_intermediate — B-spline W bases, custom families, heterogeneous slots.
• vop05_distributional_K_dharma — distributional fits with \(K > 1\) slot families and DHARMa residual diagnostics.

For the theoretical background, the entry point is v00_framework_overview; for identifiability conditions, v01_amm_identifiability; for the operational meaning of the population reference, v02_gnoseological_validity.

2. Installation

The package depends on cmdstanr for the Stan back-end. Install once per machine:

install.packages(
  "cmdstanr",
  repos = c("https://stan-dev.r-universe.dev", getOption("repos"))
)
cmdstanr::install_cmdstan()

install_cmdstan() downloads and compiles the C++ Stan toolchain in the user’s home directory. The first run takes a few minutes; subsequent runs reuse the cached install.

Then install gdpar itself (development version):

# install.packages("remotes")
remotes::install_github("IsadoreNabi/gdpar")

Load the package:

library(gdpar)

3. A synthetic dataset

The quickstart uses a deliberately simple dataset: one outcome \(y\), one covariate \(x\), no grouping, no multivariate covariate path, and a Gaussian likelihood. The data-generating process is

\[y_i = \theta_{\text{ref}} + a\,x_i + \varepsilon_i, \qquad \varepsilon_i \sim \mathcal{N}(0, \sigma^2).\]

Here \(\theta_{\text{ref}} = 1.5\) is the population reference, \(a = 0.8\) is the additive slope on \(x\), and \(\sigma = 0.3\) is the residual noise scale. We simulate \(n = 200\) observations:

set.seed(20260526)
n <- 200L
x <- rnorm(n)
theta_ref_true <- 1.5
a_true <- 0.8
sigma_true <- 0.3
y <- theta_ref_true + a_true * x + rnorm(n, sd = sigma_true)
dat <- data.frame(y = y, x = x)
head(dat)

The data frame dat has two columns: the outcome y and the covariate x. No further preparation is needed; gdpar() will centre and scale x internally for the AMM design matrices.

4. A first fit

The call is one line:

fit <- gdpar(
  formula = y ~ x,
  data    = dat,
  family  = gdpar_family("gaussian"),
  path    = "bayes"
)

A few things happen here that are worth knowing before reading the output.

• formula = y ~ x follows the standard R modelling formula. The left-hand side is the outcome; the right-hand side lists the covariate(s) that enter the AMM deviation \(a(x) + b(x)\,\theta_{\text{ref}} + W(\theta_{\text{ref}})\,x\).
• family = gdpar_family("gaussian") selects the Gaussian likelihood with identity link. Other primary families (Poisson, NB, Bernoulli) accept the same constructor.
• path = "bayes" selects Path 1 (hierarchical Bayesian via Stan), which is the only currently implemented path in gdpar 0.0.0.9001. Calls with path = "vcm" or path = "hyper" will abort with gdpar_unsupported_feature_error.
• The function runs a pre-flight diagnostic (a short NCP fit) to choose CP vs NCP parametrization per component; see vop01 for the details. Then it samples the long fit. With \(n = 200\), default chains and warmup, the wall-time is on the order of a minute on a current laptop.

The return object fit has class gdpar_fit and carries the Stan draws, the AMM design matrices, the parametrization decision, and the post-fit identifiability and validity diagnostics.

5. Reading the print

The print method shows a compact one-screen summary:

print(fit)

The default print covers, in order:

• Call and path / family / link declared by the user.
• Sample size \(n\), covariate dimension \(p\), and slot count \(K\) (here \(p = 1\), \(K = 1\)).
• Parametrization decision per component (a and W) reported by the pre-flight, plus the filter that fired (see vop01 §3).
• Sampling diagnostics: number of divergences, maximum tree-depth saturations, E-BFMI per chain. The line passes when there are zero divergences and E-BFMI is above the default threshold.
• Identifiability verdict: the joint diagnostic of (C1)-(C7) declared by gdpar_check_identifiability() on the fit (see v01 for the conditions).

A clean output ends with a line of the form Identifiability: OK | Sampling: OK. If either fails, the print directs the user to the corresponding diagnostic helper.

6. Coefficients

The coef method returns the posterior summaries:

coefs <- coef(fit)
coefs

For a \(K = 1\) fit, coef(fit) returns a single object of class gdpar_coef with named slots for theta_ref, a_coef, b_coef, W_coef, sigma_y (when applicable), plus the hyperparameter slots (mu_theta_ref, sigma_theta_ref, sigma_a, sigma_b, sigma_W). Each slot is a data frame with one row per parameter and columns mean, sd, q5, q50, q95. The recovered theta_ref should sit close to the true value \(1.5\) within the posterior credible band; the recovered a_coef close to \(0.8\); the recovered sigma_y close to \(0.3\).

For \(K > 1\), coef(fit) returns a named list of gdpar_coef objects, one per slot; see vop05 for the multi-slot case.

7. Predictions

The predict method generates posterior predictions, optionally on new data:

new_x <- data.frame(x = seq(-2, 2, length.out = 11))
preds <- predict(fit, newdata = new_x, level = 0.9)
preds

The output is a data frame with one row per newdata row and columns mean, lower, upper (and per-slot columns when \(K > 1\)). The level argument controls the credible-band width. If newdata is omitted, predict returns posterior predictions on the training data.

8. Diagnostics

The diagnostics() helper aggregates the post-fit diagnostic battery in one call:

diag <- diagnostics(fit)
print(diag)

The output reports, in order, the sampling diagnostics (already in print(fit)), the identifiability diagnostics, the gnoseological validity diagnostics ((HOM) / (REG) / (CLOS) tests; see v02), and the contraction diagnostic (\(\theta_{\text{ref}}\) posterior contraction relative to its prior; see v04 and R/contraction_diagnostic.R). Each block ends with a one-line verdict (OK / WARN / FAIL) and a pointer to the relevant vignette when the verdict is not OK.

For Path 1 Gaussian fits like this quickstart, the expected output is OK on all four blocks. For low-information fits, weak priors, or violations of structural conditions, the corresponding block flags WARN or FAIL and the user follows the vignette pointer for the remediation recipe.

9. What to read next

If the quickstart ran end-to-end and the output looks plausible, the next operational steps depend on what you want to do:

• Multiple covariates (\(p > 1\)): vop02_arbitrary_p. The formula interface extends naturally; the per-component parametrization decision is made by the multivariate pre-flight (filters 1-3, see vop01 §3 for the scalar / multivariate cross-reference).
• Grouped reference (one \(\theta_{\text{ref}}\) per group): vop03_grouped_anchors. The group argument of gdpar() enables this regime; see v01 §6.6.2 and the (C7) condition for the identifiability story.
• B-spline modulation \(W\): vop04_amm_intermediate §2.
• Heterogeneous families per slot (\(K > 1\) with mixed families): vop04_amm_intermediate §3 + vop05_distributional_K_dharma.
• Custom family: vop05_distributional_K_dharma §2.4 (lognormal_loc_scale recipe) and §6 (general gdpar_family_custom registry).
• CATE/ITE estimation: v08_cate_ite_positioning (theory) + v08b_cate_ite_bridge_implementation (T-learner bridge) + vop06_meta_learner_comparison (operational + external meta-learner comparison via gdpar_compare_meta_learners).
• Empirical Bayes for \(\theta_{\text{ref}}\): v07_eb_vs_fb (theory) + v07b_eb_multivariate (multivariate extension) + vop07_eb_workflow (operational, gdpar_eb + gdpar_compare_eb_fb).
• Asymptotic theory: v04_asymptotics_path1_bayesian for Path 1.

If the quickstart failed at any step —installation, fit, or diagnostics— the relevant vignette pointer above is the canonical reference for remediation; this quickstart deliberately does not cover error recovery, custom families, or advanced diagnostics.

10. Summary

A first gdpar fit needs three ingredients: data (a data frame), a formula, and a family. Everything else has defaults that are operational out of the box for the AMM Gaussian case. The print, coef, predict, and diagnostics methods cover the day-to-day reading of the fit; the eight operational vignettes (vop00-vop07) plus the theoretical vignettes (v00-v09 series) provide the depth for any specific use case.