# Spatial modeling using mgcv soap film smoothing

**URL:** <https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463>\
**Category:** brms\
**Created:** [August 21, 2019, 7:13pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463 "2019-08-21T19:13:37Z")\
**Posts on this page:** 16\
**Page:** 1

<div class="post-metadata">

**Author:** ![denise-colombano](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/denise-colombano/32/3917_2.png) [@denise-colombano](https://discourse.mc-stan.org/u/denise-colombano)\
**Post date:** [August 21, 2019, 7:13pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/1 "2019-08-21T19:13:37Z")

</div>

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/](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

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [August 21, 2019, 9:08pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/2 "2019-08-21T21:08:24Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![denise-colombano](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/denise-colombano/32/3917_2.png) [@denise-colombano](https://discourse.mc-stan.org/u/denise-colombano)\
**Post date:** [August 21, 2019, 9:22pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/3 "2019-08-21T21:22:47Z")

</div>

The error is as follows:

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

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [August 21, 2019, 9:24pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/4 "2019-08-21T21:24:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![denise-colombano](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/denise-colombano/32/3917_2.png) [@denise-colombano](https://discourse.mc-stan.org/u/denise-colombano)\
**Post date:** [August 21, 2019, 9:27pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/5 "2019-08-21T21:27:00Z")

</div>

Yes- absolutely - will prepare that right now thanks.

---

<div class="post-metadata">

**Author:** ![StaffanBetner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/staffanbetner/32/11056_2.png) [@StaffanBetner](https://discourse.mc-stan.org/u/StaffanBetner)\
**Post date:** [March 13, 2020, 12:42pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/6 "2020-03-13T12:42:20Z")

</div>

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

```

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [March 13, 2020, 3:10pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/7 "2020-03-13T15:10:44Z")

</div>

Thanks! I will take a look.

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [March 13, 2020, 3:24pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/8 "2020-03-13T15:24:11Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jayb](https://avatars.discourse-cdn.com/v4/letter/j/71c47a/32.png) [@jayb](https://discourse.mc-stan.org/u/jayb)\
**Post date:** [March 26, 2020, 11:16pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/9 "2020-03-26T23:16:00Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![ucfagls](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ucfagls/32/1073_2.png) [@ucfagls](https://discourse.mc-stan.org/u/ucfagls)\
**Post date:** [March 26, 2020, 11:30pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/10 "2020-03-26T23:30:28Z")

</div>

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).

---

<div class="post-metadata">

**Author:** ![paul.buerkner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/paul.buerkner/32/3303_2.png) [@paul.buerkner](https://discourse.mc-stan.org/u/paul.buerkner)\
**Post date:** [March 27, 2020, 7:06am UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/11 "2020-03-27T07:06:26Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![denise-colombano](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/denise-colombano/32/3917_2.png) [@denise-colombano](https://discourse.mc-stan.org/u/denise-colombano)\
**Post date:** [April 1, 2020, 7:06pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/12 "2020-04-01T19:06:00Z")

</div>

Thanks for looking into this- I ended up using mgcv instead of calling it through brms and the result was very slick.

---

<div class="post-metadata">

**Author:** ![StaffanBetner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/staffanbetner/32/11056_2.png) [@StaffanBetner](https://discourse.mc-stan.org/u/StaffanBetner)\
**Post date:** [April 2, 2020, 9:03am UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/13 "2020-04-02T09:03:09Z")

</div>

> [@ucfagls](#):
>
> This is a limitation of how the GAMs smooth functions are represented as random effects.

But isn’t the limitation rather in the implementation in `mgcv::smooth2random` rather than in the equivalence between smooths and random effects?

---

<div class="post-metadata">

**Author:** ![ucfagls](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/ucfagls/32/1073_2.png) [@ucfagls](https://discourse.mc-stan.org/u/ucfagls)\
**Post date:** [April 2, 2020, 4:01pm UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/14 "2020-04-02T16:01:20Z")

</div>

I’d say it is a limitation of how random effects are implemented in `(g)lmer` and how the `so` basis, in particular, is constructed. In the same way that `te()` smooths don’t work in `gamm4()` and `brm()` because the way those penalties are formed is incompatible with how random effects are represented in those functions.

Writing that got me thinking, however; I’d forgotten that Simon had implemented separate basis types for the hoop (boundary; `bs = 'sf'`) and the film inside the hoop (`bs = 'sw')`. These should have separable penalties which can be fit in `lme()` and `(g)lmer()` and hence `smooth2random()` should work for these and thence (fingers crossed) should also work in _brms_.

Here is a modified version of your reproducible example:

```
library('tibble')
library('dplyr')
library('mgcv')
library('brms')
library('gamm4')

x <- c(0, 0, 1, 1,0)
y <- c(0, 1, 1, 0, 0)
x2 <- runif(100)
y2 <- runif(100)
z <- rnorm(100)

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

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

bndry <- list(list(x = x, y = y))
gamm(z ~ s(x, y, bs = "sf", xt = list(bnd = bndry)) +
         s(x, y, bs = "sw", xt = list(bnd = bndry)),
    data = testdata, method = "REML", knots = knots_test)
# works
gamm4(z ~ s(x, y, bs = "sf", xt = list(bnd = bndry)) +
          s(x, y, bs = "sw", xt = list(bnd = bndry)),
     data = testdata, method = "REML", knots = knots_test)
# works

brm(z ~ s(x, y, bs = "sf", xt = list(bnd = bndry)) +
        s(x, y, bs = "sw", xt = list(bnd = bndry)),
    data = testdata, knots = knots_test)

```

The `brm()` call ends with the original error:

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

```

but I believe @paul.buerkner has this fixed in the GitHub version of _brms_. I haven’t had chance to test that version yet so if someone wanted to try this out I suspect it should work, at least in principle now.

---

<div class="post-metadata">

**Author:** ![StaffanBetner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/staffanbetner/32/11056_2.png) [@StaffanBetner](https://discourse.mc-stan.org/u/StaffanBetner)\
**Post date:** [March 29, 2022, 7:14am UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/15 "2022-03-29T07:14:57Z")

</div>

I can confirm that this approach works for estimation, but not for any post-processing function (e.g. `posterior_predict`, `conditional_effects`). I get an error from `tcrossprod`: “non-conformable arguments”.

Backtrace:

```nohighlight
<error/rlang_error>
Error in `tcrossprod()`:
! non-conformable arguments
---
Backtrace:
     x
  1. +-brm_model %>% predict()
  2. +-stats::predict(.)
  3. \-brms:::predict.brmsfit(.)
  4. +-rstantools::posterior_predict(...)
  5. \-brms:::posterior_predict.brmsprep(...)
  6. \-brms::get_dpar(object, dpar = dp)
  7. +-brms:::predictor(x, i = i, fprep = prep)
  8. \-brms:::predictor.bprepl(x, i = i, fprep = prep)
  9. \-brms:::predictor_sm(prep, i)
 10. \-brms:::.predictor_fe(X = Zs, b = s)
 11. \-base::tcrossprod(b, X)

```

---

<div class="post-metadata">

**Author:** ![StaffanBetner](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/staffanbetner/32/11056_2.png) [@StaffanBetner](https://discourse.mc-stan.org/u/StaffanBetner)\
**Post date:** [March 29, 2022, 7:24am UTC](https://discourse.mc-stan.org/t/spatial-modeling-using-mgcv-soap-film-smoothing/10463/16 "2022-03-29T07:24:58Z")

</div>

It seems like there is something going on when we have both smoothers in the model.

```nohighlight
brm(z ~ s(x, y, bs = "sf", xt = list(bnd = bndry)),
    data = testdata, knots = knots_test) -> brm_model1

brm_model1 %>% conditional_smooths()

```

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/a/a6fe1942945475a49795706d78b618082f4e9610.png)

```nohighlight
brm(z ~ s(x, y, bs = "sw", xt = list(bnd = bndry)),
    data = testdata, knots = knots_test) -> brm_model2

brm_model2 %>% conditional_smooths()

```

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/6/6786e4eeb75312025c93bb71c3940df41e552809.jpeg)