Global smoother or global smoother plus group-level smoothers?

Question

I am trying to analyze a dataset that contains measurements of a physiological trait that has been repeatedly measured from before sunrise till after sunset on different individuals within different species. In this case, one replicate is one day of repeated measurements from one individual within a species. There are at least four replicates within each species that approximately covers a whole day combined (some species have not been measured in the exact same timeframe). There are several groups in the dataset as there are species from different habitats (two) and different types (three) in the dataset (see at picture for an example of the dataframe []). In addition, a single species is only found in one of the groups. For example, one species does not appear in both habitats and can only be grouped within one type (the "type" is related to the type of species).

I wish to answer a couple of questions:

1. Is the measured trait under a temporal variation across all groups (i.e. across both habitats)?
2. Is the measured trait under a temporal variation within each habitat and type?

To clarify, temporal variation means that the measured trait either increases or decreases significantly throughout a day (for example increase from morning till noon and then decrease from noon till evening). To answer these questions, I thought GAMs with automatic smoothing selection would be appropriate. The "common" approach in my field appears to correspond to just fitting a GAM to the trait with time of the day as the predictor variable without group-level smoothers. This is because global effects are important, but species-specific differences are of less importance (not really of any interest). In this perspective I could fit this simple GAM across the data from both habitats to answer question 1, and then separately to the data from each habitat and type to answer question 2:

Model1 <- gam(Trait ~ s(Timeofday),
family = gaussian,
data = my_data, method = "REML")

However, I’m concerned that this approach is somewhat erroneous because in some (few) cases the diurnal variation of the measured trait differs quite substantially between species. After reading Gavin Simpson’s paper about hierarchical GAMs (https://doi.org/10.7717/peerj.6876) I came up with a set of alternative models, but I’m not really sure as to how I decide which is the appropriate:

Model_2 <- gam(Trait ~ s(Timeofday) + s(Species, bs = "re"),
family = gaussian,
data = my_data, method = "REML")

Model_3 <- gam(Trait ~ s(Timeofday) + s(Timeofday, Species, bs = "fs", m = 1),
family = gaussian,
data = my_data, method = "REML")

Model_4 <- gam(Trait ~ s(Timeofday) + s(Timeofday, by = Species, m = 1) + s(Species, bs = "re"), family = gaussian, data = my_data, method = "REML")

Moving from model 1 to model 2 does almost nothing to the AIC values. Model 3 captures significantly more deviance than model 1 and 2, accompanied by a substantially lower AIC value as well. When looking at how model 3 fits each species, the model appears to do a poor job because some species differ in how wiggly they are. Model 4 appears to solve this problem (as expected), and the AIC value is indeed lower than the AIC value of model 3. However, I’m concerned that model 4 is overfitting, and I really have no basis for arguing that species should differ in how wiggly they are. Some species are quite wiggly, and this is most likely a result of few samples rather than the actual trend of the species. I tried to fit this model instead:

Model_5 <- gam(Trait ~ s(Timeofday) + Species + s(Timeofday, Species, bs = "fs", m = 1),
family = gaussian,
data = my_data, method = "REML")

Model 5 appears to fit the data well with a much lower AIC value compared to model 1 and 2 (but similar to model 4). I’m completely new to GAMs, so I’m a little confused as to what model 5 actually does. Regardless, I would be verry thankful if someone could explain how I might decide which model is appropriate to answer these questions.

Thanks.

Answer 1

It may be easier to start with linear modeling before using additive modeling if you are not familiar with this technique. Indeed you could add a quadratic term to a linear model (Timeofday + Timeofday^2) to model a non-linear relation between Trait and Timeofday.

Moreover I thought that a smoother could only be added to a quantitative variable, so I am confused about the s(Species) term.

If indeed you are not interested in differences between species but you see very different patterns between species, perhaps it would be better to use mixed models with Species as random variable or fit one model per species.