Generalized additive model: choosing between cubic and thin-plate splines

Question

I am using the gam function (from the mgcv package) to model a continuous response (a soil nutrient) in relation to a continuous predictor plus categorical predictors plus site random effects (modest number of sites) with:

    gam(log(y) ~ s(x, bs = "cr", by = fac1) + 
              fac1 + fac2 + s(site, bs = "re", by = flag), 
        method = "REML")

The data come from a landscape study where sampling was challenging (to say the least), so the number of samples is small and the coverage of the range of $x$ values is quite patchy in the upper part of the range.

Use of a cubic spline for $x$, as above, results in quite a different shape for the fitted function compared to the default thin-plate spline for that part of the space where $x$ values are sparse. However, both models have almost identical results in terms of deviance explained and pattern of residuals. The shape of the function produced by the cubic spline is much preferred by my colleagues who conducted the study because they feel they can reconcile it with theory more readily than the alternative. However, I am wondering if:

• There are other criteria that I can investigate to choose between the models?
• Any theoretical reason to choose between using "cr" or "tp" for a univariate smooth with patchy data?

Update: plots of the alternative models

Answer 1

Quite an old question that seemed to be mostly answered in the comments already. But in any case, it seems clear that the thin plate regression splines are the clear contender here. The curves on the cubic regression splines make no sense and have a ridiculous standard error region.

Something to add beyond that...there really aren't a lot of observations here, which contributes to the problem greatly. A future design would undoubtedly be better with even a little more data sprinkled in. Its hard to understand exactly what the actual function relationship is when you have such extreme values on one end and a clearly flat association with lower values. The real relationship may actually have a much stronger exponential rise with increases in $x$. But without more data, that remains totally unclear.

A full Bayesian implementation may be better in this respect. Incorporating prior information when fitting the curves (such as those possible with brms or bamlss) would do wonders to show what a more plausible fit should be given how few datapoints there are.