General Purpose Plotting of Nonparametric Objects

Plotting is provided via plot S3 methods, which generate plots of nonparametric statistical objects such as regressions, quantile regressions, partially linear regressions, single-index models, densities and distributions, given training data and a bandwidth object. plot(...) is the supported public interface.

Usage

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

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

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

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

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

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

Arguments

Plot Object

This argument identifies the object to plot.

x: an object generated by the np package. The most convenient inputs are fitted estimator objects returned by invocations of npudens, npcdens, npreg, npconmode, npplreg, npindex, or npscoef. Bandwidth objects returned by npudensbw, npcdensbw, npregbw, npplregbw, npindexbw, or npscoefbw are also supported when plotting should construct the fitted object.

Plot Options Passed Through `...`

Named options passed via ... include:

Data and evaluation controls

data: replacement training data for the partially linear formula-bandwidth plot route, when explicit training inputs are omitted. Other plot routes do not consume this argument and reject it, including explicit NULL. Refit with the intended data before plotting; see np-data-interfaces.
grid.control: optional detailed grid controls created by np_grid_control.
neval: the number of evaluation points used for each plotted continuous margin, defaulting to 50L. Continuous numeric variables use this grid even when their observed values are integers or have few distinct values. Categorical margins use their levels. For npudens/npudist perspective surfaces with numeric ordered margins, the grid uses the ordered levels observed in the training data. Their one-dimensional categorical slices use all declared levels, including unused levels. A surface with two continuous margins has neval^2 evaluation points.
xdat, ydat, zdat: training data used when plotting bandwidth objects. xdat is the regressor/conditioning data, ydat is the response or conditional outcome data, and zdat is the auxiliary semiparametric covariate block when relevant.
xq, yq, zq: quantiles used to hold non-plotted variables fixed when constructing slices.
xtrim, ytrim, ztrim: trimming/expansion controls for the evaluation support.

Plot content controls

coef, coef.index: coefficient-function plotting controls for semiparametric families where supported.
common.scale: a logical value specifying whether multiple panels should share a common plotting scale when feasible.
gradient.order, gradients: derivative plotting controls for families that support gradients/derivatives. For npreg local-polynomial objects, gradient.order requests the continuous-predictor derivative order to generate on the plot evaluation grid; if omitted the first derivative/effect is plotted. Location-scale quantile-regression nplsqreg objects use the same coordinate-wise local-polynomial contract while preserving their existing tau-indexed output; no mixed partials are generated. When a requested continuous derivative order exceeds that coordinate's fitted degree, its values remain NA and its plot panel remains present but empty if another derivative or categorical effect is available. Plotting stops with an informative error only when no requested component is available. Ordered and unordered categorical panels remain first-order finite differences, denoted by $\Delta$ labels. For npconmode objects, gradients=TRUE plots stored class-probability gradients/effects for the selected response level, computed by npconmode(..., gradients=TRUE, level=...).
level: response level to plot for class-probability displays when supported, notably plot.conmode. If omitted, the default for npconmode objects is the base/reference response level levels(y)[1]. These plots require probabilities stored by npconmode(..., probabilities=TRUE). Currently plot.conmode displays stored object-fed probability/effect payloads only; probability/effect interval bands and surface grids are not yet implemented for this route.
layout: "auto" or TRUE requests automatic par(mfrow=...) management for multi-panel plots. "current" or FALSE preserves the user's current graphics layout.
output: one of "plot", "data", or "both". "both" is the public spelling for the historical "plot-data" behavior.
proper, proper.control, proper.method: controls for proper conditional density/distribution handling where implemented.
tau: the quantile index or indices used for quantile-regression-related plotting routes when relevant. A vector of tau values overlays the corresponding conditional quantile curves and reuses the stored bandwidth object. With errors="asymptotic" or errors="bootstrap", intervals are constructed separately for each requested tau; bootstrap work is shared across tau values within each resample. For npqreg objects whose conditional-distribution bandwidths were selected with nomad=TRUE, plotting reuses the stored LP regtype.engine and selected degree.engine metadata.

Interval and bootstrap controls

alpha: the nominal significance level used for interval/band construction. The default is 0.05, corresponding to a nominal variability level of $1-\alpha=0.95$. This nominal level describes the construction; it is not by itself a claim that smoothing bias is negligible or that the displayed bounds have literal confidence coverage for the unknown target.
B: the bootstrap replication count.
band: one of "pmzsd", "pointwise", "bonferroni", "simultaneous", or "all". The common bandwidth-engine default is "pmzsd"; direct npcopula and npconmode plots retain their method-specific "pointwise" default. The meaning of a band name depends on whether errors="bootstrap" or errors="asymptotic"; see Details.
boot.control: optional detailed bootstrap controls created by np_boot_control. This controls non-fixed bootstrap behavior, wild-bootstrap multipliers, and dependent-bootstrap block length.
bootstrap: the bootstrap resampling scheme. Regression-family plots support "wild", "inid", "fixed", and "geom"; density/distribution-family plots support "inid", "fixed", and "geom".
boxplot.outliers, factor.boxplot: boxplot-style bootstrap summary controls for factor/categorical displays when available.
center: one of "estimate" or "bias-corrected". When bootstrap errors are requested, center="bias-corrected" centers plotted uncertainty bands or bars on a bootstrap bias-corrected estimate. The package chooses the validated correction engine for the requested bootstrap route. Wild bootstrap regression routes use the standard wild-bootstrap correction. Pairwise and dependent-bootstrap routes "inid", "fixed", and "geom" use an oversmoothed pilot correction for fixed-bandwidth estimators where that route has been implemented and validated. Unsupported combinations fail clearly rather than silently using a weak correction. For one-dimensional regression curves, the ordinary estimate remains visible and the selected bias-corrected estimate is drawn as an additional center line when available. Fixed-bandwidth regression gradient plots use the same returned gradient-bias convention described below.
errors: one of "none", "bootstrap", or "asymptotic". The default is "none". Bootstrap errors can use either a symmetric normal-quantile rule based on bootstrap standard deviations or full bootstrap quantile/rank distributions, according to band. Asymptotic errors use pointwise standard errors and normal quantiles where available.

Rendering and view controls

data.overlay: logical; when TRUE, overlay raw data on regression, quantile-regression, single-index, partially linear, and smooth coefficient plot surfaces and expand plot limits to include the data. For npqreg level plots the data layer is drawn first by default, so quantile curves are overlaid on top; ordered and factor conditioning variables use the natural boxplot-style display from plot(). The option is ignored for derivative/gradient plots and coefficient plots. Overlay points default to small, light markers and can be customized with standard points() arguments passed via ... (for example pch, cex, col, and bg). Default is TRUE.
legend: legend control for legends drawn by the package, including vector-tau npqreg level plots and band="all" uncertainty legends. TRUE draws the default legend, FALSE, NULL, or NA suppresses it, a single character string supplies the legend position, and a list is merged into the default legend call for custom labels, placement, or styling. For plot routes with more than one package-drawn legend, named sublists such as legend=list(tau=list(x="topright"), bands=FALSE) may be used when a route distinguishes those legend roles.
data.rug: logical; when TRUE, draw a rug-style support cue for the plotted coordinates. For ordinary 1D base plots this behaves like a standard rug. For eligible 2D surface plots, the base renderer draws short floor-rug ticks and renderer="rgl" draws the analogous 3D floor-rug ticks on supported surface routes. In this rollout tranche, data.rug=TRUE is implemented across the current base plot families and across supported rgl surface routes; unsupported renderer/geometry combinations still fail clearly rather than being ignored silently. On mixed base plot layouts, rug marks are drawn for continuous panels and omitted for categorical panels. This option is distinct from data.overlay: data.rug shows where observations lie in the plotted covariate support, while data.overlay is the regression-only 3D response overlay used on supported surface plots. On supported regression surfaces, both options may be used together. Default is FALSE.
perspective: a logical value specifying whether bivariate continuous surfaces are displayed with persp-style plots when relevant.
phi, theta, view: view controls for bivariate continuous surface plots.
render.control: optional detailed rendering controls created by np_render_control, including band-versus-bar uncertainty display controls.
renderer: the surface-rendering backend for eligible bivariate continuous plots. "base" uses the existing base R persp path. "rgl" is currently a fail-closed backend for supported plain 2D surface routes in npreg/npregbw, npplreg/npplregbw, npscoef/npscoefbw, npudens/npudensbw, npudist/npudistbw, npcdens/npcdensbw, and npcdist/npcdistbw. In this tranche, renderer="rgl" requires perspective=TRUE, view="fixed" and output="plot"; unsupported combinations still error rather than falling back silently. The suggested packages rgl and htmlwidgets must be installed when renderer="rgl" is requested. When drawn, the rgl surface is displayed via the active rgl backend, which on this rollout path is typically the RStudio Viewer or a browser/WebGL window rather than the base graphics plot pane. The base and rgl renderers share the same viridis palette mapper and neutral-gray three-face reference-grid default. When col is omitted, the rgl surface uses a value-based hcl.colors(..., "viridis") gradient. When the inherited base-view defaults theta=0 and phi=20 are left unchanged, the rgl backend uses theta=0, phi=-70, and fov=55 to approximate the base perspective in the square interactive viewport. The widget automatically uses up to two backing-store pixels per CSS pixel on high-DPI displays. Axis titles retain rgl's native screen-facing text orientation while the scene is rotated. For the currently supported npreg, npplreg, and plain mean-surface npscoef surface routes, data.overlay=TRUE adds observed data as an interactive point cloud. For all currently supported rgl surface routes, errors != "none" adds interval surfaces; band="all" draws pointwise, simultaneous, and bonferroni band families with distinct, colorblind-safer colors when those bands are available for the route. Unsupervised and conditional density/distribution surface routes do not use observed-data point overlays in this tranche. For npscoef/npscoefbw, renderer="rgl" currently applies only to the plain mean-surface route; coefficient-mode plots remain on the base renderer. Common rgl display controls may be supplied directly in ... on supported routes, including surface arguments such as alpha, back, and front, and view arguments such as fov and zoom. Additional backend-specific arguments may be passed through using prefixed names in ..., namely rgl.persp3d.<arg>, rgl.view3d.<arg>, rgl.par3d.<arg>, rgl.grid3d.<arg>, and rgl.widget.<arg>, and rgl.legend3d.<arg>. For the new overlay layers, rgl.points3d.<arg> customizes observed-point overlays and rgl.surface3d.<arg> customizes interval surfaces. When band="all" is used on supported rgl routes, a matching construction-aware interactive legend is drawn using the same default labels as the base renderer, and rgl.legend3d.<arg> may be used to customize it, for example rgl.legend3d.bty="o" to restore a legend box or rgl.legend3d.magnify=1 to request the original legend texture resolution. By default, these legends follow the base-graphics uncertainty legend convention and are drawn without a surrounding box.

Graphical parameters

border, cex.axis, cex.lab, cex.main, cex.sub, col, lty, lwd, main, sub, type, xlab, xlim, ylab, ylim, zlab, zlim: standard graphical controls forwarded to the underlying plotting primitives when relevant.

Reproducibility controls

random.seed: an optional random seed used to make bootstrap plotting reproducible.

Additional Arguments

Further plotting controls are passed through ... to the relevant plot method.

...: additional arguments supplied to control plotting behavior or passed through to underlying plotting helpers where supported. Supported named plotting controls are summarized above; standard graphical parameters are forwarded where relevant.

Details

Documentation guide: see np.kernels for kernels and np.options for global options.

The preferred public interface is plot on fitted or bandwidth objects (e.g., plot(fit) or plot(bw)).

plot is a general purpose plotting routine for visually exploring objects generated by the np library, such as regressions, quantile regressions, partially linear regressions, single-index models, densities and distributions. There is no need to call plot directly: plotting is handled by class-specific S3 plot methods for objects generated by the np package.

Visualizing one and two dimensional datasets is a straightforward process. The default behavior of plot is to generate a standard 2D plot to visualize univariate data, and a perspective plot for bivariate data. When visualizing higher dimensional data, plot resorts to plotting a series of 1D slices of the data. For a slice along dimension $i$, all other variables at indices $j \ne i$ are held constant at the quantiles specified in the $j$th element of xq. The default is the median.

The slice itself is evaluated on a uniformly spaced sequence of $neval$ points. The interval of evaluation is determined by the training data. The default behavior is to evaluate from min(txdat[,i]) to max(txdat[,i]). The xtrim variable allows for control over this behavior. When xtrim is set, data is evaluated from the xtrim[i]th quantile of txdat[,i] to the 1.0-xtrim[i]th quantile of txdat[,i].

Furthermore, xtrim can be set to a negative value in which case it will expand the limits of the evaluation interval beyond the support of the training data, by measuring the distance between min(txdat[,i]) and the xtrim[i]th quantile of txdat[,i], and extending the support by that distance on the lower limit of the interval. plot uses an analogous procedure to extend the upper limit of the interval.

Numeric estimator output follows R's extended-real conventions: mathematically defined one-sided limits are retained as -Inf or Inf, while unavailable or indeterminate values are returned as NA_real_. Plotting does not rewrite the estimator object or values returned with output="data". Instead, a temporary graphics-only copy maps every non-finite coordinate to missing geometry. Consequently lines break at an infinite or unavailable point, interval segments at that point are omitted, and axis limits are computed from finite values. If a panel has no finite geometry, a finite fallback scale is used unless the user supplies ylim. This allows exact beta endpoint derivatives and other extended-real estimates to retain their mathematical meaning without passing invalid coordinates or limits to a graphics backend.

Plot variability-band constructions are:


  errors       band             construction
  bootstrap    pmzsd            center +/- z_(1-alpha/2) * bootstrap SD
  bootstrap    pointwise        pointwise bootstrap quantiles
  bootstrap    bonferroni       Bonferroni-adjusted bootstrap quantiles
  bootstrap    simultaneous     rank-based joint bootstrap band
  bootstrap    all              all three bootstrap quantile/rank families
  asymptotic   pmzsd/pointwise  center +/- z_(1-alpha/2) * asymptotic SE
  asymptotic   bonferroni       normal rule with alpha replaced by alpha/m
  asymptotic   simultaneous     finite-grid Gaussian rule for supported targets
  asymptotic   all              marginal families plus Gaussian where supported

where $m$ is the number of evaluation points used in the plotted curve/surface ($m=\texttt{neval}$ for univariate curves, typically $m=\texttt{neval}^2$ for full 2D perspective surfaces).

Thus errors="bootstrap", band="pmzsd" summarizes each bootstrap distribution by its standard deviation and imposes a symmetric normal-quantile shape on the displayed band. In contrast, bootstrap "pointwise" and "bonferroni" bands use the corresponding lower and upper empirical bootstrap quantiles separately and may be asymmetric at each evaluation point. Bootstrap "simultaneous" uses a joint rank construction rather than either pointwise standard deviations or pointwise tail quantiles. With band="all", these distinct bootstrap families are displayed together where supported.

For asymptotic errors, "pmzsd" and "pointwise" are two names for the same symmetric normal construction. Let $T(x)$ denote the plotted functional (mean, gradient, density, distribution, etc.) and $\widehat{se}(x)$ its asymptotic standard error. Both use $T(x)\pm z_{1-\alpha/2}\widehat{se}(x)$. "bonferroni" applies the same rule with $\alpha/m$ in place of $\alpha$.

Regression, single-index and smooth-coefficient mean and gradient curves support asymptotic Gaussian simultaneous variability bands within each panel's finite evaluation grid. This supports fixed, generalized nearest-neighbor and adaptive nearest-neighbor bandwidths with second-order Gaussian or Epanechnikov kernels and the existing LC, LL or LP specification. Numeric, unordered-factor and ordered-factor panels are supported, including mixed-data models and all-categorical smoothing under its existing local-constant specification. When there are no continuous smoothing variables, continuous-kernel settings do not restrict joint calibration. Regression, partially linear and smooth-coefficient mean surfaces support the same bands with both perspective and renderer="rgl" rendering. Gradients retain their slice plots; single-index plots retain their existing scalar-index geometry. Other continuous kernels retain a warning and unavailable simultaneous limits. With band="all", supported curves retain the pointwise and Bonferroni families and add the Gaussian simultaneous family.

For partially linear plots, a material influence reconstruction discrepancy (more than 1% of a standard error, with a roundoff allowance) signals numerical instability in the plotted fit. A warning identifies the affected points. Pointwise and Bonferroni limits are NA only at those points, and simultaneous limits are NA for the entire affected panel. The fitted curve is preserved; other panels retain their available bands.

Density and distribution plots (npudens, npudist, npcdens and npcdist) also support Gaussian simultaneous variability bands for ordinary panels and complete perspective/rgl level grids. Supported specifications include fixed, generalized nearest-neighbor and adaptive nearest-neighbor bandwidths with second-order Gaussian or Epanechnikov kernels. Conditional LC, LL and LP fits retain their selected degree and basis. Available conditional continuous derivatives and categorical first differences are supported in slice plots; higher generalized-NN derivatives retain their existing restrictions. Unconditional adaptive-NN uniform kernels are supported on regular spacing and support-face domains. Fixed continuous conditional LC specifications also support beta kernels. Proper-projected density/CDF plots and other continuous kernel specifications retain their documented simultaneous-band restrictions. Conditional quantile plots have a separate finite-grid influence construction; see npqreg for its supported level and gradient targets and conditioning assumptions.

These density-family bands use the incumbent marginal standard errors and the correlation of centered observation-level empirical influence rows. Writing $D$ for the diagonal matrix of those standard errors and $R$ for that correlation, the working covariance is $D R D$. The conditional ratio includes numerator-denominator dependence; LP uses the accepted local solve, including smooth ridge-scale and intercept compensation terms. Categorical contrasts subtract endpoint influence rows before computing covariance. Conditional NN geometry is held fixed, and unconditional generalized-NN uses its frozen-query-radius approximation; unconditional adaptive-NN includes the same-sample random-radius contribution. All three band families share the same marginal scale. Simultaneity is within the complete displayed panel/grid, not across separate panels.

The Gaussian-multiplier calibration is based on the framework developed by Victor Chernozhukov, Denis Chetverikov and Kengo Kato (2013, 2015). Related Gaussian-process and simultaneous-inference results for local regression distribution estimators are due to Matias D. Cattaneo, Michael Jansson and Xinwei Ma (2024); for local polynomial conditional density estimators, see Matias D. Cattaneo, Rajita Chandak, Michael Jansson and Xinwei Ma (2024). These works provide the related methodological foundation; the implementations here retain this package's kernel estimators, marginal standard errors and density-specific influence contracts described above. They do not implement the distinct estimators or robust bias-correction procedures in those papers. This is finite-grid variability conditional on the selected model and smoothing specification, not bias correction, bandwidth-search uncertainty or honest continuous-domain confidence coverage. Gaussian approximation can be poor with small nearest-neighbor neighborhoods or sparse evaluation points, including grids extending to sample extremes. Increasing the number of Gaussian draws reduces simulation error but does not correct these limitations. Nonregular ANN geometry, unavailable marginal uncertainty and resource limits produce a warning and NA simultaneous limits, preserving other available bands. Unexpected calculation errors and material target reconstruction failures remain errors.

The separately named plot argument gaussian.draws has a fixed default of 9999 and accepts positive integer counts subject to representability and workspace limits, without a fixed one-million upper cap. It applies only to asymptotic simultaneous/all requests; B remains bootstrap-only. The usual random.seed (default 42) controls calibration without changing the caller's random-number state. A count below ceiling(qnorm(.975)^2 * (1-alpha) / (.1^2 * alpha)) warns and is used unchanged. This is an approximate 95 percent Monte Carlo target for relative error within 10 percent in the simulated maximum's upper-tail probability. It is about 7299 draws at alpha=0.05 and 38031 at alpha=0.01. The maximum's tail probability is $\alpha$ for either a slice or a full neval^2 perspective/rgl grid; the number of points affects the maximum distribution and computational cost, but does not divide this tail probability. This recommendation does not guarantee critical-value accuracy or statistical coverage and is not a draw-count prescription in Chernozhukov, Chetverikov and Kato (2013). Explicit counts are never increased silently. Simulation uses bounded blocks and the bootstrap progress mechanism to report completed draws, elapsed time and estimated time remaining during interactive computation. It uses a conservative 64 MiB workspace budget; reduce neval or gaussian.draws if the workspace is too large. Regression, single-index and smooth-coefficient slice covariance calibration is limited to 512 grid points. Density-family plots switch to an observation-level factor above 512 points and budget up to 128 MiB for influence construction. This admission accounts for training size as well as grid size and LP workspaces: with the default 50 by 50 density-family surface, the factor/calibration allowance alone limits admission to about 1677 training observations; additional workspaces can lower that limit. Reducing neval reduces the surface grid quadratically. An over-budget request warns and retains the available marginal bands with NA simultaneous limits. Surfaces and partially linear plots instead use observation-level influence factors, avoiding a grid-by-grid covariance matrix and that grid-size limit.

Write $S$ for the training smoother and $L$ for its evaluation influence map. The covariance is $L\operatorname{diag}(z_i^2)L^\mathsf{T}$, where $z_i=((I-S)y)_i/\sqrt{\sum_j(\delta_{ij}-S_{ij})^2}$. For regression derivatives, replace $L$ by the corresponding derivative influence map; the training smoother and normalized residuals remain those of the fitted mean. Available higher-order continuous partials use the existing gradient.order contract. For categorical regression gradients, use the difference of the two endpoint influence rows. Unordered factors compare each level with the first retained level; ordered factors use adjacent-level differences, with the first forward difference repeated at the first level. This retains the covariance between endpoints, gives zero width for an unordered reference contrast, and calibrates jointly over all displayed levels. Simultaneity remains within each panel, not across separate panels. A simulated maximum of the absolute standardized joint Gaussian coordinates supplies one multiplier for the existing standard errors. Its simulated critical value is bounded by the analytical pointwise and Bonferroni critical values, ensuring that these three band widths remain ordered despite Monte Carlo variation. The raw simulated critical value is retained as critical.raw in the calibration metadata for nonzero-variance panels. For regression, single-index and smooth-coefficient plots, if the joint calculation cannot reproduce the plotted target or standard errors, has unidentified residual information, or exceeds its workspace/grid budget, the plot warns and returns unavailable (NA) simultaneous limits. With band="all", pointwise and Bonferroni intervals are retained. This can occur in sparse, locally ridged tails even after bandwidth selection. Invalid arguments, invalid covariance objects and unexpected computation errors still fail explicitly.

This Gaussian/multiplier connection follows the framework developed by Victor Chernozhukov, Denis Chetverikov and Kengo Kato (2013); see also their Gaussian comparison and anti-concentration results (2015). Conditional on the data and the fitted smoothing specification, simulating this covariance has the same Gaussian law as applying independent standard normal multipliers to the observation-level influence contributions. The covariance is assembled once per panel, so calibration requires no repeated estimator fits. These references concern the Gaussian approximation and multiplier framework; the covariance above retains this package's estimator-specific standard-error convention and does not by itself establish bias correction or honest confidence coverage.

Density-family output objects retain band.calibration and, for band="all", band.errors. With conditional gradients these are lists indexed by the explanatory-variable component. Calibration records the draw count, grid size, critical value and working dependence contract.

Single-index bands additionally condition on the fitted index coefficients; gradient panels multiply the scalar-index derivative by the corresponding fitted index coefficient and its widths by that coefficient's absolute value; negative coefficients exchange lower and upper widths for asymmetric bootstrap intervals before scaling. They retain the scalar-index regression standard-error convention. Smooth-coefficient bands describe the noniterated plotted mean or existing slope-coefficient gradient selected by coef.index, and use fitted raw training residuals in place of the normalized residuals above. These gradients hold the smoothing coordinates fixed, including when explanatory and smoothing variables overlap; they are not total derivatives through those coordinates. The separate coef=TRUE display remains excluded. A joint request reports unavailable simultaneous limits if its influence map cannot reproduce the incumbent plotted target and standard errors, including incompatible locally deficient solve configurations. These three families retain their existing pointwise intervals.

Partially linear mean plots use the complete response map of the fitted estimator, including coefficient estimation. Write $R=X-S_xX$ and $R_e=X_e-L_xX$, with each column using its own nuisance smoother, and let $K=R^+(I-S_y)$. The evaluation map is $H=L_y+R_eK$; the training map is $S=S_y+RK$. The plot covariance is $H\operatorname{diag}(z_i^2)H^\mathsf{T}$, using the complete training map in the normalized-residual formula above. Thus the shared-observation cross-covariances are retained. This corrects the previous partially linear plot helper, which added component variances separately; its pointwise and Bonferroni plot widths consequently change as well. Coefficient estimates and their reported covariance are unchanged. With raw residuals in place of normalized residuals, this response map gives exactly the conditional covariance of independent fixed-design wild-bootstrap multipliers. Partially linear curves and surfaces use influence-factor calibration. Their nuisance influence projections reuse the scalar fit's accepted local factorization or global inverse; the fitted estimator and its ridge policy are unchanged. Reconstruction differences below one percent of a plotted standard error are treated as negligible on the interval scale. Material discrepancies raise an error.

Nearest-neighbor radii are held at the fitted design for this conditional asymptotic calculation. Training residuals retain the estimator's training evaluation convention, which can differ from explicit evaluation at the same coordinates. A blocked response-map calculation preserves that distinction without retaining a full training smoother. This can add computation for NN and partially linear asymptotic bands; it does not repeat bandwidth selection or propagate its sampling variability.

Bandwidths and the polynomial specification are held fixed. This construction neither corrects smoothing bias nor guarantees finite-sample confidence coverage. It does not cover points between the grid or jointly cover different panels. Returned data include compact band.calibration metadata; band="all" also returns the three width matrices in band.errors. Surface widths follow the flattened evaluation-grid order. For a single-index gradient object, each width family is an evaluation-point by bound by coefficient array, aligned with the returned gradient columns. Factor-based calibration reports rank=NA because no eigendecomposition is computed.

On an eligible single-panel plot rendered by base graphics, a missing sub argument allows the plot method to add a concise variability subtitle after finite bounds have actually been drawn. This includes individual figures in a sequence with options(plot.par.mfrow=FALSE) when the active graphics layout is a full-frame single panel. Levels and supported gradients/effects follow the same rule. Categorical interval bars use interval rather than band terminology. The wording distinguishes symmetric $z$-times-standard-error bands from bootstrap quantile and rank constructions. For example, at the default alpha=0.05, labels include "95% symmetric variability band (normal quantile × asymptotic SE)" and "95% pointwise bootstrap variability band (quantiles)". Bias-corrected bootstrap labels also identify their center.

The automatic subtitle uses the effective subtitle typography (including cex.sub) and only appears when an approved one-line wording fits the existing plot width and bottom-margin clearance. It does not enlarge margins or reduce the plotting region; when no wording fits, it is silently omitted. Supplying sub, including sub="" or sub=NULL, suppresses the automatic label without changing that route's existing handling of the supplied subtitle.

Package-managed multi-panel base figures instead receive one shared annotation at the top of each page, using existing whitespace without changing margins, panel sizes, or fonts. Simultaneous and Bonferroni wording identifies the within-panel scope; for multiple quantile curves, the scope is within each quantile curve, not joint across quantiles. An annotation is omitted if its text does not fit or clearance from a main title cannot be established. Caller-managed multi-panel layouts and partial figure regions are not labelled with a figure-wide annotation. Bootstrap boxplots and coefficient displays of plus/minus one standard error are not labelled as nominal variability intervals. The global plot.par.mfrow option also governs automatic layout for class-probability plots.

On these eligible base plots, with band="all", default legend labels also identify the construction, for example "Pointwise (quantiles)" and "Simultaneous (rank-based)" for bootstrap bands, or "Pointwise (z × SE)" and "Bonferroni (adjusted z × SE)" for asymptotic bands. A default legend omits any family for which no finite band was drawn. Explicit user-supplied legend labels are preserved without relabelling or filtering.

These automatic annotations apply to the base renderer. The "rgl" renderer uses the same construction-aware default labels in its existing band="all" legend, but retains its sparse interactive presentation: it adds no automatic footer or subtitle, and panning, zooming, responsive sizing, and widget geometry are unchanged.

When center="bias-corrected" and output="data" (or output="both") is requested, returned plot objects that support bootstrap bias correction include bias and bias.corrected components where available. The sign convention is bias = estimate - bias.corrected, so the corrected estimate can be recovered as estimate - bias. For example, regression plot objects use mean - bias; density objects use dens - bias; and distribution objects use dist - bias. Regression-gradient plot objects use the analogous convention gbias = grad - gradient.bias.corrected, so the corrected gradient can be recovered as grad - gbias. Two-dimensional plots do not draw a second corrected surface in this tranche; use output="data" to recover the corrected surface values.

Kernel estimators are generally biased because smoothing trades variance reduction for bias. For this reason, pointwise and simultaneous error bounds centered on the fitted kernel estimate are often best interpreted as variability bounds rather than literal confidence intervals for the unknown target function. Confidence-interval interpretation requires the smoothing bias to be negligible relative to sampling variability, or to be corrected with sufficient accuracy.

When bootstrap errors are requested, center="bias-corrected" recenters the plotted bootstrap bounds on a method-specific Efron-Tibshirani bootstrap bias-corrected estimate, where available. This produces bounds that are closer in spirit to confidence intervals than variability bounds centered on the original biased estimate. Bias correction for nonparametric estimators is difficult, however, and no universal correction is optimal across estimators, bandwidths, supports, and data types.

For wild-bootstrap regression routes, center="bias-corrected" uses the standard wild-bootstrap correction, without a second oversmoothed pilot bandwidth. This follows the wild-bootstrap regression literature and the package's finite-sample evidence that an additional pilot can underadjust in this route. For supported fixed-bandwidth regression-style routes with bootstrap="inid", "fixed", or "geom", center="bias-corrected" estimates the bias center from a second pilot bootstrap matrix while retaining the ordinary bootstrap matrix for variability widths. These combinations can behave differently in finite samples, so the returned bias and bias.corrected components should be inspected when bias correction is scientifically important. For supported fixed regression gradient plots, inspect gbias, gradient.bias.corrected, or route-specific gradient-bias components such as gc1bias; different bootstrap schemes can give materially different derivative-bias corrections.

The pilot rule is family-specific. For fixed-bandwidth regression plots, the non-wild pilot uses the full oversmoothing-rate rule: $$g_j = h_j n^{2r / ((p + 2r)(p + 4r))},$$ where $h_j$ is the selected bandwidth, $r$ is the continuous kernel order, and $p$ is the number of continuous smoothing variables in that regression route. For ordinary second-order kernels with one continuous smoothing variable this reduces to $g_j = h_j n^{4/45}$. Categorical variables do not enter this continuous-rate exponent. Supported semiparametric regression and single-index plot routes use the square-root version of this rule.

For unconditional density and distribution plots with continuous variables, including mixed-data routes supported by the estimator, the smooth-bootstrap perturbation rule uses $g_j = h_j$, following the density smooth-bootstrap literature. Continuous variables are perturbed and categorical variables are resampled row-wise without perturbation. The smooth perturbation law is drawn from the selected ordinary continuous kernel when that kernel defines a probability law: second-order Gaussian, uniform, and second-order Epanechnikov kernels are supported. Epanechnikov perturbations use the package's standardized support $[-\sqrt{5},\sqrt{5}]$ before multiplication by the pilot bandwidth. Higher-order Gaussian and Epanechnikov kernels remain valid estimation kernels, but their signed equivalent-kernel forms do not define ordinary random perturbations, so center="bias-corrected" with pair/block/geometric bootstrap fails clearly for those higher-order smooth-bootstrap perturbation routes. For conditional density/distribution plots, the continuous dependent-variable side uses the same density-style rule $g_y = h_y$; categorical dependent variables are resampled row-wise without perturbation. The conditioning-variable side uses a full regression-style pilot: continuous conditioning bandwidths are enlarged, and categorical conditioning lambdas are moved toward their admissible upper bounds without crossing them. Categorical conditioning values are drawn independently given the selected donor, using the selected categorical kernel normalized over the declared possible levels for that donor. Factor levels and their ordering are retained. For conditional density this is $$g_{xj} = h_{xj} n^{2r / ((p + 2r)(p + 4r))},$$ and for conditional distribution this is $$g_{xj} = h_{xj} n^{r / ((p + r)(p + 2r))},$$ where here $p$ is the total number of continuous dependent and explanatory variables. The same conditional pilot rule is used for supported fixed-bandwidth conditional density/distribution gradient plots, including mixed-data routes. In mixed-data conditional gradients, continuous dependent variables are perturbed by the density-style rule, categorical dependent variables are resampled row-wise without perturbation, continuous conditioning variables use the conditional regression-style pilot, and categorical conditioning lambdas are moved toward their admissible upper bounds as in the corresponding level route and actually used in the categorical draws. The conditional pilot uses donor-truncated continuous probability laws when finite domain bounds are supplied, with the same normalization in its reference calculation.

Each conditional pilot resample is refitted at the original estimation bandwidths and the selected LC/LL/LP specification. The bias estimate subtracts the known conditional target of the generating pilot law from the mean pilot-bootstrap estimate, not the original fitted estimate. Conditional PDF/CDF levels, continuous derivatives and categorical first differences use this same reference contract; the endpoints of a categorical contrast share each resampled dataset. Ordinary bootstrap draws still determine variability widths. A pilot reference without support, or a derivative at a nonsmooth compact-kernel knot, fails explicitly rather than silently changing the requested centering method. These are practical, evidence-backed rules-of-thumb guided by the smoothed-bootstrap and regression two-bandwidth literatures, not claims of universal optimality.

In practice, users who want confidence-style bounds should prefer bootstrap errors with center="bias-corrected". The ordinary fitted curve or surface remains the default estimate because bias correction can increase variance and may have larger mean squared error than the uncorrected estimator.

Asymptotic standard errors are taken from fitted-object components such as merr, gerr, derr, conderr, and congerr where implemented. Ordinary regression, single-index mean/gradient and transformed-response least-squares quantile errors inherit the residual normalization described in npreg. This does not replace other families' covariance formulas.

Bootstrap resampling is conducted pairwise on $(y,X,Z)$ (i.e., by resampling rows of $(y,X)$ or $(y,X,Z)$ as appropriate). Bootstrap method support differs by estimator family:


  Regression-family (npreg/npindex/npscoef/npplreg): wild, inid, fixed, geom
  Density/distribution-family (npudens/npudist/npcdens/npcdist/npconmode): inid, fixed, geom

hence "wild" is only available for regression-family plotting. Bias-corrected bootstrap centering is currently restricted to supported fixed-bandwidth plot routes. Non-fixed, adaptive-nearest-neighbor, and generalized-nearest-neighbor bandwidth routes fail closed for center="bias-corrected". Conditional quantile npqreg bootstrap plots currently require center="estimate"; center="bias-corrected" fails clearly for this route. Copula npcopula and instrumental-variable npregiv bias-corrected bootstrap plot routes are not claimed in this release tranche and fail clearly when unsupported bias-corrected bootstrap controls are requested. Partially linear npplreg level plots support the regression-family bootstrap methods above; gradients=TRUE is not a partially linear plot route and fails clearly.

Implementation notes for speed:


  wild            : fast np*hat linear-operator bootstrap path
  inid/fixed/geom : fast direct helper path (no internal bandwidth search)

For non-fixed density/distribution bootstrap, an explicit experimental approximation is available via boot.control = np_boot_control(nonfixed = "exact") or boot.control = np_boot_control(nonfixed = "frozen"). The default "exact" route recomputes the non-fixed geometry for each resample. The experimental "frozen" route reuses the original-sample non-fixed geometry throughout the bootstrap run. This option is currently implemented only for unconditional and conditional density/distribution bootstrap routes and remains off by default. For generalized/adaptive nearest-neighbor runs, "frozen" is an approximation that can alter interval/band width by holding the original-sample nearest-neighbor geometry fixed across bootstrap resamples; "exact" remains the recommended setting for production inference. This approximation can be more noticeable for conditional density/distribution plotting than for the regression-style plot families because the conditional bootstrap paths freeze both numerator and denominator nearest-neighbor geometry before recombining them. In practice, conditional distribution bands are often closer, while conditional density bands can differ more materially from "exact" under generalized/adaptive nearest-neighbor bandwidths. For smooth coefficient plots (npscoef) under non-fixed bandwidths, "exact" can also be much more expensive than "frozen" on large jobs, because the coefficient field must be recomputed for each bootstrap resample rather than reusing the original-sample geometry. This recomputation cannot in general be avoided without a more aggressive approximation: for npscoef the local weighted systems that define the coefficient vector depend on the bootstrap resample weights/counts at each evaluation point, so unlike npplreg there is no single global coefficient vector that can be updated once per draw. inid admits general heteroskedasticity of unknown form, though it does not allow for dependence. fixed conducts Kunsch's (1988) block bootstrap for dependent data, while geom conducts Politis and Romano's (1994) stationary bootstrap.

For local polynomial conditional density/distribution plotting (npcdens/npcdist with regtype="ll" or regtype="lp") and proper=TRUE, the plotted estimate is rendered proper slice-by-slice on the fixed evaluation grid: each conditional density slice is projected to be nonnegative and to integrate to one using trapezoidal quadrature weights from the evaluation $y$-grid, while each conditional distribution slice is projected to be monotone and bounded in $[0,1]$. When errors="bootstrap", the bootstrap resample surfaces are computed first on that same fixed grid and then properized resample-by-resample using the same grid geometry before "pointwise", "bonferroni", "simultaneous", and "all" bands are constructed. Thus the bootstrap distribution used to form these bands is built from properized resample surfaces. The final lower/upper band surfaces are interval envelopes and are not themselves separately re-projected to satisfy the density/distribution shape constraints.

For consistency of the block and stationary bootstrap, the (mean) block length $b$ should grow with the sample size $n$ at an appropriate rate. If $b$ is not given, then a default growth rate of $const \times n^{1/3}$ is used. This rate is “optimal” under certain conditions (see Politis and Romano (1994) for more details). However, in general the growth rate depends on the specific properties of the DGP. A default value for $const$ ($3.15$) has been determined by a Monte Carlo simulation using a Gaussian AR(1) process (AR(1)-parameter of 0.5, 500 observations). $const$ has been chosen such that the mean square error for the bootstrap estimate of the variance of the empirical mean is minimized.

The default bootstrap replication count is B=9999, including copula plots. Class-probability plots retain their default B=399. Estimator and test bootstrap defaults are separate. For pointwise percentile intervals the bare tail-resolution reference is $B \ge \lceil 2/\alpha - 1 \rceil$. For "bonferroni", "simultaneous", and "all", the warning uses $B_{\min}=\lceil 2m/\alpha-1 \rceil$, where $m$ is the actual within-panel evaluation count, including $m=\texttt{neval}^2$ on a full perspective grid. For rank-based simultaneous intervals this is a conservative Bonferroni reference. These are bare rank-scale floors, not integrality requirements, Monte Carlo precision targets or coverage guarantees. Meeting a floor may still leave noisy endpoints; larger B can stabilize them. At alpha=0.05, a 50-point curve has a floor of 1999, while a 50 by 50 surface has a floor of 99999. Below the floor the plot warns and uses the requested count unchanged. Bootstrap standard-deviation bands ("pmzsd") do not estimate empirical tail quantiles and are not subject to this tail-resolution warning.

Typical plotting calls:


  ## Asymptotic pointwise/bonferroni intervals
  plot(fit, errors="asymptotic", band="pointwise")
  plot(fit, errors="asymptotic", band="bonferroni")

  ## Regression-family bootstrap (wild available)
  plot(fit, errors="bootstrap", bootstrap="wild")

  ## Density/distribution-family bootstrap (use inid/fixed/geom)
  plot(fit, errors="bootstrap", bootstrap="inid")

Compound public plot argument names use dots, consistent with the package interface. Retired underscore spellings are not accepted. Use errors, bootstrap, B, band, alpha, center, output, data.overlay, data.rug, layout, and the control constructors documented in np_boot_control, np_grid_control, and np_render_control.

Value

Setting output will instruct plot what data to return. Option summary:
plot: instruct plot to just plot the data and return NULL
both: instruct plot to plot the data and return the data used to generate the plots. The data will be a list of objects of the appropriate type, with one object per plot. For example, invoking plot on 3D density data will have it return a list of three npdensity objects. If biases were calculated, they are stored in a component named bias
data: instruct plot to generate data only and no plots

References

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

Cao-Abad, R. (1991), “Rate of convergence for the wild bootstrap in nonparametric regression,” The Annals of Statistics, 19, 2226-2231.

Cattaneo, Matias D., Rajita Chandak, Michael Jansson and Xinwei Ma (2024), “Boundary adaptive local polynomial conditional density estimators,” Bernoulli, 30(4), 3193-3223. doi:10.3150/23-BEJ1711 .

Cattaneo, Matias D., Michael Jansson and Xinwei Ma (2024), “Local regression distribution estimators,” Journal of Econometrics, 240(2), 105074. doi:10.1016/j.jeconom.2021.01.006 .

Chernozhukov, Victor, Denis Chetverikov and Kengo Kato (2013), “Gaussian approximations and multiplier bootstrap for maxima of sums of high-dimensional random vectors,” The Annals of Statistics, 41(6), 2786-2819. doi:10.1214/13-AOS1161 .

Chernozhukov, Victor, Denis Chetverikov and Kengo Kato (2015), “Comparison and anti-concentration bounds for maxima of Gaussian random vectors,” Probability Theory and Related Fields, 162, 47-70. doi:10.1007/s00440-014-0565-9 .

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.

Hall, P. (1992a), “Effect of bias estimation on coverage accuracy of bootstrap confidence intervals for a probability density,” The Annals of Statistics, 20, 675-694.

Hall, P. (1992b), “On bootstrap confidence intervals in nonparametric regression,” The Annals of Statistics, 20, 695-711.

Hardle, W. and J.S. Marron (1991), “Bootstrap simultaneous error bars for nonparametric regression,” The Annals of Statistics, 19, 778-796.

Kunsch, H.R. (1989), “The jackknife and the bootstrap for general stationary observations,” The Annals of Statistics, 17, 1217-1241.

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

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

Politis, D.N. and J.P. Romano (1994), “The stationary bootstrap,” Journal of the American Statistical Association, 89, 1303-1313.

Silverman, B.W. and G.A. Young (1987), “The bootstrap: to smooth or not to smooth?” Biometrika, 74, 469-479.

De Angelis, D. and G.A. Young (1992), “Smoothing the bootstrap,” International Statistical Review, 60, 45-56.

Scott, D.W. (1992), Multivariate Density Estimation. Theory, Practice and Visualization, New York: Wiley.

Silverman, B.W. (1986), Density Estimation, London: Chapman and Hall.

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.

Examples

if (FALSE) { # \dontrun{
# EXAMPLE 1: For this example, we load Giovanni Baiocchi's Italian GDP
# panel (see Italy for details), then create a data frame in which year
# is an ordered factor, GDP is continuous, compute bandwidths using
# likelihood cross-validation, then create a grid of data on which the
# density will be evaluated for plotting purposes

data("Italy")
with(Italy, {

data <- data.frame(ordered(year), gdp)

# Compute bandwidths using likelihood cross-validation (default). Note
# that this may take a minute or two depending on the speed of your
# computer...

bw <- npudensbw(dat=data)
fit <- npudens(tdat=data, bws=bw)

# You can always do things manually, as the following example demonstrates

# Create an evaluation data matrix

year.seq <- sort(unique(year))
gdp.seq <- seq(1,36,length=50)
data.eval <- expand.grid(year=year.seq,gdp=gdp.seq)

# Generate the estimated density computed for the evaluation data

fhat <- fitted(npudens(tdat = data, edat = data.eval, bws=bw))

# Coerce the data into a matrix for plotting with persp()

f <- matrix(fhat, length(unique(year)), 50)

# Next, create a 3D perspective plot of the PDF f

persp(as.integer(levels(year.seq)), gdp.seq, f, col="lightblue",
      ticktype="detailed", ylab="GDP", xlab="Year", zlab="Density",
      theta=300, phi=50)

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

Sys.sleep(5)

# However, plot simply streamlines this process and aids in the
# visualization process (<ctrl>-C will interrupt on *NIX systems, <esc>
# will interrupt on MS Windows systems).

plot(fit)

# plot also streamlines construction of variability bounds (<ctrl>-C
# will interrupt on *NIX systems, <esc> will interrupt on MS Windows
# systems)

plot(fit, errors = "asymptotic")

# EXAMPLE 2: For this example, we simulate multivariate data, and plot the
# partial regression surfaces for a locally linear estimator and its
# derivatives.

set.seed(123)

n <- 100

x1 <- runif(n)
x2 <- runif(n)
x3 <- runif(n)
x4 <- rbinom(n, 2, .3)

y <- 1 + x1 + x2 + x3 + x4 + rnorm(n)

X <- data.frame(x1, x2, x3, ordered(x4))

bw <- npregbw(xdat=X, ydat=y, regtype="ll", bwmethod="cv.aic")
fit <- npreg(bws=bw)

plot(fit)

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

Sys.sleep(5)

# Now plot the gradients...

plot(fit, gradients=TRUE)

# Plot the partial regression surfaces with bias-corrected bootstrapped
# nonparametric confidence intervals... this may take a minute or two
# depending on the speed of your computer as the bootstrapping must be
# completed prior to results being displayed...

plot(fit,
     errors="bootstrap",
     center="bias-corrected",
     band="simultaneous")

# EXAMPLE 3: This example demonstrates how to retrieve plotting data from
# plot(). When plot() is called with the arguments
# `output="both"' (or "data"), it returns plotting objects
# named r1, r2, and so on (rg1, rg2, and so on when `gradients=TRUE' is
# set).  Each plotting object's index (1,2,...) corresponds to the index
# of the explanatory data data frame xdat (and zdat if appropriate).

# Take the cps71 data by way of example. In this case, there is only one
# object returned by default, `r1', since xdat is univariate.

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

# Compute bandwidths for local linear regression using cv.aic...

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

# Generate the plot and return plotting data, and store output in
# `plot.out' (NOTE: the call to `output' is necessary).

plot.out <- plot(fit,
                 perspective=FALSE,
                 errors="bootstrap",
                 B=25,
                 output="both")

# Now grab the r1 object that plot plotted on the screen, and take
# what you need. First, take the output, lower error bound and upper
# error bound...

logwage.eval <- fitted(plot.out$r1)
logwage.se <- se(plot.out$r1)
logwage.lower.ci <- logwage.eval + logwage.se[,1]
logwage.upper.ci <- logwage.eval + logwage.se[,2]

# Next grab the x data evaluation data. xdat is a data.frame(), so we
# need to coerce it into a vector (take the `first column' of data frame
# even though there is only one column)

age.eval <- plot.out$r1$eval[,1]

# Now we could plot this if we wished, or direct it to whatever end use
# we envisioned. We plot the results using R's plot() routines...

with(cps71, plot(age, logwage, cex=0.2, xlab="Age", ylab="log(Wage)"))
lines(age.eval,logwage.eval)
lines(age.eval,logwage.lower.ci,lty=3)
lines(age.eval,logwage.upper.ci,lty=3)

# Bias-corrected bootstrap output includes the estimated bias and the
# bias-corrected center when available. The sign convention is
# bias = estimate - bias.corrected, so bias.corrected = estimate - bias.

plot.out <- plot(fit,
                 perspective=FALSE,
                 errors="bootstrap",
                 center="bias-corrected",
                 B=25,
                 output="data")

logwage.eval <- fitted(plot.out$r1)
logwage.bias <- plot.out$r1$bias
logwage.bc <- plot.out$r1$bias.corrected
stopifnot(all.equal(logwage.bc, logwage.eval - logwage.bias))

# If you wanted plot() data for gradients, you would use the argument
# `gradients=TRUE' in the call to plot() as the following
# demonstrates...

plot.out <- plot(fit,
                 perspective=FALSE,
                 errors="bootstrap",
                 B=25,
                 output="both",
                 gradients=TRUE)

# Now grab object that plot() plotted on the screen. First, take the
# output, lower error bound and upper error bound... note that gradients
# are stored in objects rg1, rg2 etc.

grad.eval <- gradients(plot.out$rg1)
grad.se <- gradients(plot.out$rg1, se = TRUE)
grad.lower.ci <- grad.eval + grad.se[,1]
grad.upper.ci <- grad.eval + grad.se[,2]

# Next grab the x evaluation data. xdat is a data.frame(), so we need to
# coerce it into a vector (take `first column' of data frame even though
# there is only one column)

age.eval <- plot.out$rg1$eval[,1]

# We plot the results using R's plot() routines...

plot(age.eval,grad.eval,cex=0.2,
     ylim=c(min(grad.lower.ci),max(grad.upper.ci)),
     xlab="Age",ylab="d log(Wage)/d Age",type="l")
lines(age.eval,grad.lower.ci,lty=3)
lines(age.eval,grad.upper.ci,lty=3)

# EXAMPLE 4: Variations on local polynomial conditional density
# estimation with proper = TRUE.

data("Italy")

Italy2 <- within(Italy, {
  year <- as.numeric(as.character(year))
})

# Plot only: make the plotted surface proper on the plot evaluation grid.

fhat <- npcdens(gdp ~ year, data = Italy2,
                regtype = "lp", degree = 3, nmulti = 1)

plot(fhat, proper = TRUE)

# Fit an object whose fitted values are themselves proper.

ctrl_fit <- list(
  mode = "slice",
  apply = "fitted",
  slice.grid.size = 101L,
  slice.extend.factor = 0.1
)

fhat_fit <- npcdens(
  gdp ~ year,
  data = Italy2,
  regtype = "lp",
  degree = 3,
  nmulti = 1,
  proper = TRUE,
  proper.control = ctrl_fit
)

fit_proper <- fitted(fhat_fit)
fit_raw <- fhat_fit$condens.raw

# Display the repaired and raw fitted values for cases where the raw
# fitted density is negative.

head(cbind(fit_proper, fit_raw)[which(fit_raw < 0), ])

# Predict on a common explicit y-grid for several years, and render
# those predictions proper.

g.grid <- seq(min(Italy2$gdp), max(Italy2$gdp), length.out = 200)

nd_grid <- expand.grid(
  gdp = g.grid,
  year = c(1955, 1975, 1995)
)

pred_grid <- predict(fhat, newdata = nd_grid, proper = TRUE)

# Predict on paired rows with different gdp grids by year, and still
# make the predictions proper via slice mode.

g1 <- seq(quantile(Italy2$gdp, 0.10),
          quantile(Italy2$gdp, 0.60), length.out = 60)
g2 <- seq(quantile(Italy2$gdp, 0.30),
          quantile(Italy2$gdp, 0.90), length.out = 35)

nd_slice <- rbind(
  data.frame(gdp = g1, year = rep(1960, length(g1))),
  data.frame(gdp = g2, year = rep(1985, length(g2)))
)

pred_slice <- predict(
  fhat,
  newdata = nd_slice,
  proper = TRUE,
  proper.control = list(mode = "slice")
)

# One object that carries properization for fitted values and for later
# predict() calls.

ctrl_both <- list(
  mode = "slice",
  apply = "both",
  slice.grid.size = 101L,
  slice.extend.factor = 0.1
)

fhat_both <- npcdens(
  gdp ~ year,
  data = Italy2,
  regtype = "lp",
  degree = 3,
  nmulti = 1,
  proper = TRUE,
  proper.control = ctrl_both
)

fit_both <- fitted(fhat_both)
pred_both <- predict(
  fhat_both,
  newdata = nd_slice,
  proper.control = ctrl_both
)
plot(fhat_both)
})
} # }