R Notebook

Now that we have started working with some spatial statistics, this next unit will help us start building towards spatial regression. Before we can get there, we need to learn one foundational topic in spatial statistics: spatial autocorrelation. One assumption of any regression model is independence. Let’s say we are working with state-level data. If each state represents an independent observation, then no one state influences any other. But is this really ever true?

In Geography, Tobler’s Law states: Everything is related to everything else, but near things are more related than distant things. We might expect, therefore, that neighboring states influence each other. This tends to be the case - if I pick any two states in the South, they are likely to be more related to each other than, say, a Northern state. So my independence assumption is violated - we often call this “Spatial Autocorrelation”. Data with high spatial autocorrelation is data where neighboring geographic units are more related than distant geographic units. In today’s unit, we will practice methods to check for spatial autocorrelation. Later on, we will work on methods to control for this problem in spatial regression.

To illustrate this, let’s load in some data. Suppose I’m a researcher studying poverty in NYC. I have income data at the neighborhood level, and I have created a Poor-index based on the ratio of a household’s income to the poverty level. I then aggregated the number of households with an index value of 1 or below to create my poor count variable. Are impoverished households clustered near each other in NYC? This will help us to learn whether we should concentrate resources in certain areas more than others. Is there spatial inequality alongside economic inequality?

#Reasearch Question - Whats the impact of Crimerate and poverty rate on Unemployment rate.

library(sf)     
library(dplyr)   
library(spData) 
library(ggplot2)
library(ggthemes)
library(spdep)
library(sf)     
library(dplyr)   
library(spData) 
library(ggplot2)
library(ggthemes)
library(spdep)
library(spatialreg)

Warning: package ‘spatialreg’ was built under R version 4.3.3Loading required package: Matrix

Attaching package: ‘spatialreg’

The following objects are masked from ‘package:spdep’:

    get.ClusterOption, get.coresOption,
    get.mcOption, get.VerboseOption,
    get.ZeroPolicyOption, set.ClusterOption,
    set.coresOption, set.mcOption,
    set.VerboseOption, set.ZeroPolicyOption

library(GWmodel)

Warning: package ‘GWmodel’ was built under R version 4.3.3Loading required package: robustbase
Warning: package ‘robustbase’ was built under R version 4.3.3Loading required package: sp
Loading required package: Rcpp
Welcome to GWmodel version 2.3-2.

library(tidyr)

Attaching package: ‘tidyr’

The following objects are masked from ‘package:Matrix’:

    expand, pack, unpack

setwd("C:/Users/Donny/OneDrive/R studio/chicago_airbnb")

Warning: The working directory was changed to C:/Users/Donny/OneDrive/R studio/chicago_airbnb inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

chicago <- st_read("airbnb_chicago 2015.shp")

Reading layer `airbnb_Chicago 2015' from data source 
  `C:\Users\Donny\OneDrive\R studio\chicago_airbnb\airbnb_Chicago 2015.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 77 features and 20 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304
Geodetic CRS:  WGS 84

ggplot()+
  geom_sf(data = chicago, aes(fill=unemployed))+
  scale_fill_steps(
    name = "Unemployment Rate",
    low = "lightsteelblue1",
    high = "steelblue1",
    n.breaks = 4,
    show.limits = T)+
  theme_void()

chicago_list <- chicago %>% 
  poly2nb(st_geometry(chicago)) %>% 
  nb2listw(zero.policy = TRUE)

moran.plot(chicago$unemployed, 
           chicago_list, 
           zero.policy = TRUE, 
           labels = F,
           xlab = 'Unemployment rate across Chicago',
           ylab = 'Lagged Unemployment rate of neighbors',
                   pch=20)

NA
NA

Since our slope is positive, we have positive autocorrelation, so theres likely to be clustering of neighborhoods with high unemployment rates being near eachother

#Now we can calculate the "Global Moran's I" for the data
#start with the weights, then give R the corresponding variable
chicago_list %>% 
  moran.test(chicago$unemployed, ., zero.policy = TRUE)

    Moran I test under randomisation

data:  chicago$unemployed  
weights: .    

Moran I statistic standard deviate = 8.035, p-value
= 4.678e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.564814885      -0.013157895       0.005174172 

Since we have a Moran’s I statistic of 0.564, and a p-value of< 0.001, our data wasn’t produced by a spatially random process. Since our Moran’s I is somewhat close to 1, this shows a slight correlaton between areas that are near one another and rejects the idea of spatial randomness or independence.

Next, we’ll take a look at Local Indicators of Spatial Autocorrelation (LISA) plots. This will help us to identify clusters of low or high values in the graph, as well as outliers. This uses the Local Moran’s I statistic. Instead of calculating one Moran’s I value for the whole plot (as is the case with the global statistic), the local version calculates a local I for each spatial unit (neighborhood)

model <- chicago %>% 
  mutate(crimerate = num_crimes/population)%>% 
  select(c(unemployed, crimerate,population, geometry)) %>% 
  filter(!population %in% "NA") %>% 
  mutate(population = as.numeric(population))
model

Simple feature collection with 77 features and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304
Geodetic CRS:  WGS 84
First 10 features:
   unemployed  crimerate population
1        18.2 0.27486567      18238
2        28.7 0.22068266       5918
3        33.9 0.61335188       2876
4        24.3 0.29258060      21929
5        15.7 0.15206547      17841
6         8.2 0.09292786      39493
7        28.6 0.46368524      11717
8         8.4 0.11802500      25681
9        23.4 0.29192164      25983
10        8.7 0.13058500      54991
                         geometry
1  MULTIPOLYGON (((-87.60914 4...
2  MULTIPOLYGON (((-87.59215 4...
3  MULTIPOLYGON (((-87.6288 41...
4  MULTIPOLYGON (((-87.60671 4...
5  MULTIPOLYGON (((-87.59215 4...
6  MULTIPOLYGON (((-87.67441 4...
7  MULTIPOLYGON (((-87.60604 4...
8  MULTIPOLYGON (((-87.58038 4...
9  MULTIPOLYGON (((-87.57714 4...
10 MULTIPOLYGON (((-87.65456 4...

lm_error <- errorsarlm(population ~ unemployed + crimerate + population,
              data = model,
              listw = chicago_list,
              zero.policy = TRUE, 
              na.action = na.omit)

Warning: the response appeared on the right-hand side and was droppedWarning: problem with term 3 in model.matrix: no columns are assignedWarning: inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
  reciprocal condition number = 4.33906e-18 - using numerical Hessian.

summary(lm1 <- lm(population ~ crimerate + unemployed, data = model))

Call:
lm(formula = population ~ crimerate + unemployed, data = model)

Residuals:
   Min     1Q Median     3Q    Max 
-31365 -16634  -3391  14017  69637 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  48588.4     5629.7   8.631  8.3e-13 ***
crimerate     7012.6    24592.7   0.285   0.7763    
unemployed    -980.4      456.7  -2.147   0.0351 *  
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 21600 on 74 degrees of freedom
Multiple R-squared:  0.09113,   Adjusted R-squared:  0.06657 
F-statistic:  3.71 on 2 and 74 DF,  p-value: 0.02914

#run the Moran's I test for regression residuals
lm.morantest(lm1, chicago_list)

    Global Moran I for regression residuals

data:  
model: lm(formula = population ~ crimerate +
unemployed, data = model)
weights: chicago_list

Moran I statistic standard deviate = 3.1539,
p-value = 0.0008056
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.195748533     -0.026580278      0.004969411 

seResiduals <- rep(0, length(model$population))

resIndex <- lm_error$residuals %>% names() %>% as.integer();

seResiduals[resIndex] <- lm_error$residuals
chicago_list %>%
  moran.test(seResiduals, ., zero.policy = TRUE) 

    Moran I test under randomisation

data:  seResiduals  
weights: .    

Moran I statistic standard deviate = -0.25804,
p-value = 0.6018
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     -0.031597282      -0.013157895       0.005106447 

Now after doing the error model, our P value is much higher and our Morans I is extremely close to 0, which together would indicate that neighboring chicago areas are independent of eachother after the transormation of the dataset.

I chose to do a error model because of the spatial autocorrelation found to be present in the original analysis. My P value was much smaller than 0.001, so it meant that the idea of independence or randmoness between neighboring areas wasn’t really a possibility, which means there was some auto correlation going on. This P value being as small as it was helped me decide that the error model would be a good fit since these models tend to be used when there is spacial correlation.

Unemployment had a negative and statisitcally significant impact on the population, while Crimerate didn’t as much. I was expecting it to be the opposite, where crimerate would have a greater impact on the population but surprisingly it was the other way around.