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

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
``````
1 Like

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

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.

1 Like

Thank you!