## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, collapse = TRUE, comment = "#>" ) have_cmdstan <- requireNamespace("cmdstanr", quietly = TRUE) && isTRUE(try(cmdstanr::cmdstan_version(), silent = TRUE) != "") have_dharma <- requireNamespace("DHARMa", quietly = TRUE) have_bayesplot <- requireNamespace("bayesplot", quietly = TRUE) # When both cmdstanr and DHARMa are available the §4.1 [dharma-api] chunk is # evaluated; it consumes a `fit_K2` that the user-facing recipe in §2.3 declares # under `eval = FALSE`. Build a minimal Gaussian K = 2 fit here so the rendered # DHARMa section is reproducible end-to-end. The setup chunk is `include = FALSE` # so the construction is silent in the rendered vignette. fit_K2 <- NULL if (have_cmdstan && have_dharma && requireNamespace("gdpar", quietly = TRUE)) { library(gdpar) fit_K2 <- tryCatch({ set.seed(2026L) .n_setup <- 60L .x1_setup <- rnorm(.n_setup) .x2_setup <- rnorm(.n_setup) .mu_setup <- 0.4 + 0.6 * (.x1_setup - mean(.x1_setup)) .ls_setup <- -0.2 + 0.4 * (.x2_setup - mean(.x2_setup)) .y_setup <- rnorm(.n_setup, .mu_setup, exp(.ls_setup)) .d_setup <- data.frame(y = .y_setup, x1 = .x1_setup, x2 = .x2_setup) gdpar( gdpar_bf(y ~ a(x1), sigma ~ a(x2)), data = .d_setup, family = gdpar_family("gaussian"), chains = 2L, iter_warmup = 200L, iter_sampling = 200L, refresh = 0L, show_messages = FALSE ) }, error = function(e) NULL) } have_fit_K2 <- !is.null(fit_K2) ## ----api-three-forms, eval=FALSE---------------------------------------------- # 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") # ) ## ----k2-gaussian, eval=FALSE-------------------------------------------------- # 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 ## ----family-custom-K-signature, eval=FALSE------------------------------------ # 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 # ) ## ----family-custom-K-lognormal, eval=FALSE------------------------------------ # 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 # ) ## ----k2-predict, eval=FALSE--------------------------------------------------- # # 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) ## ----residuals-api, eval=FALSE------------------------------------------------ # # 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)) ## ----posterior-predict-api, eval=FALSE---------------------------------------- # # 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) # } ## ----dharma-api, eval=have_fit_K2--------------------------------------------- dh <- gdpar_dharma_object(fit_K2) class(dh) DHARMa::testResiduals(dh) ## ----zinb-example, eval=FALSE------------------------------------------------- # set.seed(515L) # n <- 120L # x1 <- rnorm(n); x2 <- rnorm(n); x3 <- rnorm(n) # mu_eta <- 1.0 + 0.5 * (x1 - mean(x1)) # log_phi <- -0.3 + 0.2 * (x2 - mean(x2)) # logit_pi <- -1.0 + 0.6 * (x3 - mean(x3)) # mu_true <- exp(mu_eta) # phi_true <- exp(log_phi) # pi_true <- 1 / (1 + exp(-logit_pi)) # zero_struc <- rbinom(n, 1, pi_true) # y_count <- rnbinom(n, size = phi_true, mu = mu_true) # y <- ifelse(zero_struc == 1L, 0L, y_count) # d <- data.frame(y = y, x1 = x1, x2 = x2, x3 = x3) # # fit_zinb <- gdpar( # gdpar_bf(y ~ a(x1), phi ~ a(x2), pi ~ a(x3)), # data = d, # family = gdpar_family("zinb"), # chains = 2L, iter_warmup = 600L, iter_sampling = 600L, # refresh = 0L # ) # # # Per-slot coefficient summary # co <- coef(fit_zinb) # names(co) # co$mu # co$pi ## ----zinb-residuals, eval=FALSE----------------------------------------------- # # G2 quantile residuals — robust to mixture structure when jittering # # discrete responses is enabled (default for ZIP/ZINB/hurdle). # r_q <- residuals(fit_zinb, type = "quantile", randomize_seed = 99L) # hist(r_q, breaks = 20L, # main = "Bayesian quantile residuals — ZINB K=3", # xlab = "residual") # # # DHARMa-side diagnostics if available # if (requireNamespace("DHARMa", quietly = TRUE)) { # dh <- gdpar_dharma_object(fit_zinb) # DHARMa::testZeroInflation(dh) # } ## ----family-custom-signature, eval=FALSE-------------------------------------- # gdpar_family_custom( # name, # character scalar; must not collide with a built-in # link, # one of "identity", "log", "logit" # did_holds, # logical; explicit user declaration of D-ID # did_condition, # character scalar (NA_character_ if unconditional) # stan_loglik_block, # Stan snippet for the model block (per-observation # # target += ... ; references eta[i] and y_real[i] or # # y_int[i] per y_type) # stan_log_lik_block, # Stan snippet for generated quantities log_lik[i] # stan_y_pred_block, # Stan snippet for generated quantities y_pred[i] # y_type, # one of "real", "integer" # did_reference # citation supporting did_holds # ) ## ----family-custom-lognormal, eval=FALSE-------------------------------------- # my_family <- gdpar_family_custom( # name = "my_log_normal", # link = "log", # did_holds = TRUE, # did_condition = NA_character_, # stan_loglik_block = # "target += normal_lpdf(log(y_real[i]) | eta[i], sigma_y[1]);", # stan_log_lik_block = # "log_lik[i] = normal_lpdf(log(y_real[i]) | eta[i], sigma_y[1]);", # stan_y_pred_block = # "y_pred[i] = exp(normal_rng(eta[i], sigma_y[1]));", # y_type = "real", # did_reference = "User declaration" # )