Force GAM model fit to be monotonic and go through a fixed point (x0, y0) with R mgcv

Question

I am trying to fit a GAM model to data under two constraints simultatenously: (1) the fit is monotonic (increasing), (2) the fit goes through a fixed point, say, (x0,y0).

So far, I managed to have these two constraints work separately:

• For (1), based on mgcv::pcls() documentation examples, by using mgcv::mono.con() to get linear constraints sufficient for monotonicity, and estimate model coefs via mgcv::pcls(), using the constraints.

• For (2), based on this post, by setting the value of spline at knot location x0 to 0 + using offset term in the model formula.

However, I struggle to combine these two constraints simultaneously. I guess a way to go is mgcv::pcls(), but I could work out neither (a) doing a similar trick of setting the value of spline at knot location x0 to 0 + using offset nor (b) setting equality constraint(s) (which I think could yield my (2) constraint setup).

I also note that the approach for setting the value of spline at knot location x0 to 0 for my constrain condition (2) yields weirdly wiggly outcome (as compared to unconstrained GAM fit) -- as showed below.

Attempt so far: fit a smooth function to data under two constraints separately

Simulate some data

library(mgcv)
set.seed(1)
x <- sort(runif(100) * 4 - 1)
f <- exp(4*x)/(1+exp(4*x))
y <- f + rnorm(100) * 0.1
dat <- data.frame(x=x, y=y)

GAM unconstrained (for comparison)

k <- 13
fit0 <- gam(y ~ s(x, k = k, bs = "cr"), data = dat)
# predict from unconstrained GAM fit
newdata <- data.frame(x = seq(-1, 3, length.out = 1000))
newdata$y_pred_fit0 <- predict(fit0, newdata = newdata)

GAM constrained: (1) the fit is monotonic (increasing)

k <- 13

# Show regular spline fit (and save fitted object)
f.ug <- gam(y~s(x,k=k,bs="cr"))

# explicitly construct smooth term's design matrix
sm <- smoothCon(s(x,k=k,bs="cr"),dat,knots=NULL)[[1]]
# find linear constraints sufficient for monotonicity of a cubic regression spline
# it assumes "cr" is the basis and its knots are provided as input
F <- mono.con(sm$xp)

G <- list(
  X=sm$X,
  C=matrix(0,0,0),  # [0 x 0] matrix (no equality constraints)
  sp=f.ug$sp,       # smoothing parameter estimates (taken from unconstrained model)
  p=sm$xp,          # array of feasible initial parameter estimates
  y=y,           
  w= dat$y * 0 + 1  # weights for data    
)
G$Ain <- F$A        # matrix for the inequality constraints
G$bin <- F$b        # vector for the inequality constraints
G$S <- sm$S         # list of penalty matrices; The first parameter it penalizes is given by off[i]+1
G$off <- 0          # Offset values locating the elements of M$S in the correct location within each penalty coefficient matrix.  (Zero offset implies starting in first location)

p <- pcls(G);       # fit spline (using smoothing parameter estimates from unconstrained fit)

# predict 
newdata$y_pred_fit2 <- Predict.matrix(sm, data.frame(x = newdata$x)) %*% p
# plot
plot(y ~ x, data = dat)
lines(y_pred_fit0 ~ x, data = newdata, col = 2, lwd = 2)
lines(y_pred_fit2 ~ x, data = newdata, col = 4, lwd = 2)

Blue line: constrained; red line: unconstrained

GAM constrained: (2) fitted go through (x0,y0)=(-1, -0.1)

k <- 13

## Create a spline basis and penalty
## Make sure there is a knot at the constraint point (here: -1)
knots <- data.frame(x = seq(-1,3,length=k)) 
# explicit construction of a smooth term in a GAM
sm <- smoothCon(s(x,k=k,bs="cr"), dat, knots=knots)[[1]]

## 1st parameter is value of spline at knot location -1, set it to 0 by dropping
knot_which <- which(knots$x == -1)
X <- sm$X[, -knot_which]                  ## spline basis
S <- sm$S[[1]][-knot_which, -knot_which]  ## spline penalty
off <- dat$y * 0 + (-0.1)                   ## offset term to force curve through (x0, y0)

## fit spline constrained through (x0, y0)
gam_1 <- gam(y ~ X - 1 + offset(off), paraPen = list(X = list(S)))

# predict (add offset of -0.1)
newdata_tmp <- Predict.matrix(sm, data.frame(x = newdata$x))
newdata_tmp <- newdata_tmp[, -knot_which]    
newdata$y_pred_fit1 <- (newdata_tmp %*% coef(gam_1))[, 1] + (-0.1)

# plot
plot(y ~ x, data = dat)
lines(y_pred_fit0 ~ x, data = newdata, col = 2, lwd = 2)
lines(y_pred_fit1 ~ x, data = newdata, col = 3, lwd = 2)
# lines at cross of which the plot should go throught
abline(v=-1, col = 3); abline(h=-0.1, col = 3)

Green line: constrained; red line: unconstrained

Answer 1

I think you could augment the data vectors x and y with (x0, y0) and then put a (really) high weight on the first observation (i.e. add a weight vector to your G list).

Alternatively to the simple weighting strategy, we can write the quadratic programming problem starting from the results of the preliminary smoothing. This is illustrated in the second R-code below (in this case I used p-spline smoothers, see Eilers and Marx 1991).

Hope this helps a bit (a similar problem is discussed here).

Rcode example 1 (weight strategy)

set.seed(123)
N = 100
x <- sort(runif(N) * 4 - 1)
f <- exp(4*x)/(1+exp(4*x))
y <- f + rnorm(N) * 0.1
x = c(-1, x)
y = c(-0.1, y)
dat = data.frame(x = x, y= y)

k <- 13
fit0 <- gam(y ~ s(x, k = k, bs = "cr"), data = dat)
# predict from unconstrained GAM fit
newdata <- data.frame(x = seq(-1, 3, length.out = 1000))
newdata$y_pred_fit0 <- predict(fit0, newdata = newdata)

k <- 13

# Show regular spline fit (and save fitted object)
f.ug <- gam(y~s(x,k=k,bs="cr"))

# explicitly construct smooth term's design matrix
sm <- smoothCon(s(x,k=k,bs="cr"),dat,knots=NULL)[[1]]
# find linear constraints sufficient for monotonicity of a cubic regression spline
# it assumes "cr" is the basis and its knots are provided as input
F <- mono.con(sm$xp)

G <- list(
X=sm$X,
C=matrix(0,0,0),  # [0 x 0] matrix (no equality constraints)
sp=f.ug$sp,       # smoothing parameter estimates (taken from unconstrained model)
p=sm$xp,          # array of feasible initial parameter estimates
y=y,           
w= c(1e8, 1:N * 0 + 1)  # weights for data    
)
G$Ain <- F$A        # matrix for the inequality constraints
G$bin <- F$b        # vector for the inequality constraints
G$S <- sm$S         # list of penalty matrices; The first parameter it penalizes is given by off[i]+1
G$off <- 0          # Offset values locating the elements of M$S in the correct location within each penalty coefficient matrix.  (Zero offset implies starting in first location)

p <- pcls(G);       # fit spline (using smoothing parameter estimates from unconstrained fit)

# predict 
newdata$y_pred_fit2 <- Predict.matrix(sm, data.frame(x = newdata$x)) %*% p
# plot
plot(y ~ x, data = dat)
lines(y_pred_fit0 ~ x, data = newdata, col = 2, lwd = 2)
lines(y_pred_fit2 ~ x, data = newdata, col = 4, lwd = 2)
abline(v = -1)
abline(h = -0.1)

rm(list = ls())
library(mgcv)
library(pracma)
library(colorout)

set.seed(123)
N   = 100
x   = sort(runif(N) * 4 - 1)
f   = exp(4*x)/(1+exp(4*x))
y   = f + rnorm(N) * 0.1
x0  = -1
y0  = -0.1 
dat = data.frame(x = x, y= y)

k = 50
# Show regular spline fit (and save fitted object)
f.ug = gam(y~s(x,k=k,bs="ps"))

# explicitly construct smooth term's design matrix
sm = smoothCon(s(x,k=k,bs="ps"), dat,knots=NULL)[[1]]

# Build quadprog to estimate the coefficients 
scf = sapply(f.ug$smooth, '[[', 'S.scale')
lam = f.ug$sp / scf
Xp  = rbind(sm$X, sqrt(lam) * f.ug$smooth[[1]]$D)  
yp  = c(dat$y, rep(0, k - 2))
X0  = Predict.matrix(sm, data.frame(x = x0))
sm$deriv = 1
X1 = Predict.matrix(sm, data.frame(x = dat$x))
coef_mono = pracma::lsqlincon(Xp, yp, Aeq = X0, beq = y0, A = -X1, b = rep(0, N))

# fitted values
fit = sm$X %*% coef_mono
sm$deriv = 0
xf = seq(-1, 3, len = 1000)
Xf = Predict.matrix(sm, data.frame(x = xf))
fine_fit = Xf %*% coef_mono

# plot
par(mfrow = c(2, 1), mar = c(3,3,3,3))
plot(dat$x, dat$y, pch = 1, main= 'Data and fit')
lines(dat$x, f.ug$fitted, lwd = 2, col = 2)
lines(dat$x, fit, col = 4, lty = 1, lwd = 2)
lines(xf, fine_fit, col = 3, lwd = 2, lty = 2)
abline(h = -0.1)
abline(v = -1)

plot(dat$x, X1 %*% coef_mono, type = 'l', main = 'Derivative of the fit', lwd = 2)
abline(h = 0.0)

Answer 2

The following package seems to implement what you are looking for:

The proposed shape constrained smoothing has been incorporated into generalized additive models with a mixture of unconstrained and shape restricted smooth terms (mono-GAM). [...]

The proposed modelling approach has been implemented in an R package monogam. The model setup is the same as in mgcv(gam) with the addition of shape constrained smooths. In order to be consistent with the unconstrained GAM, the package provides key functions similar to those associated with mgcv(gam).

• Additive models with shape constraints