Conditional distribution as the inferential target

Let \(Y\) denote the response and let \(X \in \mathbb{R}^p\) denote a predictor vector. For a fixed covariate value \(x\), the conditional cumulative distribution function and conditional density are

\[ F(y \mid x) = \Pr(Y \le y \mid X = x), \qquad f(y \mid x) = \frac{\partial}{\partial y} F(y \mid x), \]

and the conditional quantile function is

\[ Q(\tau \mid x) = \inf\{y : F(y \mid x) \ge \tau\}, \qquad 0 < \tau < 1. \]

In CausalMixGPD, this conditional distribution is the primary object of inference. Scalar summaries such as means or upper quantiles are derived from the fitted distribution rather than modeled separately. That design is important because it keeps different summaries internally coherent and avoids issues such as quantile crossing that can arise when quantiles are estimated one at a time.

Bulk model: Dirichlet process mixtures

Below a tail threshold, the package uses a DPM regression model. If \(k(\cdot \mid x; \theta)\) denotes a chosen kernel family, the bulk conditional density can be written as

\[ f_{\mathrm{DP}}(y \mid x) = \int k(y \mid x; \theta) \, dH(\theta) = \sum_{j=1}^{\infty} w_j k(y \mid x; \theta_j), \]

where \(H\) denotes a random mixing distribution having a Dirichlet process prior such that the mixing proportions obey \(w_j \ge 0\) and \(\sum_{j} w_j = 1\), and the kernel parameters can be either fixed or covariate-dependent (Ferguson 1973; Sethuraman 1994; Neal 2000). The mixture formulation provides much flexibility in dealing with multimodality, skewness, nonuniform local dispersion, and potential latent grouping without requiring any prespecification of the number of mixture components.

From an implementation standpoint, the bulk model specification is largely driven by four main inputs: kernel, backend, components, and param_specs. The first two control the kernel function \(k(\cdot)\) and the underlying representation, respectively, the third is the internal truncation level, and the fourth decides whether the kernel parameters are fixed, drawn from priors, or regressed.

Tail model: generalized Pareto exceedances

When the upper tail requires explicit modeling, the package replaces the part of the density above a threshold \(u(x)\) by a GPD component. For exceedances over the threshold, define

\[ Z = Y - u(x) \quad \text{given } Y > u(x). \]

Under the peaks-over-threshold framework, the conditional excess distribution is approximated by a GPD when the threshold is sufficiently high (Balkema and Haan 1974; Pickands 1975; Coles 2001). For \(y > u(x)\), with scale parameter \(\sigma(x) > 0\) and shape parameter \(\xi\), the GPD cumulative distribution function is

\[ F_{\mathrm{GPD}}(y \mid u(x), \sigma(x), \xi) = \begin{cases} 1 - \left(1 + \xi \frac{y-u(x)}{\sigma(x)}\right)^{-1/\xi}, & \xi \ne 0, \\ 1 - \exp\left(-\frac{y-u(x)}{\sigma(x)}\right), & \xi = 0. \end{cases} \]

The support is \(y \ge u(x)\) when \(\xi \ge 0\), and \(u(x) \le y \le u(x) - \sigma(x)/\xi\) when \(\xi < 0\). The shape parameter controls the tail regime: negative values imply a finite upper endpoint, \(\xi = 0\) gives an exponential-type tail, and positive values correspond to heavy tails.

Spliced bulk-tail construction

The key statistical idea in the package is the spliced model. Let \(F_{\mathrm{DP}}(\cdot \mid x; \Theta)\) and \(f_{\mathrm{DP}}(\cdot \mid x; \Theta)\) denote the bulk CDF and density, and let \(\Phi\) collect the tail parameters. Write

\[ p_u(x; \Theta) = F_{\mathrm{DP}}(u(x) \mid x; \Theta) \]

for the bulk probability mass below the threshold. Then the spliced conditional CDF is

\[ F(y \mid x; \Theta, \Phi) = \begin{cases} F_{\mathrm{DP}}(y \mid x; \Theta), & y \le u(x), \\ p_u(x; \Theta) + \{1 - p_u(x; \Theta)\} F_{\mathrm{GPD}}(y \mid x; \Phi), & y > u(x), \end{cases} \]

and the corresponding conditional density is

\[ f(y \mid x; \Theta, \Phi) = \begin{cases} f_{\mathrm{DP}}(y \mid x; \Theta), & y \le u(x), \\ \{1 - p_u(x; \Theta)\} f_{\mathrm{GPD}}(y \mid x; \Phi), & y > u(x). \end{cases} \]

This formulation preserves the flexibility of the mixture in the central region while forcing the tail to obey an EVT-motivated form. The splice is continuous at \(u(x)\) by construction, and the tail density is normalized by the bulk survival at the threshold.

The quantile function follows the same two-part logic. For quantile levels below the threshold mass, the package inverts the bulk CDF. For levels above the threshold mass, the probability index is rescaled to the exceedance scale:

\[ \widetilde\tau(x; \Theta) = \frac{\tau - p_u(x; \Theta)}{1 - p_u(x; \Theta)}. \]

Then the upper-tail quantiles are obtained from the GPD quantile function evaluated at \(\widetilde\tau(x; \Theta)\). This distinction explains why upper quantile predictions are one of the main places where a spliced fit can differ meaningfully from a bulk-only fit.

Means, restricted means, and existence conditions

When the kernel mean exists and the tail satisfies \(\xi < 1\), the conditional mean of the spliced model is finite. In that case the one-arm interface can report ordinary means through predict(type = "mean"). When the tail is too heavy for a finite mean or when the user prefers a bounded summary, the package provides restricted means through predict(type = "rmean"), which returns

\[ E\{\min(Y, c) \mid X = x\} \]

for a user-supplied cutoff \(c\). Restricted means are often easier to interpret in heavy-tailed settings because they remain finite and are less dominated by a few very large posterior draws. The software paper emphasizes this distinction explicitly in the one-arm prediction section.

Conditional density and survival summaries

The package article notes that density and survival summaries are returned by evaluating the corresponding functional draw by draw and then summarizing over retained MCMC draws. We illustrate that at three representative predictor profiles formed from empirical quartiles.

x_new <- as.data.frame(lapply(
  onearm_dat[, c("x1", "x2", "x3")],
  stats::quantile,
  probs = c(0.25, 0.50, 0.75),
  na.rm = TRUE
))
rownames(x_new) <- c("q25", "q50", "q75")
x_new

#>              x1          x2            x3
#> q25 -0.59851252 -0.57151934 -0.7224762264
#> q50  0.02874559 -0.01658997  0.0004347721
#> q75  0.69727001  0.41546653  0.5516024827

y_grid <- seq(min(onearm_dat$y), quantile(onearm_dat$y, 0.99), length.out = 120)
pred_dens <- predict(
  fit_spliced,
  newdata = x_new,
  y = y_grid,
  type = "density",
  interval = "credible",
  level = 0.95,
  store_draws = FALSE
)
plot(pred_dens, type = "density", facet = "covariate")

y_grid <- seq(min(onearm_dat$y), quantile(onearm_dat$y, 0.99), length.out = 120)
pred_surv <- predict(
  fit_spliced,
  newdata = x_new,
  y = y_grid,
  type = "survival",
  interval = "credible",
  level = 0.95,
  store_draws = FALSE
)
plot(pred_surv)

These summaries show how the fitted conditional distribution changes across covariate profiles. In practice, that helps determine whether predictors act mostly through a location shift, through dispersion, or through upper-tail behavior.

Quantiles, means, and restricted means

The one-arm workflow is especially useful when the scientific target is not just the density itself but one or more functionals of the conditional distribution. The article highlights quantiles and mean-type summaries as central examples, with restricted means remaining available even when the ordinary mean is unstable or undefined under a heavy tail.

pred_quant <- predict(
  fit_spliced,
  newdata = x_new,
  type = "quantile",
  index = c(0.25, 0.50, 0.75, 0.90, 0.95),
  interval = "credible",
  level = 0.95,
  store_draws = FALSE
)

plot(pred_quant)

pred_mean <- predict(
  fit_spliced,
  newdata = x_new,
  type = "mean",
  interval = "hpd",
  level = 0.90,
  store_draws = FALSE
)

plot(pred_mean)

cutoff_val <- as.numeric(stats::quantile(onearm_dat$y, 0.95))
pred_rmean <- predict(
  fit_spliced,
  newdata = x_new,
  type = "rmean",
  cutoff = cutoff_val,
  interval = "hpd",
  level = 0.90,
  store_draws = FALSE
)

knitr::kable(pred_rmean$fit_df, format = "html", row.names = FALSE, digits = 4)

The restricted mean is often a good complement to upper-tail quantiles. It preserves the original response scale, is always finite for finite cutoffs, and gives a more stable summary when a few extreme draws can dominate the ordinary mean.

Model specification options

kernel selects the bulk family. Positive-support choices include "gamma", "lognormal", "invgauss", and "amoroso"; real-line kernels include "normal", "laplace", and "cauchy". The support of the response should guide this choice.

backend determines how the DPM is represented internally. The main options are "sb" for stick-breaking and "crp" for a Chinese restaurant process representation.

components controls the truncation size used for the finite approximation to the infinite mixture.

GPD is implicit in the wrapper choice: dpmgpd() activates the spliced tail, whereas dpmix() fits the bulk-only model.

Parameter and prior customization

param_specs is the central argument for advanced model control. Each kernel or tail parameter can be assigned one of three modes: a fixed value, a prior-driven random value, or a covariate-linked regression structure. This is the main route for changing whether a threshold, scale, or kernel parameter is held constant or allowed to vary with predictors.

alpha_random controls whether the DPM concentration parameter is treated as random under its prior or fixed.

epsilon controls downstream truncation and summary behavior for very small mixture components.

MCMC and monitoring options

mcmc contains the iteration controls niter, nburnin, thin, nchains, and seed.

monitor chooses how much of the posterior state is retained. The usual values are "core" and "full".

monitor_latent = TRUE requests latent allocation labels in the retained output, and monitor_v = TRUE stores stick-breaking fractions for SB fits.

Prediction options

predict() supports type = "density", "survival", "quantile", "mean", "rmean", "median", "sample", and "fit".

For density and survival summaries, the argument y supplies the evaluation grid.

For quantile summaries, the argument index supplies the probability levels.

For restricted means, the argument cutoff defines the truncation point and nsim_mean can control Monte Carlo evaluation when applicable.

interval can be set to "credible", "hpd", or NULL, and level controls the posterior interval coverage level.