Specifying generalized additive models in R (using the mgcv package)

I'm currently trying to figure out how to specify a generalized additive model in four different ways in order to investigate how my outcome variable changes surrounding different events. So what I want to model is

1. fixed effect of year from event only, with no person-specific model parameters (i.e., each person is assigned an identical curve),
2. multilevel with person-specific varying intercepts (i.e., the shape of the curve is the same for each person, but an intercept allows for different overall levels of my outcome variable),
3. multilevel with person-specific varying intercepts and slopes (i.e., in addition to the intercept, a slope parameter allows for different linear change in my outcome variable over time), and
4. multilevel with person-specific varying curves. In a next step I would like to compare these models on the basis of predictive fit and deviance explained.

So far, I got this, which (if I got it right?) does not include any specification. It just represents the estimated population-level trajectories of my outcome variable surrounding the event (as "year" was event-centered):

gam_event1 <- gam(outcome ~ s(year), data = event1.sample, method ="REML")

This is my output plotted: But how can I specify these four different models described above? Unfortunately looking at the mgcv documentation wasn't a great help to me.

I'm relatively new to R and even newer to GAM so please excuse me if this is a really stupid question! Your help is highly appreciated!

You can add the following terms to the model:

1. The smooth, year effect is as you have it in the example code you showed:

outcome ~ s(year)

2. A random intercept (assuming person is a factor coding for the individual subjects)

outcome ~ s(year) + s(person, bs = 're')

3. A population year effect, random intercept, and random slope (of year is what I understand you want)

outcome ~ year + s(person, bs = 're') + s(year, person, bs = 're')

4. Random smooths can be generated in several ways, but the simplest, which assumes that all smooths share a common smoothness parameter (which doesn't mean they share the same shape!) is via the fs basis

outcome ~ s(year, person, bs = 'fs')

If you want this model to include a population-level effect, then something like

outcome ~ s(year) + s(year, person, bs = 'fs')

is required, and you may to use select = TRUE or add m = 1 to the fs smoother to penalise any deviation from a flat function (i.e. from the population level curve).

You don't want to fit all these models separately and try to compare them with a LRT (via the anova method) as, for example comparing model 1 with model 2 will result in unreliable p-values that typically are anti-conservative. You are best served fitting model 2 and using the summary() method for a test of the random effect.

Likewise, you can fit model 4 and just look at the summary and the fitted curves, as this model includes the random intercept model and the random intercept and slope model as special cases.

These are not the only ways to specify these models in mgcv; options include using mgcv::gamm() and gamm4::gamm4(), plus there are some other tricks; see paraPen argument to gam() and ?smooth.construct.fs.smooth.spec for some ideas.

• Thank you very much! This was great help :) I used the summary() function and additionally the AIC() function to compare models and it seems to work. Your additional hints to gamm, gamm4 and paraPen are also very helpful, I'll check it out as well. – Marie B. Mar 22 '18 at 9:40