I have repeated measurements of individuals, like this

y  |  id  | time  |  V1   | V2
100    1       1     23.2   0.8
88     1       2     22.6   0.9
98     2       1     10.6   1.1
83     2       2     11.1   1.3

y is a continous outcome variable, id is the patient id, time is the time point of oberservation and V1 and V2 are covariates.

In this data example we have 2 patients that have 2 observations each, one at time 1 and another one at time 2.

My real data sets has hundreds of patients (ids) and 2-5 observations (time) for each patient.

I now know that the effect of V1 on y is non-linear and what to model this with a GAM.

In the mgcv package there is a function called gamm(). In my example, I use it like this:

m <- gamm(y ~ s(V1) + V2 + time ,family=gaussian,
          data=dat,random= ~(time|id))

Is this correct? Does this model integrate the fact that V1 can change over time for each individual?

  • 1
    $\begingroup$ I suggest to implement random effects within gam as a smoother instead, i.e., bs = "re" or even bs = "fs". Such a model is easier to handle since it doesn't consist of two models internally as gamm does. Anyway, what does "V1 can change over time for each individual" mean exactly? Is V1 perhaps strongly collinear with time? $\endgroup$
    – Roland
    Feb 11, 2019 at 13:16
  • $\begingroup$ @Roland how would I formulate this using gam ? Do I still get parameter estimates for V2 ? V1 is a variable like BMI that may be strongly correlated with its observation at the previous time point. $\endgroup$
    – spore234
    Feb 11, 2019 at 13:20
  • $\begingroup$ Can you explain what you mean by "may be strongly correlated with its observation at the previous time point"? This suggests a lagged effect, but you include V1 as a predictor of y contemporaneously (unless you lagged it yourself?). If the effect of V1 is intended to vary with time, then your model isn't accounting for that. Some further explanation would be helpful. $\endgroup$ Feb 12, 2019 at 3:00
  • $\begingroup$ @GavinSimpson let's assume V1 is BMI and I want to integrate individual BMI changes over time into the model $\endgroup$
    – spore234
    Feb 14, 2019 at 7:59
  • 1
    $\begingroup$ @spore234 the main effects and interaction would be te(BMI, time, k = c(a,b)) with suitable values for a and b if needed. This is the smooth equivalent of BMI * time. If you want to decompose this, then: s(BMI) + s(time) + ti(BMI, time). If you want the linear effect of BMI to vary smoothly with time, i.e. a varying coefficient model, then s(time, by = BMI) assuming BMI is coded as a numeric variable. It really depends on what you want to represent in the model. $\endgroup$ Feb 15, 2019 at 19:55

1 Answer 1


What @Roland is getting at is to use a random spline basis for the time by id random part. So your model would become:

m <- gam(y ~ s(V1) + V2 + s(time, id, bs = 'fs'),
         family=gaussian, data=dat, method = "REML")

This model says that the effect of time is smooth and varies by id, with a separate smooth being estimated for each id but each smooth is assumed to have the same wiggliness (a single smoothness penalty is estimated for all the time smoothers) but can differ in shape.

To estimate a separate "global" effect the model could be

m <- gam(y ~ s(V1) + V2 + s(time) + s(time, id, bs = 'fs'),
         family=gaussian, data=dat, method = "REML")

If you want similar models but where each smooth can have different wiggliness as well as shape, then the by smoothers can be used:

## without a "global" effect
m <- gam(y ~ s(V1) + V2 + s(id, bs = 're') + s(time, by = id),
         family=gaussian, data=dat, method = "REML")

## with a "global" effect
m <- gam(y ~ s(V1) + V2 +
         s(id, bs = 're') + s(time) + s(time, by = id, m = 1),
         family=gaussian, data=dat, method = "REML")

The m=1 means that the smoother uses a penalty on the squared first derivative, which penalises departure from a flat function of no effect. As this is on the subject specific smooths, the model is penalising deviations from the "global" smooth.

Some colleagues and I have described these models in some detail in a paper submitted to PeerJ, which is available as a preprint. A new version in response to reviewers comments should be up in a few days (we've submitted it to the journal).

  • 1
    $\begingroup$ Nice manuscript. Looks like a useful reference for a paper's methods section. $\endgroup$
    – Roland
    Feb 12, 2019 at 6:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.