They actually deal with non-additive "low-order" interactions quite well. In R's mgcv for example, let's say you had data from many years of temperature readings across a wide geographic area, so your data are (lat, long, year, temperature). mgcv lets you fit a model like:

  gam(temperature ~ te(long, lat) + s(year) + ti(long, lat,year))  

where you have (1) a nonlinear two-way interaction (i.e. a smooth surface) across two spatial dimensions, (2) a univariate nonlinear effect of time, and (3) a three-way nonlinear interaction, i.e. "does the pattern of temperature distributions shift over time?"

You still can't do arbitrary high-order interactions like you can get out of tree-based methods (xgboost & friends) but that's a small price to pay for valid confidence intervals and p-values. For example, the model above will give you a p-value for the ti() term, which you can use as formal statistical evidence to say -- at what level of confidence -- a spatiotemporal trend exists.

This Rmarkdown file (not rendered sadly) shows how to do this and other tricks https://github.com/eric-pedersen/mgcv-esa-workshop/blob/mast...

Hey cool. I'll take a closer look then. Thanks! I assume that there are mixed model variants out there too.