Kernel Quantile Regression with Mixed Data Types

npqreg computes a kernel quantile regression estimate of a one (1) dimensional dependent variable on \(p\)-variate explanatory data, given a set of evaluation points, training points (consisting of explanatory data and dependent data), and a bandwidth specification using the methods of Li and Racine (2008) and Li, Lin and Racine (2013). A bandwidth specification can be a condbandwidth object, or a bandwidth vector, bandwidth type and kernel type.

Usage

npqreg(bws, ...)

# S3 method for class 'formula'
npqreg(bws, data = NULL, newdata = NULL, ...,
       se = FALSE)


# S3 method for class 'condbandwidth'
npqreg(bws,
       txdat = stop("training data 'txdat' missing"),
       tydat = stop("training data 'tydat' missing"),
       exdat,
       tau = 0.5,
       gradients = FALSE,
       tol = 1.490116e-04,
       small = 1.490116e-05,
       itmax = 10000,
       ...,
       se = FALSE)

# Default S3 method
npqreg(bws, txdat, tydat, nomad = FALSE, ...,
       se = FALSE)

# S3 method for class 'qregression'
predict(object, se.fit = FALSE, ...)

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

Arguments

Data, Bandwidth Inputs And Formula Interface

These arguments identify the bandwidth specification, formula/data interface, and training data.

bws

a bandwidth specification. This can be set as a condbandwidth object returned from an invocation of npcdistbw, or as a vector of bandwidths, with each element \(i\) corresponding to the bandwidth for column \(i\) in txdat. If specified as a vector, then additional arguments will need to be supplied as necessary to specify the bandwidth type, kernel types, and so on.

data

This argument belongs to the formula interface (including supported formula-bearing bandwidth-object refits). Native calls must supply their vectors, matrices or data frames directly, or use with(data, ...). Supplying data, including explicit data = NULL, to a native route that does not consume it is an error; see np-data-interfaces.

an optional data frame, list or environment (or object coercible to a data frame by as.data.frame) containing the variables in the model. If not found in data, the variables are taken from environment(bws), typically the environment from which npcdistbw was called.

txdat

a \(p\)-variate data frame of explanatory data (training data) used to calculate the regression estimators. Defaults to the corresponding training data recoverable from the bandwidth object in the current session.

tydat

a one (1) dimensional numeric or integer vector of dependent data, each element \(i\) corresponding to each observation (row) \(i\) of txdat. Defaults to the corresponding training data recoverable from the bandwidth object in the current session.

object

an object of class "qregression" returned by npqreg.

x

an object of class "qregression" returned by npqreg.

Local-Polynomial Degree And Bandwidth Search

This argument controls the recommended automatic local-polynomial NOMAD route, which jointly selects continuous polynomial degree and bandwidths when conditional-distribution bandwidths are computed inside npqreg.

nomad

logical shortcut passed through to npcdistbw when bandwidths are computed inside npqreg. When TRUE, the conditional-distribution bandwidth route fills any missing values among regtype, search.engine, degree.select, bernstein.basis, degree.min, degree.max, degree.verify, and bwtype with the recommended automatic local-polynomial degree-and-bandwidth NOMAD preset documented in npcdistbw. Additional NOMAD tuning arguments such as nomad.nmulti may also be supplied through ...; nmulti remains the outer restart count while nomad.nmulti controls inner native crs NOMAD API multistarts within each outer restart. After fitting, inspect fit$bws$nomad.shortcut on the returned object fit to see the normalized shortcut metadata.

Evaluation Data And Returned Quantities

These arguments control where the quantile regression is evaluated and which fitted quantities are returned.

exdat

a \(p\)-variate data frame of points on which the regression will be estimated (evaluation data). By default, evaluation takes place on the data provided by txdat.

gradients

a logical value indicating that you want gradients of the conditional quantile with respect to the conditioning variables computed and returned in the resulting qregression object. Defaults to FALSE.

newdata

An optional data frame in which to look for evaluation data. If omitted, the training data are used.

se.fit

logical value. If TRUE, predict.qregression returns a list with components fit and se.fit; otherwise it returns fitted conditional quantiles.

se

a scalar logical controlling asymptotic standard-error computation. Defaults to FALSE. With TRUE, quantile standard errors are computed, and available gradient standard errors are computed when gradients=TRUE. Point quantiles and gradients do not depend on this switch. se(fit) and gradients(fit, se=TRUE) only extract stored results; they never initiate search or refitting. To add omitted uncertainty, use npqreg(bws=fit$bws, tau=fit$tau, se=TRUE) (also gradients=TRUE for gradient errors), with the original training/evaluation data if needed. The bandwidth object avoids repeating search. predict(..., se.fit=TRUE) explicitly requests its required uncertainty during the existing prediction fit.

tau

a numeric scalar or vector specifying the quantile probability or probabilities \(\tau\). Defaults to 0.5.

Quantile Solver Controls

These arguments control the one-dimensional numerical quantile extraction step.

itmax

integer maximum number of iterations allowed in the one-dimensional quantile refinement. Defaults to 10000.

small

absolute interval-width allowance, in response units, added to the relative stopping tolerance in the one-dimensional quantile refinement. Defaults to 1.490116e-05 (approximately 1000*sqrt(.Machine$double.eps)).

tol

relative tolerance on the one-dimensional quantile location refinement. Refinement stops when the bracket width is at most tol times the initial observed response range plus small, or when no distinct floating-point midpoint remains. The relative scale does not depend on the response origin. When rescaling the response, rescale small in the same units to retain the same accuracy contract. Defaults to 1.490116e-04 (approximately 10000*sqrt(.Machine$double.eps)).

Additional Arguments

Further arguments are passed to the bandwidth-selection counterpart, prediction/evaluation route, or plot route as appropriate.

...

additional arguments supplied to npcdistbw when npqreg computes bandwidths internally, or arguments needed to interpret a numeric bws vector. This is where bandwidth selection controls such as bwmethod, bwtype, and bwscaling, kernel/support controls such as cxkertype, cykertype, cxkerorder, cykerorder, cxkerbound, and cykerbound, categorical kernel controls such as uxkertype, uykertype, oxkertype, and oykertype, search controls such as nmulti and scale.factor.search.lower, and local-polynomial/NOMAD controls such as regtype, degree, bernstein.basis, degree.select, and nomad.nmulti are supplied. In predict.qregression, additional arguments are passed to npqreg for evaluation with the stored bandwidth object; common examples are newdata, native exdat, and tau. In plot.qregression, additional arguments are passed through the package plot route; common controls include tau, gradients, output, legend, and graphics arguments. See npcdistbw and plot.np for the complete bandwidth-selection and plotting argument surfaces.

Details

Conditional fitting and replay use the literal categorical kernels for both explanatory and response variables. Near-upper categorical bandwidths are not replaced by approximate constant factors. The same target is used for levels, gradients, standard errors, conditional hats, and reconstructed plot uncertainty. This conditional-family policy does not change direct regression or unconditional kernel-sum evaluation and does not modify the caller's np.largelambda option.

With a formula-derived bandwidth object, supplying only tydat replaces the evaluated response in retained-training-sample order: one row per observation retained after the original subset and missing-value selection. Original excluded rows are not restored; row names or equal lengths do not imply original-data order. The effective na.action is then applied jointly to the retained predictors and replacement response; any new omission indices refer to this retained sample. The selected bandwidths remain fixed, and the returned object retains the replacement sample for subsequent refits. Use data for full formula-data replacement, not together with a partial response replacement. Explicit native training arguments retain their separate positional contract.

For native (non-formula) objects, named newdata columns are matched to the retained variable names and reordered as needed. Extra columns are ignored; missing, duplicate or ambiguous names are errors. An unnamed matrix or vector remains positional for a single evaluation role. Explicit native evaluation arguments retain their positional contract and take precedence over newdata.

For a native bandwidth object, newdata supplies exdat; explicitly supplied native evaluation arguments take precedence. This evaluates at the requested rows without repeating bandwidth selection.

Formula variables of class ts, including lagged variables, are aligned on their common time support before fitting. Other variables must already have the length of that intersection. Supply explicitly aligned rows when using native data arguments instead of a formula.

Local-polynomial calculations inherit numerical ridging from npcdist. See that estimator's help for the solve policy and its distinct conditional uncertainty formula.

Documentation guide: see np.kernels for kernels, np.options for global options, and plot, plot.np for plotting options.

Given a conditional distribution bandwidth object, npqreg estimates the conditional distribution function \(F(y|x)\) and extracts the requested conditional quantile. For \(0 < \tau < 1\), the conditional quantile at probability \(\tau\) is $$ q_\tau(x) = \inf\{y : F(y|x) \ge \tau\}. $$ Equivalently, \(q_\tau(x)\) is a quasi-inverse of the conditional distribution in the sense of Nelsen (2006): an inverse agrees with \(F\) on the range of \(F\), while outside that range the generalized inverse is defined by the lower endpoint at which \(F\) reaches or exceeds the requested probability. Numerically, npqreg inverts the selected conditional distribution estimator represented by bws. This includes the selected bandwidth type, kernels, local-polynomial regression type, selected polynomial degree, basis, and Bernstein-basis setting inherited from npcdistbw. If the bandwidth object was selected with nomad=TRUE, the returned conditional-distribution bandwidth object is an LP object: its regtype/regtype.engine metadata identify the selected local-polynomial route and its degree/degree.engine metadata record the selected continuous-coordinate polynomial degree. npqreg, predict, and plot reuse this stored LP metadata; plotting additional tau values does not recompute or downgrade the selected degree. The one-dimensional inversion is carried out over the observed support of the dependent variable using the same selected conditional CDF estimator that is later used for quantile standard errors and gradients. The arguments tol, small, and itmax control this one-dimensional refinement. Conditional quantile inversion requires a monotone dependent-variable conditional distribution estimate; therefore npqreg supports cykerorder = 2 for this route and rejects higher-order dependent-variable CDF kernels, which can be nonmonotone.

During this inversion, an external predictor row identified by the selected CDF as having no kernel support returns NA for its quantiles and requested standard errors and gradients. Supported rows retain their original order, including repeated evaluation points and vector-tau requests. An unavailable categorical counterfactual affects only the corresponding contrast, not an otherwise supported base quantile. One warning summarizes affected external rows after the public evaluation or plot call completes. Other numerical failures, and required training or bootstrap evaluations, retain their strict error behavior.

The existing point-only shortcut for an exactly constant response is unchanged: it returns that response value without evaluating the conditional CDF. External predictor support is not checked by this shortcut; requests for standard errors or gradients still enter the conditional evaluation machinery.

Let \(f(y|x) = \partial F(y|x)/\partial y\) denote the conditional density. The asymptotic standard error of the conditional quantile is computed by the first-order delta method, $$ se\{\hat q_\tau(x)\} = \frac{se\{\hat F(\hat q_\tau(x)|x)\}} {\hat f(\hat q_\tau(x)|x)} , $$ using the selected conditional distribution standard-error machinery and the selected conditional density evaluated at the fitted quantile. This corresponds to the quantile variance expression in Li, Lin and Racine (2013).

Conditional-quantile plots support first derivatives only. The default and explicit gradient.order=1 retain that target; other derivative orders are rejected before plot fitting or drawing. Higher polynomial degrees in the fitted bandwidth object do not add higher-order quantile derivatives. This restriction does not apply to nplsqreg or conditional-distribution gradient plots.

A scalar quantile can be displayed as a surface by explicitly supplying perspective=TRUE when there are two continuous or numeric-labelled ordered conditioning variables. Omitting perspective, or supplying perspective=FALSE, retains the usual slice layout. Vector quantiles and gradient plots also retain their slice layout. Surface uncertainty is calibrated over the complete displayed grid for that one tau. Both the base renderer and renderer="rgl", view="fixed" support level surfaces and their interval meshes. Returned evaluation rows retain the trained predictor names and ordered-factor support; numeric axis coordinates are used only for drawing. The data overlay can be disabled with data.overlay=FALSE.

Local-polynomial conditional-response factors use local response moments as a working covariance approximation conditional on the accepted design and response map. At fixed nonshrinking bandwidths they need not equal the full empirical-caseweight derivative, which also perturbs the design map. Their leading interpretation requires smooth shrinking neighborhoods and nondegenerate local response variation.

Asymptotic plot(..., band="simultaneous") calibrates a joint Gaussian maximum over the retained grid within one panel and one tau; calibration is not joint across quantiles or panels. Levels use inverse-CDF influence rows at their own fitted quantiles. Continuous gradients use the selected CDF-derivative influence and the existing leading plug-in standard error. Positive-degree categorical effects use paired own-quantile endpoint influences with shared observation columns. Local-constant and degree-zero categorical effects retain their different public target, the same-root CDF contrast divided by the selected density: their joint bands include density and root-evaluation terms in a full local delta calculation, with separately computed joint-band standard errors. band="all" retains existing pointwise and Bonferroni widths even when the simultaneous calculation uses this different scale.

Local-constant and local-polynomial routes retain the fitted kernel, bandwidth type, degree and basis. Quantile inversion still requires a second-order response kernel. A valid joint grid requires finite positive slopes, interior unclamped quantiles, negligible inversion residuals, identifiable local influence and regular kernel/support geometry. The nearest-neighbor same-root calculation holds the realized neighbor counts and geometry fixed in its empirical-mass approximation, while retaining movement of a generalized-neighbor query radius when differentiating in the response coordinate. It does not claim an unconditional law for estimating the counts or radii.

These are selected-smoothing variability approximations, not smoothing-bias corrections or continuum confidence bands. An asymptotic coverage interpretation additionally requires independent observations, consistent local covariance and marginal scales, negligible smoothing/ridge bias and non-dominating influence rows. The continuous-gradient leading approximation requires shrinking explanatory bandwidth \(h_x\), with \(h_x^2/h_y\longrightarrow0\), and sufficient local derivative sample size; fixed nonshrinking smoothing retains a working approximation. A growing or data-dependent grid needs a corresponding uniform linearization and Gaussian-maximum approximation. Unsupported or nonregular grids receive a specific unavailable-band notice. gaussian.draws controls simulation precision and random.seed restores the caller's RNG state.

If gradients=TRUE, npqreg also computes gradients of the conditional quantile with respect to the conditioning variables for which gradients are defined. Differentiating \(F(q_\tau(x)|x) = \tau\) gives $$ \nabla_x q_\tau(x) = -\frac{F_x(q_\tau(x)|x)}{f(q_\tau(x)|x)}, $$ where \(F_x(y|x)\) is the derivative of the same selected conditional distribution estimator with respect to \(x\). For regtype="lc", this uses the local-constant conditional-gradient machinery; for regtype="ll" it uses the canonical local-polynomial degree-one route; and for regtype="lp" it uses the selected or supplied degree vector. The corresponding first-order gradient standard errors are computed componentwise as $$ se\{\nabla_x \hat q_\tau(x)\} = \frac{se\{\hat F_x(\hat q_\tau(x)|x)\}} {\hat f(\hat q_\tau(x)|x)} . $$

This is a plug-in estimate of the population implicit-CDF identity, using the selected density in the denominator. It is not in general the exact derivative of the fitted finite-smoothing quantile curve. In particular, bounded and beta response kernels use observation-centred CDFs and target-centred densities, so numerical differentiation of the fitted quantiles need not reproduce the reported plug-in gradients. The same selected-density convention is used for the delta-method standard errors.

For categorical conditioning variables in positive-degree local-polynomial fits (including local-linear fits), first-difference standard errors retain the paired covariance between the two category endpoints. Each endpoint uses its own fitted quantile and positive conditional-density slope in the implicit-CDF calculation. The selected basis, smoothing parameters and realized neighbour radii are held fixed. This requested-only calculation supports fixed, generalized nearest-neighbour and adaptive nearest-neighbour bandwidths; a nonpositive density slope or nonsmooth/nonfinite influence produces NA uncertainty with a notice.

For continuous conditioning variables, these gradient standard errors are componentwise leading asymptotic approximations. They propagate the conditional CDF-gradient uncertainty but do not jointly propagate fluctuations of the density denominator and estimated quantile location. Their use requires appropriate smoothness, positive-density and relative bandwidth-rate conditions. They are not a full finite-smoothing covariance or a joint covariance across coordinates or tau values.

In a heterogeneous fitted local-polynomial degree vector, a degree-zero coordinate has no available first-derivative estimate and is returned as NA; the established all-degree-zero compatibility route is unchanged. Its plot panel remains present but empty when another derivative or categorical first-difference effect is available; a request with no available component is rejected.

When npqreg is called without an explicit bws object, it first computes conditional distribution bandwidths using npcdistbw and stores them in the returned object's bws component. If a scalar tau was used initially and additional quantiles are later desired as fitted objects, reuse those selected bandwidths directly, for example npqreg(bws = fit$bws, tau = c(0.25, 0.5, 0.75)). If the goal is only to inspect additional quantiles graphically, use plot(fit, tau = c(0.25, 0.5, 0.75)); this reuses the stored bandwidths and recomputes only the one-dimensional quantile extraction step for the requested tau values. Vector-tau plots are overlaid and include a legend; use legend=FALSE, legend=NULL, or a legend=list(...) control to suppress or customize it.

The predict method follows the usual S3 newdata convention. For formula fits, supply a data frame of evaluation covariates via predict(fit, newdata=...). For non-formula fits, newdata is translated to the native evaluator argument exdat when exdat is not supplied. The native exdat argument remains available for advanced workflows, contains already transformed predictor coordinates for formula fits, and takes precedence if both newdata and exdat are supplied. If tau is omitted in predict, the fitted object's stored tau value is used.

Value

npqreg returns a qregression object. The generic functions fitted (or quantile), se, predict, and gradients extract stored estimated values, asymptotic standard errors, and gradients, or compute predictions from the retained bandwidth object. predict uses the object's stored tau value by default; supply tau= to override it. Furthermore, the functions summary and plot support objects of this type. The returned object has the following components:

eval

evaluation points

quantile

estimation of the quantile regression function (conditional quantile) at the evaluation points. If tau has length greater than one this is an evaluation-by-tau matrix.

quanterr

asymptotic standard errors of the quantile regression estimates, obtained from the conditional distribution standard error and the estimated conditional density at the fitted quantile. If tau has length greater than one this is an evaluation-by-tau matrix. This component is NULL when se=FALSE.

quantgrad

gradients of the conditional quantile with respect to the conditioning variables at each evaluation point, when gradients=TRUE. If tau has length greater than one this is an evaluation-by-gradient-by-tau array.

quantgerr

asymptotic standard errors for gradients, when gradients=TRUE and se=TRUE. If tau has length greater than one this is an evaluation-by-gradient-by-tau array. This component is NULL when se=FALSE.

se

whether asymptotic uncertainty was requested

tau

the quantile probability or probabilities computed

fit.controls

resolved tol, small, and itmax controls used for numerical quantile extraction. Prediction and quantile plots (including bootstrap fits and centers) reuse these values unless the corresponding argument is explicitly supplied. Older objects without this component use the historical defaults.

Book And Method Pointers

The conditional quantile target is the generalized inverse \(q_\tau(x)=\inf\{y:F(y\mid x)\ge \tau\}\) of the conditional distribution. The standard errors and gradients described above are first-order delta-method quantities evaluated using the same selected conditional CDF, conditional density, bandwidths, kernels, and local-polynomial degree inherited from the supplied npcdistbw object.

For book-length derivations, see Li and Racine (2007), Chapter 6 Conditional CDF and Quantile Estimation, especially Sections 6.3-6.5, and Racine (2019), Chapter 4 Conditional Probability Density and Cumulative Distribution Functions. The quasi-inverse terminology follows Nelsen (2006).

References

Aitchison, J. and C.G.G. Aitken (1976), “Multivariate binary discrimination by the kernel method,” Biometrika, 63, 413-420.

Hall, P. and J.S. Racine and Q. Li (2004), “Cross-validation and the estimation of conditional probability densities,” Journal of the American Statistical Association, 99, 1015-1026.

Koenker, R. W. and G.W. Bassett (1978), “Regression quantiles,” Econometrica, 46, 33-50.

Koenker, R. (2005), Quantile Regression, Econometric Society Monograph Series, Cambridge University Press.

Li, Q. and J.S. Racine (2007), Nonparametric Econometrics: Theory and Practice, Princeton University Press.

Li, Q. and J.S. Racine (2008), “Nonparametric estimation of conditional CDF and quantile functions with mixed categorical and continuous data,” Journal of Business and Economic Statistics, 26, 423-434.

Li, Q. and J. Lin and J.S. Racine (2013), “Optimal Bandwidth Selection for Nonparametric Conditional Distribution and Quantile Functions”, Journal of Business and Economic Statistics, 31, 57-65.

Nelsen, R.B. (2006), An Introduction to Copulas, Second Edition, Springer.

Wang, M.C. and J. van Ryzin (1981), “A class of smooth estimators for discrete distributions,” Biometrika, 68, 301-309.

Author

Tristen Hayfield tristen.hayfield@gmail.com, Jeffrey S. Racine racinej@mcmaster.ca

Usage Issues

If you are using data of mixed types, then it is advisable to use the data.frame function to construct your input data and not cbind, since cbind will typically not work as intended on mixed data types and will coerce the data to the same type.

For formula fits, character-valued missing-data policies are resolved through the stats model-frame owner; pass a caller-local policy as a function object rather than its name. With replacement training data, an explicit na.action, including NULL, takes precedence over the saved policy. Native fitting still requires complete rows, so a permissive policy does not guarantee the same retained sample as stats::model.frame.

When reusing a bandwidth object without explicit training or evaluation data, newly constructed bandwidth objects use retained training values. A fit with explicit replacement training data retains that sample in its returned bandwidth object; evaluation grids do not replace training data. Older objects without retained samples use their stored call environment when available and may require explicit training inputs after saveRDS() and readRDS(). Custom formula transformations used on new evaluation data must still be available.

See also

np.kernels, np.options, plot, plot.np, quantreg

Examples

if (FALSE) { # \dontrun{
# EXAMPLE 1 (INTERFACE=FORMULA): For this example, we compute a
# bivariate nonparametric quantile regression estimate for Giovanni
# Baiocchi's Italian income panel (see Italy for details)

data("Italy")
with(Italy, {

# Compute conditional distribution bandwidths and extract three
# conditional quantiles using the same selected bandwidths.

model.q <- npqreg(gdp~ordered(year), tau=c(0.25, 0.50, 0.75))

# Plot the overlaid quantiles.

plot(model.q)

# If a scalar tau was used first, additional quantiles can reuse the
# selected bandwidths without recomputing cross-validation. Use npqreg()
# when the additional fitted values are needed as an object, or plot()
# when graphical inspection is all that is desired.

model.med <- npqreg(gdp~ordered(year), tau=0.50)
model.q <- npqreg(bws=model.med$bws, tau=c(0.25, 0.50, 0.75))
plot(model.med, tau=c(0.25, 0.50, 0.75))

})

# EXAMPLE 1 (INTERFACE=DATA FRAME): For this example, we compute a
# bivariate nonparametric quantile regression estimate for Giovanni
# Baiocchi's Italian income panel (see Italy for details)

data("Italy")
with(Italy, {
data <- data.frame(ordered(year), gdp)

# First, compute the likelihood cross-validation bandwidths (default).
# Note - this may take a few minutes depending on the speed of your
# computer...

bw <- npcdistbw(xdat=ordered(year), ydat=gdp)

# Note - numerical search for computing the quantiles will take a
# minute or so...

model.q <- npqreg(bws=bw, tau=c(0.25, 0.50, 0.75))

plot(model.q)

})
} # }