Shape constrained (e.g., monotonic) smooths in brms

tiflo · January 3, 2019, 2:04am

I was delighted to see that brms allows shape constrained smooths as from the scam package as part of its formula syntax. However, the post-processing such as marginal_effects() does not seem to work. I found this old post by Paul https://groups.google.com/d/msg/brms-users/hYF6jVEvOvU/EKT6qaP_AwAJ, but even if I follow those steps and load the script (available at the link above), I do get the same error message. I’d appreciate any pointers you might have.

library(scam)
Predict.matrix.mpi.smooth<-function(object,data)
## prediction method function for the `mpi' smooth class
{ m <- object$m+1; # spline order, m+1=3 default for cubic spline
  q <- object$df ## +1
  Sig <- matrix(0,q,q)   # Define Matrix Sigma
  # elements of matrix Sigma for increasing smooth
  for (i in 1:q)  Sig[i,1:i] <- 1
  ## find spline basis inner knot range...
  ll <- object$knots[m+1];ul <- object$knots[length(object$knots)-m]
  m <- m + 1
  x <- data[[object$term]]
  n <- length(x)
  ind <- x<=ul & x>=ll ## data in range
  if (sum(ind)==n) { ## all in range
     X <- splines::spline.des(object$knots,x,m)$design
     X <- X[,2:(q+1)]%*%Sig ## X <- X%*%Sig 
  } else { ## some extrapolation needed 
     ## matrix mapping coefs to value and slope at end points...
     D <- splines::spline.des(object$knots,c(ll,ll,ul,ul),m,c(0,1,0,1))$design
     X <- matrix(0,n,ncol(D)) ## full predict matrix
     if (sum(ind)> 0)  X[ind,] <- splines::spline.des(object$knots,x[ind],m)$design ## interior rows
     ## Now add rows for linear extrapolation...
     ind <- x < ll 
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ll)%*%D[1:2,]
     ind <- x > ul
     if (sum(ind)>0) X[ind,] <- cbind(1,x[ind]-ul)%*%D[3:4,]
     X <- X[,2:(q+1)]%*%Sig 
  }
  X
}

library(brms)
b = brm(time ~ s(age, bs = "mpi"),
    data = kidney)
marginal_effects(b)
# yields: Error in X %*% re$trans.U : non-conformable arguments


> version
               _                           
platform       x86_64-apple-darwin15.6.0   
arch           x86_64                      
os             darwin15.6.0                
system         x86_64, darwin15.6.0        
status                                     
major          3                           
minor          5.1                         
year           2018                        
month          07                          
day            02                          
svn rev        74947                       
language       R                           
version.string R version 3.5.1 (2018-07-02)
nickname       Feather Spray

tiflo · January 6, 2019, 10:08pm

@paul.buerkner I should have tagged you in this, as you wrote the script I’m copying above. Added you now.

paul.buerkner · January 9, 2019, 1:47pm

Not sure what caused the change, but apparently, mgcv now always calls scam::Predict.matrix.mpi.smooth internally instead of the function defined in the global environment. As a result, defining the new function won’t change anything because it is not used by mgcv.

My suggestion is that you download the source of scam from CRAN, replace the function under discussion in the source code and build scam locally. Not the most straight forward solution but one that will likely work.

tiflo · January 9, 2019, 7:47pm

Thank you!

tiflo · January 24, 2026, 6:40pm

@paul.buerkner Just an update that Natalya Pya Arnqvist (the developer of the scam package; not on discourse) has kindly addressed the issue discussed here. She and I did some testing and the monotonic smooths offered by scamnow seem to play nicely with brms’s predict function. The new version of the scam package should now be available on CRAN. Thanks to Natalya!

avehtari · January 24, 2026, 7:34pm

Awesome!

paul.buerkner · January 26, 2026, 1:16pm

Great! Thank you for the update!

zacrobinson5 · July 9, 2026, 5:40pm

Pardon my ignorance if I’m missing something here, but from what I can tell the updates to the scam package allow the predict function to work, but the shape constraints are not occurring:

##load packages

pacman::p_load(tidyverse,tidybayes,brms,scam)

## simulate nonlinear data

set.seed(123)
dat = data.frame(x = 1:100) %>%
mutate(y = x^2 * 0.01 + rnorm(n = nrow(.)))

## visualize data

ggplot(data = dat,
aes(x = x,
y = y)) +
geom_point()

## fit monotonic decreasing brms model

mod_brm = brm(
data = dat,
y ~ s(x, bs = "mpd"),
seed = 123,
cores = 4,
chains = 4,
backend = "cmdstanr"
)

## fit monotonic decreasing scam model

mod_gam = scam(
data = dat,
y ~ s(x, bs = "mpd")
)

## create predictions

preds =
rbind(
data.frame(x = 1:100,
model = "brm()") %>%
mutate(y = predict(mod_brm,newdata = data.frame(x = 1:100))[,"Estimate"]),
data.frame(x = 1:100,
model = "scam()") %>%
mutate(y = predict(mod_gam,newdata = data.frame(x = 1:100)))
)

## plot predcitions

ggplot(data = preds,
aes(x = x,
y = y))+
theme_test()+
geom_point(data = dat,
color = "black")+
geom_line(aes(color = model))

tiflo · July 10, 2026, 1:08pm

Hi Zac,

in your example, you seem to be using a decreasing shape-constraint, but the data have an increasing trend. So, I think the solution the model found—a flat line—is the best possible under the constraint you imposed. Sorry, if I’m misreading this!

Florian

zacrobinson5 · July 10, 2026, 2:43pm

Hi Florian,

No worries at all. The example here is intended to show how brms is NOT enforcing the monotonic constraints that the scam package is. In this example, as you mentioned, the data are monotonically increasing, thus, a smooth term with bs = “mpd” should result in a roughly flat line at the mean. This behavior is observed from the predicted values in the scam package but not brms.

Hopefully that helps to clarify what I was trying to demonstrate - any help here would be appreciated, thanks!

tiflo · July 10, 2026, 4:30pm

Ooops, yes, I see it now (I misread the graph; shouldn’t try to answer queries before my first cup of coffee =)).

That indeed is an issue @paul.buerkner!

paul.buerkner · July 22, 2026, 12:22pm

It is an issue indeed, but I don’t quite know what do you about it. Since brms leaves all the spline details to mgcv, I have little control over what is happening inside the splines. Can you check what happens with you run your reprex with gam, gamm and gamm4 functions? Do they behave like scam or like brms?

zacrobinson5 · July 23, 2026, 2:05pm

Hey Paul, thanks for chiming in here. I have adjusted the code as you suggested. gam() doesn’t even recognize the monotonicbsarguments, and the other functions (i.e., gamm() and gamm4()) result in the same behavior as brms). Seems like the fix might need to happen inside mgcv then?

##load packages

pacman::p_load(tidyverse,tidybayes,brms,scam,mgcv,gamm4)

simulate nonlinear data

set.seed(123)
dat = data.frame(x = 1:100) %>%
mutate(y = x^2 * 0.01 + rnorm(n = nrow(.)))

visualize data

ggplot(data = dat,
aes(x = x,
y = y)) +
geom_point()

fit monotonic decreasing brms model

mod_brm = brm(
data = dat,
y ~ s(x, bs = "mpd"),
seed = 123,
cores = 4,
chains = 4,
backend = "cmdstanr"
)

fit monotonic decreasing scam model

mod_scam = scam(
data = dat,
y ~ s(x, bs = "mpd")
)

fit monotonic decreasing gam model

mod_gam = gam(

data = dat,

y ~ s(x, bs = "mpd")

)

Error in s(x, bs = "mpd") : unused argument (bs = "mpd")

fit monotonic decreasing gamm model

mod_gamm = gamm(
data = dat,
y ~ s(x, bs = "mpd")
)

fit monotonic decreasing gamm4 model

mod_gamm4 = gamm4(
data = dat,
y ~ s(x, bs = "mpd")
)

create predictions

preds =
rbind(
data.frame(x = 1:100,
model = "brm()") %>%
mutate(y = predict(mod_brm,newdata = data.frame(x = 1:100))[,"Estimate"]),
data.frame(x = 1:100,
model = "scam()") %>%
mutate(y = predict(mod_scam,newdata = data.frame(x = 1:100))),
data.frame(x = 1:100,
model = "gamm()") %>%
mutate(y = predict(mod_gamm$gam,newdata = data.frame(x = 1:100))),
data.frame(x = 1:100,
model = "gamm4()") %>%
mutate(y = predict(mod_gamm4$gam,newdata = data.frame(x = 1:100)))
)

plot predcitions

ggplot(data = preds,
aes(x = x,
y = y))+
theme_test()+
facet_wrap(~model)+
geom_point(data = dat,
color = "black")+
geom_line(aes(color = model),
linewidth = 1)

jsocolar · July 23, 2026, 6:00pm

The monotonicity constraint doesn’t arise from anything special about the basis, but rather from constraints on the coefficients. I think that the monotonicity constraint arises from an ordering constraint on coefficients of B-spline basis functions; I cannot easily say what the corresponding constraint is when the spline is expressed in random effect form. In either case, brms inherits the random-effect form of the basis and then applies coefficients in the usual way. You can see that the direction of monotonicity definitely is not being enforced from the fact that there are no constraints on the coefficients (so you could negate the all and switch the sign of the slope everywhere). Moreover, I’m pretty sure that brms isn’t enforcing monotonicity at all here, though I can’t easily check because when I try to extract predictions, linear predictors, etc from this model I get an unrelated error. Maybe a quirk of my setup, though I have no idea what’s causing it.

tiflo · July 23, 2026, 7:27pm

I’ve notified the scam developer (Natalya Pya Arnqvist) of this post to see whether she might have some insights to offer. She’s out of office at the moment, so it might be a wee moment.

Topic		Replies	Views
Monotonic/shape constrained splines still possible in brms? brms	8	2115	May 19, 2021
Monotonic splines in brms? brms specification	6	1107	September 3, 2025
Gam contrainted to pass through fixed value - add argument paraPen? Modeling brms	0	251	March 26, 2024
Spatial modeling using mgcv soap film smoothing brms	15	2654	March 29, 2022
Hierarchical brms with splines (limiting prediction below threshold) brms techniques , brms	6	1616	August 2, 2022

Shape constrained (e.g., monotonic) smooths in brms

Related topics