3
$\begingroup$

I have run two GAM models to see if there is a non-linear relationship between my response variable (blood metabolites) and time.

I'm currently accounting for individuals (ID_new) with a random smooth (bs="fs").

 mod1 <- gam(met1 ~ s(Timepoint, k=7) + s(Timepoint, ID_new, k=7, bs="fs"), data=df1, method="REML")

When I use

  summary(mod1)

  Approximate significance of smooth terms:
                      edf  Ref.df     F p-value  
s(Timepoint)         1.195   1.333 1.561  0.2221  
s(Timepoint,ID_new) 18.860 278.000 0.085  0.0828 

From my understanding this shows that there is no effect of time, i.e. no non-linear association between the response variable and time and also no effect when considering how that varies per subject with random smooth (same smoothing parameter, but different shaped smooth). Is this on the right lines and how would I interpret it if the s(Timepoint) was significant but not s(Timepoint, ID_new)?

Secondly, would it make sense to use this model allowing the smoothing parameter for each person to vary;

 mod2 <- gam(met1 ~ s(Timepoint) + s(Timepoint, by="ID_new", k=7) + s(ID_new, bs="re", data=df1, method="REML")

I ran this just to see, and still non significant effect of time as a global trend, but six of the participants had significant p values and wasn't sure how to interpret that (6 subjects have non-linear relationships when considering the metabolite over time); however, that non linearity is not captured when using the bs=fs approach.

Not fishing for significance here as a sidenote, just want to understand the summaries a little better.

Thank you

$\endgroup$

1 Answer 1

2
$\begingroup$

In my opinion you are correct to be cautious about simply using $p$-values to guide interpretation. {mgcv} has amazing functionality and the significance tests are useful, but to really make useful inferences from hierarchical GAMs requires additional, prediction-based investigations. Here I give some suggestions to help guide interpretation of these models, making use of the powerful {marginaleffects} package. You can also do similar kinds of investiations using the {gratia} package. First I simulate a "shared" nonlinear function using a squared exponential Gaussian Process:

library(dplyr)
library(marginaleffects)
library(ggplot2); theme_set(theme_classic())
library(mgcv)
set.seed(55)

# A function to simulate from a squared exponential Gaussian Process
sim_gp = function(N, c, alpha, rho){
  Sigma <- alpha ^ 2 *
    exp(-0.5 * ((outer(1:N, 1:N, "-") / rho) ^ 2)) +
    diag(1e-9, N)
  c + mvnfast::rmvn(1,
                    mu = rep(0, N),
                    sigma = Sigma)[1,]
}

# Simulate the "shared" function
N <- 50
x <- 1:N
N_groups <- 5
shared <- sim_gp(N, c = 5, alpha = 0.25, rho = 15)
plot(shared, type = 'l', lwd = 2, ylab = 'F(x)', xlab = 'x')

Next I simulate "deviation" functions for the groups, where each group's final function is a sum of the "shared" and "deviation" function

# Simulate "deviation" functions for each group and add them
# to the "shared" function. For each group, allow the deviation
# to be chosen randomly
devs <- matrix(NA, nrow = N_groups, ncol = N)
for(group in 1:N_groups){
  devs[group, ] <- shared + sim_gp(N, 
                                   c = 0, 
                                   alpha = runif(1, 0.15, 0.75), 
                                   rho = runif(1, 6, N / 2))
}

# Plot the group-level functions
cols <- c('darkred', 'darkblue', 'black',
          'darkorange', 'darkgreen', 'grey60')
plot(devs[1,], type = 'l', ylim = range(devs),
     col = cols[1], lwd = 2, ylab = 'F(x)', xlab = 'x')
for(group in 2:N_groups){
  lines(devs[group, ], col = cols[group], lwd = 2)
}

I then simulate Gaussian noise and wrangle the data into a data.frame before fitting a hierarchical GAM

dat <- do.call(rbind, lapply(1:N_groups, function(group){
  data.frame(x = x,
             grp = paste0('gr_', group),
             y = rnorm(N, mean = devs[group, ],
                       sd = 1))
})) %>%
  dplyr::mutate(grp = as.factor(grp))
dplyr::glimpse(dat)
#> Rows: 250
#> Columns: 3
#> $ x   <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,…
#> $ grp <fct> gr_1, gr_1, gr_1, gr_1, gr_1, gr_1, gr_1, gr_1, gr_1, gr_1, gr_1, …
#> $ y   <dbl> 5.665085, 5.628093, 5.429635, 4.664341, 3.444215, 5.157122, 4.8685…

# Fit a hierarchical GAM to the data
mod <- gam(y ~ s(x) + s(x, grp, bs = 'fs'),
           data = dat)
#> Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has repeated
#> 1-d smooths of same variable.

# Summary and default mgcv plot
summary(mod)
#> 
#> Family: gaussian 
#> Link function: identity 
#> 
#> Formula:
#> y ~ s(x) + s(x, grp, bs = "fs")
#> 
#> Parametric coefficients:
#>             Estimate Std. Error t value Pr(>|t|)    
#> (Intercept)   4.9490     0.0789   62.72   <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Approximate significance of smooth terms:
#>            edf Ref.df     F  p-value    
#> s(x)      1.00      1 0.223    0.637    
#> s(x,grp) 12.71     48 0.948 4.07e-06 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> R-sq.(adj) =   0.15   Deviance explained = 19.7%
#> GCV = 1.1019  Scale est. = 1.0371    n = 250
plot(mod, pages = 1, seWithMean = TRUE)

The plots from {mgcv} are fine for a quick glance, but we can't get a sense of what each group's estimated function is by simply looking at these (because an individual group's function is composed of the "shared" and "deviation" functions. {marginaleffects} has the predictions() family of functions, which will make use of predict.gam() to automate this for us. For example, here I show estimated smooths on the response scale for each group. These take into account all of the effects in the model

# Plot expectations (response scale) using marginaleffects
plot_predictions(mod, by = c('x', 'grp')) +
  ylab('F(x)')

plot_predictions(mod, by = c('x', 'grp', 'grp')) +
  ylab('F(x)')

For further investigations, we can for example ask How do these functions differ among groups?. The comparisons() family of functions can answer this question:

# How do these functions differ? Use pairwise contrasts to investigate
plot_comparisons(mod, newdata = datagrid(x = seq(1, N, by = 0.25),
                                         grp = unique),
                 variables = list(grp = 'all'),
                 by = "x") +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(y = "Estimated difference in F(x)",
       title = "Difference smooths",
       subtitle = "Pairwise contrasts between groups")

We can also target the 1st derivatives of these functions to understand where they are changing the most. This is done using the slopes() family of functions:

# Where are these functions changing the most? Calculate first derivatives
plot_slopes(mod, variables = 'x', by = c('x', 'grp', 'grp'),
            newdata = datagrid(x = seq(1, N, by = 0.25),
                               grp = unique)) +
  ylab("F'(x)") +
  geom_hline(yintercept = 0, linetype = 'dashed')

You may be interested in finding values along the predictor $x$ at which the rate of change for a group differed 'significantly' from zero. The hypotheses() family of functions can do this:

# Do these slopes differ 'significantly' from zero?
hypotheses(slopes(mod,
                  newdata = datagrid(x = seq(1, N, by = 0.25),
                                     grp = unique),
                  variables = "x",
                  by = c("x", "grp")))
#> 
#>  Term    Contrast  x  grp  Estimate Std. Error      z Pr(>|z|)   S   2.5 %
#>     x mean(dY/dX)  1 gr_1  0.017984     0.0596  0.302    0.763 0.4 -0.0989
#>     x mean(dY/dX)  1 gr_2 -0.012400     0.0596 -0.208    0.835 0.3 -0.1292
#>     x mean(dY/dX)  1 gr_3 -0.050307     0.0596 -0.844    0.399 1.3 -0.1671
#>     x mean(dY/dX)  1 gr_4 -0.057161     0.0597 -0.957    0.339 1.6 -0.1742
#>     x mean(dY/dX)  1 gr_5 -0.063904     0.0596 -1.073    0.283 1.8 -0.1807
#>  97.5 %
#>  0.1348
#>  0.1044
#>  0.0665
#>  0.0599
#>  0.0529
#> --- 975 rows omitted. See ?print.marginaleffects --- 
#>     x mean(dY/dX) 50 gr_1 -0.037518     0.0596 -0.630    0.529 0.9 -0.1543
#>     x mean(dY/dX) 50 gr_2  0.070460     0.0595  1.183    0.237 2.1 -0.0462
#>     x mean(dY/dX) 50 gr_3  0.038296     0.0595  0.643    0.520 0.9 -0.0784
#>     x mean(dY/dX) 50 gr_4  0.000656     0.0596  0.011    0.991 0.0 -0.1161
#>     x mean(dY/dX) 50 gr_5  0.043577     0.0597  0.730    0.465 1.1 -0.0734
#>  97.5 %
#>  0.0793
#>  0.1872
#>  0.1550
#>  0.1174
#>  0.1606
#> Columns: rowid, term, contrast, x, grp, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, predicted_lo, predicted_hi, predicted 
#> Type:  response

Finally, you could ask how estimated 1st derivatives differ among groups (i.e. where along the covariate $x$ was group 2 growing substantially faster than group 1, for example):

# How do these derivatives differ? Use pairwise contrasts again to investigate
plot_comparisons(mod, newdata = datagrid(x = seq(1, N, by = 0.25),
                                         grp = unique),
                 variables = list(grp = 'all'),
                 by = "x",
                 comparison = "dydx") +
  geom_hline(yintercept = 0, linetype = 'dashed') +
  labs(y = "Estimated difference in F'(x)",
       title = "Difference derivatives",
       subtitle = "Pairwise contrasts between groups")

Created on 2024-04-16 with reprex v2.1.0

$\endgroup$
3
  • $\begingroup$ Thank you very much for this, very helpful! Was wondering how you would interpret the significant p value in the summary for s(x, grp) - would this indicate that you cannot put a horizontal line through the 95 % confidence intervals for any of the groups? And even if it wasn't significant you would have proceeded with some of your exploratory analyses above? $\endgroup$
    – Adam9
    Commented Apr 16 at 21:33
  • $\begingroup$ I have 2000 response variables, and 28 different subjects, I was wondering what your thoughts would be on running a GI model, as I don't believe each individual needs to have the same smoothing parameter, and the summary from a gam(response ~ s(Time) + s((Time, by=subject) + s(subject, bs=”re”), shows a breakdown for the global trend and also each individual, which I find easier to understand rather than s(time, subject) from a GS model (bs="fs"). $\endgroup$
    – Adam9
    Commented Apr 16 at 21:43
  • $\begingroup$ 1) I would always proceed with interpretation, regardless of p-values. You can see that I simulated a true, nonlinear "global" function but the p-value for this effect is non-significant. Does that mean it isn't important? Probably not. 2) You can run this model as well, the setup of the smooths are slightly different because the deviation smooths are no longer sharing a smoothing parameter, but none of the calls to {marginaleffects} will change at all. So you can try multiple models and make useful comparisons $\endgroup$ Commented Apr 16 at 21:45

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.