# Model building in generalized additive mixed models (GAMMs)

I am currently try to set up a GAMM model that is based on data that my team and I have gathered. The data contains self-reported performance of 30 startups during 12 weeks (weekly measures) of the first phase of the Covid19 crisis. The goal of the project are to analyze

1. whether there is a linear or nonlinear overall trend in the performance (e.g., decrease, up-shaped or something more fancier)

2. whether there are random effects, e.g., different forms of change across the 30 teams

3. whether these differences can be explained by measured covariates (e.g., degree of team cohesion).

I have some knowledge in GA(M)Ms but have problems how to build up the model.

What I did until now are the following steps (done with the gam() function in mgcv)

a) I started with the time-only model (y ~ s(time) ). In one version, I tested a linear time effect, in a second, a smooth term. The AIC supported the linear model (hence research question 1 is answered: There is a linear increase in the performance

b) Then I added random intercepts s(ID, bs="re", k=30) [ID = ID of the team)

c) Then random slopes (s(tid, time, bs="re")

d) then random smooths (s(time, tid, bs="fs", xt="cr", m=1, k=5)) [where the random intercepts and slopes were excluded) The AIC supported the random smooth model.

e) Then starts my problem: When I add the predictors as a ti(time,X) product (having time and X as main effects in the model), none of the tested X's are significant. If I eliminate the random smooth component, this changes. I found no real advises in the GAMMS sources I read so far, hence any advises would be helpful (either references or direct recommendations how to set up the steps. To boil it down to one question: Do I have to keep the random smooth part, when the interaction between time and X are entered?

The final code for the model is this

model <- gam(Y ~ time + s(X, bs="cr", k=5) + ti(time, X) + s(time, ID, bs="fs", xt="cr", m=1, k=5), data=TeamData, method="REML")

Thanks a lot for your response. Three questions

1. I learned (I guess) that I have to include random effects from the beginning, right?

2. When using your approach of y ~ time + s(time, m=c(2,0)... together with all random effects, I got an error "Indexing outside the boundaries". Would it make sense to estimate two models — one with y ~ time vs. y ~ s(time) (plus the random smooth component) and compare both with the AIC? If I do that, the AICs are almost identical and together with the edf of time in the nonlinear model of 1.01, I would conclude that the trend indeed is linear. Your idea is more solid, I agree, though.

3. With regard to the X covariate, you misunderstood my goal: X is a stable team characteristic (i.e., a time-invariant predictor). I would include that after the linear time + random smooth model as a tensor product. You already solved my main problem when and whether the random effects are in the model (and stay there).

• I edited you question so that it can be seen what is code and what is not code, i.e. text. You changed that back. I would recommend formatting your question so that this is clear. Sep 5, 2020 at 17:47
• Oywey, sorry, this was not intentional. I edited a few typos and stylistic issues and did not see that you had changed that.... Sep 5, 2020 at 20:43
• No worries! If you like you can try it yourself. Code blocks can be formatted by highlighting the code block (i.e. the model specification in your example) and then click the code sample button on top of the edit menu (or by indenting every line where this is code with four spaces). Code embedded in text such as this example can be done by placing it in back ticks. You can always look in the preview and see how those changes affect the post. See also here: meta.stackexchange.com/questions/22186/… Sep 5, 2020 at 21:28
• I'm curious why you decided to use GAMMs as opposed to methods specifically built for longitudinal data. Can you explain a bit more about your decision to use GAMMs instead of something like, say, and ARIMA model? Sep 6, 2020 at 6:42
• The answer is simply that the GAMM allows me to answer all my questions. I spent some time with time series analysis, but as a psychologist experienced with multilevel modeling I find GA(M)Ms much more relatable--especially when it comes to analyzing predictors. Do I miss anything what a GA(M)M cannot do what a ARIMA model can do? Sep 6, 2020 at 7:02

Your first model, with y ~ s(time) isn't valid as it failed to account for a lot of the structure in the data. Simply identifying that on average over all the startups there is no non-linear effect isn't the same as saying that the individual effects are non-linear either.

If you want to decompose a linear and non linear effect, you can fit models of the form

~ time + s(time, m = c(2,0))


where the m specification keeps the second order penalty but removes the null space (0). This latter step means that the linear term is excluded from the basis, which is what you want if you already have a linear parametric time effect.

But I think this might be too complex to include whilst trying to identify the other structures. As the linear fit is included in s(time), I would just include it through the smooth, figure out the other structures and then once you have the rest of the model as you want it, you can replace s(time) with time + s(time, m = c(2,0)) for the explicit test of linear vs non-linear change in time.

For the first hypothesis, I would fit this model:

y ~ time + s(time, m = c(2,0)) + s(tid, bs = 're') + s(tid, time, bs = 're') +
s(time, tis, bs = 'fs', m = c(2,0))


The only overlap I believe here will be the random intercepts, but I think that is OK.

As that accounts for the clustering at the tid level

As for "e)", you are testing for temporally varying effects of X; you could add them to the model above and see what happens. If you leave off the last smooth, the fs term, and include the ti(time, X) term, that's fine, but I would keep the two re terms or perhaps only just the s(tid, bs = 're'), as that reflects group structure that you don't want being explained by the global time-vary effects of X.

I would then check the residuals and plot them against tid and time to see if there is unmodelled temporal structure which might invalidate the assumptions of the model (and thus the tests of the terms in the model).

With such limited data, it may not be possible to identify time varying effects of covariates and account for the clustering at the startup level.

It also sounds like you are testing X1 first, then removing it and adding X2, testing that and so on. If so I don't think that is a valid way to proceed.