R/monotonic.R

In robjhyndman/demography: Forecasting Mortality, Fertility, Migration and Population Data

#' Monotonic interpolating splines
#'
#' Perform cubic spline monotonic interpolation of given data points, returning either a list of points obtained by the interpolation or a function performing the interpolation. The splines are constrained to be monotonically increasing (i.e., the slope is never negative).
#'
#' These are simply wrappers to the \code{\link[stats]{splinefun}} function family from the stats package.
#'
#' @param x,y vectors giving the coordinates of the points to be interpolated. Alternatively a single plotting structure can be specified: see \code{\link[grDevices]{xy.coords}}.
#' @param n interpolation takes place at n equally spaced points spanning the interval [\code{xmin}, \code{xmax}].
#' @param xmin left-hand endpoint of the interpolation interval.
#' @param xmax right-hand endpoint of the interpolation interval.
#' @param ... Other arguments are ignored.
#'
#' @return \item{cm.spline}{returns a list containing components \code{x} and \code{y} which give the ordinates where interpolation took place and the interpolated values.}
#' \item{cm.splinefun}{returns a function which will perform cubic spline interpolation of the given data points. This is often more useful than \code{spline}.}
#'
#' @references Forsythe, G. E., Malcolm, M. A. and Moler, C. B. (1977) \emph{Computer Methods for Mathematical Computations}.
#' Hyman (1983) \emph{SIAM J. Sci. Stat. Comput.} \bold{4}(4):645-654.
#' Dougherty, Edelman and Hyman 1989 \emph{Mathematics of Computation}, \bold{52}: 471-494.
#'
#' @author Rob J Hyndman
#'
#' @examples
#' x <- seq(0, 4, l = 20)
#' y <- sort(rnorm(20))
#' plot(x, y)
#' lines(spline(x, y, n = 201), col = 2) # Not necessarily monotonic
#' lines(cm.spline(x, y, n = 201), col = 3) # Monotonic
#' @keywords smooth
#' @aliases monotonic
#' @export
cm.spline <- function(x, y = NULL, n = 3 * length(x), xmin = min(x), xmax = max(x), ...)
                      # wrapper for spline()
# Function retained for backwards compatibility
{
  stats::spline(x, y, n = n, xmin = xmin, xmax = xmax, method = "hyman")
}


#' @rdname cm.spline
#' @export
cm.splinefun <- function(x, y = NULL, ...)
                         # wrapper for splinefun()
# Function retained for backwards compatibility
{
  stats::splinefun(x, y, method = "hyman")
}


# Function to do cubic smoothing spline fit to y ~ x
# with constraint of monotonic increasing for x>b.
# Based on code provided by Simon Wood
# Last updated: 1 February 2014

smooth.monotonic <- function(x, y, b, k = -1, w = NULL, newx = x) {
  weight <- !is.null(w)
  if (k < 3 & k != -1) {
    stop("Inappropriate value of k")
  }
  # Unconstrained smooth.
  miss <- is.na(y)
  if (weight) {
    miss <- miss | w < 1e-9
  }
  yy <- y[!miss]
  xx <- x[!miss]
  if (weight) {
    w <- w[!miss]
    w <- w / sum(w) * length(w)
    f.ug <- mgcv::gam(yy ~ s(xx, k = k), weights = w)
    #        assign("w",w,pos=1)
  } else {
    f.ug <- mgcv::gam(yy ~ s(xx, k = k))
  }

  if (max(xx) <= b) {
    return(mgcv::predict.gam(f.ug, newdata = data.frame(xx = newx), se.fit = TRUE))
  }

  # Create Design matrix, constraints etc. for monotonic spline....
  mgcv::gam(yy ~ s(xx, k = k), data = data.frame(xx = xx, yy = yy), fit = FALSE) -> G
  if (weight) {
    G$w <- w
  }
  nc <- 200 # number of constraints
  xc <- seq(b, max(xx), l = nc + 1) # points at which to impose constraints
  A0 <- mgcv::predict.gam(f.ug, data.frame(xx = xc), type = "lpmatrix")
  # A0%*%p will evaluate spline at the xc points
  A1 <- mgcv::predict.gam(f.ug, data.frame(xx = xc + 1e-6), type = "lpmatrix")
  A <- (A1 - A0) / 1e-6 # approximate constraint matrix
  # (A%%p is  -ve gradient of spline at points xc)

  G$Ain <- A # constraint matrix
  G$bin <- rep(0, nc + 1) # constraint vector
  G$sp <- f.ug$sp # use smoothing parameters from un-constrained fit
  k <- G$smooth[[1]]$df + 1
  G$p <- rep(0, k)
  G$p[k] <- 0.1 # get monotonic starting parameters, by
  # setting coefficiants of polynomial part of term
  G$p[k - 1] <- -mean(0.1 * xx) # must ensure that gam side conditions are
  # met so that sum of smooth over x's is zero
  #    G$p <- rep(0,k+1)
  #    G$p[k+1] <- 0.1
  #    G$p[k] <- -mean(0.1*xx)
  G$y <- yy
  G$off <- G$off - 1 # indexing inconsistency between pcls and internal gam
  G$C <- matrix(0, 0, 0) # fixed constraint matrix (there are none)
  p <- mgcv::pcls(G) # fit spline (using s.p. from unconstrained fit)

  # now modify the gam object from unconstrained fit a little, to use it
  # for predicting and plotting constrained fit.
  f.ug$coefficients <- p
  return(mgcv::predict.gam(f.ug, newdata = data.frame(xx = newx), se.fit = TRUE))
}

smooth.monotonic.cobs <- function(x, y, b, lambda = 0, w = NULL, newx = x, nknots = 50) {
  oldwarn <- options(warn = -1)

  weight <- !is.null(w)

  miss <- is.na(y)
  if (weight) {
    miss <- miss | w < 1e-9
  }
  yy <- y[!miss]
  xx <- x[!miss]
  if (weight) {
    w <- w[!miss]
    w <- w / sum(w) * length(w)
    f.ug <- cobs::cobs(xx, yy, w = w, print.warn = FALSE, print.mesg = FALSE, lambda = lambda, nknots = nknots)
  } else {
    f.ug <- cobs::cobs(xx, yy, print.warn = FALSE, print.mesg = FALSE, lambda = lambda, nknots = nknots)
  }

  fred <- stats::predict(f.ug, interval = "conf", nz = 200)
  fit <- stats::approx(fred[, 1], fred[, 2], xout = newx)$y
  se <- stats::approx(fred[, 1], (fred[, 4] - fred[, 3]) / 2 / 1.96, xout = newx)$y

  if (max(xx) > b) {
    delta <- (max(xx) - min(xx)) / 10
    xxx <- xx[xx > (b - delta)]
    yyy <- yy[xx > (b - delta)]
    if (weight) {
      f.mono <- cobs::cobs(xxx, yyy, constraint = "increase", w = w[xx > (b - delta)], print.warn = FALSE, print.mesg = FALSE, lambda = lambda, nknots = nknots)
    } else {
      f.mono <- cobs::cobs(xxx, yyy, constraint = "increase", print.warn = FALSE, print.mesg = FALSE, lambda = lambda, nknots = nknots)
    }
    fred <- stats::predict(f.mono, interval = "conf", nz = 200)
    newfit <- stats::approx(fred[, 1], fred[, 2], xout = newx[newx > (b - delta)])$y
    newse <- stats::approx(fred[, 1], (fred[, 4] - fred[, 3]) / 2 / 1.96, xout = newx[newx > (b - delta)])$y
    preb <- sum(newx <= (b - delta))
    newfit <- c(rep(0, preb), newfit)
    newse <- c(rep(0, preb), newse)
    postb <- sum(newx > b)
    n <- length(newx)
    cc <- c(rep(0, preb), seq(0, 1, length = n - preb - postb), rep(1, postb))
    fit <- (1 - cc) * fit + cc * newfit
    se <- (1 - cc) * se + cc * newse
  }
  options(warn = oldwarn$warn)
  return(list(fit = fit, se = se))
}