R Notebook
Last week, we started talking about spatial autocorrelation and how to measure it. This week, we’ll talk more specifically about what to do if you want to run a regression model that has spatial dependence. First, we’ll start with the simplest model: a spatial error model. We use this model when we believe that our error term is spatially correlated.
#load in data and packages first
library(sf) ## Warning: package 'sf' was built under R version 4.3.2## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUElibrary(dplyr) ## Warning: package 'dplyr' was built under R version 4.3.3##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(spData) ## Warning: package 'spData' was built under R version 4.3.2## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(ggplot2)## Warning: package 'ggplot2' was built under R version 4.3.3library(ggthemes)
library(spdep)## Warning: package 'spdep' was built under R version 4.3.3library(spatialreg)## Warning: package 'spatialreg' was built under R version 4.3.3## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 4.3.2##
## 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.ZeroPolicyOptionlibrary(GWmodel)## Warning: package 'GWmodel' was built under R version 4.3.3## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 4.3.3## Loading required package: sp## Warning: package 'sp' was built under R version 4.3.2## Loading required package: Rcpp## Welcome to GWmodel version 2.3-2.library(tidyr)## Warning: package 'tidyr' was built under R version 4.3.3##
## Attaching package: 'tidyr'## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpack#Let's load in some packages and data
setwd("C:/Users/amode/OneDrive/Documents/DIDA370/chicago_airbnb")
chic <- st_read("airbnb_Chicago 2015.shp")## Reading layer `airbnb_Chicago 2015' from data source
## `C:\Users\amode\OneDrive\Documents\DIDA370\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#As always, let's plot the data first
ggplot()+
geom_sf(data = chic, aes(fill=as.numeric(poverty)))+
scale_fill_steps(
name = "Inequality",
low = "lightsteelblue1",
high = "tomato1",
n.breaks = 5,
show.limits = T)+
theme_void()
RESEARCH QUESTION: Does hardship, crime, or without high school rates
effect unemployment rates in Chicago?
#create our variables for the model:
model <- chic %>%
select(c(poverty, num_crimes, harship_in, without_hs, geometry)) %>%
filter(!poverty %in% "NA") %>%
mutate(poverty = as.numeric(poverty))
#Run a basic model
summary(lm1 <- lm(poverty ~ num_crimes + harship_in + without_hs, data = model))##
## Call:
## lm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.4884 -3.7719 -0.7289 3.5364 15.2060
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.0167604 1.4205383 4.940 4.81e-06 ***
## num_crimes 0.0002010 0.0001082 1.858 0.0672 .
## harship_in 0.5054880 0.0371429 13.609 < 2e-16 ***
## without_hs -0.5730250 0.0889068 -6.445 1.09e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.371 on 73 degrees of freedom
## Multiple R-squared: 0.7916, Adjusted R-squared: 0.783
## F-statistic: 92.43 on 3 and 73 DF, p-value: < 2.2e-16#Let's test for spatial autocorrelation in the residuals!
#Create the weights first
chic_list <- model %>%
poly2nb(st_geometry(model)) %>%
nb2listw(zero.policy = TRUE)
#run the Moran's I test for regression residuals
lm.morantest(lm1, chic_list)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model)
## weights: chic_list
##
## Moran I statistic standard deviate = 7.1458, p-value = 4.473e-13
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.462498501 -0.033672615 0.004821243The p-value in this Moran test is weakly significant, which means there is some evidence of spatial autocorrelation in the residuals
lm_error <- errorsarlm(poverty ~ num_crimes + harship_in + without_hs,
data = model,
listw = chic_list,
zero.policy = TRUE,
na.action = na.omit)
summary(lm_error)##
## Call:errorsarlm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model, listw = chic_list, na.action = na.omit, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.46844 -2.83613 -0.66768 2.38386 14.87910
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.5442e+00 2.0486e+00 2.7063 0.006803
## num_crimes 7.0210e-05 9.2334e-05 0.7604 0.447021
## harship_in 5.3690e-01 3.9277e-02 13.6695 < 2.2e-16
## without_hs -5.2320e-01 9.4885e-02 -5.5141 3.505e-08
##
## Lambda: 0.7366, LR test value: 35.225, p-value: 2.9372e-09
## Asymptotic standard error: 0.084296
## z-value: 8.7382, p-value: < 2.22e-16
## Wald statistic: 76.357, p-value: < 2.22e-16
##
## Log likelihood: -219.0248 for error model
## ML residual variance (sigma squared): 14.876, (sigma: 3.857)
## Number of observations: 77
## Number of parameters estimated: 6
## AIC: 450.05, (AIC for lm: 483.27)Since Lamda is insignificant we should not use this model
#This code just allows me to manually extract the residual values
#Normally, R wants an LM object for this test -
#Running a spatial error model does not output an lm object, which is why we need to do this
#Create a blank vector the length of the variable
seResiduals <- rep(0, length(model$poverty))
#create an index based on the order of the residuals in the lm_error object
resIndex <- lm_error$residuals %>% names() %>% as.integer();
#use the index to fill in the seResiduals object with corresponding residual values
seResiduals[resIndex] <- lm_error$residuals
#perform the Moran Test
chic_list %>%
moran.test(seResiduals, ., zero.policy = TRUE) ##
## Moran I test under randomisation
##
## data: seResiduals
## weights: .
##
## Moran I statistic standard deviate = -0.46932, p-value = 0.6806
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.04648430 -0.01315789 0.00504237The linear error model is not significant.
lm_lag <- lagsarlm(poverty ~ num_crimes + harship_in + without_hs,
data = model,
listw = chic_list,
zero.policy = TRUE,
na.action = na.omit)
summary(lm_lag)##
## Call:lagsarlm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model, listw = chic_list, na.action = na.omit, zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.2137 -3.3077 -1.0158 3.8907 17.1961
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.92318523 1.86512855 2.1034 0.03543
## num_crimes 0.00016931 0.00010147 1.6686 0.09520
## harship_in 0.44728541 0.04608668 9.7053 < 2.2e-16
## without_hs -0.49649799 0.09083357 -5.4660 4.603e-08
##
## Rho: 0.21044, LR test value: 5.0977, p-value: 0.023957
## Asymptotic standard error: 0.094539
## z-value: 2.226, p-value: 0.026015
## Wald statistic: 4.955, p-value: 0.026015
##
## Log likelihood: -234.0885 for lag model
## ML residual variance (sigma squared): 25.357, (sigma: 5.0355)
## Number of observations: 77
## Number of parameters estimated: 6
## AIC: 480.18, (AIC for lm: 483.27)
## LM test for residual autocorrelation
## test value: 26.263, p-value: 2.9796e-07The model shows a lot of potential due to its p-value of .01
seResiduals <- rep(0, length(model$poverty))
resIndex <- lm_lag$residuals %>% names() %>% as.integer();
seResiduals[resIndex] <- lm_lag$residuals
#perform the Moran Test
chic_list %>%
moran.test(seResiduals, ., zero.policy = TRUE) ##
## Moran I test under randomisation
##
## data: seResiduals
## weights: .
##
## Moran I statistic standard deviate = 4.4918, p-value = 3.532e-06
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.307990182 -0.013157895 0.005111864The p-value here is not significant, so there is no evidence of spatial autocorrelation in the residuals.
LM <- lm.LMtests(lm1, chic_list, test = "all")## Please update scripts to use lm.RStests in place of lm.LMtestsLM##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model)
## test weights: listw
##
## RSerr = 38.646, df = 1, p-value = 5.081e-10
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model)
## test weights: listw
##
## RSlag = 5.8264, df = 1, p-value = 0.01579
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model)
## test weights: listw
##
## adjRSerr = 35.293, df = 1, p-value = 2.836e-09
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model)
## test weights: listw
##
## adjRSlag = 2.4737, df = 1, p-value = 0.1158
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = poverty ~ num_crimes + harship_in + without_hs,
## data = model)
## test weights: listw
##
## SARMA = 41.12, df = 2, p-value = 1.178e-09the R Slag and the adjRSlag test are the only one that are significant. this reinforces the idea that the lag model is the best fit for this data and varibales.