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==