Difference smooths for differences of differences?

[sebsau]

Is it possible to plot difference-of-difference smooths? I.e., difference smooths that don't take smooths but difference smooths as their input?

I have a study that crosses two factors with two levels each. Let's call them factor A (with levels 1 and 2) and factor B (also with levels 1 and 2). I am interested in the spatial extent of differences between the factor levels. Specifically, whether there is an interaction of these factors, i.e., whether the spatial difference of A1-A2 is different from the spatial difference of B1-B2.

If I understand this correctly, then this could be shown with a differences-of-differences smooth. One would fit a GAM in mgcv with a smooth with a four-level factor as by-variable that combines the two original factors and has the levels A1, A2, B1, B2.

response ~ 1 + ABFACTOR + s(x, y, by = ABFACTOR) (where x and y are the spatial coordinates)

To show the interaction, or the differences of differences, one would then create two difference smooths (A1-A2 and B1-B2) and then substract them from each other.

Is this possible with the difference smooth plotting functions already implemented in gratia?

[gavinsimpson]

It isn't possible in {gratia} automagically. If you are willing to accept some assumptions (at the moment, that the posterior distribution of the parameters of the smooth/model are well approximated by a multivariate normal) then you could plug bits of existing {gratia} functionality together to achieve what you want via posterior sampling.

Computing what you want could easily be done by simply subtracting one difference from another using {dplyr} tools applied to the two differences. What is more difficult is the uncertainty - I'm not sure if the difference of the difference is a simple function of the parameters or not. But accepting some assumptions and some amount of compute time you could evaluate the difference of your two differences of smooths for a number of draws from the posterior of the differences (or the model) and generate an uncertainty estimate that way.

If you want to know more about this, a reproducible example using some simulated data would be helpful and I'd figure out what was needed to do the sampling and compute the differences