Mapping And Modelling Geographic Data In R

Session 5: Spatial and Geographically Weighted Regression Analysis

Introduction

In this session we draw together some of what we have learned about mapping and creating spatial weights in R, using it to undertake spatial forms of regression, specifically a spatially lagged error model, a lagged y model, a mixed model and Geographically Weighted Regression. The session is more an introduction to how to fit these sorts of models in R rather than an in-depth tutorial into their specification and use. Some of the models fit the (simulated) data better than others. However, that is not to say they would necessarily be preferred when using other data and in other analytical contexts. Which model to choose will depend on the purpose of the analysis and the theoretical basis for the modelling.

Please note:

This session assumes an intermediate knowledge of linear regression modelling with some understanding of how to fit and interpret these models in R. If you do not have this knowledge then as you are reading keep in mind the foundational aim of the session which is to use a spatial weights matrix to help quantify geographical patterns in data and then, if those patterns exist, to help model them in a geographical way, doing so by incorporating that matrix into the regression models.

You are strongly advised to have read https://www.dropbox.com/s/tlz7nenwkdny4ds/Harris_chapter-revised.pdf?dl=0 and https://www.dropbox.com/s/5uh9grmtilrup7b/Wiley-AAG5.pdf?dl=0 before undertaking this session.

To get started,

rm(list=ls())
# Remove eveything currently in the workspace
if(!"sp" %in% installed.packages()) install.packages("sp")
require(sp)
load("session5.RData")

This loads a workspace containing an sp object of Beijing’s districts (spatial polygons), synthetic land parcel and district data (spatial points), and a spatial weights object for the land parcels with a Gaussian decay to the 35th neighbour. The choice of 35 was determined using a correlation approach – how correlated is the value at each location with its first, second, third, …, etc. neighbour? See Session 4 for further details.

Using a Moran’s test, it is shown that the (log) of the land price values show significant spatial variation,

if(!"spdep" %in% installed.packages()) install.packages("spdep")
require(spdep)
moran.test(combined.map$LNPRICE, spknear35gaus)

… our task is to try and explain that variation through a model-based approach.

OLS regression

We begin by fitting a regression model. It is possible that the regression model will explain the patterning of the land prices: land prices could be higher the closer the land is to the Central Business District, for example, or closer to rivers. And, to some extent it does. A test of the residuals, however, reveals an apparent spatial dependencies in them that violates the assumption of independence. It isn’t especially strong but it is present.

model1 <- lm(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map)
lm.morantest(model1, listw2U(spknear35gaus))

In addition to their apparent lack of independence we may also note that the residuals appear to show evidence of heteroskedasticity – of non-constant variance (they are therefore not independent nor identically distributed). Heteroskedasticity is sometimes indicative of spatial variation – that a ‘one size fits all’ approach to modelling will not work. But it can also be because the model is incorrectly specified (missing variables perhaps, non-linear relationships, unmodelled interactions amongst the variables, etc.).

A plot to check for heteroskedasticty is,

plot(model1, which = 1)

Unfortunately, that same code won’t work for the spatial models we produce below. Instead, we can write a function to produce what we need.

hetero.plot <- function(model) {
  fm <- as.character(formula(model))
  plot(residuals(model) ~ fitted(model),
       xlab = paste("Fitted values\n",fm[2],fm[1],fm[3]),
       ylab = "Residuals")
  abline(h=0, lty="dotted") 
  lines(lowess(fitted(model), residuals(model)), col="red")     
}
# And, now to run it...
hetero.plot(model1)

There are no quick fixes for the violated assumption of independent and identically distributed errors (heteroskedasticity). The violation suggests we cannot take on trust the standard errors, t- and p-values shown under the model summary. It is likely that the standard errors for at least some of the predictor variables have been under-estimated (because if we have spatial dependencies in the residuals then they likely arise from spatial dependencies in the data which in turn mean we have less degrees of freedom than we think we have).

summary(model1)

If we have doubts about this model because of the spatial patterning of the errors (a pattern which is likely but not necessarily caused by the patterning of the Y variable) then we need to consider other approaches.

Spatial Econometric Models

Spatial Error Model

One option is to fit a spatial simultaneous autoregressive error model which decomposes the error (the residuals) into two parts: a spatially lagged component and a remaining error: y=Xβ+λWξ+ε

Fitting the model and comparing it with the standard regression model we find that two of the predictor variables (DELE and POPDEN) no longer are significant at a conventional level and that the standard errors for many have risen. The lambda (a measure of spatial autocorrelation) is significant. The model fits the data better than the previous model: the AIC score is lower and the log likelihood value greater, as is the pseudo-R2 value:

if(! "spatialreg" %in% installed.packages()) install.packages("spatialreg")
require(spatialreg)
model2 <- errorsarlm(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map, listw = spknear35gaus)
summary(model2)

The AIC is given in the model output above but we can also obtain it using,

AIC(model1)
AIC(model2)

As a rule of thumb, a reduction in the AIC score of ten or more is strong evidence in favour of the model with the lower AIC value – which here is the spatial error model. (You can this rule of thumb at https://stats.stackexchange.com/questions/282686/what-percentage-decrease-in-aic-is-significant)

We can also use the log likelihood values to compare the models:

logLik(model1)
logLik(model2)

Here, higher is better so the evidence is again in favour of the spatial error model. More formally,

LR <- 2 * (logLik(model2) - logLik(model1)) 
      # This is the likelihood ratio
dfchg <- attr(logLik(model2), "df") - attr(logLik(model1), "df")
      # This is the difference in the degrees of freedom
LR > qchisq(0.99, df = dfchg)
      # This compares the likelihood ratio with a chi-square
      # distribution

… which passes a statistical test at a 99 per cent confidence in favour of the spatial error model. (If you want to read up on this test: https://en.wikipedia.org/wiki/Likelihood-ratio_test)

Finally, the pseudo-R-squared for the spatial error model is

cor(combined.map$LNPRICE, fitted(model2))^2

which is higher than the R-squared value for the original model; that was

summary(model1)$r.squared

(The R-squared value can be calculate as the square of the correlation between the Y values and what the model predicts for them, the fitted values)

In summary, the spatial error model appears to fit the data better although the problem of heteroskedasticity remains:

hetero.plot(model2)

A spatially lagged y model

Although the spatial error model (above) fits the data better than the standard OLS model, it tells us only that there is an unexplained spatial structure to the residuals, not what caused them. It may offer better estimates of the model parameters and their statistical significance but it does not presuppose any particular spatial process generating the patterns in the land price values. A different model that explicitly tests for whether the land value at a point is functionally dependent on the values of neighbouring points is the spatially lagged y model: y=ρWy+Xβ+ε. It models an ‘overspill’ or chain effect where the outcome (the Y value) at any location is affected by the outcomes at surrounding locations.

The model is fitted in R using,

model3 <- lagsarlm(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map, spknear35gaus)
summary(model3)

The model is an improvement on the OLS model. However, there is significant autocorrelation remaining in the residuals (see the bottom lines of the model output) and it does not appear to fit the data as well as the error model: the lagged y model has a greater AIC, lower log likelihood and lower pseudo-R2, each indicating a poorer fit.

AIC(model3)
logLik(model3)
cor(combined.map$LNPRICE, fitted(model3))^2

The heteroskedasticity remains:

hetero.plot(model3)

(potentially related to the significant autocorrelation in the residuals).

Note that the beta estimates of the lagged y-model cannot be interpreted in the same way as for a standard OLS model. For example, the beta estimate of 0.004 for the POPDEN variable does not mean that if (hypothetically) we increased that variable by one unit at each location we should then expect the (log) of the land parcel to everywhere increase by 0.004 holding the other X variables constant. That is the correct interpretation for a standard (OLS) regression model and also for the spatial error model but not for where the lag of the Y variable is included as a predictor variable. The reason is because if we did raise the POPDEN value it would start something akin to a ‘chain reaction’ through the feedback of Y via the lagged Y values. The total impact is a sum of the direct effect – what is predicted to happen through the 1 unit change in POPDEN – and the indirect effect, which is that caused by ‘the chain reaction’ / feedback / ‘overspill’ in the system:

impacts(model3, listw = spknear35gaus)

In fact, an increase of 1 unit in, for example, the POPDEN variable, will play out slightly differently in different places (because they have different neighbours with different Y values and at different distances from each other). The code below, which is from Ward & Gleditsch (2008, pp.47), will calculate the impact (of a one unit change in POPDEN) at each location but below I have limited it to the first ten else it will take a long time to run.

n <- nrow(combined.map)
I <- matrix(0, nrow=n, ncol=n)
diag(I) <- 1
rho <- model3$rho
beta <- model3$coefficients["POPDEN"]
weights.matrix <- listw2mat(spknear35gaus)
results <- rep(NA, times=10)
for (i in 1:10) {
  xvector <- rep(0, times=n)
  xvector[i] <- 1
  impact <- solve(I - rho * weights.matrix) %*% xvector * beta
  results[i] <- impact[i]
}
results

Choosing between the models using Lagrange Multiplier (LM) Tests

Before fitting the spatial error and lagged y models (above), we could have looked for evidence in support of them using the function lm.LMtests(…). This tests the basic OLS specification against the more general spatial error and lagged y models. Robust tests also are given. There is evidence for both of the spatial models in favour of the simpler OLS model but it is stronger (in purely statistical terms) for the error model.

lm.LMtests(model1, spknear35gaus, test="all")

Anselin and Rey (2014, p.110) offer the following decision tree that can, in the absence of a more theoretical basis for the model choice (e.g. the type of process being modelled), be used in conjunction with the test results above to help select the model.

browseURL("anselin_rey.jpg")

Using the decision tree, there seems to be greater evidence in favour of the spatial error model than the spatially lagged y model, although potentially both are ‘beaten’ by a third model, the spatial autoregressive moving average (SARMA) model! We won’t dwell on this other than to note that it can account for spatial dependence among the error terms as well as spatial dependence in the Y variable. To fit it in R:

model4 <- spautolm(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map, spknear35gaus, family = "SMA")
summary(model4)

On balance it may fit the data less well than model 2 given its greater complexity (remember, a lower value is preferred for AIC scores, and a higher value for log likelihood ones):

fit <- matrix(nrow = 2, ncol = 2)
fit[1,] <- c(AIC(model2), AIC(model4))
fit[2,] <- c(logLik(model2), logLik(model4))
colnames(fit) <- c("Model2", "Model4")
rownames(fit) <- c("AIC", "logLik")
fit

Including lagged X variables

So far we have accommodated the spatial dependencies in the data by consideration to the error and with consideration to the dependent (y) variable. Attention now turns to the predictor variables. An extension to the lagged y model is to lag all the included x variables too (this is sometimes referred to as the Spatial Durbin model)

model5 <- lagsarlm(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map, spknear35gaus, Durbin = TRUE)
summary(model5)

What we find is that the lag of the distance to the CBD (measured on a log scale) is significant but note that the direction of the relationship is different from that for the original DCBD variable – the sign has reversed. The same has happened to the dummy variable indicating the sale of the land parcel is in the years 2008 to 2009. Taking the first case, what this may suggest is that the relationship between land price value and the (log of) distance to the CBD is not linear. Adding the square of this variable to the original OLS model improves the model fit: the difference of about 5 in the AIC scores is medium evidence in favour of the new model.

model1b <- update(model1, . ~ . + I(DCBD^2))
    # Adds the square of the DCBD variable
AIC(model1b)
AIC(model1)
LR <- 2 * (logLik(model1b) - logLik(model1)) 
dfchg <- attr(logLik(model1b), "df") - attr(logLik(model1), "df")
LR > qchisq(0.99, df = dfchg)

We could also try adding the square of the DCBD variable to the spatial error model:

model2b <- update(model2, . ~ . + I(DCBD^2))
AIC(model2b)
AIC(model2)
LR <- 2 * (logLik(model2b) - logLik(model2)) 
dfchg <- attr(logLik(model2b), "df") - attr(logLik(model2), "df")
LR > qchisq(0.99, df = dfchg)

… it isn’t a success. Perhaps adding the spatial lag of the DCBD variable will work better?

lag.DCBD <- lag.listw(spknear35gaus, combined.map$DCBD)
model2c <- update(model2, . ~ . + lag.DCBD)
AIC(model2c)
AIC(model2)
LR <- 2 * (logLik(model2b) - logLik(model2)) 
dfchg <- attr(logLik(model2c), "df") - attr(logLik(model2), "df")
LR > qchisq(0.99, df = dfchg)

No, it doesn’t!

Including lagged X variables in the OLS model

We can, if we wish, include specific lagged X variables in the OLS model. The process is to create them then include them in the model. The lag of DCBD and the lag of Y0809 are the most obvious candidates to include (from model 4, above). To create the lagged variables and add them to the model:

lag.DCBD <- lag.listw(spknear35gaus, combined.map$DCBD)
lag.Y0809 <- lag.listw(spknear35gaus, combined.map$Y0809)
model1c <- update(model1, .~. + lag.DCBD + lag.Y0809)
AIC(model1c)
AIC(model1)

The result seems to fit the data better than the original OLS model (AIC score of 2848.3 Vs 2865.7 – remember, the lower the better) but actually the lag of DCBD appears not to be significant in this model whilst significant spatial autocorrelation appears to remain:

summary(model1c)
lm.morantest(model1c, spknear35gaus)

Taking Stock

At this point, it is worth looking back and taking stock. Which of the models appears to fit the data best?

If your head is spinning a bit don’t worry because that serves to make a useful point. There are lots of models that could be fitted to the data and we have by no means exhausted the possibilities above. Therefore we need a model fitting strategy. Sometimes that will be informed by theory or by other studies. Otherwise, I would consider the following:

1. Map the Y variable and look for a spatial pattern
2. Formally test for spatial dependency in the Y variable by constructing a spatial weights matrix and using global and/or localised Moran’s values
3. for (2) you will need to give consideration to how you specify the weights matrix – you may want to try different weighting schemes. What do the tell you about the scale of the spatial pattern in the data? (Is it localised clustering - locations only related to their most immediate neighbours or…?)
4. Fit the OLS model and go through the usual process of refining that model (e.g. drop insignificant variables)
5. Test the residuals for a spatial pattern using the Moran’s test for regression residuals.
6. If there is a pattern, either choose a spatial regression model on theoretical grounds (because it captures the process you want to model) or use the Lagrange multiplier tests.
7. Consider also the possibility that the spatial pattern exists because the OLS model was mis-specified (e.g. the relationship between one or more Xs and the Y variable is non-linear).
8. Compare your spatial regression model with the results of a Geographically Weighted Regression (see below).

Geographically Weighted Regression

An alternative way of attempting to explain the spatial variation in the land value prices is to allow the effect sizes of the predictor variables to themselves vary over space. For the previous spatial regression models they were fixed – there was one estimate of, for example, the effect of distance to the CBD on land prices. But what if being close to the CBD matters more in some places than it does in others?

Geographically Weighted Regression (GWR) is a method where the estimate of βx at point location i is not simply the global estimate for all points in the study region but a local estimate based on surrounding points weighted by the inverse of their distance away (hence geographically weighted). If there are 1000 locations in a study region then GWR will fit 1000 localised models. The results of those models are then compared to see how much the regression coefficients vary across the study region.

Fitting the GWR model

The stages of fitting a Geographically Weighted Regression model are first to load the GWR library, calculate a distance matrix containing the distances between points (this is needed for the inverse distance weighting), calibrate the bandwidth for the local model fitting (or you can specify a bandwidth yourself), fit the model then look for evidence of spatial variations in the estimates.

For the first couple of stages:

if(! "GWmodel" %in% installed.packages()) install.packages("GWmodel")
require(GWmodel)
# Install and load the package
distances <- gw.dist(dp.locat=coordinates(combined.map))
# Calculate the distances between the points

Here the bandwidth decreases to zero at the 30th neighbour, encouragingly similar to the 35th neighbour value we have used for our spatial weights throughout this session. Rather less encouragingly, it seems surprisingly sensitive to the shape of the kernel (the way the weighting tapers with distance). For comparison:

bw2 <- bw.gwr(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map, adaptive=T, dMat=distances, kernel = "bisquare")
bw2

This sensitivity to the choice of kernel function is worth investigating further. However, we will skip over it (!) and use the Gaussian version to next fit the model:

gwr.model <- gwr.basic(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map, adaptive=T, dMat=distances, bw=bw, kernel = "gaussian")
gwr.model

Looking at the results we can see that the model has a better fit to the data than any other fitted thus far although this is perhaps not surprising (as a GWR is really lots of models fitted one after another across the study region it is always likely to out-perform a one size fits all approach). The summary of the GWR estimates is of how the (local) beta estimates vary across the study region - of, for example, how the effect of distance to the CBD varies from one place to another. From the inter-quartile range (shown in the output), the effect of distance to CBD on land value prices is typically found to be from -0.682 to -0.030 (in some locations, as the distance from the CBD increases, the log of the land price values drops by 0.682 but in others the drop is much less, dropping by 0.030).

The function montecarlo.gwr(…) uses a randomisation approach to undertake significance testing for the variability of the estimated local beta estimates (the regression parameters) - i.e. of how significantly they exhibit a geographical variation. The default number of simulations is 99 (which, with the actual estimates, gives 100 sets of values in total). That is quite a low number but can be used for illustrative purposes. In practice it would be better to raise it to 999, to 9999 or even more. Even with 99 it will take a little while to run.

montecarlo.gwr(LNPRICE ~ DCBD + DELE + DRIVER + DPARK + POPDEN + JOBDEN + Y0405 + Y0607 + Y0809, data=combined.map, adaptive=T, bw=bw, kernel = "gaussian")

All but one (or marginally two) of the variables appear to show significant spatial variation (p < 0.01).

The localised estimates (the result from each fitting of the regression model as it moves across the study region) are stored in the Spatial Data Frame (SDF):

names(gwr.model)
head(gwr.model$SDF)

The estimates for the effect of distance to CBD on land price, for example, are founded in the column, gwr.model\(SDF\)DCBD

head(gwr.model$SDF$DCBD)

Since each one of these estimates corresponds to a location in the original map (and in the same order), it is straightforward to map them using the native plot functions for sf objects (see https://r-spatial.github.io/sf/articles/sf5.html)

if(! "rgeos" %in% installed.packages()) install.packages("rgeos")
# This package is needed to convert the sp objects into sf
if(! "tmap" %in% installed.packages()) install.packages("tmap")
combined.map$DCBD_beta <- gwr.model$SDF$DCBD
# Create a new variable in the map data which are local GWR
# estimates for the DCBD variable
sf_map <- st_as_sf(combined.map)
sf_districts <- st_as_sf(districts)
# Converstion to simple features
plot(sf_map["DCBD_beta"], reset = FALSE, pch = 19, cex = 0.8)
# pch controls the plotting 'character' - see ?points
# cex is the character expansion (i.e. size)
# You need reset = FALSE to add further elements to the map
plot(st_geometry(sf_districts), add = TRUE)

However, not all of those regression coefficients are necessarily significant. To highlight those that that are (p < 0.05):

if(! "tidyverse" %in% installed.packages()) install.packages("tidyverse")
require(tidyverse)
sf_map$DCBD_beta_tval <- gwr.model$SDF$DCBD_TV
sf_map %>%
  filter(abs(sf_map$DCBD_beta_tval) > 1.96) -> sub.map
plot(sub.map["DCBD_beta"], reset = FALSE, pch = 19, cex = 0.8)
plot(st_geometry(sf_districts), add = TRUE)

Arguably a more stringent definition of significance should be applied as GWR undertakes a process of repeat model fitting. To adjust the p-values for multiple hypothesis tests see the help menu, ?gwr.t.adjust(gwm.Obj).

Geographically Weighted Statistics

In closing, it should be noted that GWR is an extension of a broader class of statistics – geographically weighted statistics.

For example, to calculate geographically weighted summary statistics for the land price variable:

gwstats <- gwss(combined.map, vars = "LNPRICE", adaptive=T, bw=bw,
     kernel = "gaussian")
head(gwstats$SDF)
names(gwstats$SDF)

What has been calculated is the geographically weighted mean land price around each location and the geographically weighted standard deviation amongst other values. The geographically weighted median, and interquartile range could also be calculated (see ?gwss).

The geographically weighted standard deviation, for example, shows how much the land prices vary locally around each location. We can, of course, map them:

sf_map$GWSD <- gwstats$SDF$LNPRICE_LSD
plot(sf_map["GWSD"], reset = FALSE, pch = 19, cex = 0.8)
plot(st_geometry(sf_districts), add = TRUE)

The areas shown in yellow have the most local variation in the land prices.

You can - and should! - read more about geographically weighted regression in the further reading (see below).

Conclusion

We have covered a lot of ground in this session, drawing on a range of statistical methods, many or all of which may be new to you. The important take home messages are:

1. How they incorporate spatial relationships, encoded by means of a spatial (or distance) weights matrix into their calculations
2. The process of detecting a spatial pattern in a variable of interest, fitting a non-spatial model to explain it and then trying some spatial versions instead.

It is (b) that will form the basis for your forthcoming assignment.

Futher Reading

Gollini I, Lu B, Charlton M, Brunsdon C & Harris P (2015) GWmodel: An R Package for Exploring Spatial Heterogeneity Using Geographically Weighted Models. Journal of Statistical Software, 63 (just focus on the sections about geographically weighted summary statistics and don’t worry too much about the code used to map the results as it is a bit dated)

References

Anselin L & Rey SJ (2014) Modern Spatial Econometrics in Practice: A Guide to GeoDa, GeoDaSpace and PySAL. Chicago: GeoDa Press.

Ward M & Gleditsche KS (2008) Spatial Regression Models. London: Sage.