R Notebook

run a spatial regression model

library(sf)     
library(dplyr)   
library(spData) 
library(ggplot2)
library(ggthemes)
library(spdep)
library(spatialreg)
library(GWmodel)
library(tidyr)
#Let's load in some packages and data
setwd("/Users/laurendemarco/Downloads/DIDA370/chicago_airbnb")

Warning: The working directory was changed to /Users/laurendemarco/Downloads/DIDA370/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.

chi_data <- st_read("airbnb_Chicago 2015.shp")

Reading layer `airbnb_Chicago 2015' from data source 
  `/Users/laurendemarco/Downloads/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

#start with a general map of % in poverty in chicago
ggplot()+
  geom_sf(data = chi_data, aes(fill=as.numeric(poverty)))+
  scale_fill_steps(
    name = "% in Poverty",
    low = "lightsteelblue1",
    high = "tomato1",
    n.breaks = 5,
    show.limits = T)+
  labs(title = "Percentage of Poverty in Chicago Neighborhods") +
  theme_void()

Research Question: What is the effect of crowding, unemployment & number of crimes and theft on poverty level in Chicago neighborhoods. Independent Variables: Crowded percentage, unemployed percentage, number of crimes and number of theft Dependent Variable: Poverty percentage

#create our variables for the model:
model <- chi_data %>% 
  select(c(poverty, crowded, unemployed, num_crimes, num_theft, without_hs, 
           dependency, income_pc, population, geometry)) %>% 
  filter(!poverty %in% "NA") %>% 
  mutate(poverty = as.numeric(poverty))

#Run a basic model

summary(lm1 <- lm(poverty ~ crowded + unemployed + num_crimes + num_theft, data = model))

Call:
lm(formula = poverty ~ crowded + unemployed + num_crimes + num_theft, 
    data = model)

Residuals:
     Min       1Q   Median       3Q      Max 
-13.1427  -5.3129   0.2276   4.0086  17.1818 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1.579e+00  2.298e+00  -0.687   0.4942    
crowded      6.526e-01  2.175e-01   3.000   0.0037 ** 
unemployed   1.217e+00  1.215e-01  10.016 2.76e-15 ***
num_crimes  -1.823e-05  2.336e-04  -0.078   0.9380    
num_theft    1.013e-03  8.314e-04   1.219   0.2269    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.553 on 72 degrees of freedom
Multiple R-squared:  0.694, Adjusted R-squared:  0.677 
F-statistic: 40.82 on 4 and 72 DF,  p-value: < 2.2e-16

#Let's test for spatial autocorrelation in the residuals!
chi_list <- model %>% 
  poly2nb(st_geometry(model)) %>% 
  nb2listw(zero.policy = TRUE) 

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

    Global Moran I for regression residuals

data:  
model: lm(formula = poverty ~ crowded + unemployed + num_crimes + num_theft, data = model)
weights: chi_list

Moran I statistic standard deviate = 4.8264, p-value = 6.95e-07
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.295357914     -0.037740465      0.004763105 

#interpretation of the model:
#we have a p-value < 0.05! 
#low p value, this indicates to us that we wanna reject the null hypothesis and suggest this is not just spatial randomness, but we are see sptail autocorrelation

#take another step and run error model
#this is an improvement
lm_error <- errorsarlm(poverty ~ crowded + unemployed + num_crimes + num_theft,
              data = model,
              listw = chi_list,
              zero.policy = TRUE, 
              na.action = na.omit)

summary(lm_error)

Call:errorsarlm(formula = poverty ~ crowded + unemployed + num_crimes + 
    num_theft, data = model, listw = chi_list, na.action = na.omit,     zero.policy = TRUE)

Residuals:
      Min        1Q    Median        3Q       Max 
-13.64199  -3.64360  -0.58219   3.04990  14.07546 

Type: error 
Coefficients: (asymptotic standard errors) 
               Estimate  Std. Error z value  Pr(>|z|)
(Intercept)  2.1855e-01  2.6689e+00  0.0819 0.9347374
crowded      8.4446e-01  2.2158e-01  3.8110 0.0001384
unemployed   1.1132e+00  1.2689e-01  8.7733 < 2.2e-16
num_crimes   9.8838e-05  2.1734e-04  0.4548 0.6492840
num_theft   -1.1861e-04  8.7533e-04 -0.1355 0.8922160

Lambda: 0.59497, LR test value: 15.847, p-value: 6.8681e-05
Asymptotic standard error: 0.11113
    z-value: 5.3538, p-value: 8.6115e-08
Wald statistic: 28.663, p-value: 8.6115e-08

Log likelihood: -243.5054 for error model
ML residual variance (sigma squared): 29.918, (sigma: 5.4698)
Number of observations: 77 
Number of parameters estimated: 7 
AIC: 501.01, (AIC for lm: 514.86)

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
chi_list %>%
  moran.test(seResiduals, ., zero.policy = TRUE) 

    Moran I test under randomisation

data:  seResiduals  
weights: .    

Moran I statistic standard deviate = -0.20458, p-value = 0.581
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
     -0.027859550      -0.013157895       0.005164471 

#this test is in significant, this tell us "that we effectively controlled for the spatial variation in the error term"

#compare with Lag model

lm_lag <- lagsarlm(poverty ~ crowded + unemployed + num_crimes + num_theft,
              data = model,
              listw = chi_list,
              zero.policy = TRUE, 
              na.action = na.omit)

summary(lm_lag)

Call:lagsarlm(formula = poverty ~ crowded + unemployed + num_crimes + 
    num_theft, data = model, listw = chi_list, na.action = na.omit,     zero.policy = TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-14.0479  -3.8431  -0.4222   3.4368  15.7009 

Type: lag 
Coefficients: (asymptotic standard errors) 
               Estimate  Std. Error z value  Pr(>|z|)
(Intercept) -4.7420e+00  2.3439e+00 -2.0231  0.043065
crowded      6.4373e-01  1.9849e-01  3.2431  0.001182
unemployed   1.0104e+00  1.3194e-01  7.6581 1.887e-14
num_crimes  -6.8589e-05  2.1306e-04 -0.3219  0.747516
num_theft    1.0169e-03  7.5760e-04  1.3423  0.179493

Rho: 0.30524, LR test value: 8.0045, p-value: 0.0046661
Asymptotic standard error: 0.10369
    z-value: 2.9438, p-value: 0.0032419
Wald statistic: 8.6661, p-value: 0.0032419

Log likelihood: -247.4266 for lag model
ML residual variance (sigma squared): 35.466, (sigma: 5.9554)
Number of observations: 77 
Number of parameters estimated: 7 
AIC: 508.85, (AIC for lm: 514.86)
LM test for residual autocorrelation
test value: 5.2633, p-value: 0.02178

#signifiance, but not as strong
#AIC is also lower for the error model

seResiduals <- rep(0, length(model$poverty))
resIndex <- lm_lag$residuals %>% names() %>% as.integer();
seResiduals[resIndex] <- lm_lag$residuals

#perform the Moran Test
chi_list %>%
  moran.test(seResiduals, ., zero.policy = TRUE) 

    Moran I test under randomisation

data:  seResiduals  
weights: .    

Moran I statistic standard deviate = 1.9207, p-value = 0.02739
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.124758402      -0.013157895       0.005156063 

#again, not as strong so I think spatial error is better

#can look at this LM test
LM <- lm.LMtests(lm1, chi_list, test = "all")

Please update scripts to use lm.RStests in place of lm.LMtests

LM

    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence

data:  
model: lm(formula = poverty ~ crowded + unemployed + num_crimes + num_theft, data = model)
test weights: listw

RSerr = 15.761, df = 1, p-value = 7.187e-05


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence

data:  
model: lm(formula = poverty ~ crowded + unemployed + num_crimes + num_theft, data = model)
test weights: listw

RSlag = 8.9223, df = 1, p-value = 0.002817


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence

data:  
model: lm(formula = poverty ~ crowded + unemployed + num_crimes + num_theft, data = model)
test weights: listw

adjRSerr = 6.9463, df = 1, p-value = 0.008399


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence

data:  
model: lm(formula = poverty ~ crowded + unemployed + num_crimes + num_theft, data = model)
test weights: listw

adjRSlag = 0.10772, df = 1, p-value = 0.7428


    Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial dependence

data:  
model: lm(formula = poverty ~ crowded + unemployed + num_crimes + num_theft, data = model)
test weights: listw

SARMA = 15.869, df = 2, p-value = 0.0003583

#this confirms that we should use the error model
#the Langronge Test  shows that it is the better choice because 
#adjRSerr, p-value = 0.008399 and ajRSlag, p-value = 0.7428. 
#The error model has a significant p-value while the lag model does not. 

Brief discussion of results from the error model that I decided was the best choice for my data: Crowding and unemployment are the only two significant variables from this regression model. They both had a positive impact on poverty levels, so as they increase… so do poverty levels in Chicago areas.