Does brms have issues with (perfect) linear dependencies between (smooth) covariates?

Interfaces brms
, , ,
1

Hello (my very first question in this forum :)),

I have a data-set about an experiment, where the influence of time (t) on y is not linear (its about growth). For every individual, the data is observed with and without Treatment.

The model I would like to fit is:

y ~ Treatment + (1|individual) + s(t, by = Treatment) + s(t, by = individual) [1]

as I intend to capture the individual-specific deviations from the population parameter by ‘random’-effects. Within the frequentist (GAM) framework, columns of the design matrix X need to be removed:

  • one level of Treatment from s(t, by = Treatment)
  • one level of individual from s(t, by = individual)

I wonder how the brms framework deals with this as the summary() reports estimates for s(t) for every level of Treatment (X,Y) and individual (A,B,C, and D) as this (reducing the output to the relevant part):

Smooth Terms: 

sds(stTreatmentY)  
sds(stTreatmentX) 
sds(stindividualA) 
sds(stindividualB)
sds(stindividualC)
sds(stindividualD)

Group-Level Effects: 
~individual (Number of levels: 4)    

Population-Level Effects: 

Intercept
TreatmentX
st:TreatmentY
st:TreatmentX
st:individualA
st:individualB
st:individualC
st:individualD

I fear that my model equation (as it is) makes not sense and my guess is that I am probably not the first who is looking to capture a Treatment-effect on population-level using splines with random-effects for the individual deviations?

  • Operating System: macOS Big Sur 11.4
  • brms Version: 2.15.0

Cheers

1 Like
2

Hi,
sorry for not getting to you earlier, your question is relevant and well written. If I understand you

Unfortunately, you cannot use splines as “random” or “varying” effects, because the penalized splines from mgcv which brms relies on are already implemented as a varying effect (see e.g. Random effects and penalized splines are the same thing - Higher Order Functions for background).

So if you want a smooth term that is at the same time a “random” effect, you need to work directly with non-penalized splines. That’s IMHO usually not a big deal, since you rarely want to allow the per-individual spline to be very complex (you can’t learn them reliably anyway), so a low-degree-of-freedom spline should work just well. So something like:

y ~ Treatment + s(t, by = Treatment) + (1 + bs(t, df = 2)|individual)

Could be sensible. (this requires the splines package. You may also want to use ns but I would expect the results to be quite comparable)

I’ll tag @paul.buerkner if he knows about other solutions.

EDIT: I realized probably forgot to attend to the core of your question - if you have s(t, by = Treatment), it is not a problem that you get estimates separately for both groups. Intuitively: you fit a separate splines for both groups which is well-defined and the splines from s() do not contain intercept which is the only term that could interact with the the other terms already in the model.

Note that in GAM models you do not remove some of the terms “just because”. You remove them, because they are shared and implicitly already in the model. So e.g. if you have two binary predictors A and B and have y ~ A * B you are actually fitting 2x2 separate means. The four groups are called “Intercept”, “Intercept + ATrue”, “Intercept + BTrue” and “Intercept + ATrue + BTrue + ATrue:BTrue”. So it is not that you would treat A as having just one level. What happens is that one of the levels is already present in the model thanks to the intercept, so you do not add it once more. The per-group splines are however coded so that the terms are not shared (akin to what you get with y ~ 0 + A). Hope the previous paragraph did make at least some sense :-)

Best of luck with your model!

1 Like
3

The solution @martinmodrak provided looks very sensible to me too.

1 Like
4

Dear martinmodrak,

thank you so much for the reference - that has helped me a lot. Your suggested solution makes a lot of sense.

About the intercept issue: agreed - if every intercept of each of the splines is removed it works.

Thanks again

1 Like
5

(Sorry for coming to this late, but for others also coming across this thread I thought I’d suggest an alternative.)

The “random” equivalent of a factor-smooth interaction using penalized splines could be fitted using mgcv’s “fs” basis. This basis includes the random intercept for subjects so you won’t need the (1|individual) term, but the wiggly parts of the basis will be penalized, just like the s(t, by = Treatment) term(s). All the subject-specific smooths will use the same value of the smoothing parameters (there are three smoothing parameters for the fs basis: one each for the shrinkage/penalization of the random intercepts, random slopes (the linear function in the basis, not usually penalised as it is in the null space of the wiggliness penalty), and the wigglyness.

Your model would the be:

y ~ Treatment + s(t, by = Treatment) + s(t, individual, bs = "fs")

As for efficiency, you’ll lose the ability to specify the random intercepts using brms, and instead the random intercept and slopes will use their random-effect-as-smooth equivalents, which require full penalty matrices. The limiting factor here will be the number of subjects.

1 Like