#Assignment 4- linear model is Moran's I and interpret spatial regression model 
#Research Question: What factors have the highest correlation to poverty in Chicago? 

library(sf)     

## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(dplyr)   

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(spData) 

## 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)
library(ggthemes)
library(spdep)
library(spatialreg)

## Loading 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)

## Loading required package: robustbase

## Loading required package: sp

## Loading required package: Rcpp

## Welcome to GWmodel version 2.3-1.

library(tidyr)

## 
## Attaching package: 'tidyr'

## The following objects are masked from 'package:Matrix':
## 
##     expand, pack, unpack

setwd("/Users/Aiesha/Desktop/DIDA 370/chicago_airbnb")
chicago <- st_read("airbnb_Chicago 2015.shp") 

## Reading layer `airbnb_Chicago 2015' from data source 
##   `/Users/Aiesha/Desktop/DIDA 370/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=as.numeric(poverty)))+
  scale_fill_steps(
    name = "Poverty",
    low = "lightsteelblue1",
    high = "tomato1",
    n.breaks = 5,
    show.limits = T)+
  theme_void()

#Simple Regression Model, poverty is dependent 

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

summary(lm1 <- lm(poverty ~ unemployed + without_hs + dependency + num_spots, data = model)) 

## 
## Call:
## lm(formula = poverty ~ unemployed + without_hs + dependency + 
##     num_spots, data = model)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5526  -3.7946  -0.2939   3.8109  17.0702 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.308938   5.952610   0.892  0.37544    
## unemployed   1.350643   0.135636   9.958 3.53e-15 ***
## without_hs   0.213038   0.069361   3.071  0.00301 ** 
## dependency  -0.259300   0.183618  -1.412  0.16221    
## num_spots    0.009649   0.007601   1.269  0.20839    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.451 on 72 degrees of freedom
## Multiple R-squared:  0.7034, Adjusted R-squared:  0.6869 
## F-statistic: 42.69 on 4 and 72 DF,  p-value: < 2.2e-16

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

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

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = poverty ~ unemployed + without_hs + dependency +
## num_spots, data = model)
## weights: chicago_list
## 
## Moran I statistic standard deviate = 3.7122, p-value = 0.0001027
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.216016690     -0.039309965      0.004730657

lm_error <- errorsarlm(poverty ~ unemployed + without_hs + dependency + num_spots,
              data = model,
              listw = chicago_list,
              zero.policy = TRUE, 
              na.action = na.omit)

summary(lm_error)

## 
## Call:errorsarlm(formula = poverty ~ unemployed + without_hs + dependency + 
##     num_spots, data = model, listw = chicago_list, na.action = na.omit, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -13.34566  -4.23562  -0.42275   2.75451  12.64263 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -1.3336792  6.0290520 -0.2212 0.8249299
## unemployed   1.0418742  0.1478779  7.0455 1.848e-12
## without_hs   0.3057668  0.0800513  3.8196 0.0001336
## dependency   0.0269340  0.1853243  0.1453 0.8844466
## num_spots    0.0028586  0.0072320  0.3953 0.6926451
## 
## Lambda: 0.64063, LR test value: 12.733, p-value: 0.00035935
## Asymptotic standard error: 0.10315
##     z-value: 6.2109, p-value: 5.2677e-10
## Wald statistic: 38.575, p-value: 5.2677e-10
## 
## Log likelihood: -243.8556 for error model
## ML residual variance (sigma squared): 29.672, (sigma: 5.4472)
## Number of observations: 77 
## Number of parameters estimated: 7 
## AIC: 501.71, (AIC for lm: 512.44)

# From this we can see that lambda is 0.64 which is statistically significant because it is closer to 1. This suggests that there is significant spatial autocorrelation in the residuals. In addition, the AIC value of 501.71 which is lower than the AIC for the linear model of 512.44 which means that there was in improvement in the model. 

#Spatial autocorrelation in the residuals of this model

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

## 
##  Moran I test under randomisation
## 
## data:  seResiduals  
## weights: .    
## 
## Moran I statistic standard deviate = -0.29478, p-value = 0.6159
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.034384258      -0.013157895       0.005185026

#Since the p-value is 0.6159, which is greater than 0.05, this means the Moran test is is insignificant. 

#Spatial lag model
lm_lag <- lagsarlm(poverty ~ unemployed + without_hs + dependency + num_spots,
              data = model,
              listw = chicago_list,
              zero.policy = TRUE, 
              na.action = na.omit)

summary(lm_lag)

## 
## Call:lagsarlm(formula = poverty ~ unemployed + without_hs + dependency + 
##     num_spots, data = model, listw = chicago_list, na.action = na.omit, 
##     zero.policy = TRUE)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -12.81473  -3.82159  -0.80766   3.40439  15.32690 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  0.8903049  5.6148974  0.1586  0.874015
## unemployed   1.1205339  0.1452171  7.7163 1.199e-14
## without_hs   0.2002305  0.0646313  3.0980  0.001948
## dependency  -0.1998272  0.1689150 -1.1830  0.236807
## num_spots    0.0094193  0.0069885  1.3478  0.177715
## 
## Rho: 0.27716, LR test value: 6.6178, p-value: 0.010097
## Asymptotic standard error: 0.10419
##     z-value: 2.6601, p-value: 0.007812
## Wald statistic: 7.0761, p-value: 0.007812
## 
## Log likelihood: -246.913 for lag model
## ML residual variance (sigma squared): 35.127, (sigma: 5.9268)
## Number of observations: 77 
## Number of parameters estimated: 7 
## AIC: NA (not available for weighted model), (AIC for lm: 512.44)
## LM test for residual autocorrelation
## test value: 2.0017, p-value: 0.15712

#Check for spatial autocorrelation in these residuals
seResiduals <- rep(0, length(model$poverty))
resIndex <- lm_lag$residuals %>% names() %>% as.integer();
seResiduals[resIndex] <- lm_lag$residuals

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

## 
##  Moran I test under randomisation
## 
## data:  seResiduals  
## weights: .    
## 
## Moran I statistic standard deviate = 1.2666, p-value = 0.1027
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.077755919      -0.013157895       0.005152423

LM <- lm.LMtests(lm1, chicago_list,test = "all")
LM

## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = poverty ~ unemployed + without_hs + dependency +
## num_spots, data = model)
## weights: chicago_list
## 
## LMerr = 8.4306, df = 1, p-value = 0.00369
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = poverty ~ unemployed + without_hs + dependency +
## num_spots, data = model)
## weights: chicago_list
## 
## LMlag = 7.159, df = 1, p-value = 0.007459
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = poverty ~ unemployed + without_hs + dependency +
## num_spots, data = model)
## weights: chicago_list
## 
## RLMerr = 2.2136, df = 1, p-value = 0.1368
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = poverty ~ unemployed + without_hs + dependency +
## num_spots, data = model)
## weights: chicago_list
## 
## RLMlag = 0.94203, df = 1, p-value = 0.3318
## 
## 
##  Lagrange multiplier diagnostics for spatial dependence
## 
## data:  
## model: lm(formula = poverty ~ unemployed + without_hs + dependency +
## num_spots, data = model)
## weights: chicago_list
## 
## SARMA = 9.3726, df = 2, p-value = 0.009221

#This shows that the LMerr test is the most significant because it has the lowest p-value of 0.00369 while LMlag has a p-value of 0.007459. We can also see that the p-value for RLMerr test is 0.1369 which is more significant than the RLMlag test with a p-value of 0.3318. 

#Let's pick the optimum bandwidth for our regression
bwVal <- GWmodel::bw.gwr(poverty ~ unemployed + without_hs + dependency + num_spots,
                         data = model %>% sf::as_Spatial(), 
                         approach = 'AICc', kernel = 'bisquare', 
                         adaptive = TRUE)

## Adaptive bandwidth (number of nearest neighbours): 55 AICc value: 508.0122 
## Adaptive bandwidth (number of nearest neighbours): 42 AICc value: 502.0065 
## Adaptive bandwidth (number of nearest neighbours): 33 AICc value: 502.5195 
## Adaptive bandwidth (number of nearest neighbours): 46 AICc value: 502.8133 
## Adaptive bandwidth (number of nearest neighbours): 37 AICc value: 502.3383 
## Adaptive bandwidth (number of nearest neighbours): 42 AICc value: 502.0065

bwVal

## [1] 42

lm.gwr <- gwr.basic(poverty ~ unemployed + without_hs + dependency + num_spots,
                     data = model %>% sf::as_Spatial(), 
                     bw = bwVal, kernel = "bisquare", adaptive = TRUE)

#Let's check the output of the model!
print(lm.gwr)

##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2024-04-26 18:44:18 
##    Call:
##    gwr.basic(formula = poverty ~ unemployed + without_hs + dependency + 
##     num_spots, data = model %>% sf::as_Spatial(), bw = bwVal, 
##     kernel = "bisquare", adaptive = TRUE)
## 
##    Dependent (y) variable:  poverty
##    Independent variables:  unemployed without_hs dependency num_spots
##    Number of data points: 77
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.5526  -3.7946  -0.2939   3.8109  17.0702 
## 
##    Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
##    (Intercept)  5.308938   5.952610   0.892  0.37544    
##    unemployed   1.350643   0.135636   9.958 3.53e-15 ***
##    without_hs   0.213038   0.069361   3.071  0.00301 ** 
##    dependency  -0.259300   0.183618  -1.412  0.16221    
##    num_spots    0.009649   0.007601   1.269  0.20839    
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 6.451 on 72 degrees of freedom
##    Multiple R-squared: 0.7034
##    Adjusted R-squared: 0.6869 
##    F-statistic: 42.69 on 4 and 72 DF,  p-value: < 2.2e-16 
##    ***Extra Diagnostic information
##    Residual sum of squares: 2996.493
##    Sigma(hat): 6.320858
##    AIC:  512.4438
##    AICc:  513.6438
##    BIC:  475.5694
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: bisquare 
##    Adaptive bandwidth: 42 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: Euclidean distance metric is used.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                      Min.     1st Qu.      Median     3rd Qu.    Max.
##    Intercept  -61.6494894 -23.4452369  -3.5147861   3.4656277 13.4621
##    unemployed   0.8167436   1.0086143   1.1296304   1.6670148  1.9304
##    without_hs  -0.0954411   0.0577141   0.1513663   0.2958841  0.6051
##    dependency  -0.4148439  -0.1978481   0.1693656   0.5116104  1.4683
##    num_spots   -0.0011847   0.0058865   0.0161795   0.1013088  0.3107
##    ************************Diagnostic information*************************
##    Number of data points: 77 
##    Effective number of parameters (2trace(S) - trace(S'S)): 19.63031 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 57.36969 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 502.0065 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 474.6029 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 449.7179 
##    Residual sum of squares: 1749.689 
##    R-square value:  0.8268223 
##    Adjusted R-square value:  0.7665145 
## 
##    ***********************************************************************
##    Program stops at: 2024-04-26 18:44:18

#Add the coefficient for education to the model data
model$gwr_exp_coeff <- lm.gwr$SDF$without_hs
#There's a built in function to plot this, but we need to transform our sf object into an sp object instead
model1 <- as(model, Class = "Spatial")
#plot it!
p1 <- spplot(model1, "gwr_exp_coeff")

p1

model$gwr_exp_coeff <- lm.gwr$SDF$unemployed
model1 <- as(model, Class = "Spatial")
#plot it!
p2 <- spplot(model1, "gwr_exp_coeff")

p2

library(cowplot)

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggthemes':
## 
##     theme_map

plot_grid(p1, p2, 
          labels = c('A) Without HS Diploma', 'B) Unemployed'),
          hjust = -0.3, vjust = 4, label_size = 10)

Based on the results of my model, unemployed and without HS diploma are statistically significant and have on impact on poverty. They also have some of the lowest p-values with unemployed being more statistically significant than a HS diploma.