Package {mixqr}

mixqr: An Extensible Framework for Mixtures of Quantile and Expectile Regressions

Description

An extensible expectation-maximization (EM) framework for finite mixtures of quantile regressions (clusterwise / mixture-of-experts quantile regression). A single EM substrate with an engine/extension contract carries a family of capabilities: the core free-weight mixture of Wu and Yao (2016) – a fast asymmetric-Laplace path and the nonparametric kernel-density EM with components constrained to have their tau-quantile equal to zero; expectile and M-quantile component-loss families; component-specific penalized variable selection; and joint multi-quantile estimation with a shared classification and non-crossing curves. The companion package mixqrgate adds location-varying gating.

Main entry points

• mixqr() – fit a mixture of quantile (or, via family=, expectile / M-quantile) regressions.

• mixqr_pen() – component-specific penalized variable selection.

• mixqr_nc() – joint multi-quantile estimation with non-crossing.

• mixqr_select() – choose the number of components by AIC/BIC.

• sim_mixqr2(), sim_mixqr3() – Wu & Yao simulation designs.

• register_mixqr_engine() – register a custom EM engine (the extension contract on which the capabilities are built).

Author(s)

Maintainer: Kailas Venkitasubramanian kailasv@gmail.com [copyright holder]

References

Wu, Q. and Yao, W. (2016). Mixtures of quantile regressions. Computational Statistics & Data Analysis 93, 162–176.

Hall, P. and Presnell, B. (1999). Density estimation under constraints. Journal of Computational and Graphical Statistics 8, 259–277.

Koenker, R. and Bassett, G. (1978). Regression quantiles. Econometrica 46, 33–50.

See Also

Useful links:

• https://github.com/kvenkita/mixqr

• https://kvenkita.github.io/mixqr/

• Report bugs at https://github.com/kvenkita/mixqr/issues

Number of active (non-zero) slopes per component.

Description

Number of active (non-zero) slopes per component.

Usage

active_slopes(beta, tol = 1e-08)

Construct the ALD engine.

Description

Construct the ALD engine.

Usage

ald_engine_ctor(tau, m, control = mixqr_control())

Arguments

Value

A mixqr_engine list of closures.

Initial densities for the ALD engine given responsibilities (for iteration 0).

Description

Initial densities for the ALD engine given responsibilities (for iteration 0).

Usage

ald_init_density(p, X, y, beta, tau)

Flag near-coincident component slope vectors (near-violation of Thm 2.1).

Description

Uses a scale-RELATIVE distance ⁠||b_a - b_b|| / (||b_a|| + ||b_b||)⁠ so the threshold is meaningful regardless of slope magnitude.

Usage

check_distinct_slopes(beta, tol = 0.001)

Arguments

Value

The minimum pairwise relative slope distance (invisibly); warns if ⁠< tol⁠.

Component coefficients of a mixqr fit

Description

Component coefficients of a mixqr fit

Usage

## S3 method for class 'mixqr'
coef(object, include_pi = FALSE, ...)

Arguments

Value

A ⁠(p+1) x m⁠ coefficient matrix (one column per component), or a list if include_pi = TRUE.

Confidence intervals for a mixqr fit

Description

Wald intervals from the fitted standard errors. Requires variance = "sparsity" or "stochEM"; "stochEM" is recommended (sparsity SEs are classification-conditional and understate uncertainty).

Usage

## S3 method for class 'mixqr'
confint(object, parm, level = 0.95, ...)

Arguments

Value

A matrix of lower/upper bounds.

Constrained KDE update for all components (eq. 2.5 unequal / eq. 2.8 equal).

Description

Constrained KDE update for all components (eq. 2.5 unequal / eq. 2.8 equal).

Usage

constrained_kde(
  e,
  p,
  tau,
  variant = c("unequal", "equal"),
  control = mixqr_control()
)

Arguments

Details

Exported as part of the mixqr extension API (see weighted_rq()).

Value

A length-m list of density objects, each with an eval closure, f0, grid, and a constraint_ok flag.

Coupled (cross-tau pooled) E-step -> single shared n x m posterior. p_ij proportional to pi_j times prod_l f_jl(e_ijl) raised to the weight w_l.

Description

Coupled (cross-tau pooled) E-step -> single shared n x m posterior. p_ij proportional to pi_j times prod_l f_jl(e_ijl) raised to the weight w_l.

Usage

coupled_estep(X, y, beta_arr, pi, sigma_mat, tau, wL)

Count within-component quantile crossings at the data anchors. For each observation and component, the L fitted quantiles x'beta_j(tau_l) should be nondecreasing in l; a decrease is a crossing.

Description

Count within-component quantile crossings at the data anchors. For each observation and component, the L fitted quantiles x'beta_j(tau_l) should be nondecreasing in l; a decrease is a crossing.

Usage

crossing_count(X, beta_arr)

K-fold cross-validated held-out negative predictive log-likelihood.

Description

Scores each m by the held-out mixture predictive density ⁠-sum_test log sum_j pi_j g_j(y_i - x_i' beta_j)⁠, using each engine's fitted component error densities g_j (ALD or KDE). Unlike a min-component check loss, this PENALISES complexity – extra components overfit the training density and predict held-out points worse – so it is a valid cross-m selector for either engine (the kdEM density is available even though it has no global likelihood).

Usage

cv_negloglik(
  formula,
  data,
  tau,
  m,
  engine,
  error_density,
  folds,
  nstart,
  control,
  weights
)

Draw a valid mixing-probability vector ~ N(pi_hat, multinomial cov) (eq. 3.4). Rejection-samples the first m-1 elements until a valid simplex point.

Description

Draw a valid mixing-probability vector ~ N(pi_hat, multinomial cov) (eq. 3.4). Rejection-samples the first m-1 elements until a valid simplex point.

Usage

draw_pi(pi, n, maxtry = 200L)

Objective used for multi-start root selection. ALD: log-likelihood (maximise). kdEM: negative total weighted check loss (so larger is always better).

Description

Objective used for multi-start root selection. ALD: log-likelihood (maximise). kdEM: negative total weighted check loss (so larger is always better).

Usage

em_objective(engine, state, X, y, tau, p)

Engine ethanol-combustion data (Brinkman 1981)

Description

Measurements from a single-cylinder automobile test engine burning ethanol, exhibiting two homogeneous regimes when the equivalence ratio is regressed on nitrous-oxide concentration – the engine example of Wu & Yao (2016, Fig. 5). Derived from the public-domain ethanol data redistributed in the lattice package (Brinkman 1981).

Usage

engine

Format

A data frame with 88 rows and 3 columns:

equivalence

equivalence ratio: richness of the air/ethanol mix (response).

nox

concentration of nitrous oxide (NO and NO2) in the exhaust, normalized by engine work (predictor).

compression

compression ratio of the engine.

Source

Brinkman, N. D. (1981). Ethanol fuel – a single-cylinder engine study of efficiency and exhaust emissions. SAE Transactions 90, 1410–1427. Redistributed via the lattice package ethanol dataset.

References

Wu, Q. and Yao, W. (2016). Mixtures of quantile regressions. CSDA 93, 162–176.

Examples

fit <- mixqr(equivalence ~ nox, data = engine, tau = 0.5, m = 2)
summary(fit)

Generic responsibility E-step from stored component densities. p_ij = pi_j g_j(e_ij) / sum_l pi_l g_l(e_il). (Wu & Yao eq. 2.2)

Description

Generic responsibility E-step from stored component densities. p_ij = pi_j g_j(e_ij) / sum_l pi_l g_l(e_il). (Wu & Yao eq. 2.2)

Usage

estep_from_density(state, X, y, tau, weights = NULL)

Asymmetric-normal working density with theta-expectile pinned at 0. g(e) proportional to exp(-w_theta(e) e^2 / (2 sigma^2)), w_theta(e)=theta if e>0 else 1-theta.

Description

Asymmetric-normal working density with theta-expectile pinned at 0. g(e) proportional to exp(-w_theta(e) e^2 / (2 sigma^2)), w_theta(e)=theta if e>0 else 1-theta.

Usage

expectile_density(e, sigma, tau)

Construct the expectile engine (shares the mixqr EM driver and S3 surface).

Description

Construct the expectile engine (shares the mixqr EM driver and S3 surface).

Usage

expectile_engine_ctor(tau, m, control = mixqr_control())

Nonparametric error density at 0 from (weighted) residuals.

Description

Wu & Yao (2016, p.166) read the sparsity f(0) off the kernel density estimate. Using the KDE value – rather than a parametric/ALD density – is essential for valid SEs when the true error density is non-ALD (e.g. bimodal), where the ALD f(0) can badly understate the variance.

Usage

f0_kde(e, w, h = NULL)

Fitted values, residuals and number of observations

Description

fitted() returns the classified conditional tau-quantile for each observation; residuals() the corresponding residuals (or, with type = "component", the full ⁠n x m⁠ residual matrix); nobs() the sample size.

Usage

## S3 method for class 'mixqr'
fitted(object, ...)

## S3 method for class 'mixqr'
residuals(object, type = c("hard", "component"), ...)

## S3 method for class 'mixqr'
nobs(object, ...)

Arguments

Value

A numeric vector or matrix.

Retrieve a registered mixqr engine constructor

Description

Retrieve a registered mixqr engine constructor

Usage

get_mixqr_engine(name)

Arguments

Value

The constructor function.

Silverman bandwidth for a (weighted) residual vector.

Description

Silverman bandwidth for a (weighted) residual vector.

Usage

hp_bandwidth(e, p, override = NULL)

Per-point Hall-Presnell / empirical-likelihood weights enforcing the constraint.

Description

Returns non-negative per-point weights w_i with sum(w_i) = 1 and ⁠sum(w_i v_i) = tau⁠ (so the resulting KDE has its tau-th quantile = 0), obtained by the KL-from-p tilt ⁠w_i propto p_i / (1 + lambda (v_i - tau))⁠ with lambda solving the moment constraint (Hall & Presnell 1999; Owen 2001). Returns NULL if the constraint is infeasible (tau outside the range of v).

Usage

hp_el_weights(p, v, tau)

Asymmetric Huber rho with cutoff k = c * sigma.

Description

Asymmetric Huber rho with cutoff k = c * sigma.

Usage

huber_rho(u, k)

Generate nstart initial responsibility matrices.

Description

Generate nstart initial responsibility matrices.

Usage

init_seeds(X, y, m, strategy = "ald", nstart = 20L, manual = NULL)

Arguments

Value

A list of nstart ⁠n x m⁠ responsibility matrices.

Construct the kdEM (Wu & Yao) engine.

Description

Construct the kdEM (Wu & Yao) engine.

Usage

kdem_engine_ctor(
  tau,
  m,
  control = mixqr_control(),
  variant = c("unequal", "equal")
)

Arguments

Value

A mixqr_engine list of closures.

Total weighted check loss – kdEM objective surrogate for root selection. ⁠sum_ij p_ij rho_tau(e_ij)⁠ (Wu & Yao have no likelihood; DESIGN sec. 2.4).

Description

Total weighted check loss – kdEM objective surrogate for root selection. ⁠sum_ij p_ij rho_tau(e_ij)⁠ (Wu & Yao have no likelihood; DESIGN sec. 2.4).

Usage

kdem_objective(state, X, y, tau, p)

One-hot-with-floor responsibility matrix from hard labels.

Description

One-hot-with-floor responsibility matrix from hard labels.

Usage

labels_to_p(labels, m, hot = 0.85)

List registered mixqr engines

Description

List registered mixqr engines

Usage

list_mixqr_engines()

Value

Character vector of engine names.

Log-likelihood, AIC and BIC of a mixqr fit

Description

Available for the parametric ALD engine (a genuine working likelihood); NA for the semiparametric kdEM engine, which has no global likelihood.

Usage

## S3 method for class 'mixqr'
logLik(object, ...)

## S3 method for class 'mixqr'
AIC(object, ..., k = 2)

## S3 method for class 'mixqr'
BIC(object, ...)

Arguments

Value

A "logLik" object (logLik) or a numeric criterion (AIC, BIC).

Build one constraint-preserving KDE density object from residual points.

Description

The public eval closure uses O(n + G) linear interpolation off the stored grid (DESIGN sec.2.2), not an O(n^2) kernel re-sum, so the E-step scales linearly. The grid density and f0 are computed exactly.

Usage

make_kde_density(e_pts, p_pts, h, tau, grid_n = 512L)

Fit a finite mixture of quantile regressions

Description

Estimates a finite mixture of tau-quantile regressions (clusterwise quantile regression) at a single quantile level. Two engines are available: a fast parametric asymmetric-Laplace mixture ("ald", genuine likelihood and AIC/BIC) and the kernel-density EM of Wu & Yao (2016) ("kdEM", nonparametric component error densities constrained to have their tau-th quantile = 0).

Usage

mixqr(
  formula,
  data,
  tau = 0.5,
  m = 2L,
  family = c("quantile", "expectile", "mquantile"),
  engine = "ald",
  error_density = c("unequal", "equal"),
  init = c("ald", "kmeans", "random", "manual"),
  nstart = 20L,
  control = mixqr_control(),
  weights = NULL,
  manual_init = NULL,
  variance = c("none", "sparsity", "stochEM"),
  vcontrol = mixqr_vcontrol(),
  ...
)

Arguments

Value

An object of class "mixqr".

Bias under asymmetric errors (Wu & Yao sec.6)

The semiparametric (kdEM) estimator solves estimating equations whose score I(a <= 0) - tau is not orthogonal to the nuisance tangent space of the unknown error densities, so it can be BIASED when component error densities are asymmetric and the clusters have imbalanced overlap (Wu & Yao 2016, sec.6, Fig.6). Watch fit$diagnostics$overlap (responsibility entropy; larger = more overlap = more bias risk) and cross-check against the parametric ald engine. Well-separated clusters and (near-)symmetric errors are the safe regime.

Standard errors

Sparsity SEs (variance = "sparsity", eq.3.3) are CONDITIONAL on the fitted classification and understate uncertainty (they are flagged as such by summary()). Use variance = "stochEM" (the stochastic-EM multiple-imputation estimator, Algorithm 3.1) for inference; it propagates classification and mixing uncertainty.

References

Wu, Q. and Yao, W. (2016). Mixtures of quantile regressions. CSDA 93, 162–176.

Examples

set.seed(1)
d <- sim_mixqr2(n = 300)
fit <- mixqr(y ~ x, data = d, tau = 0.5, m = 2, engine = "ald", nstart = 10)
fit
coef(fit)

Control parameters for mixqr()

Description

Control parameters for mixqr()

Usage

mixqr_control(
  tol = 1e-06,
  maxit = 500L,
  reltol = TRUE,
  bandwidth = NULL,
  kde_grid = 512L,
  label_order = "slope",
  distinct_tol = 0.001,
  trace = FALSE,
  seed = NULL
)

Arguments

Value

A list of class "mixqr_control".

Run one EM fit from a single initial responsibility matrix.

Description

Run one EM fit from a single initial responsibility matrix.

Usage

mixqr_em(
  engine,
  X,
  y,
  tau,
  m,
  p_init,
  error_density = "unequal",
  weights = NULL,
  control = mixqr_control()
)

Arguments

Value

A list with state, posterior, iter, converged, objective, objective_path.

The frozen mixqr engine contract (documentation only)

Description

A mixqr_engine is a list of closures with these signatures. The generic driver mixqr_em() owns convergence, multi-start, labelling and bookkeeping; engines supply only the statistical updates.

Details

estep(state, X, y, tau, weights)

returns an ⁠n x m⁠ responsibility matrix p whose rows sum to 1.

mstep_pi(p, weights)

returns a length-m mixing-probability vector.

mstep_beta(p, X, y, tau, weights, beta_prev)

returns a ⁠(ncol(X)) x m⁠ coefficient matrix (weighted quantile regression).

update_density(p, X, y, beta, tau, control)

returns a length-m list of density objects, each a list with at least eval (a function), f0 (density at 0) and type.

loglik(state, X, y, tau)

returns the log-likelihood or NA.

npar(state)

returns the number of free parameters.

state is a list with elements pi, beta, density, tau, m, error_density.

Largest lambda that zeros all slopes (top of the regularisation path).

Description

Largest lambda that zeros all slopes (top of the regularisation path).

Usage

mixqr_lambda_max(X, y, tau)

Fit a joint multi-tau mixture of quantile regressions (non-crossing, shared labels)

Description

Estimates a finite mixture of quantile regressions at a vector of quantile levels tau jointly, with one latent classification shared across all levels (a coupled E-step that pools the per-level component likelihoods) and guaranteed non-crossing of the within-component quantile curves. This closes the two problems Wu & Yao (2016, sec.5) leave open: within-component crossing and cross-tau classification ambiguity.

Usage

mixqr_nc(
  formula,
  data,
  m = 2L,
  tau = c(0.25, 0.5, 0.75),
  noncrossing = c("rearrange", "none"),
  class_coupling = c("pool", "median"),
  nstart = 10L,
  control = mixqr_control(),
  weights = NULL
)

Arguments

Value

An object of class c("mixqr_multitau", "mixqr") with a coefficient array coefficients (⁠[p+1, m, L]⁠), the shared posterior/classification, mix_prop, the tau grid, and a crossing diagnostic (raw crossings, 0 after rearrangement).

References

Wu, Q. and Yao, W. (2016). Mixtures of quantile regressions. CSDA 93, 162–176. Chernozhukov, V., Fernandez-Val, I. and Galichon, A. (2010). Quantile and probability curves without crossing. Econometrica 78, 1093–1125.

Examples

d <- sim_mixqr_cross(n = 160, seed = 1)
fit <- mixqr_nc(y ~ x, data = d, m = 2, tau = c(0.1, 0.25, 0.5, 0.75, 0.9))
fit$crossing      # raw crossings -> 0 after rearrangement
predict(fit)[1, , ]   # non-crossing quantiles for observation 1

Parameter names for the stacked beta-by-component then pi (first m-1) vector.

Description

Parameter names for the stacked beta-by-component then pi (first m-1) vector.

Usage

mixqr_par_names(fit)

Fit a penalized (variable-selecting) mixture of quantile regressions

Description

Component-specific penalized variable selection for the mixqr() model: each latent component gets its own sparse slope vector, a covariate active in one quantile-regression cluster and dropped in another. The weighted check-loss M-step gains a SCAD / adaptive-LASSO / LASSO / MCP penalty (intercept free); the penalty strength lambda is selected by a mixture BIC over a path, and components whose mixing weight falls below pi_min are pruned.

Usage

mixqr_pen(
  formula,
  data,
  tau = 0.5,
  m = 2L,
  penalty = c("SCAD", "aLASSO", "LASSO", "MCP"),
  lambda = NULL,
  nlambda = 40L,
  a = 3.7,
  alasso_gamma = 1,
  penalty_factor = NULL,
  pi_min = 0.02,
  nstart = 10L,
  control = mixqr_control(),
  weights = NULL
)

Arguments

Value

A "mixqr" object with an extra ⁠$selection⁠ slot (active covariates per component, chosen lambda, the BIC path, and df).

References

Khalili, A. and Chen, J. (2007). Variable selection in finite mixture of regression models. JASA 102, 1025–1038. Sherwood, B., Li, S. and Maidman, A. (2025). rqPen. R Journal.

Examples

set.seed(1)
n <- 250; p <- 8; x <- matrix(rnorm(n * p), n)
colnames(x) <- paste0("x", 1:p)
z <- rbinom(n, 1, 0.5)
y <- ifelse(z == 0, 1 + 2 * x[, 1], -1 + 2 * x[, 2]) + rnorm(n)
d <- data.frame(y = y, x)
fit <- mixqr_pen(y ~ ., data = d, tau = 0.5, m = 2, penalty = "SCAD")
selectedVars(fit)

Select the number of mixture components

Description

Fits mixqr over a range of component counts and scores them by an information criterion or cross-validated check loss.

Usage

mixqr_select(
  formula,
  data,
  tau = 0.5,
  m = 1:4,
  criterion = c("BIC", "AIC", "cv"),
  engine = "ald",
  error_density = c("unequal", "equal"),
  folds = 5L,
  weights = NULL,
  nstart = 20L,
  control = mixqr_control()
)

Arguments

Details

For criterion %in% c("AIC", "BIC") the score uses the parametric ALD working-likelihood (the semiparametric kdEM engine has no global likelihood; Wu & Yao 2016, p. 164). Note that AIC/BIC for the number of mixture components is heuristic: testing m vs m+1 places a component on the parameter-space boundary (pi_j = 0), so the usual penalty asymptotics do not hold (McLachlan & Peel 2000, ch. 6; Chen, Chen & Kalbfleisch 2001). The "cv" criterion (cross-validated held-out check loss) avoids the likelihood entirely and works for either engine.

Value

A list with table (data frame of scores by m), best (chosen m), criterion, and fit (the chosen model refit on all data).

Control parameters for stochastic-EM variance estimation (Algorithm 3.1)

Description

Control parameters for stochastic-EM variance estimation (Algorithm 3.1)

Usage

mixqr_vcontrol(B = 500L, burnin = 50L, thin = 1L, seed = NULL)

Arguments

Value

A list of class "mixqr_vcontrol".

Analytically-normed asymmetric-Huber working density factory for one component. The norming constant of exp(-w_tau(u) rho_H^k(u) / sigma) is computed in closed form (an asymmetric half-Gaussian core on the interval from -k to k plus exact exponential tails), which is exact for any tau and faster than grid quadrature.

Description

Analytically-normed asymmetric-Huber working density factory for one component. The norming constant of exp(-w_tau(u) rho_H^k(u) / sigma) is computed in closed form (an asymmetric half-Gaussian core on the interval from -k to k plus exact exponential tails), which is exact for any tau and faster than grid quadrature.

Usage

mquantile_density_obj(ej, pj, tau, cc)

Construct the M-quantile engine. control$mq_c sets the cutoff (default 1.345).

Description

Construct the M-quantile engine. control$mq_c sets the cutoff (default 1.345).

Usage

mquantile_engine_ctor(tau, m, control = mixqr_control())

Expectile M-step over all components -> (ncol(X) x m) coefficient matrix.

Description

Expectile M-step over all components -> (ncol(X) x m) coefficient matrix.

Usage

mstep_beta_expectile(p, X, y, tau, weights, beta_prev = NULL)

M-quantile M-step over all components.

Description

M-quantile M-step over all components.

Usage

mstep_beta_mquantile(p, X, y, tau, weights, beta_prev = NULL, cc = 1.345)

Penalized M-step over all components -> (ncol(X) x m) coefficient matrix.

Description

Penalized M-step over all components -> (ncol(X) x m) coefficient matrix.

Usage

mstep_beta_penalized(
  p,
  X,
  y,
  tau,
  weights,
  beta_prev = NULL,
  penalty = "SCAD",
  lambda = 0.1,
  a = 3.7,
  penalty_factor = NULL,
  gamma = 1
)

Weighted-QR M-step over all components -> (ncol(X) x m) coefficient matrix.

Description

Weighted-QR M-step over all components -> (ncol(X) x m) coefficient matrix.

Usage

mstep_beta_qr(p, X, y, tau, weights, beta_prev = NULL)

Weighted mixing-probability M-step (eq. 2.3).

Description

Weighted mixing-probability M-step (eq. 2.3).

Usage

mstep_pi_weighted(p, weights = NULL)

Single multivariate-normal draw via eigen-decomposition (no MASS).

Description

Single multivariate-normal draw via eigen-decomposition (no MASS).

Usage

mvrnorm_one(mu, Sigma)

Component ordering permutation for a label-switching constraint.

Description

Component ordering permutation for a label-switching constraint.

Usage

order_components(beta, label_order = "slope")

Arguments

Value

An integer permutation of the components.

Penalized engine at a fixed lambda (reuses ALD density, E-step, pi-update).

Description

Penalized engine at a fixed lambda (reuses ALD density, E-step, pi-update).

Usage

penalized_engine_ctor(
  tau,
  m,
  control,
  penalty,
  lambda,
  a = 3.7,
  penalty_factor = NULL,
  gamma = 1
)

Penalized weighted-QR M-step for one component via rqPen (single lambda). Returns a length-(p+1) coefficient vector (intercept first, unpenalised).

Description

Penalized weighted-QR M-step for one component via rqPen (single lambda). Returns a length-(p+1) coefficient vector (intercept first, unpenalised).

Usage

penalized_rq(
  X,
  y,
  tau,
  w,
  penalty,
  lambda,
  a = 3.7,
  penalty_factor = NULL,
  beta_prev = NULL,
  gamma = 1
)

Multinomial covariance block for pi[1:m-1].

Description

Multinomial covariance block for pi[1:m-1].

Usage

pi_cov(pi, n)

Plot a mixqr fit

Description

Draws the data coloured by hard classification with the component quantile-regression lines, and a panel of the estimated component error densities (mirroring Wu & Yao Figs. 1, 3–5).

Usage

## S3 method for class 'mixqr'
plot(x, which = c("both", "fit", "density"), ...)

Arguments

Value

Invisibly x.

Plot the non-crossing component quantile curves of a multi-tau fit

Description

Overlays, for each component, the (rearranged, non-crossing) fitted quantile curves at every tau against the first covariate.

Usage

## S3 method for class 'mixqr_multitau'
plot(x, ...)

Arguments

Value

x, invisibly.

Predict from a mixqr fit

Description

Predict from a mixqr fit

Usage

## S3 method for class 'mixqr'
predict(
  object,
  newdata = NULL,
  type = c("quantile", "quantile_byclass", "posterior", "class", "density"),
  ...
)

Arguments

Value

A matrix, vector, or list depending on type.

Classifying genuinely new x

"posterior", "class", and "quantile_byclass" require the RESPONSE in newdata, because component membership is inferred from the residual densities. This core therefore cannot assign a class to a covariate vector x alone; covariate-based gating (membership as a function of x) is provided by the location-varying-gating extension (QMM sub-project 03).

Predict non-crossing component quantiles from a multi-tau fit

Description

Predict non-crossing component quantiles from a multi-tau fit

Usage

## S3 method for class 'mixqr_multitau'
predict(object, newdata = NULL, ...)

Arguments

Value

An ⁠[n, m, L]⁠ array of fitted quantiles. With noncrossing = "rearrange" each component's L quantiles are sorted to be nondecreasing in tau (Chernozhukov et al. 2010), guaranteeing no crossing.

Register a mixqr EM engine

Description

An engine is a constructor ⁠function(tau, m, control)⁠ returning a list of closures with the FROZEN signatures documented in mixqr_engine_contract(). Sub-projects 03–06 of the QMM suite ship engines this way without forking the driver.

Usage

register_mixqr_engine(name, constructor)

Arguments

Value

Invisibly TRUE.

Permute mixture components to satisfy a label-switching constraint.

Description

Permute mixture components to satisfy a label-switching constraint.

Usage

relabel(state, label_order = "slope")

Arguments

Value

The state with components reordered.

Moore-Penrose-safe inverse via SVD, returning the numerical rank.

Description

Moore-Penrose-safe inverse via SVD, returning the numerical rank.

Usage

safe_solve_rank(A, tol = 1e-10)

k-means seed on standardised (x, y).

Description

k-means seed on standardised (x, y).

Usage

seed_kmeans(X, y, m)

Random seed (row-normalised uniform); never the degenerate constant 1/m.

Description

Random seed (row-normalised uniform); never the degenerate constant 1/m.

Usage

seed_random(n, m)

Active (selected) covariates per component of a penalized mixqr fit

Description

Active (selected) covariates per component of a penalized mixqr fit

Usage

selectedVars(object)

Arguments

Value

A named list (one entry per component) of selected covariate names.

Simulate the Wu & Yao two-component design (eq. 4.1-4.2)

Description

⁠Y = 10 - 10x + e⁠ (component 1) or ⁠Y = -10 + 10x + e⁠ (component 2), with ⁠pi = (0.5, 0.5)⁠, x ~ U(0, 1), and a bimodal error ⁠e ~ 0.5 N(-1, 1^2) + 0.5 N(2, 2^2)⁠ (median 0; median regression, tau = 0.5).

Usage

sim_mixqr2(n = 300L, pi1 = 0.5)

Arguments

Value

A data frame with columns x, y, z (true component label).

Simulate the Wu & Yao three-component design (eq. 4.3-4.4)

Description

Three components with ⁠pi = (1/3, 1/3, 1/3)⁠, covariates ⁠x1, x2 ~ U(0, 1)⁠, means ⁠-20 x1 - 20 x2⁠, 0, ⁠20 x1 + 20 x2⁠, and shifted-lognormal errors exp(N(1, sigma^2)) - exp(1) with ⁠sigma^2 in {1, 0.5, 0.25}⁠ (each median 0).

Usage

sim_mixqr3(n = 300L)

Arguments

Value

A data frame with columns x1, x2, y, z.

Simulate a two-component mixture whose per-tau fits cross (for O4 validation)

Description

One component is well-behaved; the other has a near-flat slope and strong heteroscedasticity (error scale growing across the design), so its independently fitted per-tau quantile lines cross in finite samples – exactly the within-component crossing Wu & Yao (2016, sec.5) flag as common "when the sample size is small." Use to demonstrate that noncrossing = "none" crosses while "rearrange" does not. (A valid location-scale model has non-crossing true quantiles; crossing here is the finite-sample estimation artefact the method repairs.)

Usage

sim_mixqr_cross(n = 160, seed = NULL)

Arguments

Value

A data frame with columns y, x, and the true label z.

Back-compatible name: the two-constant weight pair (see two_constant_weights()).

Description

Back-compatible name: the two-constant weight pair (see two_constant_weights()).

Usage

solve_hp_weights(A, B, C, D, tau)

Sparsity variance for one component (Koenker-Bassett / eq. 3.3). ⁠V(beta_j) = tau(1-tau)/f0^2 * solve(X' W X)⁠, W = diag(responsibility*weights).

Description

Sparsity variance for one component (Koenker-Bassett / eq. 3.3). ⁠V(beta_j) = tau(1-tau)/f0^2 * solve(X' W X)⁠, W = diag(responsibility*weights).

Usage

sparsity_beta_var(X, wj, f0, tau, label = NULL)

Stochastic-EM multiple-imputation variance (Wu & Yao Algorithm 3.1).

Description

Faithful P-step: draws beta ~ N(beta_hat, V_beta), pi ~ rejection-sampled normal (eq. 3.4), and the error density from a residual/case bootstrap, before refreshing the responsibilities. ⁠Var(theta) = V_W + (1 + 1/B) V_B⁠ (Little & Rubin 2002).

Usage

stoch_em(fit, X, y, weights, vcontrol = mixqr_vcontrol())

Summarize a mixqr fit

Description

Prints, per component, the coefficient estimates with standard errors, z-values and p-values (when a variance method was used), the mixing probabilities, and the separability / overlap diagnostics. Sparsity standard errors are flagged as classification-conditional.

Usage

## S3 method for class 'mixqr'
summary(object, ...)

Arguments

Value

An object of class "summary.mixqr".

Two-constant Hall-Presnell weights (Wu & Yao sec.3 simplification).

Description

Solve ⁠[A B; C D] [w1; w2] = [1; tau]⁠. Returns the pair and a feasibility flag (ok = TRUE only if both sides are populated, the system is well-conditioned, and both weights are non-negative).

Usage

two_constant_weights(A, B, C, D, tau)

Variance-covariance matrix of a mixqr fit

Description

The estimated parameter covariance (requires a fitted variance method).

Usage

## S3 method for class 'mixqr'
vcov(object, ...)

Arguments

Value

The variance-covariance matrix.

Sparsity standard errors for a fitted mixqr object (eq. 3.3).

Description

Conditional on the fitted classification; under-supported components trigger a warning. The chosen variance = "stochEM" is recommended for valid inference.

Usage

vcov_sparsity(fit, X, weights)

Weighted asymmetric-least-squares fit for one component (IRLS). Minimises sum_i w_i |theta - I(e_i<0)| (y_i - x_i'b)^2 (Newey-Powell eq. 1).

Description

Weighted asymmetric-least-squares fit for one component (IRLS). Minimises sum_i w_i |theta - I(e_i<0)| (y_i - x_i'b)^2 (Newey-Powell eq. 1).

Usage

weighted_expectile_ls(
  X,
  y,
  tau,
  w,
  beta_prev = NULL,
  maxit = 100L,
  tol = 1e-08
)

Weighted asymmetric-Huber fit for one component (IRLS).

Description

Weighted asymmetric-Huber fit for one component (IRLS).

Usage

weighted_mquantile(
  X,
  y,
  tau,
  w,
  cc = 1.345,
  beta_prev = NULL,
  maxit = 100L,
  tol = 1e-08
)

Weighted quantile regression for one mixture component

Description

Solves ⁠argmin_b sum_i w_i rho_tau(y_i - x_i' b)⁠ via quantreg::rq.wfit().

Usage

weighted_rq(X, y, tau, w, beta_prev = NULL)

Arguments

Details

Exported as part of the mixqr extension API: building blocks that companion packages (e.g. location-varying gating, non-crossing) reuse to build custom EM engines on the same component machinery.

Value

A coefficient vector of length ncol(X).