# mgcv gam: imposing constraints on coefficients

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?

• Have you tried running the regression without any constraints to see if the coefficient by chance does lie between -1 and 1? ... and I think you can ignore the bot, your question is quite clear. Commented Feb 22, 2023 at 16:43
• To add to @jbowman's suggestion, if the coefficient is not in that range, you can run two other models: y - x2 ~ s(x1, ...) and y + x2 ~ s(x1, ...). Pick the better of the two.
– whuber
Commented Feb 22, 2023 at 16:52
• @jbowman I run many regressions (for many stocks), and some are in the range and others are not Commented Feb 22, 2023 at 19:05
• @whuber are you sure that if the unconstrained coefficient was e.g. >1, then constrained would be 1? It seems like likelihood functions may be quite non-linear, so within the [-1, 1] range the optimal value could be 0.2 or whatever Commented Feb 22, 2023 at 19:07
• He is suggesting you check both endpoints. If the optimal value is $0.2$, then that's what you'll get with your initial estimate, and you'll know you don't need to do anything about the constraint, as it's satisfied. Since your model is linear in $x_2$, and, at least as written, the likelihood function is the default Gaussian, it will not be multimodal. Commented Feb 22, 2023 at 20:11

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. • This is really clever - why hasn't it been accepted as the right answer? – Bob Commented Apr 27 at 3:17 • @Bob—Thanks! Sometimes, posters figure things out for themselves and don't return to the site to see if there's an answer. The fact that the OP had posted an answer themselves a few hours earlier than this one (see below) indicates that's probably what happened here. Commented Apr 27 at 15:59 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
}

• It is unclear that this accomplishes what you want or why it would. Could you explain it or at least show an example in which it clearly succeeds in implementing a constraint?
– whuber
Commented Feb 23, 2023 at 15:22