Understanding effect plots of GAM with categorical interactions and differing range of predictor variable

Question

I'm trying to model complex relationships of a continuous predictor variable that has interactions with two categorical variables. I have skew in my samples as in that I have differing predictor ranges for each of the categorical combinations, as illustrated in this plot:

If I now model this with mgcv::gam as

mod<-gam(y ~ s(x, by=interaction(Type,ID)) + ID + Type, data=dat, method='REML')
plot(mod, pages=1, scale=0, shade=T)

the resulting effect plots look like this:

So, for each categorical combination, the effect plot shows a fit and rug over the entire pooled predictor range, e.g:

I'm wondering if I can even do this or if there is something I need to include in my model to control for this. Possibly, I just have to fit numerous single models? Because as far as I understand it, right now, it's fitting relationships over data that doesn't always exist, probably resulting in a spurious fit and significance? For example, x:A:3 is highly significant, but if I mentally take out that decrease at the beginning where no data exists for this combination, I see a mostly flat line...

In addition, do I need to be concerned over very differing sample sizes per categorical combination?

I know for sheer visualization, I can just extract the underlaying data used to create the effect plots and re-plot them while excluding data outside the range for each combination. But I want to make sure my model approach is even correct and the results valid. This is a simplified version for illustration purposes of a more complex model that includes additional predictors. In addition, while ID is a classic blocking variable, the effects of the different Types may possibly influence each other, making me think I really shouldn't fit individual models for each separate Type:ID combination, but then the problem of differing predictor ranges remains.

Answer 1

Your uncertainties are growing in over ranges where the splines are extrapolating to new data, so in the present model configuration the question is whether you are comfortable with the way this extrapolation happens (it goes linearly in whatever direction the spline was "wiggling" at the boundary). But you may also benefit from a hierarchical setup, where you could perhaps fit higher-level smooths that these lower level smooths are then regularized toward. For example, you could consider something like:

mod <- gam(y ~ s(x) +                              # shared smooth
            s(x, ID, bs = 'fs') +                  # avg. ID smooths
            s(x, Type, bs = 'fs') +                # avg. Type smooths
            s(x, interaction(Type,ID), bs = 'fs'), # lower-level smooths
           data = dat, method = 'REML')

The fs basis acts as "random" smooths here (which also implictly include random intercepts), so you are in effect fitting a "shared" smooth of x to learn the average shape from all the data at once. You are then allowing each level of ID to have its own "average" smooth, which acts as a deviation from the shared smooth. The rest of the model proceeds down the hierarchy. This would perhaps give you more efficient use the data structure. You may need to alter the types of penalties you use or the k values you impose to ensure you retain identifiability (see guidance from Pedersen et al on hierarchical GAMs). But overall it is good that you are investigating where these estimates come from when extrapolating to new data, it is often not well recognised that these functions extrapolate in this way.