Monotonic/shape constrained splines still possible in brms?

Interfaces brms
1

Hoping that @paul.buerkner or anyone else might be able to weigh in on whether it is (still?) possible to include monotonic spline terms in models fit with brms. I saw a thread from Jan '19 indicating that it is workable, though possibly involving some modifications to the scam pacakge. However, when I run the minimal working example posted on that thread (here) after modifying the scam source code as suggested, the resulting smooths are decidedly not monotonic. Of course, alternative suggestions for incorporating monotonic smooth of continuous covariates would be also welcome. Thanks!

# minimal example
library(brms)
library(scam) # N.B. this is compiled from source with modifications per referenced post
b = brm(time ~ s(age, bs = "mpi"), data = kidney)
marginal_smooths(b, spaghetti = T, nsamples = 10)

# system info
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          6.2                         
year           2019                        
month          12                          
day            12                          
svn rev        77560                       
language       R                           
version.string R version 3.6.2 (2019-12-12)
nickname       Dark and Stormy Night 

# brms version: 2.10.0
3 Likes
2

I don’t know to be honest. I am also confused why scam doesn’t naturally play well with the standard post-processing of mgcv that I apply in brms but I am sure the authors have good reasons for this other behavior. It may be worth asking the package maintainer of scam about this. Right now, brms works well only with smooths that follow the standard mgcv logic so can be reliability post-processed in brms.

2 Likes
3

Just a quick followup: I followed up with the author of the scam package. She confirmed that the fitting and postprocessing procedures in scam differ from the those in mgcv due to the shape constraints.

4

Thanks for following up on this. Perhaps there will be some way to eventually have penalized monotone splines (in whatever form) in brms. I am open for suggestions.

5

Monotonic splines in mcgv require “only” sampling fullfilling a constraint matrix C.

library(mgcv)
library(MASS)
library(rstan)
# x <- seq(0,1,length.out=10)
# y <- +0.2 * x + x^2*0.3 + 2
x <- 1:10
y=c(100,41,22,10,6,7,2,1,3,1)

sigma <- 0.05
y.noise <- y + rnorm(length(x), 0,sigma)
df <- data.frame(x=x, y=y.noise)

## Set up the size of the basis functions/number of knots
k <- 5
## This fits the unconstrained model but gets us smoothness parameters that
## that we will need later
unc <- gam(y ~ s(x, k = k, bs = "cr"), data = df)

## This creates the cubic spline basais functions of `x`
## It returns an object containing the penalty matrix for the spline
## among other things; see ?smooth.construct for description of each
## element in the returned object
sm <- smoothCon(s(x, k = k, bs = "cr"), df, knots = NULL)[[1]]

## This gets the constraint matrix and constraint vector that imposes
## linear constraints to enforce montonicity on a cubic regression spline
## the key thing you need to change is `up`.
## `up = TRUE` == increasing function
## `up = FALSE` == decreasing function (as per your example)
## `xp` is a vector of knot locations that we get back from smoothCon

K <- 5
df.sm=data.frame(x=x)
sm <- smoothCon(s(x, k = K, bs = "cr"), df.sm, knots = NULL
  ,  absorb.cons=FALSE  # no mean in spline. add. parameter in Stan
  ,  diagonal.penalty = TRUE  # diagonal smooth matrix
)[[1]]

Nzero <- sum(diag(sm$S[[1]]) == 0.0)

F <- mono.con(sm$xp, up = FALSE)   # get constraints: up = FALSE == Decreasing constraint!
C<-F$A

The Question is now how can we sample parameters fulfilling the constraints in matrix C?
Maybe we can extend @bgoodri Constrained multinormal?

3 Likes
6

Dear @paul.buerkner,

Could you perhaps clarify for me at which stage fitting a scam monotonic spline with brms fails?

If I try something like this:
b1 <- brm(y~ s(x, bs = “mpi”), data = dat1)

It generates samples and fitted() and predict() work, but the result is not monotonic.

But, is this just the fitted and predict methods missing something? Because if I try to predict from the lpmatrix and coefficients of a model fit with scam, it also does not create the monotonic response.

I have a use-case where there will only ever be 1 continuous predictor variable. So if the model is actually fitting correctly, then I was hoping to be able to hack together a predict method for this simple case.

extract_draws() and posterior_samples() both seem to work, and I can find what I think is the lpmatrix. But I still only seem to be able to get the unconstrained result.

Thanks for any help you can give.

Andrew.

Example here:

# Reproducible example

library(mgcv)
library(scam)
library(brms)

# setup some non monotonic data

set.seed(202000613)

dat1 <- data.frame(
  x = seq(0, 0.75*pi, length.out = 10) 
)
dat1$y <- sin(dat1$x) + rnorm(10, 0, 0.05)
plot(y~x, data  = dat1)

# fit with unconstrained spline for comparison
g1 <- gam(y ~ s(x), data = dat1)

plot(dat1$x, fitted(g1), type = "l")
points(y~x, data = dat1)


# Fit with constrained spline.
s1 <- scam(y ~ s(x, bs = "mpi"), data = dat1)

plot(dat1$x, fitted(s1), type = "l")
points(y~x, data = dat1)


# Now with brms

b1 <- brm(y ~ s(x, bs = "mpi"), data = dat1, save_all_pars = TRUE)

# not monotonic
plot(dat1$x, predict(b1)[,1], type = "l")
points(y~x, data = dat1)


# BUT
# neither is this if I try to predict "manually" from the standard 
# coefficients from the scam model

m <- predict(s1, type = "lpmatrix")
cfs <- coef(s1)

dat1$scam.pred.1 <- m %*% cfs
plot(scam.pred.1~x, data = dat1)


# I have to do this
dat1$scam.pred.2 <- m %*% s1$beta.t

plot(scam.pred.2~x, data = dat1, type = "l", ylim = c(0,1))
points(y~x, data = dat1)

# So maybe I can find the equivalent of beta.t in the brms output?
7

I am sorry, I don’t know. s() splines in brms work because it talks to mgcv in some standardized way, but this way does not work with scam provided splines for reasons the authors of scam and mgcv will know better.

8

Thanks for answering. I’ll try to dig into it a bit more myself and maybe contact the scam package author.

9

Dear Andrew,

even though s() is built internally with the mpi constraint (through reparametrization), it is not considered within predict (if predict.scam is not used).

Easiest way would be, that scam already internally reparametrize such that the usual lpmatrix %*% coef() method works as desired and yields predictions with the desired shape constraint.

If you check the output of coef(), the coefficients should (in the case of mpi) be already a monotone sequence, which is not the case.

JFM

1 Like