14
$\begingroup$

I have some data which I would like to smooth so that the smoothed points are monotonically decreasing. My data sharply decreases and then begins to plateau. Here's an example using R

df <- data.frame(x=1:10, y=c(100,41,22,10,6,7,2,1,3,1))
ggplot(df, aes(x=x, y=y))+geom_line()

plot of data to smooth

What's a good smoothing technique I could use? Also, it'd be nice if I can force the 1st smoothed point to be close to my observed point.

$\endgroup$
  • 1
    $\begingroup$ I notice your example values are integer. Are your real values counts? If they were, then (while this is no guarantee of monotonicity, for data like these it will generally give it anyway), something like this might be useful: plot(y~x,data=df); f=fitted( glm( y~ns(x,df=4), data=df,family=quasipoisson)); lines(df$x,f) $\endgroup$ – Glen_b -Reinstate Monica Feb 21 '16 at 2:58
  • $\begingroup$ Similar Q with answer: stats.stackexchange.com/questions/206073/… $\endgroup$ – kjetil b halvorsen Jan 31 '19 at 15:47
18
$\begingroup$

You can do this using penalised splines with monotonicity constraints via the mono.con() and pcls() functions in the mgcv package. There's a little fiddling about to do because these functions aren't as user friendly as gam(), but the steps are shown below, based mostly on the example from ?pcls, modified to suit the sample data you gave:

df <- data.frame(x=1:10, y=c(100,41,22,10,6,7,2,1,3,1))

## Set up the size of the basis functions/number of knots
k <- 5
## This fits the unconstrained model but gets us smoothness parameters that
## that we will need later
unc <- gam(y ~ s(x, k = k, bs = "cr"), data = df)

## This creates the cubic spline basis functions of `x`
## It returns an object containing the penalty matrix for the spline
## among other things; see ?smooth.construct for description of each
## element in the returned object
sm <- smoothCon(s(x, k = k, bs = "cr"), df, knots = NULL)[[1]]

## This gets the constraint matrix and constraint vector that imposes
## linear constraints to enforce montonicity on a cubic regression spline
## the key thing you need to change is `up`.
## `up = TRUE` == increasing function
## `up = FALSE` == decreasing function (as per your example)
## `xp` is a vector of knot locations that we get back from smoothCon
F <- mono.con(sm$xp, up = FALSE)   # get constraints: up = FALSE == Decreasing constraint!

Now we need to fill in the object that gets passed to pcls() containing details of the penalised constrained model we want to fit

## Fill in G, the object pcsl needs to fit; this is just what `pcls` says it needs:
## X is the model matrix (of the basis functions)
## C is the identifiability constraints - no constraints needed here
##   for the single smooth
## sp are the smoothness parameters from the unconstrained GAM
## p/xp are the knot locations again, but negated for a decreasing function
## y is the response data
## w are weights and this is fancy code for a vector of 1s of length(y)
G <- list(X = sm$X, C = matrix(0,0,0), sp = unc$sp,
          p = -sm$xp, # note the - here! This is for decreasing fits!
      y = df$y,
          w = df$y*0+1)
G$Ain <- F$A    # the monotonicity constraint matrix
G$bin <- F$b    # the monotonicity constraint vector, both from mono.con
G$S <- sm$S     # the penalty matrix for the cubic spline
G$off <- 0      # location of offsets in the penalty matrix

Now we can finally do the fitting

## Do the constrained fit 
p <- pcls(G)  # fit spline (using s.p. from unconstrained fit)

p contains a vector of coefficients for the basis functions corresponding to the spline. To visualize the fitted spline, we can predict from the model at 100 locations over the range of x. We do 100 values so as to get a nice smooth line on the plot.

## predict at 100 locations over range of x - get a smooth line on the plot
newx <- with(df, data.frame(x = seq(min(x), max(x), length = 100)))

To generate predicted values we use Predict.matrix(), which generates a matrix such that when multiple by coefficients p yields predicted values from the fitted model:

fv <- Predict.matrix(sm, newx) %*% p
newx <- transform(newx, yhat = fv[,1])

plot(y ~ x, data = df, pch = 16)
lines(yhat ~ x, data = newx, col = "red")

This produces:

enter image description here

I'll leave it up to you to get the data into a tidy form for plotting with ggplot...

You can force a closer fit (to partially answer your question about having the smoother fit the first data point) by increasing the dimension of the basis function of x. For example, setting k equal to 8 (k <- 8) and rerunning the code above we get

enter image description here

You can't push k much higher for these data, and you have to be careful about over fitting; all pcls() is doing is solving the penalised least squares problem given the constraints and the supplied basis functions, it's not performing smoothness selection for you - not that I know of...)

If you want interpolation, then see the base R function ?splinefun which has Hermite splines and cubic splines with monotonicty constraints. In this case you can't use this however as the data are not strictly monotonic.

$\endgroup$
  • $\begingroup$ Thanks. I'm sure your solution is an appropriate one, but it's so complex and obfuscated that I just can't use it. splinefun was my initial thought as well (I am interpolating) but spline(x=df$x, y=df$y, n=nrow(df), method="monoH.FC") and spline(x=df$x, y=df$y, n=nrow(df), method="hyman") both raise errors $\endgroup$ – Ben Feb 19 '16 at 21:01
  • 1
    $\begingroup$ If you just try, I'm sure you can use it; I have little idea of what's going on under the hood here but I worked it out, and I have indicated the places you'd need to change things. Assuming you know some R of course. Most of the detail is implementational that you can ignore if all you want to do it fit a monotonically constrained spline. Would you like me to annotate the code a bit more to highlight more about what each step is doing? The reference in ?mono.con has further detail on the method. $\endgroup$ – Reinstate Monica - G. Simpson Feb 19 '16 at 21:10
  • $\begingroup$ As for why splinefun is raising an error; I've just realised, you can fit an monotonic spline that interpolates data that is not itself monotonic. The observation at x = 6 is at greater y that the observations at x = 5. You'll just have to ignore that part of the answer :-) $\endgroup$ – Reinstate Monica - G. Simpson Feb 19 '16 at 21:12
  • $\begingroup$ Got it. And no need - I'm a pretty experienced R user. I just like to understand the math behind what I use and this solution appears to have quite a lot going on under the hood. Thanks again for your help. $\endgroup$ – Ben Feb 19 '16 at 21:16
  • $\begingroup$ I've added some notes to explain what each things is or does; the main point to note is that the monotonicity constraints are being imposed by a specific set of inequality constraints that mono.con returns for a cubic spline. ?pcls has examples for thin plate splines and additive models which are less user-friendly than the above, but which might expose some more of the math if you are familiar with the maths for those types of spline (I'm not that familiar myself). $\endgroup$ – Reinstate Monica - G. Simpson Feb 19 '16 at 21:42
13
$\begingroup$

The recent scam package by Natalya Pya and based on the paper "Shape constrained additive models" by Pya & Wood (2015) can make part of the process mentioned in Gavin's excellent answer much easier.

library(scam)
con <- scam(y ~ s(x, k = k, bs = "mpd"), data = df)
plot(con)

There are a number of bs functions you can use - in the above I used mpd for "monotonic decreasing P-spline" but it also has functions that enforce convexity or concavity either separately or alongside the monotonic constraints.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.