7 Spatial Regression

7.1 Introduction

Aims

The aims of this practical are to:

1. Understand how spatial regression is different from classical (aspatial) regression and why it matters.
2. Understand the differences between common forms of spatial regression.
3. Use spatial regression to reveal relationships in a sociological context.

Application

To that end, we’ll apply aspatial and spatial regression to predict life expectancy at the MSOA-level (Middle Layer Super Output Area) in the North West of England.

Data

The data for this practical are sourced from the Department of Health and Social Care (DHSC, Fingertips) and the Office for National Statistics (ONS), and are available in [data/dhsc] and [data/ons/msoa] respectively. Key attributes have been retained and are described below:

• ONS MSOA-level geometries and lookup table [Source]
• ONS Economic Activity Status⁷³ [Source]
• ONS Highest level of qualification⁷⁴ [Source]
• ONS Disability⁷⁵ [Source]
• ONS Household income⁷⁶ [Source]
• DHSC Index of Multiple Deprivation⁷⁷ [Source]
• DHSC Fuel poverty⁷⁸ [Source]
• DHSC Life Expectancy⁷⁹ [Source]

Tools

Ordinary Least Squares (OLS) (Spatial Statistics), Spatial Autoregression (Spatial Statistics)

---

7.2 Practical

7.2.1 Why spatial regression?

When analysing data, one of the most important questions we can ask is:

Is there an association or relationship between the variables of interest?

This question is foundational, found in the fields of epidemiology⁸⁰, ecology⁸¹, and earth science⁸², and everything in between. By answering these questions, we can better understand the mechanisms which link the variables of interest, make predictions about the future, and inform decision-making.

In Practical 4, we introduced a measure of association, Pearson’s \(r\) (Pearson, 1895), with values ranging from \(-1\) (perfect negative correlation) to \(+1\) (perfect positive correlation). As a refresher, association refers to any form of dependency between variables, whereas correlation is a specific form of association, focusing on the degree of linearity.

While measures of association are useful, they primarily provide information on the strength and direction of the correlation, but little about the nature of the relationship between the variables i.e., how much does \(Y\) change as \(X\) changes? We can address this by modelling the relationship between \(X\) and \(Y\) and producing an equation which represents this relationship. In doing so, we can develop a deeper understanding of how the variables are related, moving beyond just strength and direction, and can use our equation to make predictions i.e., what would happen to \(Y\) if we doubled \(X\)?

This is known as regression.

There are many types of regression methods, and perhaps the best known and most commonly used is ordinary least squares (OLS) regression. This approach is used to model the linear relationship between a dependent variable \(Y\) and an independent variable (or variables) \(X\), where the outcome is a best-fit regression line, as shown below.

This is produced by minimising the sum of the squares of the residuals⁸³ between the observed dependent variable \(Y\) and the predicted values of \(Y\), which are produced based on the regression equation. Other regression lines, as illustrated below, might reduce the residuals for some of the observations, but the total error would be larger.

The regression line itself can be defined as follows:

\[Y=\beta_{0}+\beta_{1}X_{1}+\beta_{2}X_{2}+\cdots+\beta_{i}X_{i}+\epsilon\]

where \(Y\) is the dependent variable of interest, \(X_{i}\) is an independent variable where \(\beta_{i}\) is the corresponding slope coefficient, \(\beta_{0}\) is where the regression line intercepts the Y-axis, and \(\epsilon\) is the random error (residual) between the observed dependent variable \(Y\) and the regression line.

Let’s use a simple example to illustrate, using the ecological data from Getzin et al. (2012) from Practical 4, which consists of paired measurements of the median gap shape complexity index (GSCI) \(X\) and species richness (SR) \(Y\) for a forest patch (Hainich) in central Germany. Here we have a single dependent variable we want to predict (SR) and a single independent variable (GSCI) which we believe to be related to SR in some way.

Here is a plot of the paired measurements, redrawn from Getzin et al. (2012):

You’ll remember that Pearson’s \(r=0.86\). As discussed above, this tells us about the strength and direction of the correlation i.e., there is strong positive correlation between median gap shape complexity index \(X\) and species richness \(Y\). However, we know little about the nature of the relationship, and could not use this measure of association for prediction e.g.,

What species richness would we expect for a median GSCI value of \(2.3\)?

This could be addressed via OLS regression, which for this dataset would produce the following best-fit regression line:

With only a single independent variable, the regression line is defined as:

\[Y=\beta_{0}+\beta_{1}X_{1}+\epsilon\]

Based on the minimising the sum of the squares of the residuals, the best-fit regression line is defined by a slope coefficient \(\beta_{1}\) of 14.150 and an intercept \(\beta_{0}\) of 5.595:

\[\text{SR} = 5.595 + 14.150\cdot\text{GSCI}\]

Prediction is now possible, as we can substitute in our median GSCI value of interest⁸⁴:

\[\text{SR} = 5.595 + 14.150\cdot2.3 = 38.14\]

Alongside the regression equation, it is also typical to calculate the coefficient of determination \(\text{R}^2\), which measures the proportion of variance in the dependent variable that can be explained by the independent variables in the model. This typically ranges from \(0\) i.e., none of the variance in the dependent variable is explained by the independent variables, to \(1\) i.e., the variance in the dependent variable is perfectly explained by the independent variables.

For this dataset \(\text{R}R^2=0.74\), which indicates that 74% of the variance in species richness can be explained by median GSCI⁸⁵.

However, there is a problem.

What property of the data are we not including in our regression?

Take some time to reflect on this question, thinking back to earlier parts of the unit, and then read the following discussion.

Discussion

We know from our analysis of correlation that Pearson’s \(r\) and other similar correlation metrics, such as Spearman’s rank correlation coefficient \(ρ\), are aspatial i.e., they do not take into account the spatial structure of the data. As a result, \(r\) and \(p\) provide only a partial perspective on the degree of correlation.

We also know from our exploration of spatial autocorrelation that we can often reject the hypothesis of spatial randomness, indicating the presence of clustering or dispersion of the data values.

In turn, we introduced a measure of bivariate spatial association Lee’s \(L\) (Lee, 2001), which incorporated both the numeric similarity of paired observations and the spatial structure of the data.

In much the same way, OLS regression is also aspatial, in that it incorporates the values for \(Y\) and \(X\) (e.g., species richness, GSCI), but with no explicit consideration of their geography.

This is important for a few key reasons.

First, the values for a dependent variable of interest \(Y\) often arise from spatial processes which need to be explicitly accounted for. As an example, in today’s practical we are going to model life expectancy for Middle Layer Super Output Areas (MSOAs), which comprise approximately 2,000 to 6,000 households, with a resident population of 5,000 to 15,000 persons. Life expectancy for each MSOA is undoubtedly influenced by the characteristics of the MSOA and its population e.g., level of income, deprivation, and fuel poverty. However, life expectancy might also be influenced by the characteristics of neighbouring MSOAs or nearby sites e.g., proximity to major roads, being downwind of polluting factories, access to nearby green- and blue-space, or the level of income in surrounding MSOAs, which might dictate the resources available for healthcare. Incorporating the geography of the MSOA, and its spatial dependence with its neighbours, may result in better predictions.

Incorporating geography can also be useful to minimise model error. When evaluating a regression, we would typically investigate the distribution of residuals and calculate summary statistics, such as RMSE, as introduced last week. The numerical distribution of residuals could be evaluated using a histogram, as shown below:

While this provides important information on the whether the model is systematically underpredicting \(Y\) (distribution \(A\)) or over-predicting \(Y\) (distribution \(C\)), it provides no information on where that systematic error is located. For example, perhaps our model systematically over-predicts life expectancy in certain regions, and under-predicts life expectancy elsewhere, indicating clustering of the model residuals i.e., positive spatial autocorrelation, as illustrated below (e.g., over-prediction in parts of southern Europe, under-prediction in northern Europe). Spatial structure of the errors might indicate the exclusion of important spatial variables. Moreover, if we can model and account for spatial clustering of the residuals (if present), we can produce models which minimise error.

Finally, and from a statistical perspective, one of the key assumptions of many tests, including OLS, is the independence of observations assumption. This states that each observation in your dataset (e.g., a species richness value) should not be influenced by or be related to any other observation in the dataset. Given Tobler’s Law (Tobler, 1970), and the evidence of spatial clustering and dispersion in many datasets, as revealed using Moran’s \(I\) and \(I_{i}\) (Moran, 1948; Anselin, 1995a), we know that observations are often spatially dependent e.g., the Republican vote share in a particular US county is related to the vote share of its neighbours. In these situations, the independence of observations assumption is violated, and the OLS regression coefficients and errors may be biased.

---

Based on the above reasoning, spatial regression methods are often more appropriate than OLS and similar aspatial regression methods when applied to spatial data, in that they incorporate the spatial structure and spatial dependence of the data, are a better representation of the spatial processes which produce the variables of interest, and can often minimise model error.

7.2.2 Pre-processing

To begin:

Open ArcGIS, initialise a new project in the practical-7 directory, establish a connection to data and then load the MSOA geometries for England and Wales, stored in MSOA_2021_EW_BGC_V3.shp.

Next, load and inspect the MSOA lookup table MSOA_2021_lookup_table.csv which contains additional attributes of interest e.g., region name [RGN22NM].

For our analysis, we are going to focus on the North West of England, which includes the counties of Cheshire, Cumbria, Greater Manchester, Lancashire, and Merseyside.

Use a dynamic table join via Joins and Relatives → Add Join to link the geometries and lookup table via the MSOA21CD attribute. When complete, use Select by Attributes to filter to the North West (RGN22NM = 'North West', \(n=932\)), and then export the selection to your project directory with a suitable name e.g., MSOA_2021_nw. To simplify the export, remove all the joined fields (Export Features → Fields → Field Map → Remove all below GlobalID).

When complete, remove the national MSOA geometries and lookup table.

The aim of our analysis is to understand and predict life expectancy at the MSOA-level i.e., our dependent variable. Our independent variables, which we expect to be related to life expectancy in some way, include economic activity, education, income, deprivation, fuel poverty and disability.

Add the datasets containing the dependent and independent variables (\(n=7\)) to ArcGIS. These can be found in [data/dhsc] and [data/ons/msoa].

Our next step is to join the important attributes to our MSOA geometry layer. However, we are using the well-established shapefile format, which limits field names to 10 characters. As we’ll working with a large number of fields, where the field name contains important information e.g., economically_inactive_and_ftstudent, we want to preserve these names, rather than see them truncated. This could be solved by assigning longer aliases to each field name, but a quicker workaround is simply to switch to a more modern format.

Use Export Features to save the MSOA geometries as a Feature Class in practical-7.gdb. When complete, remove the earlier version from ArcGIS, which you can check via Properties → Source → Data Source.

When complete:

Use Join Field to produce permanent table joins between the MSOA geometries and the standalone tables, joining based on MSOA21CD and msoa_code. Use “Transfer Fields” to join the important variables for each table, as listed below:

• [dhsc-life-expectancy.csv] = life_expectancy_male and life_expectancy_female
• [ons-household-income.csv] = total_annual_income
• [ons-economic-activity-status.csv] = employed_not_ftstudent, economically_inactive_not_ftstudent and economic_status_total
• [ons-highest-level-qualification.csv] = no_qualifications, level_four and qualification_total
• [ons-disability.csv] = disabled_perc
• [dhsc-imd.csv] = imd_decile
• [dhsc-fuel-poverty.csv] = fuel_poverty_perc

As MSOAs have varying populations, ranging from approximately 5,000 to 15,000 residents, we should take care to use proportional variables (e.g., disabled_perc, fuel_poverty_perc) or those that are not population-dependent (e.g., imd_decile, total_annual_income). For the economic activity status and qualification datasets, we can calculate proportional variables using the total number of individuals with an economic or qualification status (economic_status_total, qualification_total), focusing on the end-member classes i.e., economically active vs. inactive, no qualifications vs. level-four qualifications.

Use Calculate Field to produce proportional variables (“Field Type” = Double) representing the proportion of the population employed (employed_not_ftstudent_perc) and not economically active (economically_inactive_not_ftstudent_perc), and the proportion with no qualifications (no_qualifications_perc) and level-four qualifications (level_four_perc) e.g., employed_not_ftstudent_perc = (!employed_not_ftstudent! / !economic_status_total!) * 100

As we have permanently joined these attributes to our geometry layer:

Remove the standalone tables.

7.2.3 Data exploration

Before we model the relationship between our dependent and independent variables, we should always conduct an exploratory data analysis.

What is the average, minimum and maximum male and female life expectancy?

Modify the symbology to investigate the spatial patterning of life expectancy.

Are there are any clusters of high or low life expectancy? Why might this be?

There are lots of options for interesting exploratory analysis. For example, you could investigate the differences in life_expectancy_male and life_expectancy_female using Primary symbology → Bivariate Colours, or by calculating the difference using Add Field.

You could also use Create Chart → Scatter Plot to investigate the bivariate relationships between our dependent variable e.g., life_expectancy_male, and the independent variables e.g., total_annual_income, no_qualifications, disabled_perc, imd_decile, fuel_poverty.

It is worth emphasising that the plots and statistics (\(\text{R}^2\), coefficients \(B_{0}\) and \(B_{1}\)) are produced using OLS i.e., they are not taking geography into account. As we know from earlier weeks of the unit, including spatial correlation, there are often alternative methods that are more appropriate e.g., Lee’s \(L\).

The regressions are also bivariate, in that they are assessing the relationship between the dependent variable \(Y\) and a single independent variable \(X\). As life expectancy is likely influenced by a wide range of factors, a mutlivariate analysis, incorporating multiple independent variables, would typically produce the best model.

Building on that, we also need to remember key concepts, such as the designation of the spatial weights matrix and the importance of spatial autocorrelation, which we should always assess.

Use the Neighborhood Explorer to visualise the connections that underlie the spatial weights matrix.

Which spatial weights method is most appropriate for our data? Contiguity? Distance-based?

Interpretation

When selecting a spatial weights method, we need to consider both theory and the resulting spatial weights structure. The latter is relatively straightforward. For example, we would typically want to avoid any isolated spatial units i.e., those with no neighbours. We could also use row standardisation to account for uneven spatial data, where there may be significant differences in the number of neighbours between spatial units. An alternative approach would be to used a fixed number of neighbours \(K\) for each spatial unit.

These practical considerations are important, but more so is the theory of the nature of the interaction between spatial units for the specific geographic process.

For example, if we were focused on a geographic process which was dependent on contact between individuals e.g., the spread of disease, or the expansion of a species, then a contiguity-based approach would be most appropriate, as this is a representation of the nature of the interaction i.e., spatial units must be contiguous for disease to spread⁸⁶, or for a species to move from one spatial unit to another.

Alternatively, if the geographic process is dependent on accessibility, then a distance-based approach might be more appropriate. In the context of life expectancy, we might expect this to be influenced by the characteristics of nearby areas, such as the provision of hospitals, the proximity of sources of pollution, or the availability of high paying jobs within a certain commutable distance. These characteristics would affect life expectancy even if the spatial units were non-contiguous.

The geographic process could also be dependent on structural differences across the study area. Continuing with our life expectancy example, which we are analysing at the MSOA-level, there might be differences in how the health service is funded and managed at a higher spatial level i.e., Cheshire, Cumbria, Greater Manchester, Lancashire, and Merseyside, resulting in a different quality of care for different counties. In this situation, we might want to use block weights, which define neighbours based on the higher spatial level i.e., all MSOAs within Greater Manchester would be defined as neighbours, and so on.

To summarise, the choice of spatial weights method is critical and should be dependent on the nature of the interaction between spatial units for a specific geographic process. Different geographic processes necessitate different methods, even when the same spatial units are used (e.g., MSOA).

Assess for global and local spatial autocorrelation for some of the attributes (e.g., life_expectancy_male and life_expectancy_female) using Moran’s \(I\) (Spatial Autocorrelation (Global Moran’s I)) and local Moran’s \(I_i\) (Cluster and Outlier Analysis (Anselin Local Moran’s I)⁸⁷.

Is there evidence of spatial autocorrelation?

7.2.4 Ordinary least squares

As a first foray into regression, we will begin by running and critically engaging with Ordinary Least Squares.

Before we can do so, the ArcGIS implementation of Ordinary Least Squares (OLS) requires a unique ID field in our Attribute Table, which we don’t have, so we’ll create one:

Use Calculate Field to produce a unique_id field (“Field Type” = Long integer (32-bit)), using the SequentialNumber() expression, which can be found under “Helpers”.

Now we can run Ordinary Least Squares (OLS), using unique_id as the “Unique ID Field” and either life_expectancy_male or life_expectancy_female as the “Dependent Variable”. For simplicity, we’ll begin by using a single independent (Explanatory) variable: total_annual_income. Save the “Output Feature Class” with a suitable name, and if desired, generate an “Output Report File”, saving to the practical-7 directory.

Before we jump to interpretation, it is important to assess any warnings that appear, which in this case is as follows:

WARNING 000851: Use the Spatial Autocorrelation (Morans I) Tool to ensure residuals are not spatially autocorrelated.

This is very good advice and is described more fully here.

For now we’ll ignore this and inspect the results under “Messages”, which reveals \(\text{R}^2=0.62\) (Multiple R-Squared), matching the Scatter Chart presented above, as well as coefficients \(\beta_0\approx62.86\) and \(\beta_1\approx0.0003\), and p-values for the independent variable, alongside a wealth of information in the PDF.

What is the numeric distribution of the model residuals?

As discussed in the introduction, we can also inspect the spatial patterning of the model residuals to identify whether the model is systematically over-predicting or under-predicting life expectancy in different regions⁸⁸.

How does model performance vary across the study region? Is there evidence of clustering of the residuals?

A visual inspection of the residuals can be useful, but a better approach is to evaluate autcorrelation statistically.

Check for global and local spatial autocorrelation of the OLS residuals. Make sure to use a consistent spatial weights method, reflecting the geographic process of interest.

Is there evidence of global spatial autocorrelation?

What might explain the pattern of high (HH) and low (LL) clusters, or outliers (LH, HL)?

Interpretation

Based on Moran’s \(I\) there is clear evidence of positive global spatial autocorrelation (clustering, \(I>0\)), with HH clusters in parts of Cumbria and Lancashire, alongside LH outliers (e.g., Kendal, Morecambe). Conversely, there are LL clusters in the major cities of Manchester and Liverpool, with HL outliers in (presumably) more affluent areas (e.g., Didsbury, Chorlton).

As outlined in the esri documentation, statistically significant spatial autocorrelation of the regression residuals indicates misspecification of the OLS model, as well as a violation of OLS assumptions, including the independence of errors assumption i.e., the residual for one observation should not be related to or be predictable by the residuals of other observations.

In this situation, OLS is an inappropriate method and we should either use a spatial regression method which incorporates spatial dependence, or use methods which incorporate spatial non-stationarity⁸⁹.

In our first model, we have used a single independent variable (total_annual_income). As life expectancy is influenced by a range of factors, we might be able to explain significantly more variability by including a range of variables in a multivariate approach. If so, we need to be aware of issues of multicollinearity i.e., when independent variables are highly correlated with each other, see Kim (2019).

Re-run Ordinary Least Squares (OLS), including additional independent variables e.g., disabled_perc, imd_decile, and fuel_poverty_perc, and inspect the results table. You may also want to include economic and education attributes, and if so, make sure to use the proportional variables you calculated earlier.

How much additional variability in life expectancy can be explained with our new model?

Answer

When modelling life_expectancy_male, and including the full range of independent variables, \(R^2\approx0.82\), indicating 82% of variability in life expectancy can be explained.

There are lots of tools available to diagnose multicollinearity, one of which is the Variance Inflation Factor (VIF), which can be found in the “Summary of OLS Results” table. As a rule-of-thumb, multicollinearity is present when the VIF is higher than 5 to 10 (Kim, 2019).

Is there evidence of multicollinearity?

Answer

There is evidence of multicollinearity, as imd_decile and no_qualifications_perc have \(\text{VIF}>10\). This is not to suggest that these variables are unimportant, but it may be that their effects are being captured by other variables.

If multicollinearity is present, a typical strategy is to remove highly correlated variables.

Are all the independent variables significant in the model? i.e., \(p<0.05\).

Answer

Fuel poverty is classed as not statistically significant, with \(p\approx0.68\). As above, this does indicate that fuel poverty is unimportant, see Robinson et al. (2018) and Robinson et al. (2019), but that its effects is being captured by other variables e.g., deprivation.

In summary, OLS is inappropriate for our analysis of life expectancy, as statistically significant spatial autocorrelation of the regression residuals indicates misspecification of the OLS model, violation of OLS assumptions, or that additional important spatial factors are missing from the model. More fundamentally, variability in life expectancy likely represents the outcome of spatial processes. Incorporating geography into our modelling approach may produce better predictions and minimise error.

However, it is worth emphasising that OLS isn’t necessarily inappropriate for analysis of spatial data. For example, if the observations are independent, there is no clear spatial pattern to the model residuals, as evidenced using Moran’s \(I\) and \(I_i\), and no evidence of spatial non-stationarity… use OLS! However, the results of OLS should always be evaluated critically. In many cases, alternative techniques are more appropriate.

7.2.5 Spatial regression

An alternative to OLS and similar aspatial regression methods is spatial regression, which incorporates space explicitly in the regression framework.

There are many approaches to achieve this, and today we are going to focus on two main regression models: the spatial error model (SEM) and the spatial lag model (SLM).

Both approaches address the impact of spatial dependence, in which the configuration and spatial arrangement of the data is important and therefore needs to be considered.

How do we know that spatial dependence is important?

Answer

This is evidenced by the spatial autocorrelation present in the dependent and independent variables, and the OLS residuals.

While both model types address spatial dependence, SEM and SLM are typically selected based on different theoretical arguments and use different mathematical approaches.

The spatial lag model is defined as follows:

\[Y=\rho WY+X\beta+\epsilon\]

This is very similar to a standard OLS formulation, which could be written simply as \(Y=X\beta+\epsilon\)⁹⁰, but with an additional component \(\rho WY\), which is a spatially lagged \(Y\) variable i.e., we are including the average neighbourhood value for each observation as an additional explanatory variable.

The spatial lag model therefore includes the spatial lag of the dependent variable, with the calculation approach defined in Practical 3. In this case, the spatial lag is of life expectancy where the coefficient \(\rho\) is a single value which represents the strength and direction of the spatial lag influence, much like a coefficient \(\beta\). As \(\rho\) increases, neighbourhood effects become increasingly important.

This approach therefore accounts for spillover effects, in that the dependent variable at a specific location \(Y_i\) may be influenced by the corresponding dependent variable of its neighbours \(Y_j\), which we define based on the spatial weights matrix \(W\) i.e., life expectancy for an MSOA is influenced by life expectancy in neighbouring MSOAs, which in turn, influences life expectancy at other MSOAs, and so on.

The following example outlines a situation in which diffusion of behaviour could lead to spatial dependence:

“… the behaviour is likely to be highly social in nature, and understanding the interactions between interdependent units is critical to understanding the behaviour in question. For example, citizens may discuss politics across adjoining neighbours such that an increase in support for a candidate in one neighbourhood directly leads to an increase in support for the candidate in adjoining neighbourhoods”⁹¹

Darmofal (2015b)

In this situation, a spatial lag model would be appropriate, as it mathematically represents the underlying spatial process, in which the outcome (dependent) variable is influenced by the neighbourhood outcomes.

Unlike SLM, which includes spatial dependence as a explanatory variable, the spatial error model treats residual autocorrelation as a nuisance, rather than being the result of important spatial processes, as summarised below:

“spatial dependence observed in our data does not reflect a truly spatial process, but merely the geographical clustering of the sources of the behaviour of interest. For example, citizens in adjoining neighbourhoods may favour the same (political) candidate not because they talk to their neighbors, but because citizens with similar incomes tend to cluster geographically, and income also predicts vote choice. Such spatial dependence can be termed attributional dependence, as neighboring units have shared attributes that produce the clustering of behaviors”

Darmofal (2015b)

In turn, the spatial error model assumes that the observed spatial dependence arises due to the geographical clustering of attributes⁹², rather than being the result of spatial diffusion. This might also be due to other spatial variables that we haven’t included in the analysis. Based on this reasoning, SEM assumes that the error of an observed value can be related to the errors of its neighbours, and is defined as follows:

\[Y=X\beta+u\]

As above, this is similar to OLS, where the dependent variable \(Y\) is predicted by the independent variables \(X\) and its coefficients \(\beta\). This prediction is modified by a residual term \(u\), which is defined as:

\[u=\lambda Wu+\epsilon\]

Unlike a typical OLS, where the error term \(\epsilon\) must meet certain mathematical assumptions⁹³, here the error term \(u\) is modelled directly, where \(Wu\) is the spatial lag of the model errors and \(\lambda\) is a coefficient which represents the strength of the correlation between an observations error and that of its neighbours.

The mathematics of this are quite complicated, see Darmofal (2015a), so the important thing to understand is that unlike OLS, which selects coefficients \(\beta\) to produce the best fit between \(Y\) and \(X\) (minimising the sum of the squares of the residuals), a spatial error model finds the values for \(\beta\) and \(\lambda\) which best explain the observed data while accounting for the spatial correlation in the model errors. Unlike OLS, which treats the unexplained errors \(u\) as independent, SEMs allow an observations errors to be partly related to the errors of its neighbours, where the strength of this relationship is controlled by \(\lambda\). The values for \(\beta\) and \(\lambda\) are typically determined based on a maximum likelihood approach, which we’ll discuss more fully in Practical 10.

To summarise, both spatial lag and spatial errors models address spatial dependence. Spatial lag models assume that this dependence arises because of diffusion across spatial units, in which the outcome (dependent) variable is influenced by the outcome of its neighbours. This is addressed by including a spatially lagged outcome variable \(\rho WY\).

In contrast, spatial error models assume that spatial dependence arises because of attribute similarity, rather than true spatial processes. This is addressed by including a spatially lagged error variable \(\lambda Wu+\epsilon\), which represents the relationship between an observations error and the neighbourhood errors.

Thus, the difference between SLM and SEM is both theoretical and mathematical.

7.2.6 Spatial error model

To run a spatial error model in ArcGIS:

Open Spatial Autoregression⁹⁴, selecting Spatial Error Model (SEM) for “Model Type” and using life_expectancy as the “Dependent Variable”. To simplify comparison with OLS, use total_annual_income as the single independent variable, and for now, use Contiguity Edges Corners as the spatial weights method.

This should run successfully⁹⁵, and produce a range of output tables and a feature layer:

We’ll focus on a few key statistics to start, including those in the “Summary of SAR Results” and the “Model Diagnostics” tables.

What is the model \(\text{R}^2\)?

Answer

The coefficient of determination for this model is \(\text{R}^2\approx0.62\). Being more precise, the reported value is actually a Pseudo R2, which is described here. We won’t go into the mathematics of how this is calculated, but the important thing to note is that unlike an OLS \(\text{R}^2\), which can be interpreted quite simply: “The proportion of the total variation in \(Y\) that is explained by the independent variables \(X\)”, we need to be more cautious with the pseudo \(\text{R}^2\) values reported for spatial regressions.

While pseudo \(\text{R}^2\) is a measure of model fit, the actual interpretation is more complex, something like: “The variation in \(Y\) that can be explained by both the independent variables \(X\) and the spatial dependence of the errors”. In turn, we can’t say that approximately 62% of variability in life expectancy is explained by income, because this is ignoring the modelling of the spatial dependence of the errors via \(u=\lambda Wu+\epsilon\).

Is income positively or negatively correlated with life expectancy?

Answer

Life expectancy is still positively correlated with total annual income, with \(\beta\approx0.0003\).

The lag residual coefficient (\(\lambda\)) for my model is \(\approx0.48\). This can range from \(-1\) to \(+1\) and represents the strength of the correlation between an observations error and that of its neighbours.

Are errors for each observation positively or negatively correlated with that of its neighbours?

Answer

The \(\lambda\) value is positive (\(\approx0.48\)) indicating spatial clustering of errors i.e., high errors are found in neighbourhoods of high errors. A negative \(\lambda\) would indicate dispersion i.e., a high error may be found in a neighbourhood of low errors, or vice versa.

The tool also produces a Moran’s Scatter Plot of Residuals, described here, and which we introduced in Practical 3. This shows the residuals on the x-axis (the difference between the observations and the predicted values) and the spatial lag of the residuals on the y-axis. If there are evenly distributed across the four quadrants, there is no discernible spatial autocorrelation.

We can test for this directly using Moran’s \(I\).

Evaluate the spatial patterning of the SEM residuals using the global spatial autocorrelation tool.

How has the spatial autocorrelation of the residuals changed compared to the OLS model?

Answer

Previously, the distribution of OLS residuals was clustered, but with SEM we have now switched to a dispersed pattern which suggests we’re overcorrecting for spatial dependence (overfitting). This may be due to misspecification of the spatial weights matrix.

The use of contiguity-based approach may be less suitable for our analysis, because the processes which govern health (i.e., life_expectancy) operate over broader geographic scales, rather than only over contiguous boundaries e.g., air pollution over kilometre-scales, access to hospitals, commuting patterns, employment opportunities.

It is worth noting that the p-value for Moran’s \(I\) is \(>0.05\), as shown below, so while there is evidence of dispersion, it may not be statistically significant:

We can evaluate how the spatial patterning of the SEM residuals changes with different spatial weights:

Re-run a SEM using the Spatial Autoregression tool, but this time using K nearest neighbours as the neighbourhood type (\(K=10\)). When complete, assess for autocorrelation using Global Moran’s I and Local Moran’s I, making sure to use \(K=10\) for consistency.

How has the spatial patterning of the residuals changed?

Answer

Based on this spatial weights method, the distribution of residuals is now random and highly non-significant (\(p\approx0.79\)), as assessed using global Moran’s \(I\). Local Moran’s \(I_i\) reveals no statistically signficant hotspots or outliers.

The results above indicate that the spatial error model has successfully accounted for spatial dependence, as there is no statistically significant spatial autocorrelation of the regression residuals. In turn, we can be much more confident in the model coefficients and their use for prediction.

In this case, we can see that total annual income is positively correlated with life expectancy (\(\beta\approx0.0003\)), with a baseline life expectancy of approximately 63 years (\(\beta_0\approx63.28\)). For every one unit increase in income (£1), life expectancy would increase by \(\approx0.0003\) years, or more meaningfully, for every additional £10,000 increase in income, life expectancy would increase by \(\approx3\) years. An interesting but sobering statistic…

Have a go at including additional independent variables in the spatial error model.

Which variables are positively and negatively correlated with life expectancy, when spatial dependence is accounted for?

Interpretation

When the full suite of independent variables are included⁹⁶, most are statistically significant (i.e., \(p<0.05\)), excluding fuel_poverty_perc and employed_not_ftstudent_perc, indicating that these factors do not influence life expectancy, or are being accounted for by one or more of the other variables:

Disability and IMD decile are negatively correlated with life expectancy, whereas the percentage who are economically inactive is positively correlated, which is a little surprising. Both education attributes are negatively correlated, but the strength of this correlation is much higher for those with no qualifications.

The key thing to emphasise here is that by incorporating spatial dependence in the error term, we can be more confident in the calculated model coefficients, compared to those produced by OLS.

7.2.7 Spatial lag model

We can also evaluate the results of a spatial lag model:

Run Spatial Autoregression, but this time select the Spatial Lag Model, using the K-nearest neighbour approach (\(K=10\)) for consistency with the previous model. For simplicity, again use total_annual_income as the single independent variable.

Here are some of the key results:

As we discussed above, the spatial lag model assumes that spatial dependence arises because of diffusion across spatial units, in which the outcome (dependent) variable is influenced by the outcome of its neighbours. This is addressed by including a spatially lagged outcome variable \(\rho\), which is included above as Lag Y (rho).

As with the lag residual coefficient (\(\lambda\)) for the SEM model, this can range from \(-1\) to \(+1\) and represents the strength of the correlation between an observed dependent variable and that of its neighbours.

Here, the value is negative at \(\rho\approx-0.1\). This is indicative of weak negative spatial dependence i.e., when observed life expectancy is low, neighbourhood life expectancy is slightly higher, when we have controlled for the independent variables in the model⁹⁷. This weak spatial dependence is also evidenced by the probability value of \(p\approx0.117\), which we would classify as not statistically significant. Based on this data, we might already begin to be cautious about the use of the spatial lag model, given only weak evidence for spatial dependence of the observations.

We also need to be cautious when interpreting the coefficients of the spatial lag model. While these are included in the table above, they are not directly comparable to those produced by OLS or SEM, which can be interpreted quite simply: “How much does the dependent variable change for a one unit change in the independent variable, when all other independent variables are held constant?” This is not possible for SLM because of spatial spillover effects, and the fact that our dependent variable \(Y\) is influenced by the surrounding neighbourhood values.

For example, let’s imagine a spatial unit \(A\). With OLS, the dependent variable at that location \(Y_A\) is simply a function of the independent variable at that location \(X_A\) e.g., life expectancy at \(A\) is predicted by total annual income at \(A\):

For spatial regression, we also consider the impact of the neighbourhood, and let’s imagine \(A\) has three neighbours \(B,C\) and \(D\), as visualised below. In a spatial lag model, life expectancy at \(A\) is predicted by total annual income at \(A\) but also life expectancy at \(B,C\) and \(D\). As life expectancy at \(B\) is also influenced by total annual income at \(B\) (and the same for \(C\) and \(D\)), the independent variable total annual income is present in multiple locations! It affects life expectancy at \(A\) directly, via \(X_A\), and indirectly, via life expectancy at \(B,C,D\):

This allows us to distinguish two effects of an independent variable:

• direct i.e., the local effect for that spatial unit.
• indirect i.e., the average neighbourhood effect.

The value listed in the coefficient table above (\(\beta\approx0.0003\)) is the direct effect only, and we can find the indirect and the total effects (i.e., direct + indirect) in the coefficients summary at the bottom:

In this example, the indirect effects are very small at \(-0.000032\), even within the context of the small direct effect of \(0.000335\). Based on this interpretation, we can conclude that spatial spillover effects seem relatively unimportant for the current model setup, but this is not always the case, with a more striking example in the esri documentation, where the magnitude of the indirect effects exceeds the direct effects.

As a final consideration, it is worth evaluating the spatial patterning of the SLM residuals:

Assess for autocorrelation using Global Moran’s I, making sure to use \(K=10\) for consistency with the SLM model.

What is your interpretation of Moran’s \(I\)?

Interpretation

Statistically significant global spatial autocorrelation is still present in the model residuals, as evidenced in the Moran plot (below) and when tested using Global Moran’s \(I\) i.e., \(I>0\), \(p<0.05\).

7.2.8 Model selection

You should now have a good understanding of both the spatial error and the spatial lag models and how to evaluate them. This raises the question:

Which model should we select?

Model selection should always be based on both theory and statistical reasoning.

For the latter, we can be quite critical of the spatial lag model for this dataset and for this model specification, given the spatial autocorrelation of the model residuals. This is one of the key limitations of the OLS model which we were trying to address! Moreover, the spatial lag coefficient \(\rho\) was weak, negative and not statistically significant (\(p\approx0.117\)). By comparison, the spatial error model exhibited spatially random model residuals (\(I\approx-0.005\), \(p\approx0.79\)), and a highly significant \(\lambda\approx0.48\), indicating spatial clustering of errors.

From a theoretical perspective, we can also be cautious about whether the observed spatial dependence is a result of diffusion of behaviour, as assumed in the spatial lag model. In the political example presented by Darmofal (2015b) above, there is a clear spatial process by which political opinion and voting intention in one spatial unit might be influenced by political opinion in a neighbouring spatial unit i.e., interactions and discussion between people from different spatial units. This is less clear for our analysis of life expectancy, where it is difficult to envisage how my life expectancy is directly influenced by the life expectancy of my neighbours. Instead, spatial dependence in life expectancy is more easily explained by attribute similarity, and the co-patterning of income, education, employment, deprivation, and so on.

We also have other tools at our disposal to select the best model, one of which is the Lagrange Multiplier (LM), developed by Bera and Yoon (1993), and applied to spatial regression by Anselin et al. (1996a), with the process outlined well in the esri documentation.

We won’t go into the mathematics of the LM tests here, which you are welcome to explore in Anselin et al. (1996b). As a simple overview, the LM test begins with an OLS model and then evaluates whether incorporating the spatial lag of the outcome (SLM) or incorporating the spatial lag of the error (SEM) improves the model performance. These are known as the \(\text{LM}_{\text{lag}}\) and \(\text{LM}_{\text{error}}\) tests respectively. These circumstances are not mutually exclusive, so robust versions of these tests evaluate whether including the spatial lag of the outcome improves the model when the spatial lag of the error has already been accounted for (\(\text{Robust}\ \text{LM}_{\text{lag}}\)) and vice versa (\(\text{Robust}\ \text{LM}_{\text{error}}\)).

Based on the results of these four tests, and following the process described by Anselin and Rey (2014):

• If both \(LM_{lag}\) and \(LM_{error}\) are not significant (e.g., \(p>0.05\)), then we can simply use OLS i.e., incorporating spatial dependence does not improve the model.
• If only one is significant, that is the model of choice.
• If both are significant, then more robust test statistics are calculated as described above (\(\text{Robust}\ \text{LM}_{\text{lag}}\), \(\text{Robust}\ \text{LM}_{\text{error}}\)) and the process is repeated. If only one of the robust tests is significant, that is the model of choice.

Inspect the output from the previous Spatial Autoregression model, focusing on the LM Test Results section.

Based on the results above and the workflow of Anselin and Rey (2014), what is your preferred model?

Interpretation

Here, \(\text{LM}_{\text{error}}\) is significant at \(p<0.05\) whereas \(\text{LM}_{\text{lag}}\) is not statistically significant with \(p\approx0.69\). Based on the Lagrange Multiplier tests, we would reject SLM in favour of SEM, which accords with our theoretical and statistical interpretation above.

In future, you can allow the Spatial Autoregression tool to select the best model for you automatically based on the LM test, using “Model Type” = Auto-detect, although this should not take precedence over theoretical reasoning, and your own judgement of the output statistics.

As a final consideration, we could also run a combined autoregressive model (SAC) via the Spatial Autoregression tool. This is simply a combination of SEM and SLM, incorporating the spatial dependence of the errors as well as the spatial lag of the dependent variable.

Why might this not be a suitable choice in our situation?

Interpretation

Our theoretical interpretation above, and the results of the Lagrange Multiplier test, indicate that a spatial lag model is not a significant improvement compared to a standard OLS. As such, and with model parsimony in mind, there is no clear benefit to including SLM in a combined approach, when SEM is a marked improvement on its own.

7.2.9 Summary

In today’s practical, we have discussed the problems with aspatial regression methods when applied to spatial data. Mathematically these methods don’t account for the spatial processes which often influence the dependent variable of interest, and don’t account for spatial dependence. This is important because the resulting OLS models are often misspecified or invalid when key assumptions are not met (e.g., independence of observations, independence of error). More fundamentally, as GI Scientists we are working with spatial data. Aspatial methods ignore a key component of that data i.e., their geography. Incorporating this information into the regression framework can result in unbiased coefficients, can minimise errors, and can allow us to make more accurate predictions.

We introduced two common approaches to incorporate spatial dependence, namely the spatial error model (SEM) and the spatial lag model (SLM). These differ mathematically and theoretically, and should be selected based on the topic (e.g., is spatial dependence due to diffusion of behaviour or attribute similarity?) and the resulting model performance (e.g., statistical significance and the Lagrange Multiplier).

While these are the most commonly used spatial regression models, and are easily run in ArcGIS via the Spatial Autoregression tool, it is worth noting that there are alternatives. Examples include the SLX model, which includes spatially-lagged independent variables (unlike SLM, which includes a spatially-lagged dependent variable), and the Spatial Durbin and Spatial Durbin Error models, which combine the SLX model with SLM and SEM respectively.

…and with that, congratulations! You have completed the practical and should now be cautious about using OLS and similar aspatial methods for spatial data analysis, and should be able to implement and evaluate spatial regression methods in ArcGIS.

Finished!

7.3 Extra

In this practical we have developed global models to predict life expectancy. However, we should recognise that the model performance and the identified relationships (i.e., coefficients \(\beta\)) are specific to these data and the study area i.e., the North West of England. This raises the question:

How generalisable are our models to new locations and datasets?

To test this, we can generate new spatial regression models for a different regions of England and Wales e.g., London, and assess how model performance varies. The instructions above provide further information and context, but as a brief refresher:

Load the England and Wales MSOA geometries [MSOA_2021_EW_BGC_V3.shp] and lookup table [MSOA_2021_lookup_table.csv] and join via MSOA21CD. When complete, Select by Attributes using RGN22NM = 'London' (\(n=1002\)). Export the Selection to MSOA_2021_london, saving to practical-7.gdb to facilitate use of full field names.

Next, load the datasets containing the dependent and independent variables from [data/dhsc] and [data/ons/msoa], and join the important attributes to the London MSOA geometries (listed above). Note if you wish to include education, economic activity or disability status variables, these data would need to be obtained from the ONS. By comparison, the DHSC life expectancy, fuel poverty and IMD datasets contain values for all MSOAs in England and Wales, so no data wrangling required!

To simplify the regional comparison, you may wish to use a single independent variable e.g., the effects of deprivation (imd_decile) on life expectancy.

Run Spatial Autoregression using the same dependent and independent variables, model type (e.g., SEM) and spatial weights method (e.g., \(K=10\)), but for the North West and London regions separately.

For example, here are the coefficients for a spatial error model, where deprivation is the sole predictor of life expectancy, for the North West:

and for London:

How well do the models explain life expectancy in the two regions?

Answer

In the North West, \(\approx0.78 \%\) of variability in life expectancy can be explained by the model, as a function of deprivation and the modelling of the spatial dependence of the errors. By comparison, \(\approx0.58 \%\) of variability can explained in London. While these estimates are derived from the same modelling approaches, we should be cautious of comparing them directly, as these are pseudo \(\text{R}^2\) values, which are described here and discussed above.

How different is the baseline life expectancy between the two regions?

Answer

This is expressed by the intercept term \(B_0\), which is \(\approx82.9\) years in the North West and \(\approx84.4\) years in London. Note that this is the intercept term when deprivation is the only predictor, and this will vary as additional predictors are added.

How different is the new model? Does this change your interpretation of the links between MSOA-level attributes and life expectancy?

Answer

Both models show that life expectancy is negatively correlated with deprivation. The estimated coefficient is slightly larger in London at \(\approx-0.23\) compared to \(\approx-0.19\) in the North West, potentially revealing a larger effect of deprivation in London. Either way, deprivation shows a similar negative effect in both regions, which lends support to the generalisability of our model.

7.4 Resources

• Anselin, L. (1988). Spatial econometrics: methods and models. Kluwer Academic Publishers google schola, 2, pp.283-291.
• Anselin, L. (2002). Under the hood issues in the specification and interpretation of spatial regression models. Agricultural economics, 27(3), pp.247-267.
• Darmofal, D. (2015) Spatial Lag and Spatial Error Models. In: Spatial Analysis for the Social Sciences. Analytical Methods for Social Research. Cambridge University Press:96-118.
• Gibbons, S. and Overman, H.G. (2012). Mostly pointless spatial econometrics?. Journal of regional Science, 52(2), pp.172-191.