mgcv gam: imposing constraints on coefficients

Question

I want to use pcls to constrain beta coefficient on a linear term in my gam. Specifically, I have a model: gam(y ~ s(x1, ...) + x2, data=dat), and I want to make sure that x2's coefficient lies between -1 and 1. How can I do it?

Answer 1

The code below should do the job for you. It uses BIC as its model selection criterion; you can substitute AIC or whatever else you want, of course.

library(mgcv)
library(data.table)

# Sample data that (almost certainly) violates the constraint
fitting_data <- data.table(x1 = rnorm(100),
                           x2 = rnorm(100),
                           x3 = rnorm(100))
fitting_data[, y := 2*x1 + x2 + sin(x3) + rnorm(100)]

G_gam = bam(y ~ x1 + x2 + s(x3),  data = fitting_data)

if (coef(G_gam)["x1"] < -1 | coef(G_gam)["x1"] > 1) {
  # G_gam_p1 -> the coefficient of x1 = 1
  # G_gam_m1 -> the coefficient of x1 = -1
  G_gam_p1 <- bam((y-x1) ~ x2 + s(x3), data=fitting_data)
  G_gam_m1 <- bam((y+x1) ~ x2 + s(x3), data=fitting_data)
  
  if (BIC(G_gam_p1) < BIC(G_gam_m1)) {
    final_model <- G_gam_p1
  } else {
    final_model <- G_gam_m1
  }

} else {
  final_model <- G_gam
}

Running this with the sample data above results in:

> final_model$formula
(y - x1) ~ x2 + s(x3)

from which it is easily deduced that the coefficient on x1 has been set to 1, which is as hoped for given that the actual coefficient is 2.

Answer 2

This is the solution that I arrived at:

  G_gam = bam(
    y ~ x1 + x2 + s(x3, k=knots),
    data = fitting_data
  )
  gam_coef = coef(G_gam)
  beta = gam_coef[1]
  if (beta < beta_min_max[1] | beta > beta_min_max[2]) {
    x_matrix = predict(G_gam, fitting_data, type="lpmatrix")
    M <- list(
      X=x_matrix,
      p=some_new_starting_coefficients,
      off=array(0,0),
      sp=array(0,0), # penalties
      Ain=matrix(0,2,knots+2),
      bin=c(beta_min_max[1], -beta_min_max[2]),
      C=matrix(0,0,0), # Matrix containing any linear equality constraints on the problem
      y=fitting_data$y,
      w=fitting_data$y*0+1)
    M$Ain[1,] <- c(0,1,rep(0, knots)) # greater than constraint
    M$Ain[2,] <- c(0,-1,rep(0, knots)) # less than constraint
    M$p <- pcls(M)
    
    new_X_matrix <- predict(G_gam,newdata=new_data,type="lpmatrix")
    new_prediction = new_X_matrix %*% M$p
  }