Random Effects Eigenvector Spatial Filtering

• Random Effects Eigenvector Spatial Filtering is a framework that models spatially varying coefficients as random effects using positive Moran eigenvectors.
• It integrates eigenvector spatial filtering with mixed-effects modeling to flexibly control spatial scales and effectively manage spatial autocorrelation.
• The method achieves computational efficiency through the Nyström extension and reduced matrix operations, enabling robust analysis of large datasets.

Random Effects Eigenvector Spatial Filtering (RE-ESF) is a statistical modeling framework for capturing structured spatial variation in regression coefficients through random-effects expansions on Moran eigenvectors. The method embeds spatially varying coefficients (SVCs) into a mixed-effects regression, providing coefficient-specific spatial scale control, computational scalability for large datasets, and interpretability in terms of spatial autocorrelation as captured by the Moran coefficient. RE-ESF combines the advantages of eigenvector spatial filtering, flexible bandwidths, and random-effects modeling, delivering accurate and stable inference for spatially heterogeneous processes and covariate effects (Murakami et al., 2017, Murakami et al., 2016, Murakami et al., 2017).

1. Mathematical Formulation

At its core, RE-ESF models a spatially varying coefficient regression: $$y_i = \sum_{k=1}^K X_{i,k}\,B_k(s_i) + \varepsilon_i, \quad \varepsilon_i \sim N(0,\sigma^2)$$ where $y_i$ is the response, $X_{i,k}$ the value of covariate $k$ at location $s_i$, and $B_k(s_i)$ the spatially varying coefficient. RE-ESF represents each coefficient surface as: $$B_k(s_i) = \beta_{k0} + \sum_{\ell=1}^L e_\ell(s_i)\,y_{k,\ell}$$ with $$\mathbf y_k = (y_{k,1},\ldots,y_{k,L})' \sim N\big(\mathbf 0, \sigma_{y,k}^2\,A(a_k)\big)$$.

Here, the eigenvectors $e_\ell$ correspond to positive-eigenvalued bases of the centered connectivity matrix $MCM$, which capture spatial autocorrelation structure as measured by the Moran coefficient: $$MC(e_\ell) = \frac{N}{1'C1} \lambda_\ell, \qquad \lambda_\ell > 0$$ The diagonal matrix $A(a_k)$ parameterizes shrinkage by spatial scale: $$[A(a_k)]_{\ell\ell} = \left(\frac{\lambda_\ell}{\lambda_1}\right)^{a_k}$$ Larger $a_k$ values induce smooth, large-scale spatial variation; small $a_k$ values retain fine-scale heterogeneity (Murakami et al., 2017, Murakami et al., 2016).

The full model in matrix notation: $$\mathbf y = X\boldsymbol\beta + \tilde{E}\,\mathbf u + \boldsymbol\varepsilon,\quad \mathbf u \sim N(\mathbf 0, \Sigma_u)$$ with $$\Sigma_u = \mathrm{blockdiag}(\sigma_{y,1}^2\,A(a_1),\dots,\sigma_{y,K}^2\,A(a_K))$$ and $\tilde{E} = [X_1\circ E, \dots, X_K\circ E]$, where $\circ$ is the Hadamard product.

2. Eigenvector Selection and Spatial Basis Construction

RE-ESF employs the real symmetric Moran operator $MCM$, where $M = I_N - \frac{1}{N} \mathbf{1} \mathbf{1}'$ is the centering matrix, and $C$ is a non-negative, symmetric spatial weights matrix, e.g., $C_{ij} = \exp(-d(s_i, s_j)/b)$.

Eigen-decomposition yields a set of orthonormal eigenvectors and corresponding eigenvalues: $$M C M = E_{\mathrm{full}} \Lambda_{\mathrm{full}} E_{\mathrm{full}}'$$ RE-ESF uses the subset $E = [e_1, \ldots, e_L]$ with strictly positive eigenvalues. Unlike fixed-effects ESF, which typically selects a subset via stepwise methods, RE-ESF utilizes all $L$ positive eigenvectors, regulating their influence via shrinkage in the random-effects prior (Murakami et al., 2017, Murakami et al., 2017).

For large datasets, computation can be accelerated using the Nyström extension to approximate the leading eigenvectors (see Section 4). Empirically, retaining $L \geq 200$ yields negligible loss of statistical efficiency for datasets with up to hundreds of thousands of observations (Murakami et al., 2017).

3. Random-Effects Structure and Scale Parameters

Each coefficient process $k$ in RE-ESF possesses two variance components: $\sigma_{y,k}^2$ (overall marginal variance) and $a_k$ (controls spatial scale of smoothness). The prior covariance for the random coefficients is diagonal in the eigenbasis, imposing stronger shrinkage on high-frequency (low-eigenvalue) directions when $a_k$ is large: $$[A(a_k)]_{\ell\ell} \to 0 \text{ as } a_k \gg 1 \text{ for } \lambda_\ell \ll \lambda_1$$ allowing each regression coefficient to have a distinct degree of spatial smoothness. The interpretation is direct: $a_k$ regulates the spatial scale in the random-effects expansion of $B_k(\cdot)$, analogous to a bandwidth in kernel SVC models but estimated directly from data via REML (Murakami et al., 2017, Murakami et al., 2016).

4. Estimation Algorithm and Computational Complexity

RE-ESF estimation proceeds by maximizing the restricted (type-II) maximum likelihood (REML) of the mixed model. The log-profile restricted likelihood is: $$\ell_{\mathrm{R}}(\Theta) = -\frac{1}{2} \left\{\log |V(\Theta)| + \log |X'V(\Theta)^{-1}X| + (\mathbf y - X\hat\beta)'V(\Theta)^{-1}(\mathbf y - X\hat\beta)\right\}$$ where $$V(\Theta) = \sigma^2 I_N + \tilde{E} \Sigma_u \tilde{E}'$$ and $\hat\beta = (X'V^{-1}X)^{-1} X'V^{-1} \mathbf y$.

Computational workload is dominated by:

• Eigen-decomposition: $O(N^3)$ for full computation, but with the Nyström extension, this reduces to handling $L \times L$ blocks with L ≪ N, e.g., $L \sim 200$ (Murakami et al., 2017).
• Precomputing cross-products: $X'X$, $X'\tilde{E}$, $\tilde{E}'\tilde{E}$ (all $(K+KL) \times (K+KL)$).
• Optimization: Maximizing over $2K + 1$ parameters, requiring $O((K+KL)^3)$ per likelihood evaluation—crucially, this is independent of $N$.

For massive datasets, the two main algorithmic accelerations are:

1. Nyström Extension: Approximates the leading Moran eigenvectors using a small set of spatial "knots." The error in eigenvalues and eigenvectors is $O(n/L)$, empirically negligible for $L \geq 200$.
2. Small-Matrix Tricks: After projecting data onto the reduced eigenbasis, all subsequent calculations involve only $(K+L) \times (K+L)$ matrices, allowing $n$ to be arbitrarily large without impacting memory or CPU bottleneck (Murakami et al., 2017).

5. Monte Carlo Evaluation and Empirical Findings

Monte Carlo simulation studies have systematically benchmarked RE-ESF alongside geographically weighted regression (GWR), flexible bandwidth GWR (FB-GWR), standard ESF, and local coefficients regression. Key aspects of these designs include:

• Data-generating processes with true coefficient surfaces based on spatial moving averages at small or large scale.
• Covariates with varying spatial dependency (controlled by $r_x$ and $b_x$), closely mimicking collinearity and spatial confounding scenarios.

Performance was measured via RMSE, MAE, bias, effective degrees of freedom ($p^* = \mathrm{tr}(H)$), and CPU time. Findings (Murakami et al., 2017, Murakami et al., 2017, Murakami et al., 2016):

• Accuracy: RE-ESF achieves the lowest RMSE for significant, fine-scale SVCs, particularly when covariates are also spatially structured. It is consistently on par with or superior to FB-GWR and substantially outperforms standard GWR and ESF, especially in the presence of local scale variation and spatial confounding.
• Complexity Control: The effective degrees of freedom $p^*$ in RE-ESF remain stable under fine-scale processes due to adaptive eigenbasis shrinkage, mitigating overfitting seen in fixed-bandwidth or non-regularized models.
• Computational Efficiency: For $N=400$, RE-ESF required 3.38 s (vs. 1.50–10.31 s for flexible GWRs, 65.41 s for ESF). For large $N$, the fast implementation with $L=200$ handles $n>10^5$ in tens of seconds or less (Murakami et al., 2017).

Model	RMSE (β₁, fine-scale)	p* Stability	CPU Time (N=400)
RE-ESF	Lowest	Stable	3.38 s
FB-GWR/FB-GWRa	Slightly higher	Variable	1.50–10.31 s
ESF (no reg.)	High	Unstable	65.41 s
GWR	Poor (for local β₁)	Low	0.13–0.18 s

6. Practical Implementation, Best Practices, and Limitations

The spmoran R package implements fast RE-ESF, enabling practical analysis on large spatial datasets (Murakami et al., 2017). Standard workflow:

1. Construct Moran eigenvectors (via Nyström extension): meig <- meigen(coords, k=200)
2. Fit model: resf <- resf(y~x1+x2, data, meig)
3. Examine results: summary(resf) (β̂, γ̂, σ̂², α̂)
4. Prediction for new locations: predict(resf, newdata)

Recommendations include:

• Use $L \geq 200$.
• Only positive spatial dependence (λ_ℓ > 0) is modeled; negative spatial autocorrelation is not captured.
• Range parameter $r$ in the spatial kernel must be set or estimated a priori.
• Residual Moran's I should be checked post-fit to verify adequate spatial filtering.
• RE-ESF is robust to the choice of kernel (exponential, Gaussian, spherical).

Limitations: Negative spatial dependence cannot be modeled; estimation of the kernel range is not automatic within the baseline framework; accuracy may degrade if too few eigenvectors are chosen in highly complex spatial fields.

7. Interpretive Context and Applicability

RE-ESF provides a unified, scalable, mixed-effects formulation for spatially varying coefficient modeling, exploiting the geometric and statistical properties of Moran eigenbases for spatial filtering. The hyperparameters ($\sigma_{y,k}^2$, $a_k$) have interpretable roles as coefficient-specific spatial variance and scale, replacing the need for ad hoc bandwidth selection. The framework is particularly advantageous in contexts with heterogeneous spatial processes, severe spatial confounding, and large sample size. Empirical work in land value hedonic modeling and simulation confirms its superior stability and interpretability relative to GWR/ESF, as well as computational tractability for large $N$ (Murakami et al., 2017, Murakami et al., 2016).

RE-ESF sits at the intersection of spatial random-effects, low-rank spatial regression, and eigenvector filtering, and provides an extensible platform for further methodological development in spatial statistics.