Accommodating Space-Time Scaling Issues in GAM-Based Varying Coefficient Models

Increasingly the space-time data we analyse and use are not collected under some grand over-arching experimental design, nor for the purposes we indend to use them for. As such, the big idea behind stgam workflows is that it is naive to construct models that make assumptions about the presence and nature of data space-time dependencies, whether for the purposes of prediction or inference (process understanding). Instead these need to be explicitly examined and the most appropriate model form determined. This position is in contrast to a classic statistics perspective where data are considered to be a realisation of carefully considered data collection activity, constructed in such a way as to allow specific hypothesis to be tested.

In this context, most widely used approaches for capturing spatial and temporal dependencies in data and process heterogeneity are flawed because they assume the presence space-time interactions and dependencies. Examples include the alignment of lagged responses to nearby lagged variables in autoregressive moving average models, and existing GAM-based approaches that consider variable interaction over space but with only the aim of model selection and penalization and not process inference [10].

GAMs can be specified with smooths or splines to model non-linear relationships. These are constructed from basis functions that can include single or multiple predictor variables. If a predictor variable is included in a smooth with geographic location (X and Y) then non-linear relationships over space can be modelled - an SVC model. If the smooth is specified with geographic location and time of observation (X, Y and T) then an STVC model is specified. To illustrate this, consider each predictor variable in a SVC or STVC model. It can be specified in 3 different ways:

1. i.

   It is excluded from the model.

2. ii.

   It is included in the SVC / STVC model as a standard parametric response (as in an OLS regression).

3. iii.

   It is included in the SVC / STVC model in a smooth with location (X and Y) but not time.

There are a further 3 ways that each covariate can be specified in an STVC model:

1. iv.

   It is included in the STVC model in a smooth with time (T) but not location.

2. v.

   It is included in the STVC model in a smooth with location and time (X, Y and T).

3. vi.

   It is included in in the STVC model in 2 separate smooths, one with location (X and Y) and the other with time (T).

The intercept is treated in a similar way, but without it being absent. Thus for any SVC with $k$ predictor variables there are $2\times 3^{k+1}$ potential models and for any STVC there are $5\times 6^{k+1}$ potential models.

The stgam package [8] was created to provide a wrapper around GAM functionality in the mgcv package [17], and to allow the user to investigate the different ways each predictor variable could be specified within the GP smooth (as described above), and thus model selection and specification. Its workflow determines the most plausible model given the data and it includes functions that i) create multiple SVC and STVC models defined in different ways as described above, ii) determine the probability of each being the correct model given the data, iii) combine multiple plausible models and, iv) generate spatially and temporally varying coefficient estimates.

stgam is underpinned by 2 core concepts. First, the need to test for the presence and nature of any spatial and temporal dependences and thus to determine whether to include each predictor variable in the model, including within GAM smooths, thereby informing SVC / STVC model specification. In stgam and related initial work [3], this is done by creating multiple models with each predictor variable specified in different ways, as described above. The probability of each model being the best model given the data is then determined and the most probable model is selected. The second core concept in stgam is to evaluate the probability of each potential model being the correct model given the data. This is determined from the model Bayesian Information Criterion (BIC) value [16]. It has been shown elsewhere that the marginal posterior probability of observing the data $D$ given the model $M_i$ – $\text{Pr}(D|M_i)$ – can be approximated using BIC [3], and used to derive the probability of any individual model $M_i$ being the correct model given the data. Of course this is under the assumption that the true model, although unknown, is one of the potential models being evaluated. If multiple models are highly probable, then model averaging can be performed with model weighting from the BIC-derived posterior beliefs in each model.

However, since the release of the first version of the stgam package and paper submission to conferences and peer reviewed journals, some issues with the approach to STVC model construction have been identified. These derive from the way that the approach for SVC modelling GAM with GP smooth was extended to STVC models and the use of GP basis from the mgcv package. Whilst mgcv GP bases may be appropriate for spatial problems and SVC model construction, they are not for STVC models. This is because the mgcv GP basis results in spatial models with a potentially suboptimal estimate of the length scale of the correlation basis functions used in the GP basis, and thus temporal models that are erroneous for all time series except those with the longest of length scales.

The substantive problem with the mgcv GP basis is that it only fits within the penalized spline class of models that can be estimated by mgcv if the basis function parameters are specified. The critical parameter here is the correlation function length scale or range parameter, $\rho$ which may be specified as spherical, power exponential, or as one of three forms of Matérn correlation with $\kappa$ equal to 1.5, 2.5, or 3.5. The $\rho$ parameter determines the distance at which the correlation function falls below some value. If $\rho$ is not supplied, then mgcv fits a smooth using Matérn correlation functions with $\kappa =1.5$ and $\rho ={\max }_{ij}\parallel {𝐱}_i-x_j\parallel$ as the basis functions, for any pair of points $x_i$ and $x_j$, following [14]. If an order penalty is specified then one of the correlation functions is implemented but again using Kammann and Wand’s recommendation for the length scale, the maximum distance attained by any pair of points $x_i$ and $x_j$.

Evidently for SVCs (and similar spatial analyses) this is unlikely to be problematic except in the presence of long or short spatial dependencies. However, for temporal and spatiotemporal data, this specification of length scale in this way is problematic: it implies that the correlation between pairs of points will only fall below some small value when those points are separated by an amount of time equivalent to the time series itself. A space-time GP smooth specified in this way will be isotropic, with similar levels of non-linearity in all dimensions, such that the spatial heterogeneity is equal to temporal heterogeneity. This is clearly not correct but the current implementation of stgam has no way to control this, because the functions that create different models hardcode a GP basis in the smooths, and there is currently no option to specify an order penalty or $\rho$. The result is GP smooths specified with pure Kammann and Wand basis functions, despite their unsuitability for time series or space-time problems.

The proposed high level solutions to these problems are 1) to allow users to specify the length scales of the different components, and 2) to fit the spatio-temporal terms via a tensor product smooth of a marginal spatial smooth and a marginal time smooth, both specified with a GP basis, and effectively allowing space-time to be treated in a three-dimensional way. This can be done by specifying a 2D spatial GP smooth for the first margin of the tensor product, and a 1D temporal GP smooth for the second margin. These two bases will smoothly interact creating the desired spatio-temporal varying term due to the tensor product construction.

In mgcv syntax, the change in how the GP smooths for each predictor variable var are specified as follows:

The tensor product smooth requires the marginal basis dimensions to be specified (d above). Here the space-time smooth contains 3 covariates composed of a tensor product of a 2-D thin plate regression spline basis for location, and 1-D basis for time.

However, the correlation function length scale or range parameter, $\rho$ also needs to be specified as, both in terms of its form (spherical, power exponential, etc) and distance. In the mcgv GAM implementation this is done for each of the margins through the m argument passed to the tensor product smooth. This expands the specification of the tensor product smooth with GPs to:

Here rho1 is the length range for the spatial marginal smooth and rho2 is that for the temporal marginal. Both are specified with a power of 3, indicating the distance decay of the range. For the spatial margin intervals of 100km for a case study in a large country (USA, China, Brazil, etc) could be explored, 10 km for a small country (UK, Germany, Vietnam, etc) or 1km for a local one. These may be the equivalent of predictor variable bandwidths in MGWR. For the temporal dimension, intervals of 1, 2, 3 etc years could be explored, depending on how time was recorded in the data. The revised stgam package will support investigation of different forms and length scales for spatial and temporal margins.

The ability to vary the scales and lengths ranges in the tensor smooths as above allows each each predictor variable to be specified in a way that treats implicitly space-time as three-dimensional. Therefore a second set of investigations will be supported in the revised stgm package to allow users to determine which of the possible six forms described above is the appropriate way to include each predictor variable in SVC and STVC models. For example an alternative to the mgcv tensor product smooth with GPs, and illustrated above is one that has separate 2-D spatial smooths and 1-D temporal ones:

Essentially, this replaces the investigation of different model forms using mgcv smooths specified (s() in mgcv), with ones specified with mgcv tensor product smooths (te() in mgcv). The analysis undertaken in this short paper describes initial work investigating how both of these considerations (space-time lengths and predictor variable form) can be optimised and potentially included in a revised stgam package.

The stgam package includes two datasets describing annual economic productivity for the 48 contiguous US states (with Washington DC merged into Maryland), 1970 to 1985 from the plm R package [9] and the spatial dataset of the 48 contiguous US states spData package [1]. The productivity data contains variables describing Private Capital stock (privC), Unemployment % (unemp) and Public capital investment (pubC). The Unemployment variable over time is shown in Figure 1, with 1986 not plotted for aesthetic reasons. The aim of the analysis was

1. 1.

   to determine optimal $\rho$ values for the correlation function length scales, $\rho$, a 2D spatial margin and a 1D temporal margin (${\rho }_s$ and ${\rho }_t$, rho1 and rho2 respectively in the above code snippet).

2. 2.

   to then determine the most appropriate STVC model of Private Capital stock (privC), with Unemployment and Public capital investment as predictor variables, specify the tensor product smooths in different ways as described in Section 2.1.

In both cases, $\rho$ and model form were evaluated from the model BIC.