Random effects in GAM

Question

I'm modelling some biological function, outcome, in patients over 8 hours. Over time, I measure two additional covariates x and w. Shown below is a sample of the data. Only show the first three measurements per patient.

I'm interested in estimating the effects of x and w and their interaction on outcome. I plan to use a gam to take care of the time effect.

My question is: do I need to model the effects of x and w as a random effect this they change in time and I have multiple observations of these covariates per patient? If so, what is the appropriate way to specify this model with a gam?

Currently, I do something like

gam(outcome ~ s(time, by=patient, bs = 'gp') + x*w, data= data)

I can provide more details as necessary, including more data.

---

Data

tibble::tribble(
     ~patient,     ~time,      ~x,         ~w,          ~outcome,
            1,     3.245,     87.738,    38.571,         19.006,
            1,     8.245,     80.104,    39.958,         19.441,
            1,    13.245,     76.206,    40.912,         19.046,
            2,     5.335,     80.731,    56.822,         32.689,
            2,    10.335,     83.174,     52.52,         32.654,
            2,    15.335,     81.507,    52.374,         32.589,
            3,     4.965,     68.584,    111.69,         37.621,
            3,     9.965,     68.544,   111.751,         37.762,
            3,    14.965,     68.979,   112.118,         37.451,
            4,     20.38,     108.19,    78.406,         44.792,
            4,     25.38,    111.694,    81.583,         49.085,
            4,     30.38,    111.312,    80.914,         44.741,
            5,     6.295,     72.692,    52.059,         42.737,
            5,    11.295,     75.155,     49.55,         42.788,
            5,    23.533,     78.084,    49.131,         41.934,
            6,     7.075,     79.788,    69.537,         25.219,
            6,    12.075,     79.652,    68.692,         25.427,
            6,    17.075,     78.746,    69.277,         25.635,
            7,     3.335,      80.44,    34.839,         61.974,
            7,     8.335,     79.318,    34.325,         56.651,
            7,    13.335,     79.624,    34.346,          53.32,
            8,      9.92,     96.634,    64.871,         44.481,
            8,     14.92,    100.979,    67.194,         43.645,
            8,     19.92,     97.952,    67.456,         43.934,
            9,    74.635,      98.08,    83.723,         39.156,
            9,    79.635,     97.779,    85.049,         38.958,
            9,    84.635,     96.881,    88.736,         38.911,
           10,    70.155,     73.583,    53.296,         28.519,
           10,    75.155,     73.332,    52.935,         28.705,
           10,    80.155,     73.174,    51.781,         29.405,
           11,      3.17,     80.425,    45.579,         27.271,
           11,      8.17,     86.442,    39.839,         27.095,
           11,     13.17,     80.377,    42.182,         28.844,
           12,      3.23,     53.399,    75.199,         26.986,
           12,      8.23,     64.986,    63.395,         17.803,
           12,     13.23,     65.306,    68.421,         15.915
     )

Answer 1

I'm reasonably sure you don't want to use the gp basis for this; when used with the defaults, mgcv fits a smooth with a Gaussian Process basis comprising correlation functions as the basis functions where the effective range parameter is set to essentially be the largest observable time difference between any two observations. If you must use the gp basis then you probably want to profile over the effective range parameter, fitting a new GAM for a grid of values over the set of values for the effective range parameter you want to optimise for.

Anyway, you really don't need to do this.

You are fitting one option for group-level splines (the by variable smoother) but you could take an approach that is slightly more a "random effect" than this by using the fs basis. The main difference is that in the fs basis, all the splines (one per patient) have the same smoothness parameter (so the same wiggliness) though the don not need to have the same shape. The by variant estimates a separate smoothness parameter per patient, which essentially means you are fitting an entirely separate smoother for each patient. Both are acceptable, there's just a lot more parameters (length(levels(patient)) - 1 more) to estimate for the by variant.

If you think the effect of x (or w) changes over time then you should include terms in the model for that. There are a number of ways in which this can be done, including

1. a varying coefficient model-like approach, wherein the effect of x is linear but varying smoothly with some other covariate, in your case time. Think of this as a linear regression of y ~ x where the slope of the linear regression line changes smoothly over time.

   For this kind of term, use s(time, by = x) and so on in place of the x and w parametric terms.

   If you also want the interaction to change smoothly, then I think you need to create a new variable xw which is the product of x and w, and then use xw as the variables passed to by.

2. A full smooth effect of x on outcome that smoothly interacts with time. These are produced via tensor products and the te() function. A bivariate interaction would be

   + te(x, time)
   

   and the full blown three-way would be

   + te(x, w, time)
   

   but to smooth in high dimensions like this you really need a lot of data.

To give more specifics you'll need to explain more about how you view the effect of x, w and x:w changing over time.