# Modelling spatiotemporal data with GAMs

I am currently modelling a data set of house prices, where some of the important covariates are the house coordinates (x,y) and the date of the house sale, t. In doing so, two questions have popped up:

Q1:

In the hunt for the best model, I have tried various options for modelling the spatiotemporal dependence in the data set. My first approach has been to model this dependence using the default thin-plate splines provided in mgcv coupled with a rather large number of basis functions

Model 1: $$s(x,y,t, k = 1000)$$

Model 1 appears to do reasonably well. However, after reading a bit more on the topic, it seems that in cases where covariates are not in the same units, a better choice would be to use tensor splines. Hence, as a second approach I have tried a model of the form:

Model 2: $$te(x,y,t, k = 10)$$

I have purposedly used k=10 (which to my understanding is equivalent to 10^3 = 1000 basis functions in total) to make the number of basis functions for model 1 and model 2 comparable. To my surprise, model 2 performs a lot worse than model 1 as measured by the gcv score. One possible explanation I could come up with, is that the choice of k = 10 forces an equal number of basis functions in spatial and temporal domain, which might not be optimal - typically the temporal dependence is much simpler than the spatial one. I have therefore tried a third model of the form:

Model 3: $$te(x,y,t,k=c(16,16,4))$$

This model does perform slightly better than model 2 but still significantly worse than model 1.

Could anyone elaborate on what is going on here? I suspect that my assumptions about basis function numbers might be wrong. Also, does anyone know how mgcv estimates the number of basis functions to use in each direction for the default thin-plate splines s(x,y,z..)?

Q2:

Ultimately, I am interested in a model that models the house prise as a sum of the non-interacting spatial and temporal parts as well as their interaction, i.e.

$$P(x,y,t) = s(x,y) + s(t) + ti(x,y,t)$$

I am however unsure, whether the model specification above is correct. More specifically, I am worried that interaction between x- and y is contained both in s(x,y) as well as the ti(x,y,t). Could anyone clarify this for me?

Thanks a lot!

The model s(x,y,t, k = 1000) is almost surely incorrect. This isotropic thin plate spline assumes that there is the same wiggliness is all dimensions, and you are fitting a spatiotemporal smooth so it is quite unlikely that you have the same wiggliness in space as you do in time, regardless of whether the covariates are in different units or not.

The model te(x,y,t, k = 10) is also potentially suboptimal if we want something a bit like the s() version for space. You can fit a tensor product of a 2d isotropic thin plate spline and a 1d smooth of time. You do this via the d argument, so I would probably try something like

z ~ te(x, y, t, d = c(2,1), bs = c(200, 10), bs = c("tp", "cr"))


for example to get an isotropic smooth of the spatial coordinates. Depending on your domain, look at a Duchon spline in place of the thin plate spline (bs = "ds") but read it's help page well - it can be useful if you want the spline to tend towards a flat function as it moves away from the support of the data (so no wild extrapolation). You can also increase k for the spatial smooth as necessary.

One reason the te() might be doing worse (in your example) is that you're not enforcing the isotropy or viewing two of the dimensions in the smooth as representing a single spatial smooth.

However, you probably shouldn't be fitting these models with GCV; unless you have a very good reason not to, always fit with method = "REML" or method = "ML". GCV is known to undersmooth.

Also, does anyone know how mgcv estimates the number of basis functions to use in each direction for the default thin-plate splines s(x,y,z..)?

It doesn't work like that. The TPRS is generating 1000 3D basis functions in your example, by truncating the full rank TPRS (which would be a radial basis function at each unique data combination) via an eigendecomposition. For example, the basis functions generated via

s(x, z, k = 25)


are shown in the figure below after application of the identifiability constraints. Note that these basis functions are 2D — yours would be 3D but that's much harder to visualise.

To get the decomposed interaction version of your model, you need to create a ti() that contains both the marginal smooths that you have specified separately. This also hints at how to set up a 2D spatial smooth interacting with a temporal smooth via te().

y ~ s(x,z, bs = "tp", k = 100) + s(t, bs = "ps", k = 10)


then we need to tell mgcv to generate the same bases as margins for the pure interaction tensor product smooth

y ~ s(x,z, bs = "tp", k = 100) + s(t, bs = "ps", k = 10) +
ti(x, z, t, bs = c("tp", "ps"), d = c(2, 1), k = c(100, 10))


for example. Note how we have to repeat exactly the information on the respective bases. The k on the ti() doesn't have to be the same as the combination that you used in the "main" effects part, and doing as I did here will result in lots of basis functions to it's going to be a slow model to fit, but this choice was purely for illustration. Typically, one would use fewer basis functions for the ti() marginals as you's already got the main effects modelled through the two s() terms. My choice of "ps" for the temporal smooth was also arbitrary — if you are using s() for the "main" effect smooths, you need to be aware that the default basis for s() is not the same as for ti() (or any tensor product) as the latter use bs = "cr". So it is always advisable to specify explicitly what bases you want for both the "main" effect smooths and the ti() interaction. Here I just chose something different from either of the defaults for s() or ti(); "ps".

More specifically, I am worried that interaction between x- and y is contained both in s(x,y) as well as the ti(x,y,t). Could anyone clarify this for me?

The way the ti() smooth is formed, results in the main effects of the marginal smooths being parameterised out of the basis. This way, the main effects of space are only included via s(x, z), and the same for t.

• Thanks a lot for the very helpful answer. Following your suggestions I switched to using REML and fitted a smooth using the d-argument as: z ~ te(x,y,t, d = c(2,1), k = c(1000, 5), bs = c("tp", "cr")) I now get a better-performing tensor-product smooth but it still does not beat the original s(x,y,t) (as measured by RMSE), which I would have assumed based on your comment that the latter wrongfully assumes equal wigglyness in the temporal and spatial domain. Do you have an idea why it does not perform better? Also, what is the reason for using cubic splines in the time-domain? Feb 28 at 14:55