Section 11.3 Looking at the residuals and testing for spatial autocorrelation in regression
Residuals measure the distance between our observed \(Y\) values, \(Y\text{,}\) and the predicted \(Y\) values, \(\hat{Y}\text{.}\) So in essence they are deviations of observed reality from your model. Your regression line (or hyperplane in the case of multiple regression) is optimised to be the one that best represents your data if the assumptions of the model are met. Residuals are very helpful to diagnose, then, whether your model is a good representation of reality or not. Most diagnostics of the assumptions for regression rely on exploring the residuals.
Now that we have fitted the model, we can extract the residuals. If you look at the "fit_1" object in your RStudio environment or if you run the
str() function to look inside this object, you will see that this object is a list with different elements, one of which is the residuals. An element of this object then includes the residual for each of your observations (the difference between the observed value and the value predicted by your model). We can extract the residuals using the residuals() function and add them to our spatial data set. Let’s fit a model with all the key predictors available in the dataset and inspect the residuals.
# write model into object to save typing later
eq1_90 <- HR90 ~ RD90 + PS90 + MA90 + DV90 + UE90 + SOUTH_f
fit3_90 <- lm(eq1_90, data=ncovr)
summary(fit3_90)
ncovr$res_fit3 <- residuals(fit3_90)
## ## Call: ## lm(formula = eq1_90, data = ncovr) ## ## Residuals: ## Min 1Q Median 3Q Max ## -17.64 -2.61 -0.70 1.65 68.51 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 6.5165 1.0240 6.36 2.3e-10 *** ## RD90 3.8723 0.1427 27.13 < 2e-16 *** ## PS90 1.3529 0.1003 13.49 < 2e-16 *** ## MA90 -0.1011 0.0274 -3.69 0.00023 *** ## DV90 0.5829 0.0545 10.69 < 2e-16 *** ## UE90 -0.3059 0.0409 -7.47 1.0e-13 *** ## SOUTH_f1 2.1941 0.2205 9.95 < 2e-16 *** ## --- ## Signif. codes: ## 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 4.99 on 3078 degrees of freedom ## Multiple R-squared: 0.436, Adjusted R-squared: 0.435 ## F-statistic: 397 on 6 and 3078 DF, p-value: <2e-16
If you now look at the dataset, you will see that there is a new variable with the residuals. In those cases where the residual is negative, this is telling us that the observed value is lower than the predicted (that is, our model is overpredicting the level of homicide for that observation); when the residual is positive, the observed value is higher than the predicted (that is, our model is underpredicting the level of homicide for that observation).
We could also extract the predicted values if we wanted. We would use the
fitted() function.
ncovr$fitted_fit3 <- fitted(fit3_90)
Now look, for example, at the second county in the dataset: the county of Ferry in Washington. It had a homicide rate in 1990 of 15.89. This is the observed value. If we look at the new column we have created ("fitted_fit1"), our model predicts a homicide rate of 2.42. That is, knowing the level of resource deprivation and all the other explanatory variables (and assuming they are related in an additive manner to the outcome), we are predicting a homicide rate of 2.42. Now, this is lower than the observed value; so our model is underpredicting the level of homicide in this case. If you observed the residual, you will see that it has a value of 11.56, which is simply the difference between the observed (15.89) and the predicted value (2.42). This is normal. The world is too complex to be encapsulated by our models, and their predictions will depart from the observed world.
With spatial data one useful thing to do is to look at any spatial patterning in the distribution of the residuals. Notice that the residuals are the difference between the observed values for homicide and the predicted values for homicide, so you want your residual to not display any spatial patterning. If, on the other hand, your model displays a patterning in the areas of the study region where it predicts badly, then you may have a problem. This is telling you that your model is not a good representation of the social phenomena you are studying across the full study area: there is systematically more distortion in some areas than in others.
We are going to produce a choropleth map for the residuals, but we will use a common classification method we haven’t covered yet: standard deviations. Standard deviation is a statistical technique type of map based on how much the data differs from the mean. You measure the mean and standard deviation for your data. Then, each standard deviation becomes a class in your choropleth maps. In order to do that, we will compute the mean and the standard deviation for the variable we want to plot and break the variable according to these values. The following code creates a new variable in which we will express the residuals in terms of standard deviations away from the mean. So, for each observation, we subtract the mean and divide by the standard deviation. Remember, this is exactly what the
scale function does:
# because scale is made for matrices, we just need to get the first column using [,1]
ncovr$sd_breaks <- scale(ncovr$res_fit3)[,1]
# this is equal to (ncovr$res_fit1 - mean(ncovr$res_fit1)) / sd(ncovr$res_fit1)
summary(ncovr$sd_breaks)
## Min. 1st Qu. Median Mean 3rd Qu. Max. ## -3.537 -0.524 -0.140 0.000 0.331 13.741
Next we use a new style,
fixed, within the tm_fill function. When we break the variable into classes using the fixed argument we need to specify the boundaries of the classes. We do this using the breaks argument. In this case we are going to ask R to create 7 classes based on standard deviations away from the mean. Remember that a value of 1 would be 1 standard deviation higher than the mean, and -1 respectively lower. If we assume normal distribution, then 68% of all counties should lie within the middle band from -1 to +1.
my_breaks <- c(-14,-3,-2,-1,1,2,3,14)
tm_shape(ncovr) +
tm_fill("sd_breaks", title = "Residuals", style = "fixed",
breaks = my_breaks, palette = "-RdBu") +
tm_borders(alpha = 0.1) +
tm_layout(main.title = "Residuals", main.title.size = 0.7 ,
legend.position = c("right", "bottom"), legend.title.size = 0.8)
Notice the spatial patterning of areas of over-prediction (negative residuals, or blue tones) and under-prediction (positive residuals, or reddish tones). This visual inspection of the residuals is suggesting that spatial autocorrelation in the residuals may be present here. This, however, would require a more formal test.
Remember from Chapter 9 that in order to do this first we need to turn our
sf object into a sp class object and then create the spatial weight matrix. We cannot examine spatial dependence without making assumptions about the structure of the spatial relations among the units in our study. For this illustration we are going to follow [187] and use 10 nearest neighbours to define the matrix to codify these spatial relations. If the code below and what it does is not clear to you, revise the notes from Chapter 9, when we first introduced it.
# We coerce the sf object into a new sp object
ncovr_sp <- as(ncovr, "Spatial")
# Then we create a list of neighbours using 10 nearest neighbours criteria
xy <- coordinates(ncovr_sp)
nb_k10 <- knn2nb(knearneigh(xy, k = 10))
rwm <- nb2mat(nb_k10, style='W')
rwm <- mat2listw(rwm, style='W')
We obtain Moran’s test for regression residuals using the function
lm.morantest() as below. It is important to realize that Moran’s I test statistic for residual spatial autocorrelation takes into account the fact that the variable under consideration is a residual, computed from a regression. The usual Moran’s I test statistic does not. It is therefore incorrect to simply apply a Moran’s I test to the residuals from the regression without correcting for the fact that these are residuals.
lm.morantest(fit3_90, rwm, alternative="two.sided")
## ## Global Moran I for regression residuals ## ## data: ## model: lm(formula = eq1_90, data = ncovr) ## weights: rwm ## ## Moran I statistic standard deviate = 12, p-value ## <2e-16 ## alternative hypothesis: two.sided ## sample estimates: ## Observed Moran I Expectation Variance ## 9.544e-02 -1.382e-03 6.029e-05
You will notice that we obtain a statistically significant value for Moran’s I. If we diagnose that spatial autocorrelation is an issue, that is, that the errors (residuals) are related systematically among themselves, then we likely have a problem, and we may need to use a more appropriate approach: a spatial regression model. Before going down that road, it is important to check this autocorrelation is not the result of an incorrectly specified model. It may appear when estimating a non-linear relationship with a linear model. In this scenario, [196] suggests checking the remainders from the OLS estimated model on variable logarithms. It is also worth considering it may be the result of other forms of model mis-specification, such as omission of covariates or wrongful inclusion of predictors.
