Resource for understanding sources of uncertainty in smooth_estimates()

[Aariq]

Is there a reference or blog post out there to help me better understand what the unconditional and overall_uncertainty arguments do? I'm hoping to use overlap with confidence intervals to identify parts of a tensor product smooth that are significantly different from the model intercept, but I don't quite understand how the settings for these arguments change the interpretation of the CIs.

[Aariq]

Ok, reading Marra and Wood 2006, and I think I understand what problem overall_uncertainty fixes. When the true shape of a relationship is close to a straight line (or a plane for a 2d smooth), it will tend to get penalized to a straight line (or plane). In that case, the CIs will approach 0 at some point because of identifiability constraints and be narrower than they should be. overall_uncertainty = TRUE fixes this by "applying identifiability constraints to all other model components, but allowing the component of interest to 'carry the intercept'".

[gavinsimpson]

You (now) have the relevant reference for overall_uncertainty. More generally, when the bias in the estimated smooth is large relative to something (the variance? I forget now) then the theory of Nychka (1998), which is expanded upon by Marra and Wood, breaks down. So any situation where there is appreciable bias in the estimated smooth leads to poor coverage in the credible intervals. Such a situation would be where you have oversmoother, such as where you failed to set the basis dimension sufficiently large so that the basis might feasibly contain the true function or a close approximation to it. Another issue this tackles is the bow-tie intervals when smoothness selection penalizes the term back to a linear function. In such cases the interval is impossibly of 0 width because of the identifiabilty constraint that means the function must sum to 0 over the range of the covariate; hence where it crosses the 0 line there is no uncertainty. This isn't a reasonable reflection of the uncertainty so setting overall_uncertainty = TRUE pulls in the uncertainty in the intercept to each smooth in turn resulting in more reasonable intervals.

Importantly, this is all done post hoc by processing the model rather than something done a priori before fitting; the usual identifiability constraints are still used when fitting.

For unconditional see section 4 of Wood et al (2016) and also the comment/discussion paper in section 3 of Greven & Scheipl (2016) which is a discussion paper on Wood et al (2016). Basically it just widens the intervals to allow for the extra uncertainty in the model due to the selection of the smoothness parameter(s) λ_j. With unconditional = FALSE, the intervals are formed assuming that we knew the value(s) of the smoothness parameter(s). In practice we don't know these values and estimate them from data so these intervals yield coverage that is conditional upon the values of the smoothness parameters. Setting unconditional = TRUE uses extensions of the theory developed by Nychka and Marra and Wood to provide a simple correction to the Bayesian covariance matrix of the model parameters for the fact that we didn't know the value(s) of the smoothness parameter(s) and in fact had to select them during fitting.

All of this and more is briefly reviewed in Simon's recent (2020) invited paper, to which there is also a discussion.

[gavinsimpson]

For your application, you might want to consider simultaneous intervals instead of the Bayesian credible intervals. In some senses the Bayesian credible intervals are still pointwise (although they have an average across the function interpretation rather than true pointwise interpretation when viewed in this frequentist manner). This means that the coverage at a single point is approximately correct, but as you start to look at more points, as you plan to do, the coverage over that set of points (or the whole function) ends up being a lot less when combined than the notional coverage you used for each single point (0.95 say).

A simultaneous interval is one where the interval (we hope) covers 0.95 (say) of all possible smooths give the estimates of the smooth and it's uncertainty. There's some background on the idea in Simpson (2018) and general code in the supplements to do this for smooths. {gratia} contains a more recent implementation, with the confint() methods for GAM(M)s allowing simultaneous intervals; see ?confint.gam

[Aariq]

Thank you for the super detailed answers! It does sound like simultaneous intervals are more of what I want, so I'll look into that. I tried out gratia::confint(), but ran into some trouble—I'll post an issue with a reprex.