# How to smooth data and force monotonicity

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()


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.

• 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) Commented Feb 21, 2016 at 2:58 • Similar Q with answer: stats.stackexchange.com/questions/206073/… Commented Jan 31, 2019 at 15:47 ## 2 Answers 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:

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

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.

• 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
– Ben
Commented Feb 19, 2016 at 21:01
• 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. Commented Feb 19, 2016 at 21:10
• 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 :-) Commented Feb 19, 2016 at 21:12
• 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.
– Ben
Commented Feb 19, 2016 at 21:16
• Also, you might look at the smooth.monotone function in the fda package. From what I saw it looks like you need to jump through similar hoops as with mgcv like setting up your basis functions and generating an initial model to fit from, but you might find that function easier to use of more informative. Commented Feb 19, 2016 at 21:53

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.

• amazing thanks - hope that's not a scam! Commented May 12, 2021 at 9:58