Spatial modeling using mgcv soap film smoothing

Hi Paul,

I’m interested in reproducing Gavin Simpson’s example on soap film smoothers and lake bathymetries (URL: https://www.fromthebottomoftheheap.net/2016/03/27/soap-film-smoothers/) using brms.

His model structure is as follows:

library(mgcv)
m2 <- gam(-depth ~ s(os_x, os_y, bs = "so", xt = list(bnd = bound)),
          data = depth, method = "REML", knots = knots)

where bound is a list of spatial points (x,y) representing the lake perimeter to constrain the predictions and knots is a data frame of spatial points (x,y) representing a sequence of grid points that fall within that boundary.

What would be a comparable model structure using brmsformula syntax? I haven’t been able to figure it out on my own yet and am looking for your insight.

Thanks a ton,
Denise

  • Operating System: Windows
  • brms Version: 2.9.0

It should work with the exact same syntax. What is the error message when trying that?

The error is as follows:

Error in smooth.construct.so.smooth.spec(object, dk$data, dk$knots) :
need at least one interior knot

Hmm. Can you provide a minimal reproducible example for me to try out?

Yes- absolutely - will prepare that right now thanks.

Here is a reproducible example:

library(tidyverse)
library(mgcv)
library(brms)

x <- c(0,  0,  1, 1,0)
y <- c(0, 1, 1, 0, 0)

x2 <- runif(100)
y2 <- runif(100)

z = rnorm(100)

seq(0.1,0.9, length=5)
expand.grid(seq(0.1,0.9, length=5), seq(0.1,0.9, length=5)) %>% rename(x = Var1, y=Var2) -> knots_test
tibble(z = z, x= x2, y=y2) -> testdata

gam(formula = z ~ s(x, y, bs = "so", xt = list(bnd = list(list(x=x, y=y)))),
    data = testdata, method = "REML",
    knots = knots_test)
# works

brm(formula = z ~ s(x, y, bs = "so", xt = list(bnd = list(list(x=x, y=y)))),
    data = testdata,
    knots = knots_test)

# need at least one interior knot

Thanks! I will take a look.

1 Like

Fixed the error you see, but now we get

Error in smooth2random.mgcv.smooth(sm, names(data), type = 2) : 
  Can not convert this smooth class to a random effect 

which is no longer a bug but a limitation of brms’ spline functionality.

2 Likes

I am also very much interested in this and having the same problem, thanks paul

This is a limitation of how the GAMs smooth functions are represented as random effects. It’s the same reason why you can only use t2() for tensor product smooths with brm(), gamm4::gamm4() (and hence lmer() or glmer() in the lme package), whereas you can use te(), ti() and t2() smooths in gam(), bam() etc.

Is there a particular reason you wanted the fully Bayesian fit that brm() would give you? GAMs fitted using gam() and bam() are (empircal) Bayesian models and you can get at the posterior distribution of the model and sample from it in a similar sense to predict() and fitted() applied to a brm() model. The obvious difference is that mgcv doesn’t do the actual sampling and is using improper priors, and there isn’t a prior on the smoothness parameter(s) in the model (but you can to some extent account for this when sampling from the posterior).

Unfortunately there is nothing I can do about it at this point.