conreg {cobs}

Convex / Concave Regression

Description

Compute a univariate concave or convex regression, i.e., for given vectors, x,y,w in R^n, where x has to be strictly sorted (x_1 < x_2 < \ldots < x_n), compute an n-vector m minimizing the weighted sum of squares \sum_{i=1}^n {w_i (y_i - m_i)^2} under the constraints

(m_i - m_{i-1})/(x_i - x_{i-1}) \ge (m_{i+1} - m_i)/(x_{i+1} - x_i),

for 1 \le i \le n and m_0 := m_{n+1} := -\infty, for concavity. For convexity (convex=TRUE), replace \ge by \le and -\infty by +\infty.

Usage

conreg(x, y = NULL, w = NULL, convex = FALSE,
       method = c("Duembgen06_R", "SR"),
       tol = c(1e-10, 1e-7), maxit = c(500, 20),
       adjTol = TRUE, verbose = FALSE)

Arguments

x, y

numeric vectors giving the values of the predictor and response variable. Alternatively a single “plotting” structure (two-column matrix / y-values only / list, etc) can be specified: see xy.coords.

w

for method "Duembgen06_R" only: optional vector of weights of the same length as x; defaults to all 1.

convex

logical indicating if convex or concave regression is desired.

method

a character string indicating the method used,

"Duembgen06_R"

is an active set method written by Lutz Duembgen (University of Berne, CH) in Matlab in 2006 and translated to R by Martin Maechler.

"SR"

is an R interface to the C code of a Support Reduction algorithm written by Piet Groeneboom (TU Delft, NL) and donated to the cobs package in July 2012.

tol

convergence tolerance(s); do not make this too small!

maxit

maximal number of (outer and inner) iterations of knot selection.

adjTol

(for "Duembgen06_R" only:) logical indicating if the convergence test tolerance is to be adjusted (increased) in some cases.

verbose

logical or integer indicating if (and how much) knot placement and fitting iterations should be “reported”.

Details

Both algorithms need some numerical tolerances because of rounding errors in computation of finite difference ratios. The active-set "Duembgen06_R" method notably has two different such tolerances which were both 1e-7= 10^{7} up to March 2016.

The two default tolerances (and the exact convergence checks) may change in the future, possibly to become more adaptive.

Value

an object of class conreg which is basically a list with components

x

sorted (and possibly aggregated) abscissa values x.

y

corresponding y values.

w

corresponding weights, only for "Duembgen06_R".

yf

corresponding fitted values.

convex

logical indicating if a convex or a concave fit has been computed.

iKnots

integer vector giving indices of the knots, i.e. locations where the fitted curve has kinks. Formally, these are exactly those indices where the constraint is fulfilled strictly, i.e., those i where

(m_i - m_{i-1})/(x_i-x_{i-1}) > (m_{i+1} - m_i)/(x_{i+1}-x_i).

call

the call to conreg() used.

iter

integer (vector of length one or two) with the number of iterations used (in the outer and inner loop for "Duembgen06_R").

Note that there are several methods defined for conreg objects, see predict.conreg or methods(class = "conreg").

Notably print and plot; also predict, residuals, fitted, knots.

Also, interpSplineCon() to construct a more smooth (cubic) spline, and isIsplineCon() which checks if the int is strictly concave or convex the same as the conreg() result from which it was constructed.

Author(s)

Lutz Duembgen programmed the original Matlab code in July 2006; Martin Maechler ported it to R, tested, catch infinite loops, added more options, improved tolerance, etc; from April 27, 2007.

See Also

isoreg for isotone (monotone) regression; CRAN packages ftnonpar, cobs, logcondens.

Examples

## Generated data :
N <- 100
f <- function(X) 4*X*(1 - X)

xx <- seq(0,1, length=501)# for plotting true f()
set.seed(1)# -> conreg does not give convex cubic

x <- sort(runif(N))
y <- f(x) + 0.2 * rnorm(N)
plot(x,y, cex = 0.6)
lines(xx, f(xx), col = "blue", lty=2)
rc <- conreg(x,y)
lines(rc, col = 2, force.iSpl = TRUE)
 # 'force.iSpl': force the drawing of the cubic spline through the kinks
title("Concave Regression in R")

y2 <- y


## Trivial cases work too:
(r.1 <- conreg(1,7))
(r.2 <- conreg(1:2,7:6))
(r.3  <- conreg(1:3,c(4:5,1)))