The semiparametric spatial autoregressive model

In a very general form, the semiparametric spatial autoregressive model reads as:

$$ y_{it}=\alpha y_{it-1}+\rho \sum_{j=1}^N w_{ij,N} y_{jt} + \pi \sum_{j=1}^N w_{ij,N} y_{jt-1}+ \\ \sum_{k=1}^K \beta^*_k x^*_{k,it} + \sum_{k=1}^K \gamma^*_k \sum_{j=1}^N w_{ij,N} x^*_{k,jt}+ \sum_{\delta=1}^\Delta g_\delta(x_{\delta, it}) + \sum_{\delta=1}^\Delta h_\delta\left(\sum_{j=1}^N w_{ij,N} x_{\delta,jt}\right) + \\ \widetilde{ f}(s_{1i},s_{2i},\tau_t)+\epsilon_{it} $$

$$\epsilon_{it}=\theta \sum_{j=1}^N w_{ij,N}\epsilon_{jt}+\phi \epsilon_{it-1}+u_{it}$$ u_it ∼ i.i.d.(0, σ_u²)

where (y_it, x_1, it^*, ..., x_K, it^*, x_1, it, ...,,x_Δ, it) are data observed for a sample of spatial units i (i = 1, ..., N) in each time period t = 1, ..., T. The terms (s_1i, s_2i) indicate the spatial coordinates (latitude and longitude) of the spatial unit i (either a spatial point or the centroid of a spatial polygon), and W_ij is an element of a row-standardized spatial weights matrix. The terms included into the model are:

• y_it − 1 = the time lag of y_it;

• $\sum_{j=1}^N w_{ij,N} y_{it}$ = the contemporaneous spatial lag of y_it;

• $\sum_{j=1}^N w_{ij,N} y_{it-1}$ = the time lag of the spatial lag of y_it;

• α, ρ, π, β_k^*, and γ_k^* = fixed parameters;

• $\sum_{k=1}^K \beta^*_k x^*_{k,it}$ and $\sum_{k=1}^K \gamma^*_k \sum_{j=1}^N w_{ij,N} x^*_{k,jt}$ = parametric linear terms of some covariates x_k, it and of their spatial lags;

• g_δ(.) and h_γ(.) = nonparametric smooth functions of other covariates and of their spatial lags (they can also accommodate varying coefficient terms, smooth interaction between covariates, factor-by-curve intercations, and so on);

• f̃(s_1i, s_2i, τ_t) = an unknown nonparametric spatio-temporal trend;

• ϵ_it = an idiosyncratic error term that can follow a spatial autoregressive process $\epsilon_{it}=\theta \sum_{j=1}^N w_{ij,N}\epsilon_{it}+u_{it}$ with u_it ∼ N(0, σ²) (SEM), a time series autoregressive AR(1) process, i.e., ϵ_it = ϕϵ_it − 1 + u_it with u_it ∼ N(0, σ²), or both.

Obviously, this very general specification is hardly suitable in real data applications. However, it is worth noticing that it nests most of the spatial additive models which can be used in practice. For example, a more suitable semiparametric static model reads as:

$$y_{it}=\rho \sum_{j=1}^N w_{ij,N} y_{jt} + \sum_{k=1}^K \beta^*_k x^*_{k,it} + \sum_{\delta=1}^\Delta g_\delta(x_{\delta, it}) + \widetilde{ f}(s_{1i},s_{2i},\tau_t)+\epsilon_{it}$$

ϵ_it = ϕϵ_it − 1 + u_it u_it ∼ i.i.d.(0, σ_u²) This semiparametric SAR model turns out to be extremely useful to capture interactive spatial and temporal unobserved heterogeneity when the last one is smoothly distributed over space and time (Mı́nguez, Basile, and Durbán 2020). The dynamic extension (including y_it − 1 and $\sum_{j=1}^N w_{ij,N} y_{it-1}$) is also very promising and merits further theoretical investigation. Finally, the following semiparametric SAR model is very useful for modeling cross-setional spatial data taking into account nonlinearities, spatial dependence and spatial heterogeneity:

$$y_{i}=\rho \sum_{j=1}^N w_{ij,N} y_{j} + \sum_{k=1}^K \beta^*_k x^*_{k,i} + \sum_{\delta=1}^\Delta g_\delta(x_{\delta, i}) + \widetilde{ f}(s_{1i},s_{2i})+\epsilon_{i}$$

ϵ_i ∼ i.i.d.(0, σ_ϵ²) In many situations the spatio-temporal trend to be estimated can be complex, and the use of a single multidimensional smooth function may not be flexible enough to capture the structure in the data. To solve this problem, an ANOVA-type decomposition of f̃(s_1i, s_2i, τ_t) can be used, where spatial and temporal main effects, and second- and third-order interactions between them can be identified:

$$\widetilde{ f}(s_{1i},s_{2i},\tau_t)=f_1(s_{1i})+f_2(s_{2i})+f_{\tau}(\tau_t)+ \\ f_{1,2}(s_{1i},s_{2i})+f_{1,\tau}(s_{1i},\tau_t)+f_{2,\tau}+(s_{2i},\tau_t)+f_{1,2,\tau}(s_{1i},s_{2i},\tau_t)$$

First, the geoadditive terms given by f₁(s_1i), f₂(s_2i), f_1, 2(s_1i, s_2i) work as control functions to filter the spatial trend out of the residuals, and transfer it to the mean response in a model specification. Thus, they make it possible to capture the shape of the spatial distribution of y_it, conditional on the determinants included in the model. These control functions also isolate stochastic spatial dependence in the residuals, that is, spatially autocorrelated unobserved heterogeneity. Thus, they can be regarded as an alternative to the use of individual regional dummies to capture unobserved heterogeneity, as long as such heterogeneity is smoothly distributed over space. Regional dummies peak at significantly higher and lower levels of the mean response variable. If these peaks are smoothly distributed over a two-dimensional surface (i.e., if unobserved heterogeneity is spatially autocorrelated), the smooth spatial trend is able to capture them.

Second, the smooth time trend, f_τ(τ_t), and the smooth interactions between space and time - f_1, τ(s_1i, τ_t), f_2, τ, (s_2i, τ_t), f_1, 2, τ(s_1i, s_2i, τ_t) - work as control functions to capture the heterogeneous effect of common shocks. Thus, conditional on a smooth distribution of the spatio-temporal heterogeneity, the PS-ANOVA-SAR (SDM, SEM, SLX) model works as an alternative to the models proposed by Bai and Li (2013), Shi and Lee (2018), Pesaran and Tosetti (2011), Bailey, Holly, and Pesaran (2016) and Vega and Elhorst (2016) based on extensions of common factor models to accommodate both strong cross-sectional dependence (through the estimation of the spatio-temporal trend) and weak cross-sectional dependence (through the estimation of spatial autoregressive parameters).

Furthermore, this framework is also flexible enough to control for the linear and nonlinear functional relationships between the dependent variable and the covariates as well as the heterogeneous effects of these regressors across space. The model inherits all the good properties of penalized regression splines, such as coping with missing observations by appropriately weighting them, and straightforward interpolation of the smooth functions.

Computational approach to estimation

All the non-parametric terms are modeled using Penalized-splines (P-Splines, P. H. Eilers and Marx (1996)). This methodology assumes that the unknown function to be estimated is smooth, and can be represented as a linear combination of basis functions,

f(x) = ∑_jθ_jB_j(x), a popular election is the use of cubic B-spline basis (De Boor 1977). In the case of multidimensional functions, the basis is calculated as tensor products of the marginal basis functions in each dimension. The smoohness of the curve/surface is controlled by adding to the likelihood a penalty tunned by a smoothing parameter, λ (the number of smoothing parameters used is equal to the dimension of the function to be estimated). The penalty can take many forms, depending of the prior knowledge on the curve/surface. The most common choice is to use second order differences on adjacent B-spline coefficients (see Mı́nguez, Basile, and Durbán 2020). In the case of a univariate function, f(x), the penalty is:

Penalty = λ∑_j(θ_j − 2θ_j − 1 + θ_j − 2)².

As we mentioned above, the model is estimated via REML or ML. In order to be able to estimate simultaneously all parameters in the model (including the smothing parameters) we make use of the fact that a P-spline can be reparameterized as a mixed model, and so, all the methodology debeloped in this context can be use for estimation and inference. Furthermore, the mixed model setting allows us to impose the necessary constraints so that all terms in the model are identifiable, as well as, to add spatial (or other) random effects if they are necessary.

To speed up computations, we use a modification of the SOP (Separation of Anisotropic Penalties) algorithm derived by Rodriguez-Alvarez et al. (2015) (for variance components estimation). Also, the use of nested B-spline bases (Lee, Durban, and Eilers 2013) for the interaction components of the spatio-temporal trend contributes to the efficiency of the fitting procedure without compromising the goodness of fit of the model. See Mı́nguez, Basile, and Durbán (2020) for more details.

Installation

This package is available in GitHub (https://github.com/rominsal/pspatreg) and can be installed in the usual way.¹

The function pspatfit()

The main function in the pspatreg package is pspatfit(), which estimates spatio-temporal pernalized spline spatial regression models using either the Restricted Maximum Likelihood (REML) method or the Maximum Likelihood (ML) method. In its generic form pspatfit() appears as:

pspatfit(formula, data, na.action, listw = NULL,type = "sim", method = "eigen", Durbin = NULL, zero.policy = NULL, interval = NULL,trs = NULL, cor = "none", dynamic = FALSE, control = list())

The function pspatfit() returns a list of quantities of class pspat, including coefficients of the parametric terms and their standard errors, estimated coefficients corresponding to random effects in mixed model and their standard errors, equivalent degrees of freedom, residuals, fitted values, etc. A wide range of standard methods is also available for the pspat objects, including print(), summary(), coef(), vcov(), anova(), fitted(), residuals(), and plot().

formula <- y ~ x1 + x2 + pspl(x3, nknots = 15) + pspl(x4, nknots = 20) +
                  pspt(long, lat, year, nknots = c(18,18,8),
                       psanova = TRUE, 
                       nest_sp1 = c(1, 2, 3), 
                       nest_sp2 = c(1, 2, 3),
                       nest_time = c(1, 2, 2))

Plotting smooth terms

Plotting the estimated non-parametric smooth terms represents an important step in semiparametric regression analyses. First, the function fit_terms() computes estimated non-parametric smooth terms. Then, the functions plot_sp2d() and plot_sp3d() are used to plot and map spatial and spatio-temporal trends, respectively, while plot_sptime() is used to plot the time trend for PS-ANOVA models in 3d; finally, and plot_terms() is used to plot smooth non-parametric terms.

Computing marginal impacts

In the case of a semiparametric model without the spatial lag of the dependent variable (PS model), if all regressors are manipulated independently of the errors, ĝ_δ(x_δ, it) can be interpreted as the conditional expectation of y given x_δ (net of the effect of the other regressors). @blu.pow.03 use the term Average Structural Function (ASF) with reference to these functions. Instead, in PS-SAR, PS-SDM or in PS-SARAR model, when ρ is different from zero, the estimated smooth functions cannot be interpreted as ASF. Taking advantage of the results obtained for parametric SAR, we can compute the total smooth effect (total–ASF) of x_δ as:
ĝ_δ^T(x_δ) = Σ_q[I_n − ρ̂W_n]_ij⁻¹b_δq(x_δ)β̂_δq

where b_δq(x_δ) are P-spline basis functions, and β̂_δq the corresponding estimated parameters.

We can also compute direct and indirect (or spillover) impacts of smooth terms in the PS-SAR case as:

ĝ_δ^D(x_δ) = Σ_q[I_n − ρ̂W_n]_ii⁻¹b_δq(x_k)β̂_δq

ĝ_δ^I(x_δ) = ĝ_δ^T(x_δ) − ĝ_δ^D(x_δ)

Similar expressions can be provided for the direct, indirect and total impacts of the PS-SDM.

The function impactspar() computes direct, indirect and total impacts for continuous parametric covariates using the standard procedure for their computation (LeSage and Pace 2009).

The function impactsnopar() computes direct, indirect and total impacts functions for continuous non-parametric covariates, while the function plot_impactsnopar() is used to plot these impacts functions. It is worth noticing that total, direct and indirect impacts are never smooth over the domain of the variable x_δ due to the presence of the spatial multiplier matrix in the algorithm for their computation. Indeed, a wiggly profile of direct, indirect and total impacts would appear even if the model were linear. Therefore, in the spirit of the semiparametric approach, we included the possibility of applying a spline smoother to obtain smooth curves (using the argument smooth=TRUE in the function plot_impactsnopar()).