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

Interfaces brms
1

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
2

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

3

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
4

Thank you!