How to appropriately model complex structured, multivariate data [closed]

Question

This is an extension of this question I posted earlier. The general data I introduced there is the same.

The structure/nature of my data is:

• y ... my response variable of interest
• x ... my main predictor variable and the only variable of real interest for me
• Type ... different categorical levels for x
• ID ... my data was collected at different sites resulting in possibly different baselines, but there also may be an interaction with Type per site and/or a different relationship for y~x per site; sites are independent from each other
• Covariate1, Covariate2 ... additional variables that may impact y and that may have an interaction with Type and/or ID

Challenges:

• the complexity of the data
• possible non-overlap of x ranges per Type/ID

My main questions I'm trying to answer with my model:

• Is there an overall effect of x on y?
• Does the effect of x on y differ by Type and if so, how?
• Does the effect of x on y differ by ID and if so, how?

My questions I need help with:

1. Is it even statistically valid to model an overall y~x effect if the x ranges do not all overlap among sites (see scatterplot below) and sites are completely independent? Wouldn't I be fitting that smooth pooled over data that I possibly cannot actually pool? Therefore, should I even include the term for a shared smooth in the model? If not, I would not be able to answer my first question of an overall effect directly, but would still be able to answer it indirectly from the results for questions 2 and 3.

b) Or should I really just fit separate models per site? This would probably reduce some of the complexity I need to consider and might answer in part my other questions.

2. How do I best include my Covariates? To what level do I need to break up my model terms/interactions, especially if I'm not really interested in the effects and only include them to control for any possible variability due to those variables? E.g. my Covariates 1 and 2 may have interactions with both Type and Site as well as Type:Site, but I'm not really interested in the exact relationships of them. I also include Type, Site, and Type:Site as blocking variables, because there may just be inherent differences I need to account for. If I don't break up my Covariates by all of this, would the blocking variables capture any of this variability for which I don't include separate model terms?

Instead of a GAM, should I maybe be fitting a GAMM? I am not very experienced with mixed models, though.

3. I'm still trying to understand the difference between modeling smooths as interactions with by and as random bs='fs' smooths, in part because it's giving me quite different results both for the smoothers as well as the parametric terms.

   b) I also don't know how to test for significance if modeling this as random smooths, but if I model it as an interaction, my problem of possibly getting significances when there shouldn't be any because the smooth is interpolated for areas where no data exists, remains.

   c) Further, knowing there might be an interaction for Type:Site, would I really need to include another smooth for that interaction or would including the smooths for Type and for Site as well as the interaction as a blocking variable be sufficient?

   d) For my Covariates, is it sufficient to include them as random bs='fs' smooths knowing they probably have significant interactions by Type,Site, and Type:Site, when I'm really not interested in the nature of those effects? (see question 2)

For illustration purposes, I subset my data to include only 2 IDs and 2 Types, but I have more levels for both.

Exploring the variables to consider in my model and how:

Model:

I'm fitting GAMs with mgcv::gam. I'm restricting k to 4 because I do not expect highly complex relationships. I'm not 100% sure this is statistically appropriate?

a) without random smooths for x:

mod1 <- gam(y ~ s(x, k=4) +
              s(x, by=id, k=4) +
              s(x, by=type, k=4) +
              s(Covariate1, interaction(type,id), bs='fs', k=4) + 
              s(Covariate2, interaction(type,id), bs='fs', k=4) + 
                           type + id + interaction(type,id), 
              data=df, method='REML', select=TRUE)

b) with random smooths for x:

mod2 <- gam(y ~ s(x, k=4) +
            s(x, id, bs='fs', k=4) +
            s(x, type, bs='fs', k=4) + 
            s(Covariate1, interaction(type,id), bs='fs', k=4) + 
            s(Covariate2, interaction(type,id), bs='fs', k=4) + 
            type + id + interaction(type,id), 
          data=df, method='REML', select=TRUE)

Or would I possibly need to do something like:

mod3 <- gam(y ~ s(x, k=4) +
            s(x, by=id, k=4) +
            s(x, id, bs='fs', k=4) +
            s(x, by=type, k=4) +
            s(x, type, bs='fs', k=4) + 
            s(Covariate1, interaction(type,id), bs='fs', k=4) + 
            s(Covariate2, interaction(type,id), bs='fs', k=4) + 
            type + id + interaction(type,id), 
          data=df, method='REML', select=TRUE)

or

mod4 <- gam(y ~ s(x, k=4) +
            s(x, by=id, k=4) +
            s(x, by=type, k=4) +
            s(x, interaction(type,id), bs='fs', k=4) + 
            s(Covariate1, interaction(type,id), bs='fs', k=4) + 
            s(Covariate2, interaction(type,id), bs='fs', k=4) + 
            type + id + interaction(type,id), 
          data=df, method='REML', select=TRUE)

Also, I understand some of the effects might appear very small (this is actually part of my hypothesis), but at least some of them are still significant, and the qualitative relationship is still important within the larger context of my analysis.