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

*

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


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


*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")
Follow up questions
Thanks a lot for your response. Three questions

*

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


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


*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).
 A: 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.
