Kernel Regression with Mixed Data Types

npreg computes a kernel 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 method of Racine and Li (2004) and Li and Racine (2004). A bandwidth specification can be a rbandwidth object, or a bandwidth vector, bandwidth type and kernel type.

Usage

npreg(bws, ...)

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

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

# S3 method for class 'rbandwidth'
npreg(bws,
      txdat = stop("training data 'txdat' missing"),
      tydat = stop("training data 'tydat' missing"),
      exdat,
      eydat,
      gradient.order = 1L,
      gradients = FALSE,
      residuals = FALSE,
      se = FALSE,
      ...)

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 rbandwidth object returned from an invocation of npregbw, 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 npregbw 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. A factor response is also supported when the bandwidth object retains the fitted factor-response coding. A bandwidth object fitted to a numeric response cannot supply that coding and is rejected when reused with a factor response. Defaults to the corresponding training data recoverable from the bandwidth object in the current session.

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 these are computed inside npreg.

nomad

logical shortcut passed through to npregbw when bandwidths are computed inside npreg. When TRUE, the regression 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 npregbw. 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 regression is evaluated and which fitted quantities are returned.

se

a logical value indicating whether the residual-normalized sandwich standard errors described below are computed and retained. Defaults to FALSE. With se = FALSE, fitted values and any requested gradients are still computed, but the returned merr and gerr components are NULL; se() and gradients(..., se = TRUE) then fail with a message directing the user to refit or re-evaluate with se = TRUE. This option controls fit/evaluation output only and is not forwarded to bandwidth selection.

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.

eydat

a one (1) dimensional numeric or integer vector of the true values of the dependent variable. Optional, and used only to calculate the true errors. When evaluation predictors are supplied without evaluation outcomes, goodness-of-fit scores (including R2 and MSE) are unavailable and reported as NA; no additional training fit is performed for these diagnostics.

gradient.order

for regtype="lp" with gradients=TRUE, a positive integer (or integer vector with one entry per continuous predictor) specifying derivative order(s). Defaults to 1L. Orders exceeding degree for a variable are returned as NA, apart from the established all-degree-zero first-derivative compatibility route. Higher-order derivatives are available for continuous predictors when the requested order does not exceed the corresponding local-polynomial degree.

gradients

a logical value indicating that you want gradients computed and returned in the resulting npregression object. Defaults to FALSE. For regtype="lp", derivative components that are not defined (requested order exceeds the variable-specific polynomial degree) are returned as NA. For ordered and unordered predictors the gradient columns are first-order effects/contrasts where supported; gradient.order controls derivative order for continuous predictors. Plotting retains one panel per requested predictor: an unavailable continuous derivative is shown as an empty, correctly labelled panel when another continuous derivative or categorical effect is available. A request for which no component is available is rejected with the affected orders and degrees reported.

newdata

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

residuals

a logical value indicating that you want residuals computed and returned in the resulting npregression object. Defaults to FALSE.

y.eval

If newdata contains dependent data and y.eval = TRUE, np will compute goodness of fit statistics on these data and return them. Defaults to FALSE.

Additional Arguments

Further arguments are passed to npregbw when bandwidths are computed internally, or used to interpret a numeric bws vector.

...

additional arguments supplied to npregbw when npreg computes bandwidths internally, or arguments needed to interpret a numeric bws vector. This is where bandwidth selection controls such as bwmethod, bwtype, bwscaling, kernel/support controls such as ckertype, ckerorder, and ckerbound, categorical kernel controls such as ukertype and okertype, search controls such as nmulti and scale.factor.search.lower, and local-polynomial/NOMAD controls such as regtype, degree, basis, bernstein.basis, degree.select, and nomad.nmulti are supplied. See npregbw for the complete bandwidth-selection argument surface.

Details

With generalized nearest-neighbor bandwidths, query identity matters when selecting a radius: a training observation excludes its own occurrence, whereas explicit evaluation data remain external queries even if their values equal training observations. Consequently, predict(fit, newdata = training.data) can differ from fitted(fit). This does not omit the training contribution from the fitted estimator. See np.kernels for the radius convention.

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.

Documentation guide: see npregbw for bandwidth selection and search controls, np.kernels for kernels, np.options for global options, and plot, plot.np for plotting options.

When bws is omitted, the formula and default methods call npregbw first and pass bandwidth-selection arguments from ... to that call. When bws is already an rbandwidth object, npreg estimates with the stored bandwidth metadata in that object.

Argument groups for bandwidth selection are documented on npregbw. The most common workflow is to choose data and bandwidth inputs first, then bandwidth criterion and representation, then kernel/support controls, and finally local-polynomial/NOMAD controls when using polynomial-adaptive fits.

For S3 plotting help, see plot.np. You can list available plot methods with methods("plot").

Typical usages are (see below for a complete list of options and also the examples at the end of this help file)

    Usage 1: first compute the bandwidth object via npregbw and then
    compute the conditional mean:

    bw <- npregbw(y~x)
    ghat <- npreg(bw)

    Usage 2: alternatively, compute the bandwidth object indirectly:

    ghat <- npreg(y~x)

    Usage 3: modify the default kernel and order:

    ghat <- npreg(y~x, ckertype="epanechnikov", ckerorder=4)

    Usage 4: use the data frame interface rather than the formula
    interface:

    ghat <- npreg(tydat=y, txdat=x, ckertype="epanechnikov", ckerorder=4)
  

npreg implements a variety of methods for regression on multivariate (\(p\)-variate) data, the types of which are possibly continuous and/or discrete (unordered, ordered). The approach is based on Li and Racine (2003) who employ ‘generalized product kernels’ that admit a mix of continuous and discrete data types.

Three classes of kernel estimators for the continuous data types are available: fixed, adaptive nearest-neighbor, and generalized nearest-neighbor. Adaptive nearest-neighbor bandwidths change with each sample realization in the set, \(x_i\), when estimating the density at the point \(x\). Generalized nearest-neighbor bandwidths change with the point at which the density is estimated, \(x\). Fixed bandwidths are constant over the support of \(x\).

Data contained in the data frame txdat may be a mix of continuous (default), unordered discrete (to be specified in the data frame txdat using factor), and ordered discrete (to be specified in the data frame txdat using ordered). Data can be entered in an arbitrary order and data types will be detected automatically by the routine (see np for details).

A variety of kernels may be specified by the user. Kernels implemented for continuous data types include the second, fourth, sixth, and eighth order Gaussian and Epanechnikov kernels, and the uniform kernel. Beta associated kernels of orders 2, 4, 6, and 8 are also available for local-constant and local-polynomial regression. Unordered discrete data types use a variation on Aitchison and Aitken's (1976) kernel, while ordered data types use a variation of the Wang and van Ryzin (1981) kernel.

For ckertype="beta", suppose continuous predictor \(j\) has declared finite support \([a_j,b_j]\), width \(L_j=b_j-a_j\), public metric bandwidth \(h_j\), evaluation coordinate \(x_j\), and observation \(X_{ij}\). Define $$t_j=(x_j-a_j)/L_j,\quad z_{ij}=(X_{ij}-a_j)/L_j,\quad \tau_j=(L_j/h_j)^2.$$ The local-constant product weight is $$K_i(x)=\prod_j \left\{\frac{1}{L_j} \frac{z_{ij}^{t_j\tau_j}(1-z_{ij})^{(1-t_j)\tau_j}} {B(1+t_j\tau_j,1+(1-t_j)\tau_j)}\right\},$$ and the fitted conditional mean is the Nadaraya–Watson ratio $$\widehat m(x)=\frac{\sum_i K_i(x)Y_i}{\sum_i K_i(x)}.$$ This uses Chen's evaluation-centred beta associated-kernel contribution inside the established np local-constant ratio. It should not be confused with the Gasser–Muller beta regression smoother studied by Chen (2000). The support scale is already present in each contribution, and the implementation applies neither the legacy \(1/h_j\) factor nor the generic bounded-kernel normalizer. Although that normalizer cancels from ordinary local-constant numerator/denominator ratios, the beta kernel is not merely such a correction: its shape parameters change with the evaluation point.

At orders 4, 6, and 8, each displayed order-2 coordinate factor is replaced by the signed Igarashi–Kakizawa multiscale combination documented in np.kernels. The resulting product weights and denominator may be signed; no absolute-value, order-2 substitution, or positivity repair is applied.

Products and ratios are evaluated in the log domain, so representable beta regression results are retained even when every unscaled product would underflow to zero. In the local-constant case, \(\ell_i(x)=K_i(x)/\sum_r K_r(x)\). The reported standard error uses the residual-normalized sandwich formula below, also for beta kernels; it is not a local residual-variance estimate multiplied by an effective-sample factor.

The beta regression route requires ckerorder in c(2,4,6,8) and finite fixed or sample-range bounds for each continuous predictor. It supports local-constant, local-linear and local-polynomial fits, including mixed continuous and categorical predictors. Local-constant first derivatives requested with gradients=TRUE are analytic C implementations; their standard errors use the derivative of the normalized weights. At exact support endpoints, normalized numerator/denominator jumps are cancelled first. A genuine estimator jump is returned as a signed infinity with gerr=NA; an indeterminate numerical cancellation is returned as NA with a warning. This extended-real gradient is retained in fitted and plot-data objects. Plot methods omit the non-finite coordinate only from their temporary drawing geometry, so an endpoint singularity breaks the plotted line without changing the returned estimate. Manual and cross-validated bandwidth objects support fixed, generalized-NN, and adaptive-NN modes. An evaluation point with no supported observation weight also produces an explicit error.

Standard Errors

Ordinary regression standard errors use a residual-normalized sandwich for the actual fitted linear map, conditional on the selected bandwidths, polynomial specification and realized predictor geometry. Let \(S\) map training responses to their fitted values, let \(e=(I-S)Y\), and define $$q_i=\sum_j(\delta_{ij}-S_{ij})^2.$$ For \(q_i>0\), the normalized training residual is \(z_i=e_i/\sqrt{q_i}\). This is the pre-smoothing residual normalization in Ruppert, Wand, Holst, and Hössjer (1997, Section 1.1, equation (3)), not their separate variance-function smoothing procedure. If \(\ell_i(x)\) is the influence of response \(Y_i\) on the fitted mean at \(x\), the reported standard error is $$\operatorname{se}\{\widehat m(x)\} =\left\{\sum_i \ell_i(x)^2 z_i^2\right\}^{1/2}.$$ A derivative uses its corresponding influence row. A categorical finite difference uses the difference between the two endpoint influence rows, retaining their covariance. Its estimate and error use the same paired influence calculation rather than subtracting nearly equal fitted values. Requesting or omitting standard errors does not change the computation of a requested categorical effect.

For generalized nearest-neighbor bandwidths, categorical effects at fitted training rows retain that row's training identity when selecting the radius for both endpoints, for all regression types, polynomial degrees and bases. The identity is excluded only from radius selection, not from the fitted response weights. Explicit evaluation data use external-query radii, even when their coordinates equal training observations. Consequently fitted-row and explicitly evaluated effects can differ under this bandwidth type.

The training smoother and influence rows use the actual fitted basis and accepted regularized solve where needed. The residual-row norm is not generally \(1-S_{ii}\); that identity requires a symmetric, idempotent projection. Nor is \(q_i\) necessarily at most one. No additional degrees-of-freedom multiplier or positive variance floor is applied. The corresponding influence representation is used without retaining a full training-by-training matrix.

Treating the smoother as fixed, under independent homoskedastic errors and a mean reproduced by that smoother, \(E(e_i^2\mid X)=\sigma^2q_i\). The normalization corrects fitted residual noise contraction in that setting, not smoothing bias. With a general mean, squared residual bias also contributes. Under heteroskedasticity, \(e_i^2/q_i\) mixes donor variances and need not be unbiased for the variance of observation \(i\). Finite-sample calibration under arbitrary bias, heteroskedasticity or model selection is not guaranteed.

A zero residual-row norm leaves no residual information for that donor. Requested errors depending on such unidentified information are NA, with a warning; finite point estimates remain available. A separately certified zero influence does not require that donor's variance. Small positive residual norms are not replaced by zero or by a numerical variance floor.

This convention applies to fixed, generalized-NN and adaptive-NN bandwidths and beta kernels on supported regression routes. External evaluation errors use the training smoother's normalized residuals, not residuals centred on each external evaluation mean. Bandwidth- and degree-selection uncertainty is not incorporated. se(), derivative-error accessors, prediction standard errors and asymptotic plot intervals consume the corresponding retained errors. Pivotal npsigtest statistics use these gradient errors, recomputed from each bootstrap response when bootstrapping.

Numerical Ridging

Local-linear and local-polynomial moment solves use numerical ridging when needed, motivated by Seifert and Gasser (2000). The package uses bounded implementation-specific ridge policies, not that paper's complete data-adaptive ridge selector. Ridging stabilizes the fit; it is distinct from the residual normalization above, which uses the actual accepted smoother, including any ridge adjustment. A local system with fewer nonzero-weight donor rows than basis coefficients is structurally rank deficient and uses this same ridge policy for fitting and related hat operators. Bootstrap multiplicities do not create additional independent donor rows. This check does not assert full rank when the donor count reaches the basis dimension.

When bandwidths are obtained with regtype="lp", C-level npreg supports heterogeneous continuous polynomial degrees via degree. The basis selector currently supports basis="glp", "additive", and "tensor". For continuous predictors with degree vector \(d\), additive basis size is \(1+\sum_j d_j\), tensor basis size is \(\prod_j (d_j+1)\), and GLP uses admissible multi-indices \(\alpha\) with \(\alpha_j \le d_j\) and \(0<\sum_j \alpha_j \le \max_j d_j\) plus an intercept. Thus, for \(d=(2,1)\), GLP contains \(1,x_1,x_1^2,x_2,x_1x_2\); the additive basis omits \(x_1x_2\), while the tensor basis additionally contains \(x_1^2x_2\).

The optional flag bernstein.basis controls basis construction: FALSE (default) uses raw local-polynomial powers. These can be numerically unstable, particularly above degree one or when predictor offsets are large relative to their spread. The option bernstein.basis=TRUE is provided to improve numerical conditioning; it does not remedy genuine collinearity. Centring/scaling continuous predictors is another explicit user choice. Existing conditioning warnings and errors assess the selected basis; numerical rank deficiency need not imply genuinely redundant predictors. No automatic basis switch is made. For basis="glp", the historical selector TRUE uses a degree-graded, orthonormal shifted-Legendre representation. On the unit-scaled support \(u\in[0,1]\), its order-\(k\) factor is \(\sqrt{2k+1}P_k(2u-1)\), where \(P_k\) is a Legendre polynomial. This is an exact degree-\(k\) linear combination of Bernstein polynomials and, because the GLP multi-index set is downward closed, spans exactly the same heterogeneous polynomial space as the raw representation. It is evaluated by a deterministic Legendre recurrence, not by sample-dependent Gram–Schmidt or QR orthogonalization. Additive and tensor bases retain their fixed-coordinate-degree Bernstein representation. The homogeneous degree-0 and degree-1 cases remain equivalent to lc and ll, respectively. For continuous predictors, gradients=TRUE returns the derivative order(s) requested by gradient.order; the associated gerr entries are standard errors for those same derivative orders. Both raw and shifted-Legendre representations report fitted values and derivatives on the original response/predictor scale. In mixed-data GLP settings, unordered and ordered predictors contribute first-order effects/contrasts rather than higher-order derivatives. When npregbw(..., regtype="lp") is used with degree.select="manual", the degree vector remains fixed user input. When degree.select != "manual", npregbw can jointly select polynomial degree and bandwidth using either the cached cell-search backend or the direct search.engine="nomad"/"nomad+powell" route described in npregbw; the latter follows Hall and Racine (2015). For practitioners who want that recommended route without spelling out all LP tuning arguments, npreg(..., nomad=TRUE) and npregbw(..., nomad=TRUE) expand missing settings to the same documented automatic-LP NOMAD preset. Explicit incompatible settings fail fast rather than being silently rewritten. The direct NOMAD backend is provided by the suggested package crs, so install crs before using search.engine="nomad", "nomad+powell", or nomad=TRUE. For bernstein.basis=TRUE, evaluation points for continuous predictors must lie within training support; use bernstein.basis=FALSE for extrapolation. For regtype="ll" and regtype="lp", the training continuous design is checked for rank deficiency and extreme condition number before estimation proceeds.

The use of compactly supported kernels or the occurrence of small bandwidths can lead to numerical problems for the local linear estimator when computing the locally weighted least squares solution. To overcome this problem we rely on a form or ‘ridging’ proposed by Cheng, Hall, and Titterington (1997), modified so that we solve the problem pointwise rather than globally (i.e. only when it is needed).

Value

npreg returns a npregression object. The generic functions fitted, residuals, se, predict, and gradients, extract (or generate) estimated values, residuals, asymptotic standard errors on estimates, predictions, and gradients, respectively, from the returned object. Furthermore, the functions summary and plot support objects of this type. The returned object has the following components:

With se.fit=TRUE, predict() returns fit, se.fit, df, and residual.scale. The last field uses training MSE in its historical units, not its square root or an evaluation MSE. It reuses a known training fit's MSE without an additional fit. It is NA when training MSE is unavailable, including an object created at external evaluation points or scored against explicit eydat, and external predictions with explicit training-data or model-frame overrides. A prediction evaluated on replacement training data can use that new training fit's MSE. Evaluation outcomes do not substitute for training residual information.

eval

evaluation points

mean

estimates of the regression function (conditional mean) at the evaluation points. An evaluation row whose required kernel normalization is exactly zero or non-finite is returned as NA; other rows remain available.

merr

standard errors of the regression function estimates, or NULL when se = FALSE. Undefined normalization rows are returned as NA. Requested errors depending on unidentified training residual information are also NA; finite point estimates are retained. See the Standard Errors subsection.

grad

estimates of the gradients at each evaluation point. Undefined normalization rows are returned as NA.

gerr

standard errors of the gradient estimates, or NULL when se = FALSE. Undefined normalization rows are returned as NA. Requested errors depending on unidentified training residual information are also NA, separately from certified zero effects. See the Standard Errors subsection.

resid

if residuals = TRUE, in-sample or out-of-sample residuals where appropriate (or possible)

R2

coefficient of determination (Doksum and Samarov (1995))

MSE

mean squared error

MAE

mean absolute error

MAPE

mean absolute percentage error

CORR

absolute value of Pearson's correlation coefficient

SIGN

fraction of observations where fitted and observed values agree in sign

Book And Method Pointers

The regression target is the conditional mean \(m(x)=E[Y\mid X=x]\). Local-constant estimation can be written schematically as \(\hat m(x)=\sum_i W_i(x)Y_i\), where the weights are normalized generalized product-kernel weights over the continuous, unordered, and ordered components of \(X\). Local-polynomial regression replaces this scalar weighted average by a weighted least squares fit in powers \(b(X_i)\) of the original continuous coordinates, with kernel weights centered at \(x\). If \(\widehat\beta_x\) denotes that local weighted fit, the reported level and derivative are respectively \(b(x)^T\widehat\beta_x\) and \(D^\gamma b(x)^T\widehat\beta_x\). For the downward-closed GLP space this uncentered representation is numerically equivalent, up to floating-point solution paths, to fitting powers of \(X_i-x\) and reading the corresponding intercept and derivative coefficients.

For book-length derivations, see Li and Racine (2007), Chapter 2 Regression, especially Sections 2.1, 2.2, 2.4, and 2.5, and Chapter 4 Kernel Estimation with Mixed Data, especially Sections 4.2 and 4.4. The later workflow treatment is Racine (2019), Chapter 6 Conditional Mean Function Estimation, including the local polynomial and mixed-data marginal-effect discussion.

References

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

Cheng, M.-Y. and P. Hall and D.M. Titterington (1997), “On the shrinkage of local linear curve estimators,” Statistics and Computing, 7, 11-17.

Chen, S.X. (1999), “Beta kernel estimators for density functions,” Computational Statistics and Data Analysis, 31, 131-145.

Chen, S.X. (2000), “Beta kernel smoothers for regression curves,” Statistica Sinica, 10, 73-91.

Igarashi, G. and Kakizawa, Y. (2020), “Higher-order bias corrections for asymmetric kernel density estimators,” Journal of Nonparametric Statistics, 32(3), 716–752. doi:10.1080/10485252.2020.1770754.

Cohen, A. and G. Migliorati (2018), “Multivariate approximation in downward closed polynomial spaces,” in R. Cools and D. Nuyens (eds.), Contemporary Computational Mathematics, Springer, 233-282.

Dung, D. and T. Ullrich (2011), “Whitney type inequalities for local anisotropic polynomial approximation,” Journal of Approximation Theory, 163, 1590-1605.

Fan, J. and I. Gijbels (1996), Local Polynomial Modelling and Its Applications, Chapman and Hall.

Doksum, K. and A. Samarov (1995), “Nonparametric estimation of global functionals and a measure of the explanatory power of covariates in regression,” The Annals of Statistics, 23 1443-1473.

Hall, P. and Q. Li and J.S. Racine (2007), “Nonparametric estimation of regression functions in the presence of irrelevant regressors,” The Review of Economics and Statistics, 89, 784-789.

Hall, P. and J.S. Racine (2015), “Infinite Order Cross-Validated Local Polynomial Regression,” Journal of Econometrics, 185, 510-525.

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

Li, Q. and J.S. Racine (2004), “Cross-validated local linear nonparametric regression,” Statistica Sinica, 14, 485-512.

Pagan, A. and A. Ullah (1999), Nonparametric Econometrics, Cambridge University Press.

Racine, J.S. and Q. Li (2004), “Nonparametric estimation of regression functions with both categorical and continuous data,” Journal of Econometrics, 119, 99-130.

Ruppert, D. and M.P. Wand (1994), “Multivariate locally weighted least squares regression,” The Annals of Statistics, 22, 1346-1370.

Ruppert, D., M.P. Wand, U. Holst, and O. Hössjer (1997), “Local Polynomial Variance-Function Estimation,” Technometrics, 39(3), 262-273. doi:10.1080/00401706.1997.10485117 .

Seifert, B. and T. Gasser (2000), “Data Adaptive Ridging in Local Polynomial Regression,” Journal of Computational and Graphical Statistics, 9(2), 338-360. doi:10.1080/10618600.2000.10474884 .

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. With explicit replacement training data, an explicit na.action overrides the retained policy, including NULL; any incomplete rows subsequently excluded by native fitting are recorded in the effective training sample. User formula and custom-function environments remain available and can retain unrelated caller objects in serialized files. See npregbw for storage and missing-value policy details.

See also

np.kernels, np.options, plot, plot.np loess

Examples

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

data("Italy")
with(Italy, {

# First, compute the least-squares cross-validated bandwidths for the
# local constant estimator (default).

bw <- npregbw(formula=gdp~ordered(year))

# Now take these bandwidths and fit the model and gradients

model <- npreg(bws = bw, gradients = TRUE)

summary(model)

# Use plot() to visualize the regression function, add bootstrap
# error bars, and overlay the data on the same plot.

# Note - this may take a minute or two depending on the speed of your
# computer due to bootstrapping being conducted (<ctrl>-C will
# interrupt). Note - nothing will appear in the graphics window until
# all computations are completed (if you use
# errors="asymptotic" the figure will instantly appear).

if (interactive()) plot(bw, errors="bootstrap")
points(ordered(year), gdp, cex=.2, col="red")

})

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

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

data("Italy")
with(Italy, {

# First, compute the least-squares cross-validated bandwidths for the
# local constant estimator (default).

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

# Now take these bandwidths and fit the model and gradients

model <- npreg(bws = bw, gradients = TRUE)

summary(model)

# Use plot() to visualize the regression function, add bootstrap
# error bars, and overlay the data on the same plot.

# Note - this may take a minute or two depending on the speed of your
# computer due to bootstrapping being conducted (<ctrl>-C will
# interrupt). Note - nothing will appear in the graphics window until
# all computations are completed (if you use
# errors="asymptotic" the figure will instantly appear).

if (interactive()) plot(bw, errors="bootstrap")
points(ordered(year), gdp, cex=.2, col="red")

})

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# EXAMPLE 2 (INTERFACE=FORMULA): For this example, we compute a local
# linear fit using the AIC_c bandwidth selection criterion. We then plot
# the estimator and its gradient using asymptotic standard errors.

data("cps71", package = "np")

bw <- npregbw(logwage~age, regtype="ll", bwmethod="cv.aic", data = cps71)

# Next, plot the regression function...

if (interactive()) plot(bw, errors="asymptotic")
with(cps71, points(age, logwage, cex=.2, col="red"))

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# Next, plot the derivative...

if (interactive()) plot(bw, errors="asymptotic", gradient=TRUE)

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# EXAMPLE 2 (INTERFACE=DATA FRAME): For this example, we compute a local
# linear fit using the AIC_c bandwidth selection criterion. We then plot
# the estimator and its gradient using asymptotic standard errors.

data("cps71", package = "np")

bw <- npregbw(xdat=cps71$age, ydat=cps71$logwage, regtype="ll", bwmethod="cv.aic")

# Next, plot the regression function...

if (interactive()) plot(bw, errors="asymptotic")
with(cps71, points(age, logwage, cex=.2, col="red"))

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# Next, plot the derivative...

if (interactive()) plot(bw, errors="asymptotic", gradient=TRUE)

# Sleep for 5 seconds so that we can examine the output...

if (interactive()) Sys.sleep(5)

# EXAMPLE 3 (INTERFACE=FORMULA): For this example, we replicate the
# nonparametric regression in Maasoumi, Racine, and Stengos
# (2007) (see oecdpanel for details). Note that X is multivariate
# containing a mix of unordered, ordered, and continuous data types. Note
# - this may take a few minutes depending on the speed of your computer.

data("oecdpanel")
with(oecdpanel, {

bw <- npregbw(formula=growth~
              factor(oecd)+
              factor(year)+
              initgdp+
              popgro+
              inv+
              humancap)

if (interactive()) plot(bw, errors="asymptotic")

})

# EXAMPLE 3 (INTERFACE=DATA FRAME): For this example, we replicate the
# nonparametric regression in Maasoumi, Racine, and Stengos
# (2007) (see oecdpanel for details). Note that X is multivariate
# containing a mix of unordered, ordered, and continuous data types. Note
# - this may take a few minutes depending on the speed of your computer.

data("oecdpanel")
with(oecdpanel, {

y <- growth
X <- data.frame(factor(oecd), factor(year), initgdp, popgro, inv, humancap)

bw <- npregbw(xdat=X, ydat=y)

if (interactive()) plot(bw, errors="asymptotic")

})


# EXAMPLE 4 (INTERFACE=FORMULA): Experimental data - the effect of
# vitamin C on tooth growth in guinea pigs
#
# Description:
#
#     The response is the length of odontoblasts (teeth) in each of 10
#     guinea pigs at each of three dose levels of Vitamin C (0.5, 1, and
#     2 mg) with each of two delivery methods (orange juice or ascorbic
#     acid).
#
# Usage:
#
#     ToothGrowth
#
# Format:
#
#     A data frame with 60 observations on 3 variables.
#
#       [,1]  len   numeric  Tooth length
#       [,2]  supp  factor   Supplement type (VC or OJ).
#       [,3]  dose  numeric  Dose in milligrams.

library("datasets")
with(ToothGrowth, {

# Note - in this example, there are six cells.

bw <- npregbw(formula=len~factor(supp)+ordered(dose))

# Now plot the partial regression surfaces with bootstrapped
# nonparametric confidence bounds

if (interactive()) plot(bw, errors="bootstrap", band="simultaneous")

})

# EXAMPLE 4 (INTERFACE=DATA FRAME): Experimental data - the effect of
# vitamin C on tooth growth in guinea pigs
#
# Description:
#
#     The response is the length of odontoblasts (teeth) in each of 10
#     guinea pigs at each of three dose levels of Vitamin C (0.5, 1, and
#     2 mg) with each of two delivery methods (orange juice or ascorbic
#     acid).
#
# Usage:
#
#     ToothGrowth
#
# Format:
#
#     A data frame with 60 observations on 3 variables.
#
#       [,1]  len   numeric  Tooth length
#       [,2]  supp  factor   Supplement type (VC or OJ).
#       [,3]  dose  numeric  Dose in milligrams.

library("datasets")
with(ToothGrowth, {

# Note - in this example, there are six cells.

y <- len
X <- data.frame(supp=factor(supp), dose=ordered(dose))

bw <- npregbw(X, y)

# Now plot the partial regression surfaces with bootstrapped
# nonparametric confidence bounds

if (interactive()) plot(bw, errors="bootstrap", band="simultaneous")

})

# EXAMPLE 5 (INTERFACE=FORMULA): a pretty 2-d smoothing example adapted
# from the R mgcv library which was written by Simon N. Wood.

set.seed(12345)

# This function generates a smooth nonlinear DGP

dgp.func <- function(x, z, sx=0.3, sz=0.4)
  { (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
                 0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
  }

# Generate 500 observations, compute the true DGP (i.e., no noise),
# then a noisy sample

n <- 500

x <- runif(n)
z <- runif(n)

xs <- seq(0, 1, length=30)
zs <- seq(0, 1, length=30)

X.eval <- data.frame(x=rep(xs, 30), z=rep(zs, rep(30, 30)))

dgp <- matrix(dgp.func(X.eval$x, X.eval$z), 30, 30)

y <- dgp.func(x, z)+rnorm(n)*0.1

# Prepare the screen for output... first, plot the true DGP

split.screen(c(2, 1))

screen(1)

persp(xs, zs, dgp, xlab="x1", ylab="x2", zlab="y", main="True DGP")

# Next, compute a local linear fit and plot that

bw <- npregbw(formula=y~x+z, regtype="ll", bwmethod="cv.aic")
fit <- fitted(npreg(bws=bw, newdata=X.eval))
fit.mat <- matrix(fit, 30, 30)

screen(2)

persp(xs, zs, fit.mat, xlab="x1", ylab="x2", zlab="y",
      main="Local linear estimate")

# EXAMPLE 5 (INTERFACE=DATA FRAME): a pretty 2-d smoothing example
# adapted from the R mgcv library which was written by Simon N. Wood.

set.seed(12345)

# This function generates a smooth nonlinear DGP

dgp.func <- function(x, z, sx=0.3, sz=0.4)
  { (pi**sx*sz)*(1.2*exp(-(x-0.2)^2/sx^2-(z-0.3)^2/sz^2)+
                 0.8*exp(-(x-0.7)^2/sx^2-(z-0.8)^2/sz^2))
  }

# Generate 500 observations, compute the true DGP (i.e., no noise),
# then a noisy sample

n <- 500

x <- runif(n)
z <- runif(n)

xs <- seq(0, 1, length=30)
zs <- seq(0, 1, length=30)

X <- data.frame(x, z)
X.eval <- data.frame(x=rep(xs, 30), z=rep(zs, rep(30, 30)))

dgp <- matrix(dgp.func(X.eval$x, X.eval$z), 30, 30)

y <- dgp.func(x, z)+rnorm(n)*0.1

# Prepare the screen for output... first, plot the true DGP

split.screen(c(2, 1))

screen(1)

persp(xs, zs, dgp, xlab="x1", ylab="x2", zlab="y", main="True DGP")

# Next, compute a local linear fit and plot that

bw <- npregbw(xdat=X, ydat=y, regtype="ll", bwmethod="cv.aic")
fit <- fitted(npreg(exdat=X.eval, bws=bw))
fit.mat <- matrix(fit, 30, 30)

screen(2)

persp(xs, zs, fit.mat, xlab="x1", ylab="x2", zlab="y",
      main="Local linear estimate")
} # }