Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
62,967 characters · 13 sections · 67 citation commands
Moran's $I$ Lasso for models with spatially correlated data
\doublespacing
In conventional spatial economic modeling, the researcher is required to specify (i) a spatial weights matrix (SWM)\footnote{A spatial weights matrix is an $n\times n$ matrix describing the pair-wise relationships between the $n$ cross-sectional units.} and (ii) which parts of the model are spatially correlated. Standard specifications typically include one or more spatial lags of the dependent, exogenous and/or error term kel17. Historically, applied researchers have generally focused more on specifying the SWM rather than the empirical spatial structure lsK_sens14 and a standard robustness check in the applied spatial economic literature tests whether the estimates are sensitive to different SWMs. When estimates are found to be sensitive to the choice of SWM, researchers have attributed this sensitivity to the choice of SWM. However, as lsK_sens14 shows, estimates should not be overly sensitive to the choice of SWM, as long as they are reasonably well correlated. This implies that the sensitivity many researchers observe is driven by misspecification of the spatial economic model rather than the choice of SWM. lsK_sens14 thus argue that researchers should focus on specifying the spatial model rather than finding an ideal SWM.
The Eigenvector Spatial Filtering (ESF) approach of grif2000, grif03 uses a subset of eigenvectors from the SWM as controls to filter out terms involving the SWM in the underlying model. ESF has recently started receiving substantial attention from applied economic researchers.\footnote{Some examples include PGTN11_ESF, RSGN_JRS12ESF, esfJAE13, CS14ESF, ESF_sjpe18, bg_bcon08esf, gp11ESF.} ESF's main advantage over conventional maximum likelihood (ML) and generalised method of moments (GMM) is precisely that researchers need not specify spatial correlation explicitly in the model, or estimate the corresponding spatial parameters. In the context of ESF these are instead viewed as nuisance parameters. The fact that ESF is agnostic to the underlying spatial process is desirable for applied researchers who simply wish to obtain unbiased parameter estimates in the presence of cross-sectional dependence in the covariates. This is because it is easier to establish the presence of an underlying spatial process via a test for spatial correlation than it is to determine the exact specification of this process.
The critical challenge for ESF is that the spectral decomposition of the $n\times n$ SWM yields $n$ eigenvectors and if all are included in the model, it becomes high-dimensional and estimation by Ordinary Least Squares (OLS) is infeasible.\footnote{A high-dimensional model is defined as a model with more parameters to estimate than observations, leading to a rank-deficient Gram matrix.} grif03 argues that only a subset of eigenvectors is necessary to eliminate the cross-sectional dependence in the dependent variable. The key question becomes identifying which subset of eigenvectors is required, which we refer to as the ESF eigenvector selection problem. Several solutions to this selection problem have been proposed, such as several stepwise greedy algorithms where eigenvectors are iteratively added until some user-specified threshold is reached grif2000,grif03,tiefgrif07. These stepwise greedy algorithms are simply heuristic approximations to the full ESF selection problem, thus, they are necessarily sub-optimal. Under the assumption of sparsity (i.e. most eigenvector coefficients are zero) seya15 proposes using an $\ell_1$-penalised regression, e.g. Lasso. Given that Lasso estimates are ultimately determined by a tuning parameter, this turns the eigenvector selection problem into a tuning parameter calibration problem. seya15 propose estimating the tuning parameter using conventional $K$-fold cross-validation (CV) with prediction accuracy as the loss function. However, the existing theoretical results on CV-Lasso assume the cross-sectional units are independent czc20. This is hard to justify in the context of ESF, where the eigenvectors are derived from a matrix that encodes cross-sectional dependence. Additionally, the goal of ESF is to eliminate spatial correlation patterns, not improve prediction accuracy. There is therefor no guarantee that running CV with a prediction accuracy loss will yield consistent eigenvector selection.
We propose an alternative procedure for choosing the ESF Lasso tuning parameter, called Moran's $I$ Lasso (Mi-Lasso), which directly uses information about the level of correlation in the residuals provided by the Morans $I$ statistic moran50 to develop a point estimate for the Lasso tuning parameter. The intuition behind Mi-Lasso is that when the spatial correlation in the residuals is low, only a small set of eigenvectors will be necessary, so a high level of regularisation is required, and vice versa for a high level of residual spatial correlation. Mi-Lasso has several advantages; the method is (i) intuitive, (ii) theoretically grounded, and (iii) substantially faster than Lasso with $K$-fold cross-validation (CV) or the stepwise iterative greedy algorithms suggested in the literature.\footnote{Mi-Lasso only requires estimating a single point on a Lasso path, unlike $K$-fold cross validation which requires estimating $K$ paths. The larger $K$, the more computationally demanding the procedure.}
We establish the theoretical properties of Mi-Lasso by formalising the implicit ESF assumption that the terms which include the SWM can be approximated by a subset of eigenvectors. Under some standard spatial regularity conditions, we then derive non-asymptotic bounds for the coefficients of the eigenvectors and also assess the additional conditions required for Mi-Lasso to yield consistent eigenvector selection. In addition, given that the spectral decomposition of a square matrix power $\bm{A}^p$ always produces the same eigenvectors $\forall \; p \in \mathbb{Z}^+ $, we show that ESF also handles the case where the unknown spatial process possesses higher-order lags of the SWM. Simulations confirm that Mi-Lasso performs well for a range of levels of spatial correlation and when the data-generating process includes higher-order lags. Regarding computational time, Mi-Lasso is at least an order of magnitude faster than CV-Lasso.\footnote{Our setup explores sample sizes up to $n=10^4$, at which point the forward stepwise procedures become infeasible.}
Finally, we examine the practical performance of Mi-Lasso with an empirical application using the Boston Housing Dataset. We find that Mi-Lasso selects more than triple the number of eigenvectors compared to existing procedures. However, Mi-Lasso gives a bitter fit of the data in terms of adjusted $R^2$ and has substantially fewer insignificant eigenvectors than other selection procedure considered. Mi-Lasso is also over 60 times faster than the alternative selection procedures for this application.
The rest of this paper is organised as follows, Section (ref) describes the underlying model. Section (ref) discusses the statistical aspects of ESF and looks at existing methods for the ESF eigenvector selection problem. Section (ref) presents the Mi-Lasso procedure and derives several theoretical results. Section (ref) provides a Monte Carlo study comparing Mi-Lasso to the main existing selection procedures. Section (ref) tests the proposed method in an empirical application on house prices. Finally, Section (ref) offers our concluding remarks.
Consider the following equation, where the endogenous $n\times 1$ vector $\bm{y}$ is specified as a function of an $n\times k$ matrix of exogenous regressors $\bm{X}$ and follows some spatial process:
where $\bm{\beta}_0$ is the $k \times 1$ parameter vector of interest and $f(\bm{W,y,X,r})$ is a linear-in-parameter function of an $n\times n$ SWM of known constants $\bm{W}$,\footnote{We allow for the $\bm{W}$ to be normalised by a scalar factor as it allows for the recovery of the original autoregressive parameters kelpr10 and maintains symmetry.} $\bm{y}$, $\bm{X}$ and an $n\times 1$ vector $\bm{r}$. One example of such a model is:
where $\bm{\psi}_0$, $\rho_{i,0}$'s and $\delta_0$ describe the degree of spatial correlation in each of the $k$ exogenous variables, the dependent variable and error term. Note, simpler spatial models can be recovered by setting the spatial parameters $\rho_{i,0}$'s, $\delta_0$, and/or $\bm{\psi}_0$ equal to zero, and most spatial models set $p=1$. If the DGP of $\bm{y}$ is (ref) and (ref) then the reduced form for $\bm{y}$ is
if both $\bm{S}_1\equiv(\bm{I}-\sum^p_{i=1}\bm{W}^i\rho_{i,0})$ and $\bm{S}_2\equiv(\bm{I}-\delta_0\bm{W})$ are non-singular.
The SWM $\bm{W}$, with typical element $w_{ij}$, describes the spatial or socio-economic relationship between the cross-sectional units. When $w_{ij}\neq 0$, there is a meaningful interaction of units $j$ on unit $i$. In such cases, unit $j$ is often referred to as a neighbour of unit $i$. These interactions can stem from various sources, such as spillovers, externalities, geographic location, regulations, technology, government policy, or government expenditure. We further assume $\min_i\sum_{j=1}^nw_{ij}>0$ with probability 1, $w_{ii}=0$ by construction and $w_{ij}=w_{ji}$. The variables $\bm{Wr}$, $\bm{WX}$ and $\bm{W}^i\bm{y}$ are typically referred to as first order spatial lags of $\bm{r}$ and $\bm{X}$ and $i$th order spatial lags of $\bm{y}$.
Let $N$ denote the set of observations $N_n=N=\{1,\ldots,n\}$. All variables are normalised, as the transformed model is estimated by a Lasso-based procedure. For reasons of generality, we allow the elements of $\bm{u}_n$, $\bm{y}_n$, $\bm{W}_n$ and $\bm{X}_n$ to be dependent on $n$, that is to form triangular arrays, however, to simplify the notation we omit the $n$ index. Our analysis is conditioned on realised values of $\bm{X}$ and $\bm{W}$. We consider higher-order spatial lags only as powers of the SWM $\bm{W}$ and we allow the number of lags $p$ to be unknown.\footnote{More recent papers studying the estimation of higher-order spatial models, have generalised the concept of a higher-order spatial lag to allow for $p$ different weights matrices, thus, replacing $\bm{W}^i$ with $\bm{W}_i$ in (ref). Powers of $\bm{W}$ are viewed as a special case. Some examples are ll10_HOgmm, be13_HOgmm, gr15_HOpinc, gr18_HOml, g19_HOstoW, bde20_HOpanel, hlx21_HObay, g18_ho, g21_HOeff, gq22_spectest.} Even if $p$ is known, the estimation of such a model is non-trivial, as shown by b85_RCR. When the SWM is binary, powers of the SWM can result in the presence of circular and redundant routes. Proper higher-order spatial lags need to have these circular and redundant routes eliminated.\footnote{Both bk92_algRCR and as96_algRCR introduced algorithms to construct proper higher-order spatial lags.}
We now make the following assumptions about variables in Equation (ref)
Assumption (ref).1-(ref).3 are standard assumptions in the spatial econometrics literature kelpru98,kp99, lee04. Assumption (ref).1 (a) is required for the spectral decomposition. Assumption (ref).1 and (b) is necessary to limit the degree of dependence in $\bm{y}$. Given Assumption (ref).1 (a) if the true model is (ref)-(ref) and $\bm{W}$ is normalised by the largest eigenvalue then invertibility of $\bm{S}_1$ and $\bm{S}_2$ holds if $\sum^p_{i=1}|\rho_{i,0}|<1$ and $|\delta_0|<1$. Assumption (ref).2 ensure that the Gram matrix $\bm{X}'\bm{X}/n$ is invertible. Assumption (ref).3 requires the errors to be sub-Gaussian, this assumption allows us to derive a probability for the Lasso tuning parameter dominating the noise of the model. The finite fourth moment is needed for the selection consistency proof.
We now show how eigenvectors from a spectral decomposition of $\bm{W}$ can be used to spatially filter the model described in Section (ref). As $\bm{W}$ is a real and symmetric matrix (by Assumption (ref).1 (a)) the spectral decomposition of $\bm{W}$ is given by
where $\bm{E}$ is an $n\times n$ matrix of the $n$ eigenvectors $\bm{e}_{i\in N}$ and $\bm{\Lambda}$ is a $n\times n$ diagonal matrix of the $n$ eigenvalues ($\lambda_{i\in N}$) from $\bm{W}$. It is also important to note that the matrix of eigenvectors $\bm{E}$ of $\bm{W}$ is also the matrix of the eigenvectors of $\bm{W}^i \;\; \forall \; i\in \mathbb{Z}^+$. The proof is very simple, multiplying (ref) by $\bm{W}$ and using the orthogonal nature of the eigenvectors $\bm{E}$ to substitute $\bm{E}'\bm{E} = \bm{EE}' = \bm{I}$:
Recursive application of (ref) for any $p$ results in $\bm{W}^p=\bm{E}\bm{\Lambda}^p\bm{E}'$.
The intuition behind ESF is to use individual eigenvectors $\bm{e}_{i\in N}$ as explanatory variables to proxy for $f(\bm{W,y,X,r})$, yielding a high dimensional reduced form model:
where $\bm{E\gamma}_0$ can be viewed as a linear approximation of $f(\bm{W,y,X,r})$. The key problem with (ref) is that it is ghigh-dimensional and cannot be estimated consistently by OLS as the assumption that the regressor matrix $\bm{G}=[\bm{X},\bm{E}]$ has full column rank is violated.\footnote{This is because of $\operatorname*{rank}(\bm{G})=\operatorname*{rank}(\bm{G}'\bm{G})\leq \min(n,(n+k))$.} This implies a rank-deficient Gram matrix $\bm{\bm{G}'\bm{G}}/n$ with zero-valued eigenvalues. To handle this problem, we make the following assumptions:
Assumption (ref).1 is a weak sparsity assumption, and Assumption (ref).2 is required for the ESF approximation to be valid. While strong and untestable, they formalise the intuition of grif2000,grif03, who argue only a specific subset of eigenvectors ($\bm{E}_\Omega$) are related to the dependent variable $\bm{y}$ and will have non-zero coefficients. These assumptions imply (ref) can be reduced to the following low-dimensional equation, where $\bm{\Upsilon}_0=[\bm{\beta}_0,\bm{\gamma}_\Omega]'$ and $\bm{G}_\Omega=[\bm{X},\bm{E}_{\Omega}]$.
In principle, (ref) can be estimated by OLS. However, as $\bm{E}_{\Omega}$ is unknown, this is infeasible in practice. Thus, we now have a selection problem.
The ESF method of grif2000 is based on the Moran's $I$ statistic for spatial autocorrelation moran50. The test statistic for the Moran's $I$ ($m$) on the regression residual $\bm{M_Xy}=\hat{\bm{u}}$ of $\bm{y}=\bm{X\beta}+\bm{u}$ where $\bm{M_X}=\bm{I}-\bm{X}(\bm{X}'\bm{X})\bm{X}'$ is given by:
where $\bm{W}$ a $n\times n$ real symmetric SWM.\footnote{The assumption of symmetry of the elements of $\bm{W}$ is maintained w.l.o.g. since $\hat{\bm{u}}'\bm{W}\hat{\bm{u}} = \hat{\bm{u}}'[(\bm{W}+\bm{W}')/2]\hat{\bm{u}}$ kelpru01mI.} Substituting (ref) in (ref):
dejong84 showed that range of $m$ is determined by the maximum and minimum eigenvalues of $\bm{M_XWM_X}$. tiefboot95 showed that each of the $n$ eigenvalues of this expression represents a distinct $m$ values and all other possible $m$ values are just linear combinations of these $n$ values boot00.
It is important to note that the numerator of $m$ includes $\bm{E}'\hat{\bm{u}}$, which given the orthogonality of eigenvectors is the OLS coefficient estimate from a regression of $\hat{\bm{u}}$ on $\bm{E}$.\footnote{In other words, $\bm{E}'\hat{\bm{u}} = (\bm{E}'\bm{E})^{-1}\bm{E}'\hat{\bm{u}}$} grif03 argues that each of the $n$ eigenvectors represents mutually orthogonal spatial patterns and only a subset of eigenvectors will be relevant to the model, i.e., in a regression framework only a subset of eigenvectors will have non-zero coefficients.
Running the ESF method requires identifying $\bm{E}_\Omega$, the relevant set of eigenvectors. The first type of procedures proposed were forward stepwise greedy algorithms where eigenvectors are iteratively added until some user-specified threshold is reached grif2000,grif03,tiefgrif07, md19. grif03 proposed iteratively adding eigenvectors in a greedy manner to the base regression
until the spatial correlation in the OLS residual $\hat{\bm{u}}$ falls below a pre-specified level. Selection criteria based on alternative statistics such as the adjusted-$R^2$, the Akaike Information Criterion or Bayesian Information Criterion have also been suggested tiefgrif07, md19. tiefgrif07 specifically suggest using the standardised Moran's $I$ as the criterion for the greedy algorithm, based on its power against a wide array of autoregressive models and residual distributions ar91 and the fact it can be used for small samples kel17. The standardised version of Moran's $I$ statistic ($Z$) on the residual $\hat{\bm{u}}$ is:\footnote{Note the matrix $\bm{X}$ in the orthogonal projection matrix $\bm{M_X}$ may also include the selected eigenvectors in tiefgrif07 procedure.}
with
and
The greedy algorithm iterates over the candidate set of eigenvectors $\bm{E}_{c}$, searching for the eigenvector that minimizes $Z$. The selected eigenvector $\bm{e}_{i\in N}$ is then removed from $\bm{E}_{c}$ and added to the design matrix of (ref), and the residuals $\hat{\bm{u}}$ of this updated regression are tested to check if $|Z|< \epsilon$, where $\epsilon$ is a pre-specified threshold level of $Z$, which they suggest should be dependent on the sample size $n$.\footnote{tiefgrif07 suggest if $n<50$ then $\epsilon\approx1.0$ and if $n\approx 500$ then $\epsilon\approx0.1$.} If the condition is satisfied the iterations stop, if not the algorithm continues searching in the remaining candidate eigenvector set $\bm{E}_{c}$, with this iterative process continuing until $|Z|< \epsilon$.
grif03 argues that the candidate eigenvectors $\bm{E}_{c}$ form a subset $\bm{E}_{c} \subseteq \bm{E}$ of the full set of eigenvectors, based on several criteria. First, if $\bm{y}$ exhibits positive global spatial autocorrelation then $\bm{E}_{c}$ should be restricted to those eigenvectors with associated positive eigenvalues, as these are associated with at least weak positive spatial autocorrelation. Second, eigenvectors with small eigenvalues should be excluded from $\bm{E}_{c}$, suggesting a minimum threshold eigenvalue of 0.25, which is related to only approximately 5% of the variation attributed to spatial correlation in the dependent variable.
These forward stepwise procedures, through intuitive, have several key disadvantages. First, a lot of parameters are left to the user's discretion, such as, which statistic or information criterion to use, what threshold $\epsilon$ to use, which eigenvectors to include in the initial $\bm{E}_{c}$, and in which order to add the eigenvectors. Second, these greedy algorithms could also be at risk of data mining, with estimated models falling victim to over-fitting. Third, all these approaches are heuristics that aim to simplify the original, and infeasible, subset sum problem; therefore the solutions they obtain will be sub-optimal, with no guarantee they are close to the optimal one. Finally, these sequential methods carry a large computational burden, which becomes more acute when $n$ is large. This can be mitigated by limiting $\bm{E}_{c}$ with the rules of thumb mentioned above, but again with no guarantee these rules will consistently recover $\bm{E}_\Omega$.
This motivates seya15 to propose using Lasso t96, which shrinks many of the coefficients to zero, and can thus be used for variable selection tibs09. seya15 use Lasso under the assumption the parameter vector $\bm{\gamma}_0$ is sparse and the matrix of regressors $\bm{X}$ has full column rank, so that only the $\bm{\gamma}$ vector is penalised. The resulting Lasso estimator is:
where $\theta>0$ is the Lasso regularization or tuning parameter. Equation (ref) defines a family of estimators indexed by the tuning parameter $\theta$, a hyperparameter that ultimately determines which eigenvectors the Lasso selects.
seya15 proposed using $k$-fold cross-validation (CV) combined with the Brent algorithm brent73 to estimate $\hat{\theta}$, with prediction accuracy as the loss function. The Brent algorithm is a root-finding algorithm that allows for the optimisation to be non-convex: the algorithm first tries inverse quadratic interpolation in an attempt to achieve faster convergence which works well if the optimisation is convex. If it is non-convex and inverse quadratic interpolation fails, (slower) linear interpolation is used instead. CV using the Brent algorithm is the most time-consuming part of the seya15 Lasso procedure. Because the theoretical results on CV-Lasso hinge on the assumption that the cross-sectional units are independent czc20, it is hard to justify their validity for ESF, where eigenvectors are derived from a matrix that encodes cross-sectional dependence.\footnote{CV procedures do exist for cross-sectionally dependent data but they need to be carefully designed, for example see llz20.}
Some other methods have also been proposed. lpz13sea suggest simply including the first $j$ eigenvectors (sorted by eigenvalue magnitude) where $j$ is simply based on the sample size. Given this fixed rule, lpz13sea finds the quality of the ESF approximation is sensitive to the underlying spatial processes. ch16ESFsel argue more eigenvectors are needed when the level of spatial correlation is high compared to when the level of spatial correlation is low, thus, simple rules based on for example sample size may result in a sub-optimal set of eigenvectors being selected. ch16ESFsel instead develop the following eigenvector selection rule via simulation:
where $n_{pos}$ denotes the number of eigenvectors that exhibit positive spatial correlation (eigenvectors with positive eigenvalues). Equation (ref) was generated from a limited simulation that assumed the DGP has just spatial autoregressive disturbances, ch16ESFsel do not evaluate how their rule performs when the DGP follows some other spatial process.
The Lasso estimates are ultimately determined by tuning parameter $\theta$. Supposing $\theta=0$, the Lasso solution reduces to the OLS solution, whereas with a sufficiently large $\theta$ the penalised parameter vector is shrunk to zero (no eigenvectors selected). More moderate values of $\theta$ will result in some parameters being shrunk towards zero and some to precisely zero. As outlined above, the goal of ESF is to eliminate spatial correlation patterns in a linear regression framework. Information about these patterns will be contained in the regression residuals $\hat{\bm{u}}$, and we propose using these to determine a point estimate for $\theta$.
It seems reasonable to assume that when the level of spatial correlation in the residuals is low, only a small set of eigenvectors is necessary. Thus, a high level of regularization (value of $\theta$) is required. In contrast, when the level of spatial correlation is high, a large set of eigenvectors will be necessary. Thus, a low level of regularization (value of $\theta$) is required. Following tiefgrif07 we propose using the standardised Moran's $I$ (ref) to measure the spatial correlation of the residuals due to the previously mentioned properties. As $Z$ takes on large values when the correlation is high and small values when the correlation is low, we propose using the inverse of the square of $Z$ from the residuals of (ref) as a point estimate of $\theta$,
The square is chosen to ensure the tuning parameter is always positive.\footnote{A positive tuning parameter is necessary to ensure Lasso gives a unique solution.} The proposed estimator is called Moran $I$ Lasso (Mi-Lasso) and is outlined in Algorithm (ref).
As Lasso is a shrinkage estimator, it induces a downward bias on the estimated non-zero coefficients. Post-Lasso (pLasso) uses the Lasso estimator as selection procedure (assuming Lasso selects the correct variables), and then OLS is applied to the model selected by Lasso, straightforwardly providing unbiased estimates and standard errors.\footnote{For formal results on Post-Lasso, see bell13.} The Morans' $I$ Post-Lasso (Mi-pLasso) estimator is defined as:
To focus the theoretical analysis on the parameter vector $\bm{\gamma}$, we use the Frisch-Waugh-Lowell (FWL) partial regression theorem to partial out the $\bm{X}$ matrix. tibt11 and Lassofwl show that the FWL theorem could be used in a low-dimensional Lasso setting. Lemma (ref) shows that the FWL theorem can also be applied to the high-dimensional case of Mi-Lasso.
The proof is provided in appendix (ref).
We now introduce the following additional notation in the design. Without loss of generality, let $\bm{C}_{\Omega\Omega}=n^{-1}\tilde{\bm{E}}_{\Omega}'\tilde{\bm{E}}_{\Omega}$, $\bm{C}_{\Omega\grave{\Omega}}=n^{-1}\tilde{\bm{E}}_{\Omega}'\tilde{\bm{E}}_{\grave{\Omega}}$, $\bm{C}_{\grave{\Omega}\Omega}=n^{-1}\tilde{\bm{E}}_{\grave{\Omega}}'\tilde{\bm{E}}_{\Omega}$ and $\bm{C}_{\grave{\Omega}\grave{\Omega}}=n^{-1}\tilde{\bm{E}}_{\grave{\Omega}}'\tilde{\bm{E}}_{\grave{\Omega}}$ where $\tilde{\bm{E}}_{\Omega}$ is an $n\times s$ matrix with columns corresponding to the active set $\Omega$. $\grave{\Omega}$ is the complement set and the $n\times q$ matrix $\tilde{\bm{E}}_{\grave{\Omega}}$ is defined accordingly with $q_n=q=s-n$. Now the (re-scaled) Gram matrix $\bm{C}_n=\bm{C}=n^{-1}\tilde{\bm{E}}'\tilde{\bm{E}}$ can be expressed in block-wise form as:
Similarly we define $\bm{\gamma}=[\bm{\gamma}_{\Omega},\bm{\gamma}_{\grave{\Omega}}]'=[\gamma_1,\ldots,\gamma_{s},\gamma_{s+1},\ldots,\gamma_{n}]'$.
This section produces performance bounds for the Mi-Lasso estimates of $\bm{\gamma}$. Given the high-dimensional structure of ESF, the Gram matrix $\bm{G}'\bm{G}/n$ is singular. This implies its minimum eigenvalue will be zero. However, as shown by lasdig09 for the case of Lasso, the following restricted eigenvalue (RE) condition only requires the appropriate sub-matrix of the Gram matrix to have positive and finite eigenvalues.
Assumption (ref) requires that $\bm{\Delta}$ lies within the restricted set (ref). As $\bm{\Delta}$ is the difference between the estimate $\tilde{\bm{\gamma}}$ and the true parameter $\bm{\gamma}_0$, the restricted eigenvalue bounds the minimum change in the prediction norm from a deviation $\bm{\Delta}$ within the restricted set $\mathcal{C}(\Omega,\bar{b})$ relative to the norm of the deviation on the true support $\bm{\Delta}_{\Omega}$.
By combining Assumptions (ref) and (ref) with the RE condition, and treating $\bm{X}$ and $\bm{E}$ as constants (realisations) we can now establish the $\ell_1$ and $\ell_2$ parameter norm bounds and the $\ell_2$ prediction norm bound for the Mi-Lasso estimates of $\bm{\gamma}$.
The proof is provided in appendix (ref).
The three convergence rates presented in Theorem (ref) depend on the number of eigenvectors with non-zero coefficients, the sample size, and $Z$. They also require that the tuning parameter dominates the noise of the model. By assuming the errors are sub-Gaussian (Assumption (ref).3) we prove the probability of this event occurring goes to one as $n\to \infty$ (see proof for further details).
Corollary (ref) is satisfied if $Z=O_p(1)$,, which is reasonable as $Z$ is a measure of correlation, and $s$ grows at a rate slower than $n$, which is satisfied by Assumption (ref).4 below.
This section shows the conditions required for Mi-Lasso to consistently select the non-zero and zero elements in $\bm{\gamma}$. Following zhao06, we say that $\tilde{\bm{\gamma}} =_s \bm{\gamma}_0$ if and only if $\operatorname*{sign}(\tilde{\bm{\gamma}}) = \operatorname*{sign}(\bm{\gamma}_0)$ where $\operatorname*{sign}(\cdot)$ maps positive entry to 1, negative entry to -1 and zero to zero. We now define selection consistency for Mi-Lasso as
The following assumptions are required to prove sign consistency of Mi-Lasso.
Assumption (ref).1 is a normalisation of the transformed eigenvectors. Assumption (ref).2 bounds the eigenvalue of the eigenvectors with non-zero coefficients from below, so the inverse of $\bm{C}_{\Omega\Omega}$ is well behaved. Assumption (ref).3 and Assumption (ref).4 are important as they ensure convergence in the high dimensional space as $n\to\infty$. Assumption (ref).3 ensure there is a difference of size $n^{c_2}$ between the decay rate of $\bm{\gamma}_{\Omega}$ and $\sqrt{n}$, preventing the estimates from being dominated by the disturbance terms, which aggregate at a rate of $n^{-1/2}$. Assumption (ref).4 is a sparsity assumption that requires the square root of the size of the true model $\sqrt{s}$ to increase at a slower rate than the rate difference, preventing the Lasso estimation bias from dominating the model parameters. Assumption (ref).5 (assuming $\bm{C}_{\Omega\Omega}$ is invertible) is the Irrepresentable Condition (IC), which is the necessary condition for the consistency of Mi-Lasso selection, the inequality holds element-wise. The IC requires the correlation between the relevant and irrelevant eigenvectors to be zero or weak. In the Mi-Lasso framework, this is likely to be satisfied as the columns of $\bm{E}$ are mutually orthogonal. The columns of $\tilde{\bm{E}}$ may not be, however, as the eigenvectors are projected into the column space of $\bm{X}$. Unfortunately, in practice, the IC is impossible to verify as we do not know the true parameter vector $\bm{\gamma}_0$.
The following proposition places a lower bound on the probability of Mi-Lasso picking the true model, which quantitatively relates to the probability of Lasso selecting the correct model. Proposition (ref) is a modification of Proposition 1 in zhao06.
The proof is provided in appendix (ref).
Proposition (ref) shows that the measure of spatial correlation $Z$ determines the size of the trade-off between events $A$ and $B$. A higher level of spatial correlation will lead to larger $A$ but smaller $B$; this makes Mi-Lasso more likely to select irrelevant eigenvectors. In contrast, a larger $\nu_i$ has no impact on $A$ but leads to a larger $B$. So when IC holds with a large $\nu_i$, Mi-Lasso is more likely to select the correct model.
The proof is provided in appendix (ref).
Theorem (ref) shows that Mi-Lasso is consistent in selecting the true model if the 4$th$ moment of the errors is finite (Assumptions (ref).3), Assumptions (ref)-(ref) hold and the difference between $c_2$ and $c_1$ is 0.5. The greatest difference (between $c_2$ and $c_1$) for which Mi-Lasso is consistent is 0.5, smaller differences can also yield consistency, but this would require higher order moments of the errors to be finite. For example, if we assume the 6th or 8th moment is finite, the difference would need to be $1/3$ or 0.25 for Mi-Lasso to be consistent (see proof for further details).
To evaluate the finite sample performance of Mi-Lasso and compare it to the main existing selection procedures, we conduct two Monte Carlo exercises where the DGP is,
In both simulations, we set the `true' parameter value of $\beta = 1$ and $\psi=0.9$. The elements of $\bm{W}$, denoted $w_{ij}$, are independent draws from a Bernoulli distribution with success probability $\mu/n$ for some constant $\mu < \infty$, $w_{ii}=0$ and $w_{ij}=w_{ji}$. By construction, $\mu$ is the expected number of links for each unit, and we set $\mu\in\{4,8,12\}$. Each $\bm{W}$ is normalised by the maximal of the row (or column) sum. Sample sizes considered are $n\in\{100,250,500\}$, and we run 1000 replications.
In setup A, we set $p=1$ so we can evaluate how the method performs with different levels of spatial correlation $\rho_1\in \{0.3,0.4, 0.5, 0.6, 0.7, 0.8,0.9\}$. We consider only positive spatial correlation as this is the most common setting. In setup B, we set $p=3$ to evaluate the performance of ESF in the presence of higher-order spatial lags. In both setups the estimators compared are:\footnote{An oracle estimator is not possible hare as this requires knowledge of $\operatorname*{supp}(\bm{\gamma}_0)$, which is unknown.}
Figures (ref), (ref) and (ref) show the bias, MSE, and the number of selected eigenvectors for setup A,\footnote{Figures (ref) and (ref) are provided in appendix (ref)} revealing the different selection behaviours of these estimators. For CV-Lasso the number of selected increases very slightly as the levels of spatial correlation in the dependent variable increases, and this pattern is consistent across different sample sizes and $\mu$. In contrast, FstepZ selects more eigenvectors when the spatial correlation level is low than high for small sample sizes and the largest set of eigenvector when the level of spatial correlation in the dependent variable is small. Mi-Lasso behaviour is as expected from the intuition of the procedure, selecting a small set of eigenvectors when the level of spatial correlation is low and a large set when the level is high.
The Lasso estimators generally have a smaller bias and larger MSE than their post-Lasso (pLasso) counterparts. When the level of spatial correlation is high Mi-Lasso has the best performance in terms of bias and performs comparably to the other estimator in terms of MSE. Mi-pLasso has the smallest MSE when then level of spatial correlation is high and comparably well when the level of spatial correlation is low. Notably, FstepZ has the largest MSE when the sample size is 100 all levels of spatial correlation and $\mu$ considered and when the level of spatial correlation is low for other sample sizes. FstepZ performance in terms of bias and MSE improves as the sample size increases and the SWM becomes more dense. Generally, in terms of bias, the estimators diverge as the level of spatial correlation increases, this is because the bias is determined by an interaction between the level of $\rho_1$ and the structure of the SWM. Thus, for a given SWM, the larger $\rho_1$ the larger the bias, so mistakes/variation in selection can have a larger effect.
For setup B, we set $\rho_1=0.6$, $\rho_2=0.4$ $\rho_3=0.5$. Table (ref) shows the bias, MSE, and the number of selected eigenvectors for setup B. This table confirms that ESF can work well in the presence of higher-order spatial lags. This table shows that Mi-Lasso selects less eigenvectors as the density of the SWM ($\mu$) increases. Mi-Lasso and Mi-pLasso always has a smaller bias and a comparable MSE than CV-Lasso and CV-pLasso. FstepZ generally performs better in terms of both bias and MSE as the sample size increases.
Finally, Table (ref) shows the computational times of the different estimators used in the simulations. These results show Mi-Lasso is the fastest procedure, CV-Lasso is the second fastest, and FstespZ is the slowest procedure for a given sample size. Comparing Mi-Lasso to CV-Lasso, we find Mi-Lasso is up to 37 times faster. The most substantial computational gains are found when the sample size is 1000, but even when the sample size is very large (10,000), Mi-Lasso reamins 19 times faster than CV-Lasso, with FstepZ becoming unfeasible.
We now compare the ESF selection procedures using the Boston Housing Dataset, which was first used by hr78_bh to evaluate the relationship between house prices and demand for clean air. gp96_sbh later revisited the dataset when they noted the high spatial correlation in the dataset and proposed estimating a spatial error model instead. However, as there is no guarantee theirs is the correct specification, and given that the researcher is only concerned with the direct effect, ESF is an appropriate methodology allowing to simply control for the spatial effects.
The dataset includes 508 census tracts (spatial units). Table (ref) describes the variables used in the analysis. The eigenvectors are from a binary SWM where the tracts are connected if they share a border, and SWM is normalised by the maximal of the row (or column) sum. The following basic model (excluding the eigenvectors) is:
Table (ref) shows the parameter estimates (excluding eigenvectors) for OLS, which ignores the spatial correlation, Mi-pLasso, CV-pLasso, and FstepZ. These results show that some of the OLS estimates are biased by spatial dependence. For example, age had a positive (but insignificant) coefficient when the spatial dependence is ignored, but in the filtered estimates, the coefficient is negative and significant as expected; the coefficient on $nox^2$, $dis$, and $rm$ also have a downward bias. Additionally, the filtered estimates also give a substantially better fit of house prices, with Mi-pLasso having an adjusted R$^{2}$ of 0.978, implying an almost perfect fit of the data. Mi-pLasso standard errors are generally the same or smaller than the other estimator.
Table (ref) shows the computational times, the number of selected eigenvectors, and their significance levels, for the three ESF estimators. There is substantial variation in the number of selected eigenvectors between the procedures. Mi-Lasso selected over four and three times more eigenvectors than CV-Lasso and FstepZ. However, despite selecting substantially more eigenvectors for Mi-Lasso, only 0.5 percent of selected eigenvectors are insignificant compared to 28 percent and 21 percent for FstepZ and CV-Lasso. Mi-Lasso has more eigenvectors with coefficients significant at the 0.1 percent level than FstepZ or CV-Lasso selected in total, implying these techniques may be under-selecting in this case. Mi-Lasso is also over 65 times faster than both FstepZ and CV-Lasso.\footnote{The code to replicate the results in the section can be found in the attached filed `boston_comp.R'.}
In this paper we have formalised the ESF assumptions and evaluated the existing solutions to the ESF eigenvector selection problem. Our analysis of existing procedures has shown that a dominant selection procedure currently does not exist. The forward-iterative procedures with a user-defined cut-off and eigenvector inclusion criterion can be viewed as ad hoc and are slow, especially as the sample size increases. seya15 proposed using Lasso with prediction accuracy CV to estimate the tuning parameter. However, as ESF aims to reduce bias on $\bm{\beta}_0$ rather than improve prediction accuracy, it is unclear if this is the best way to estimate the tuning parameter. Additionally, CV-based Lasso procedure is also slow, especially when $n$ is large.
We have proposed an alternative Lasso-based procedure called Morans’ $I$ Lasso (Mi-Lasso) that uses information about the level of spatial correlation in the naïve regression residuals to determine a point estimate for the Lasso tuning parameter instead of using CV. The key benefits of Mi-Lasso are that it is intuitive, theoretically grounded, and substantially faster than seya15 CV Lasso or stepwise procedures and can thus be implemented on large data sets. We have derived performance bounds for the Mi-Lasso estimates of the eigenvectors coefficients and shown the conditions necessary for the estimator to provide consistent eigenvector selection. Our simulation results confirm the estimator performs well in terms of bias and MSE compared to existing selection procedures for a range of levels of spatial correlation and in an empirical application on house prices. Additionally, we have shown using a property of the spectral decomposition and a simulations experiment, that ESF is robust to the presence of an unknown number of higher-order spatial lags in underlying DGP.
A key limitation of the ESF literature is that there are no results on constructing robust standard errors. As all the proposed procedures can be viewed as post-model selection estimators. Thus, all the corresponding estimators suffer from the corresponding post-model selection inference problem lp08. Given the spatial dependence in the model, debiasing techniques such as Double Lasso bch14 or Partial Lasso CHS15 will not work well. A promising avenue of future research in the ESF literature is to extend Mi-Lasso (and other procedures), so standard errors robust to selection mistakes and the spatial dependence in the model can be calculated.
\noindentConflict of Interest Statement: the authors declare no conflicts of interest