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

setwd("C:/Users/matth/OneDrive/Desktop/Spatial R Class/chicago_airbnb")
chicago <-st_read("airbnb_Chicago 2015.shp")
Reading layer `airbnb_Chicago 2015' from data source 
  `C:\Users\matth\OneDrive\Desktop\Spatial R Class\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
chicago <-na.omit(chicago)
chicago_model <- chicago %>% 
  select(c(price_pp, poverty, crowded, unemployed))
  

summary(lm1 <- lm(price_pp ~ poverty + crowded + unemployed, data = chicago_model))

Call:
lm(formula = price_pp ~ poverty + crowded + unemployed, data = chicago_model)

Residuals:
    Min      1Q  Median      3Q     Max 
-59.057 -16.455  -3.644   7.625  82.727 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 102.0765     9.2667  11.015 3.11e-16 ***
poverty       1.4311     0.5796   2.469 0.016316 *  
crowded      -3.3507     1.0295  -3.255 0.001840 ** 
unemployed   -2.9318     0.8092  -3.623 0.000589 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 29.91 on 62 degrees of freedom
Multiple R-squared:  0.2797,    Adjusted R-squared:  0.2448 
F-statistic: 8.024 on 3 and 62 DF,  p-value: 0.0001339
chicago_list <- chicago_model %>% 
  poly2nb(st_geometry(chicago_model)) %>% 
  nb2listw(zero.policy = TRUE) 

lm.morantest(lm1, chicago_list)

    Global Moran I for regression residuals

data:  
model: lm(formula = price_pp ~ poverty + crowded + unemployed, data =
chicago_model)
weights: chicago_list

Moran I statistic standard deviate = 2.2369, p-value = 0.01265
alternative hypothesis: greater
sample estimates:
Observed Moran I      Expectation         Variance 
     0.138346946     -0.040104638      0.006364367 

#The p-value is under 0.05 therefore this test is significant, and because the Moran
#amount is positive that means that the data tends towards clustering.

lm_error <- errorsarlm(price_pp ~ poverty + crowded + unemployed,
              data = chicago_model,
              listw = chicago_list,
              zero.policy = TRUE, 
              na.action = na.omit)

summary(lm_error)

Call:errorsarlm(formula = price_pp ~ poverty + crowded + unemployed, 
    data = chicago_model, listw = chicago_list, na.action = na.omit, 
    zero.policy = TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-45.1416 -15.7919  -2.4674  13.3003  68.0245 

Type: error 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 89.70343   13.44721  6.6708 2.544e-11
poverty     -0.52949    0.60723 -0.8720    0.3832
crowded     -1.59032    1.20273 -1.3223    0.1861
unemployed  -0.11050    0.88622 -0.1247    0.9008

Lambda: 0.6484, LR test value: 8.2611, p-value: 0.0040504
Asymptotic standard error: 0.10434
    z-value: 6.2142, p-value: 5.1599e-10
Wald statistic: 38.616, p-value: 5.1599e-10

Log likelihood: -311.731 for error model
ML residual variance (sigma squared): 656.61, (sigma: 25.624)
Number of observations: 66 
Number of parameters estimated: 6 
AIC: 635.46, (AIC for lm: 641.72)
lm_lag <- lagsarlm(price_pp ~ poverty + crowded + unemployed,
              data = chicago_model,
              listw = chicago_list,
              zero.policy = TRUE, 
              na.action = na.omit)

summary(lm_lag)

Call:lagsarlm(formula = price_pp ~ poverty + crowded + unemployed, 
    data = chicago_model, listw = chicago_list, na.action = na.omit, 
    zero.policy = TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-50.7371 -13.0582  -4.8799  12.6640  76.1426 

Type: lag 
Coefficients: (asymptotic standard errors) 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept) 52.76053   13.04946  4.0431 5.274e-05
poverty      0.38924    0.49950  0.7793   0.43583
crowded     -2.08203    0.90641 -2.2970   0.02162
unemployed  -1.22398    0.71397 -1.7143   0.08647

Rho: 0.54132, LR test value: 11.375, p-value: 0.00074442
Asymptotic standard error: 0.11609
    z-value: 4.6629, p-value: 3.1171e-06
Wald statistic: 21.743, p-value: 3.1171e-06

Log likelihood: -310.1741 for lag model
ML residual variance (sigma squared): 653.46, (sigma: 25.563)
Number of observations: 66 
Number of parameters estimated: 6 
AIC: 632.35, (AIC for lm: 641.72)
LM test for residual autocorrelation
test value: 2.0126, p-value: 0.156
LM <- lm.LMtests(lm1, chicago_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 = price_pp ~ poverty + crowded + unemployed, data =
chicago_model)
test weights: listw

RSerr = 2.5872, df = 1, p-value = 0.1077


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

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

RSlag = 8.8233, df = 1, p-value = 0.002974


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

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

adjRSerr = 7.3177, df = 1, p-value = 0.006828


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

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

adjRSlag = 13.554, df = 1, p-value = 0.0002318


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

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

SARMA = 16.141, df = 2, p-value = 0.0003126
bwVal <- GWmodel::bw.gwr(price_pp ~ poverty + crowded + unemployed,
                         data = chicago_model %>% sf::as_Spatial(), 
                         approach = 'AICc', kernel = 'bisquare', 
                         adaptive = TRUE)
Adaptive bandwidth (number of nearest neighbours): 48 AICc value: 638.017 
Adaptive bandwidth (number of nearest neighbours): 38 AICc value: 638.5672 
Adaptive bandwidth (number of nearest neighbours): 56 AICc value: 638.2774 
Adaptive bandwidth (number of nearest neighbours): 45 AICc value: 637.7149 
Adaptive bandwidth (number of nearest neighbours): 41 AICc value: 638.1598 
Adaptive bandwidth (number of nearest neighbours): 45 AICc value: 637.7149 
bwVal
[1] 45
lm.gwr <- gwr.basic(price_pp ~ poverty + crowded + unemployed,
                     data = chicago_model %>% sf::as_Spatial(), 
                     bw = bwVal, kernel = "bisquare", adaptive = TRUE)


print(lm.gwr)
   ***********************************************************************
   *                       Package   GWmodel                             *
   ***********************************************************************
   Program starts at: 2024-04-17 10:18:24.162757 
   Call:
   gwr.basic(formula = price_pp ~ poverty + crowded + unemployed, 
    data = chicago_model %>% sf::as_Spatial(), bw = bwVal, kernel = "bisquare", 
    adaptive = TRUE)

   Dependent (y) variable:  price_pp
   Independent variables:  poverty crowded unemployed
   Number of data points: 66
   ***********************************************************************
   *                    Results of Global Regression                     *
   ***********************************************************************

   Call:
    lm(formula = formula, data = data)

   Residuals:
    Min      1Q  Median      3Q     Max 
-59.057 -16.455  -3.644   7.625  82.727 

   Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
   (Intercept) 102.0765     9.2667  11.015 3.11e-16 ***
   poverty       1.4311     0.5796   2.469 0.016316 *  
   crowded      -3.3507     1.0295  -3.255 0.001840 ** 
   unemployed   -2.9318     0.8092  -3.623 0.000589 ***

   ---Significance stars
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
   Residual standard error: 29.91 on 62 degrees of freedom
   Multiple R-squared: 0.2797
   Adjusted R-squared: 0.2448 
   F-statistic: 8.024 on 3 and 62 DF,  p-value: 0.0001339 
   ***Extra Diagnostic information
   Residual sum of squares: 55455.72
   Sigma(hat): 29.4363
   AIC:  641.7231
   AICc:  642.7231
   BIC:  607.6197
   ***********************************************************************
   *          Results of Geographically Weighted Regression              *
   ***********************************************************************

   *********************Model calibration information*********************
   Kernel function: bisquare 
   Adaptive bandwidth: 45 (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   60.79119  96.96510 108.35591 120.33745 147.9672
   poverty     -0.15189   0.77460   1.18388   1.57743   2.1210
   crowded     -6.39058  -5.50934  -3.91818  -2.69459  -2.2522
   unemployed  -3.08011  -2.70235  -2.56686  -2.26710   0.6350
   ************************Diagnostic information*************************
   Number of data points: 66 
   Effective number of parameters (2trace(S) - trace(S'S)): 13.29385 
   Effective degrees of freedom (n-2trace(S) + trace(S'S)): 52.70615 
   AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 637.7149 
   AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 619.7853 
   BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 587.3656 
   Residual sum of squares: 39456.4 
   R-square value:  0.487489 
   Adjusted R-square value:  0.3557204 

   ***********************************************************************
   Program stops at: 2024-04-17 10:18:24.231626 
chicago_model$gwr_exp_coeff <- lm.gwr$SDF$poverty

chicago_model1 <- as(chicago_model, Class = "Spatial")

p1 <- spplot(chicago_model1, "gwr_exp_coeff")

p1

chicago_model$gwr_exp_coeff <- lm.gwr$SDF$unemployed
chicago_model1 <- as(chicago_model, Class = "Spatial")
#plot it!
p2 <- spplot(chicago_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) Poverty', 'B) Unemployed'),
          hjust = -0.3, vjust = 4, label_size = 10)

#The lag model was chosen because the Lagrange Multiplier test told us that the lag model was more significant than the error model.


#Crowdedness has a significant and negative impact on price per person for the airbnbs in Chicago. As the coeffecient is negative and the p value in significant in the lag model.
LS0tDQp0aXRsZTogIlIgTm90ZWJvb2siDQpvdXRwdXQ6IGh0bWxfbm90ZWJvb2sNCi0tLQ0KDQpgYGB7cn0NCmxpYnJhcnkoc2YpICAgICANCmxpYnJhcnkoZHBseXIpICAgDQpsaWJyYXJ5KHNwRGF0YSkgDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KGdndGhlbWVzKQ0KbGlicmFyeShzcGRlcCkNCmxpYnJhcnkoc3BhdGlhbHJlZykNCmxpYnJhcnkoR1dtb2RlbCkNCmxpYnJhcnkodGlkeXIpDQoNCnNldHdkKCJDOi9Vc2Vycy9tYXR0aC9PbmVEcml2ZS9EZXNrdG9wL1NwYXRpYWwgUiBDbGFzcy9jaGljYWdvX2FpcmJuYiIpDQpjaGljYWdvIDwtc3RfcmVhZCgiYWlyYm5iX0NoaWNhZ28gMjAxNS5zaHAiKQ0KY2hpY2FnbyA8LW5hLm9taXQoY2hpY2FnbykNCmBgYA0KYGBge3J9DQoNCiN0aGUgcmVzZWFyY2ggcXVlc3Rpb24gd2lsbCBiZSB0byBmaW5kIHRoZSBlZmZlY3Qgb2YgY3Jvd2RlZG5lc3Mgb24gcHJpY2UgcGVyIHBlcnNvbiBpbiBjaGljYWdvIGFpci1ibmJzDQoNCg0KDQpjaGljYWdvX21vZGVsIDwtIGNoaWNhZ28gJT4lIA0KICBzZWxlY3QoYyhwcmljZV9wcCwgcG92ZXJ0eSwgY3Jvd2RlZCwgdW5lbXBsb3llZCkpDQogIA0KDQpzdW1tYXJ5KGxtMSA8LSBsbShwcmljZV9wcCB+IHBvdmVydHkgKyBjcm93ZGVkICsgdW5lbXBsb3llZCwgZGF0YSA9IGNoaWNhZ29fbW9kZWwpKQ0KYGBgDQpgYGB7cn0NCmNoaWNhZ29fbGlzdCA8LSBjaGljYWdvX21vZGVsICU+JSANCiAgcG9seTJuYihzdF9nZW9tZXRyeShjaGljYWdvX21vZGVsKSkgJT4lIA0KICBuYjJsaXN0dyh6ZXJvLnBvbGljeSA9IFRSVUUpIA0KDQpsbS5tb3JhbnRlc3QobG0xLCBjaGljYWdvX2xpc3QpDQpgYGANCmBgYHtyfQ0KDQojVGhlIHAtdmFsdWUgaXMgdW5kZXIgMC4wNSB0aGVyZWZvcmUgdGhpcyB0ZXN0IGlzIHNpZ25pZmljYW50LCBhbmQgYmVjYXVzZSB0aGUgTW9yYW4NCiNhbW91bnQgaXMgcG9zaXRpdmUgdGhhdCBtZWFucyB0aGF0IHRoZSBkYXRhIHRlbmRzIHRvd2FyZHMgY2x1c3RlcmluZy4NCg0KbG1fZXJyb3IgPC0gZXJyb3JzYXJsbShwcmljZV9wcCB+IHBvdmVydHkgKyBjcm93ZGVkICsgdW5lbXBsb3llZCwNCiAgICAgICAgICAgICAgZGF0YSA9IGNoaWNhZ29fbW9kZWwsDQogICAgICAgICAgICAgIGxpc3R3ID0gY2hpY2Fnb19saXN0LA0KICAgICAgICAgICAgICB6ZXJvLnBvbGljeSA9IFRSVUUsIA0KICAgICAgICAgICAgICBuYS5hY3Rpb24gPSBuYS5vbWl0KQ0KDQpzdW1tYXJ5KGxtX2Vycm9yKQ0KYGBgDQpgYGB7cn0NCmxtX2xhZyA8LSBsYWdzYXJsbShwcmljZV9wcCB+IHBvdmVydHkgKyBjcm93ZGVkICsgdW5lbXBsb3llZCwNCiAgICAgICAgICAgICAgZGF0YSA9IGNoaWNhZ29fbW9kZWwsDQogICAgICAgICAgICAgIGxpc3R3ID0gY2hpY2Fnb19saXN0LA0KICAgICAgICAgICAgICB6ZXJvLnBvbGljeSA9IFRSVUUsIA0KICAgICAgICAgICAgICBuYS5hY3Rpb24gPSBuYS5vbWl0KQ0KDQpzdW1tYXJ5KGxtX2xhZykNCmBgYA0KYGBge3J9DQpMTSA8LSBsbS5MTXRlc3RzKGxtMSwgY2hpY2Fnb19saXN0LCB0ZXN0ID0gImFsbCIpDQpMTQ0KYGBgDQpgYGB7cn0NCmJ3VmFsIDwtIEdXbW9kZWw6OmJ3Lmd3cihwcmljZV9wcCB+IHBvdmVydHkgKyBjcm93ZGVkICsgdW5lbXBsb3llZCwNCiAgICAgICAgICAgICAgICAgICAgICAgICBkYXRhID0gY2hpY2Fnb19tb2RlbCAlPiUgc2Y6OmFzX1NwYXRpYWwoKSwgDQogICAgICAgICAgICAgICAgICAgICAgICAgYXBwcm9hY2ggPSAnQUlDYycsIGtlcm5lbCA9ICdiaXNxdWFyZScsIA0KICAgICAgICAgICAgICAgICAgICAgICAgIGFkYXB0aXZlID0gVFJVRSkNCg0KYndWYWwNCmBgYA0KYGBge3J9DQpsbS5nd3IgPC0gZ3dyLmJhc2ljKHByaWNlX3BwIH4gcG92ZXJ0eSArIGNyb3dkZWQgKyB1bmVtcGxveWVkLA0KICAgICAgICAgICAgICAgICAgICAgZGF0YSA9IGNoaWNhZ29fbW9kZWwgJT4lIHNmOjphc19TcGF0aWFsKCksIA0KICAgICAgICAgICAgICAgICAgICAgYncgPSBid1ZhbCwga2VybmVsID0gImJpc3F1YXJlIiwgYWRhcHRpdmUgPSBUUlVFKQ0KDQoNCnByaW50KGxtLmd3cikNCmBgYA0KYGBge3J9DQpjaGljYWdvX21vZGVsJGd3cl9leHBfY29lZmYgPC0gbG0uZ3dyJFNERiRwb3ZlcnR5DQoNCmNoaWNhZ29fbW9kZWwxIDwtIGFzKGNoaWNhZ29fbW9kZWwsIENsYXNzID0gIlNwYXRpYWwiKQ0KDQpwMSA8LSBzcHBsb3QoY2hpY2Fnb19tb2RlbDEsICJnd3JfZXhwX2NvZWZmIikNCg0KcDENCmBgYA0KYGBge3J9DQpjaGljYWdvX21vZGVsJGd3cl9leHBfY29lZmYgPC0gbG0uZ3dyJFNERiR1bmVtcGxveWVkDQpjaGljYWdvX21vZGVsMSA8LSBhcyhjaGljYWdvX21vZGVsLCBDbGFzcyA9ICJTcGF0aWFsIikNCiNwbG90IGl0IQ0KcDIgPC0gc3BwbG90KGNoaWNhZ29fbW9kZWwxLCAiZ3dyX2V4cF9jb2VmZiIpDQoNCnAyDQpgYGANCmBgYHtyfQ0KbGlicmFyeShjb3dwbG90KQ0KDQpwbG90X2dyaWQocDEsIHAyLCANCiAgICAgICAgICBsYWJlbHMgPSBjKCdBKSBQb3ZlcnR5JywgJ0IpIFVuZW1wbG95ZWQnKSwNCiAgICAgICAgICBoanVzdCA9IC0wLjMsIHZqdXN0ID0gNCwgbGFiZWxfc2l6ZSA9IDEwKQ0KYGBgDQpgYGB7cn0NCiNUaGUgbGFnIG1vZGVsIHdhcyBjaG9zZW4gYmVjYXVzZSB0aGUgTGFncmFuZ2UgTXVsdGlwbGllciB0ZXN0IHRvbGQgdXMgdGhhdCB0aGUgbGFnIG1vZGVsIHdhcyBtb3JlIHNpZ25pZmljYW50IHRoYW4gdGhlIGVycm9yIG1vZGVsLg0KDQoNCiNDcm93ZGVkbmVzcyBoYXMgYSBzaWduaWZpY2FudCBhbmQgbmVnYXRpdmUgaW1wYWN0IG9uIHByaWNlIHBlciBwZXJzb24gZm9yIHRoZSBhaXJibmJzIGluIENoaWNhZ28uIEFzIHRoZSBjb2VmZmVjaWVudCBpcyBuZWdhdGl2ZSBhbmQgdGhlIHAgdmFsdWUgaW4gc2lnbmlmaWNhbnQgaW4gdGhlIGxhZyBtb2RlbC4NCg0KYGBgDQoNCg==