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?

  • 2
    $\begingroup$ 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. $\endgroup$
    – jbowman
    Commented Feb 22, 2023 at 16:43
  • 2
    $\begingroup$ 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. $\endgroup$
    – whuber
    Commented Feb 22, 2023 at 16:52
  • $\begingroup$ @jbowman I run many regressions (for many stocks), and some are in the range and others are not $\endgroup$
    – OlgaPp
    Commented Feb 22, 2023 at 19:05
  • $\begingroup$ @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 $\endgroup$
    – OlgaPp
    Commented Feb 22, 2023 at 19:07
  • 1
    $\begingroup$ 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. $\endgroup$
    – jbowman
    Commented Feb 22, 2023 at 20:11

2 Answers 2


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.


# 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.

  • $\begingroup$ This is really clever - why hasn't it been accepted as the right answer? $\endgroup$
    – Bob
    Commented Apr 27 at 3:17
  • $\begingroup$ @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. $\endgroup$
    – jbowman
    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(
      sp=array(0,0), # penalties
      bin=c(beta_min_max[1], -beta_min_max[2]),
      C=matrix(0,0,0), # Matrix containing any linear equality constraints on the problem
    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
  • 1
    $\begingroup$ 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? $\endgroup$
    – whuber
    Commented Feb 23, 2023 at 15:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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