Additive model with random smooth for participant

Interfaces brms
, , ,
1

Hi All,
I am trying to model some pupillometry data that I collected from one of my studies.
I have already modelled the data using mgcv and I woudl like to convert the model to brms.
I am new to both brms and mgcv so please forgive me if this is a stupid question.

From what I could read I think the model specification should be the same but I am not totally sure about the random effect part.
What I am trying to do is model the pupil data to have:

  • Interaction between Event (variable of interest; factorial with 2 levels) and Trial number
  • Smooth term on time (ms) for each Event (factorial with 2 levels)
  • Random smooth for each Subject

the model I run on mgcv is:

m1 <- bam(mean_pupil ~ Events*Trials  
          
          + s(Time,  by = Events, k=20)

          + s(Time, Subject, bs = 'fs', m=1),
          data = prediction, discrete = TRUE,
          family = scat)

my interpretation on how I should port the model the brms is simply:

mod2 = brm(mean_pupil ~  Events*Trials  

          + s(Time,  by = Events, k=20)
          
          + s(Time, Subject, bs = 'fs', m=1) , data =  prediction,
    family = student, seed = 17, cores= 4,
    iter = 2000, warmup = 1000, thin = 10,
    control = list(adapt_delta = 0.99))

Is the model conversion correct? Is this the way to specify a random smooth for each participant also in brms?

Also I am wondering if there is a way to speed up the model since when I try to run the model it takes around 16h (the dataframe has just 29768 obs)?

In case someone is interested I built the models thanks to this two helpful tutorials:

Thank you a lot !

Operating System : windows 10
brms Version : 2.18.0
2

I can’t offer specific advice on your model, but a couple thoughts about run time: first, thinning is generally not needed in stan (the thin=10 part of your brm code), and unless you are hitting divergent transitions otherwise, you can probably reduce or get rid of the adapt_delta value. The former will allow you to run fewer iterations and the latter will speed up each iteration. If you have access to the appropriate hardware, multithreading will also speed things up.

1 Like
3

@TommasoGhilardi hello, this all looks reasonable to me. One feature that I think you’re missing is a group-level intercept by subject i.e. +(1|Subject) , which allows the subjects’ functions to be shifted vertically. In my experience the shared-wiggliness random smooths tend to have difficulty capturing that alone. Of course you can also use additional by = smooths if you wanted more flexibility between subjects, but it could be a lot more expensive.

I’m not surprised that this model takes a fair bit of computation time, HGAM tend to be quite slow even with smaller datasets and this is pretty large. As @Vanessa_Brown said, this will be more efficient if you drop the thinning and keep adapt_delta lower, but in my experience HGAM tends to have divergences without adapt_delta kept quite high. Depending on your CPU, you could run more chains and reduce their length, presuming that the warmup is sufficient.

The linked paper is great, I’ll hang on to that one. I really appreciated this one: Hierarchical generalized additive models in ecology: an introduction with mgcv [PeerJ].

2 Likes
4

@Vanessa_Brown Thank you for your suggestions!! I was indeed hitting some divergent transitions that’s why I modified the adapt_delta. I now have included multithreading and this seems to have sped up the process even if not too much.

@AWoodward Thank you so much for your suggestions too!! I now included the random intercept, I just have to wait for the model to run… indeed very slow:)
I am very happy that you found the paper interesting! The one you linked is great. Very clear, it will help me a lot!