Regression with rbbnp

This article is a guide to conditional-mean kernel regression with bias_bound_regression().

Estimating a conditional expectation

bias_bound_regression() estimates \(E[Y \mid X=x]\) with bias-aware pointwise confidence intervals:

bias_bound_regression(
  x,                    # Predictor
  y,                    # Response
  eval = NULL,          # Evaluation points
  bw = "cv",            # Positive number, "cv", or "silverman"
  conf_level = 0.95,
  kernel = "schennach",
  control = NULL        # Named list of advanced settings
)

A basic fit, on data with a quadratic mean and Gaussian noise:

# Generate data: y = f(x) + noise
set.seed(42)
x_data <- runif(250) + runif(250)
y_data <- 2 * x_data - x_data^2 + rnorm(250, sd = 0.3)

# Estimate conditional expectation E[Y|X]
fit <- bias_bound_regression(
  x_data,
  y_data,
  bw = 0.1,
  kernel = "schennach"
)

# View results
fit
#> Bias-Bounded Conditional Expectation Estimation
#> 
#> Call:
#> bias_bound_regression(x = x_data, y = y_data, bw = 0.1, kernel = "schennach")
#> 
#> Sample size: n = 250
#> Bandwidth:   h = 0.1000 (user-specified) 
#> Kernel:      Schennach2004 
#> 
#> Bias bound parameters:
#>   A = 5.2202, r = 2.0000, B = 0.8378
#>   bias bounds: b1x = 0.1418, byx = 0.1709
#> 
#> Evaluation points: 100 (range: [0.0912, 1.9622])
#> Fitted values: E[Y|X] range [-0.1503, 0.9822]
#> Confidence level: 95%
#> 
#> Use summary() for detailed statistics
#> Use plot() to visualize results
#> Use fitted() to extract fitted values

The fit is an S3 object of class bbnp_regression. Use coef() for the key parameters and summary() for a fuller report:

# Key parameters
coef(fit)
#>         A         r         B         h 
#> 5.2201734 2.0000000 0.8377779 0.1000000

# Detailed summary
summary(fit)
#> Summary: Bias-Bounded Conditional Expectation Estimation
#> ============================================================
#> 
#> Call:
#> bias_bound_regression(x = x_data, y = y_data, bw = 0.1, kernel = "schennach")
#> 
#> Sample Information:
#>   Sample size (n):  250
#>   Bandwidth (h):    0.1000
#>   Kernel function:  Schennach2004
#> 
#> Bias Bound Parameters:
#>   A (amplitude):    5.2202
#>   r (decay rate):   2.0000
#>   B (Y bound):      0.8378
#>   b1x (bias f(x)):  0.1418
#>   byx (bias f_YX):  0.1709
#>   Xi interval:      [2.5062, 5.5672]
#> 
#> Range Estimation:
#>   Fitted values (E[Y|X]):
#>     min  Q1.25%  median    mean  Q3.75%     max 
#> -0.1503  0.4949  0.7871  0.6883  0.9281  0.9822 
#> 
#>   Marginal density f(x):
#>    min   mean    max 
#> 0.0474 0.5239 1.0028 
#> 
#>   Standard errors:
#>    min   mean    max 
#> 0.0268 0.0518 0.1654

Visualizing the fit

The default plot shows the estimated mean, its pointwise confidence-interval ribbon, and the original data:

plot(fit)

Element	Description
Estimate (line)	Estimated \(\hat{E}[Y \mid X=x]\)
95% CI (band)	Confidence interval
Points	Original data \((X_i, Y_i)\)

Because the plot returns a ggplot object, you can overlay the true mean for comparison:

# True function for comparison
true_fn <- function(x) 2 * x - x^2

plot(fit) +
  stat_function(fun = true_fn, aes(color = "True E[Y|X]"),
                linetype = "dashed", linewidth = 1) +
  scale_color_manual(values = c("Estimate" = "#08306B", "True E[Y|X]" = "red")) +
  labs(color = NULL) +
  theme(legend.position = "top")

Fitted values and confidence intervals

fitted() returns the estimated mean at the evaluation points:

# Get fitted values at evaluation points
predictions <- fitted(fit)
head(predictions, 10)
#>  [1] 0.2841416 0.3028584 0.3206543 0.3379210 0.3549087 0.3718462 0.3889071
#>  [8] 0.4062110 0.4238766 0.4419718

# Evaluation points
x_points <- fit$x
head(x_points)
#> [1] 0.0911572 0.1100561 0.1289551 0.1478540 0.1667529 0.1856519

In regions where the estimated marginal density \(\hat f(x)\) is very close to zero, the confidence interval may become unbounded (it can contain -Inf or Inf). This happens because the conditional mean estimator is a ratio involving \(1/\hat f(x)\). The interval is finite at interior points, where the density is well away from zero:

# Extract confidence intervals
ci <- confint(fit)

# Show the interval at interior points, where the marginal density is well away
# from zero so the ratio-based interval is finite
ok <- which(is.finite(ci[, "lower"]) & is.finite(ci[, "upper"]))
data.frame(
  x = fit$x[ok],
  estimate = fitted(fit)[ok],
  lower = ci[ok, "lower"],
  upper = ci[ok, "upper"]
)[1:6, ]
#>           x  estimate       lower     upper
#> 1 0.1100561 0.3028584 -16.5471610 28.347218
#> 2 0.1289551 0.3206543  -4.6771700  8.826715
#> 3 0.1478540 0.3379210  -2.5923395  5.435265
#> 4 0.1667529 0.3549087  -1.7184410  4.040969
#> 5 0.1856519 0.3718462  -1.2338997  3.290189
#> 6 0.2045508 0.3889071  -0.9227982  2.827247

Choosing the bandwidth

Set bw to "cv" for cross-validation or "silverman" for the faster rule-of-thumb selector. Regression cross-validation minimizes leave-one-out squared Nadaraya–Watson prediction error, so it uses both x and y (it is not density cross-validation applied only to x):

fit_cv <- bias_bound_regression(x_data, y_data, bw = "cv")
cat("CV bandwidth:", coef(fit_cv)["h"])
#> CV bandwidth: 0.1402307

Silverman’s predictor-spread rule is a faster alternative. It uses the same kernel-specific rescaling documented in the density vignette and, unlike regression CV, does not use y:

fit_silv <- bias_bound_regression(x_data, y_data, bw = "silverman")
cat("Silverman bandwidth:", coef(fit_silv)["h"])
#> Silverman bandwidth: 0.1428276

Worked examples

Built-in sample data

# Load the reproducible example
data(two_fold_uniform_example)
head(two_fold_uniform_example)
#>           x        y
#> 1 1.0859914 2.828175
#> 2 1.6292714 2.014698
#> 3 1.2303920 3.350186
#> 4 1.0686912 1.953532
#> 5 0.8441481 1.391845
#> 6 0.8789330 1.954778

# Estimate regression
fit_real <- bias_bound_regression(
  x = two_fold_uniform_example$x,
  y = two_fold_uniform_example$y,
  bw = 0.1,
  kernel = "schennach"
)

# Visualize
plot(fit_real) + ggtitle("Regression on Sample Data")

Comparing kernels

fit_sch <- bias_bound_regression(
  x_data, y_data, bw = 0.1, kernel = "schennach"
)
fit_norm <- bias_bound_regression(
  x_data, y_data, bw = 0.1, kernel = "normal"
)

grid.arrange(
  plot(fit_sch) + ggtitle("Schennach2004 Kernel"),
  plot(fit_norm) + ggtitle("Normal Kernel"),
  ncol = 1
)

Heteroscedastic errors

The method handles errors whose variance changes with X:

# Generate heteroscedastic data
set.seed(123)
x_het <- runif(250) + runif(250)
y_het <- x_het^2 + rnorm(250) * x_het  # Variance increases with x

fit_het <- bias_bound_regression(x_het, y_het, bw = 0.1)
plot(fit_het) + ggtitle("Heteroscedastic Data")

A polynomial mean

# Quadratic relationship
set.seed(456)
x_poly <- runif(250) + runif(250)
y_poly <- -x_poly^2 + 3 * x_poly + rnorm(250, sd = 0.5)

fit_poly <- bias_bound_regression(x_poly, y_poly, bw = 0.1)

# Compare with true function
true_poly <- function(x) -x^2 + 3*x

plot(fit_poly) +
  stat_function(fun = true_poly, aes(color = "True: -x^2 + 3x"),
                linetype = "dashed", linewidth = 1) +
  scale_color_manual(values = c("Estimate" = "#08306B", "True: -x^2 + 3x" = "red")) +
  labs(color = NULL) +
  theme(legend.position = "top")

S3 methods

# All available methods for bbnp_regression objects
print(fit)       # Concise summary
summary(fit)     # Detailed statistics
plot(fit)        # Visualization
coef(fit)        # Parameters (A, r, B, h)
confint(fit)     # Confidence intervals
fitted(fit)      # Fitted values