Assignment 4
This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.
#Load in data and packages first
library(sf) ## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUElibrary(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, unionlibrary(spData)
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.ZeroPolicyOptionlibrary(GWmodel)## Loading required package: robustbase## Loading required package: sp## Loading required package: Rcpp## Welcome to GWmodel version 2.3-2.library(tidyr)##
## Attaching package: 'tidyr'## The following objects are masked from 'package:Matrix':
##
## expand, pack, unpacklibrary(Gini)
setwd("/Users/jasmeetkaur/Downloads/DIDA 370/chicago_airbnb")
chicago <- st_read("airbnb_Chicago 2015.shp")## Reading layer `airbnb_Chicago 2015' from data source
## `/Users/jasmeetkaur/Downloads/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#First plot the data (variable = poverty)
ggplot()+
geom_sf(data = chicago, aes(fill=poverty))+
scale_fill_steps(
name = "Poverty Rates in Chicago",
low = "lightsteelblue1",
high = "tomato1",
n.breaks = 5,
show.limits = T)+
theme_void()
#Create a simple regression model
model <- chicago %>%
mutate(pov = poverty/population,
unemp = unemployed/population,
income = income_pc/population,
crimes = num_crimes/population) %>%
select(c(AREAID, response_r, accept_r, rev_rating, price_pp, room_type, num_spots, poverty, crowded,
dependency, without_hs, unemployed, income_pc, harship_in,
num_crimes, num_theft, population, geometry)) %>%
filter(!AREAID %in% "NA") %>%
mutate(AREAID = as.numeric(AREAID))
#Run a basic model
summary(lm1 <- lm(AREAID ~ unemployed + poverty + income_pc + num_crimes, data = model))##
## Call:
## lm(formula = AREAID ~ unemployed + poverty + income_pc + num_crimes,
## data = model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.524 -16.047 -1.211 14.762 47.773
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.9442052 10.6215451 4.137 9.41e-05 ***
## unemployed 1.6595643 0.5516481 3.008 0.00362 **
## poverty -0.8454322 0.3387132 -2.496 0.01485 *
## income_pc -0.0004019 0.0002061 -1.950 0.05503 .
## num_crimes -0.0002593 0.0004219 -0.615 0.54078
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19.75 on 72 degrees of freedom
## Multiple R-squared: 0.2616, Adjusted R-squared: 0.2206
## F-statistic: 6.378 on 4 and 72 DF, p-value: 0.0001887#Test for spatial autocorrelation in the residuals
#Specifically test autocorrelation in error term (Moran's I test)
#Create the weights first
#Need weight objects
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 = AREAID ~ unemployed + poverty + income_pc +
## num_crimes, data = model)
## weights: chicago_list
##
## Moran I statistic standard deviate = 6.8295, p-value = 4.261e-12
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.433478458 -0.037033315 0.004746396#General Notes:
#When it's greater than zero: positive
#When it's negative: negative
#p-value: reject null hypothesis (many many zeros, so small p-value; and because less than 5%, we can reject null (not seeing random spatial)
#Results:
#Moran's I: 0.4
#Spatial autocorrelation: positive#Run a spatial error model to control for it
#Need to try a different model now
lm_error <- errorsarlm(AREAID ~ unemployed + poverty + income_pc + num_crimes,
data = model,
listw = chicago_list,
zero.policy = TRUE,
na.action = na.omit)
summary(lm_error)##
## Call:errorsarlm(formula = AREAID ~ unemployed + poverty + income_pc +
## num_crimes, data = model, listw = chicago_list, na.action = na.omit,
## zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -38.0602 -6.1797 -1.0921 6.4077 64.6911
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 52.14895899 12.08221294 4.3162 1.588e-05
## unemployed 0.04140664 0.46867436 0.0883 0.9296
## poverty -0.26363545 0.29687065 -0.8880 0.3745
## income_pc -0.00025164 0.00021740 -1.1575 0.2471
## num_crimes -0.00009081 0.00034231 -0.2653 0.7908
##
## Lambda: 0.72113, LR test value: 34.553, p-value: 4.1474e-09
## Asymptotic standard error: 0.087547
## z-value: 8.237, p-value: 2.2204e-16
## Wald statistic: 67.848, p-value: 2.2204e-16
##
## Log likelihood: -319.1022 for error model
## ML residual variance (sigma squared): 201.84, (sigma: 14.207)
## Number of observations: 77
## Number of parameters estimated: 7
## AIC: 652.2, (AIC for lm: 684.76)#Results:
#p-value: 0.37
#Very close to 0
#Weakly significant
#Null hypothesis: there's no effect of poverty on the residual (the coefficient for poverty = 0 in the population)
#H1: there IS an effect of poverty on the dependent variable (not equal to 0)
#Fail to reject null (there's not enough evidence to conclude that there's a statistical significance of the effect of poverty on the dependent variable)
#Therefore, random #Results:
#Lambda = 0.72 which is statistically significant. This suggests that there is significant spatial autocorrelation in the. residuals. #Here, we can see an AIC value of 652.2, compared to an AIC for the linear model of 684.76. Because the error model's AIC is lower, this suggests an improvement in the model.
#Test the spatial autocorrelation in the residuals of this model
#Create a blank vector the length of the variable
seResiduals <- rep(0, length(model$AREAID))
#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 = -1.6914, p-value = 0.9546
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.128741113 -0.013157895 0.004669992#Residuals should be a random cloud of points in space (we don't want correlation)
#Results:
#We can see that our Moran test is insignificant since our p-value is greater than our significance level of 0.5
#Fail to reject null hypothesis, therefore there's not enough evidence to conclude that there's significant spatial autocorrelation in the residuals #Test whether one or the other is a good fit with the Lagrange Multiplier test
LM <- lm.LMtests(lm1, chicago_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 = AREAID ~ unemployed + poverty + income_pc +
## num_crimes, data = model)
## test weights: listw
##
## RSerr = 33.948, df = 1, p-value = 5.66e-09
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = AREAID ~ unemployed + poverty + income_pc +
## num_crimes, data = model)
## test weights: listw
##
## RSlag = 49.292, df = 1, p-value = 2.205e-12
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = AREAID ~ unemployed + poverty + income_pc +
## num_crimes, data = model)
## test weights: listw
##
## adjRSerr = 3.8327, df = 1, p-value = 0.05026
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = AREAID ~ unemployed + poverty + income_pc +
## num_crimes, data = model)
## test weights: listw
##
## adjRSlag = 19.177, df = 1, p-value = 1.192e-05
##
##
## Rao's score (a.k.a Lagrange multiplier) diagnostics for spatial
## dependence
##
## data:
## model: lm(formula = AREAID ~ unemployed + poverty + income_pc +
## num_crimes, data = model)
## test weights: listw
##
## SARMA = 53.125, df = 2, p-value = 2.911e-12#Pick the optimum bandwidth for the regression
bwVal <- GWmodel::bw.gwr(AREAID ~ unemployed + poverty + income_pc + num_crimes,
data = model %>% sf::as_Spatial(),
approach = 'AICc', kernel = 'bisquare',
adaptive = TRUE)## Adaptive bandwidth (number of nearest neighbours): 55 AICc value: 635.9855
## Adaptive bandwidth (number of nearest neighbours): 42 AICc value: 633.3209
## Adaptive bandwidth (number of nearest neighbours): 33 AICc value: 635.5891
## Adaptive bandwidth (number of nearest neighbours): 46 AICc value: 632.8731
## Adaptive bandwidth (number of nearest neighbours): 50 AICc value: 633.652
## Adaptive bandwidth (number of nearest neighbours): 44 AICc value: 633.0971
## Adaptive bandwidth (number of nearest neighbours): 47 AICc value: 632.778
## Adaptive bandwidth (number of nearest neighbours): 48 AICc value: 633.2201
## Adaptive bandwidth (number of nearest neighbours): 46 AICc value: 632.8731
## Adaptive bandwidth (number of nearest neighbours): 47 AICc value: 632.778bwVal## [1] 47#Now run the model based on this
lm.gwr <- gwr.basic(AREAID ~ unemployed + poverty + income_pc + num_crimes,
data = model %>% sf::as_Spatial(),
bw = bwVal, kernel = "bisquare", adaptive = TRUE)
#Check the output of the model
print(lm.gwr)## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2024-05-01 10:26:56.310643
## Call:
## gwr.basic(formula = AREAID ~ unemployed + poverty + income_pc +
## num_crimes, data = model %>% sf::as_Spatial(), bw = bwVal,
## kernel = "bisquare", adaptive = TRUE)
##
## Dependent (y) variable: AREAID
## Independent variables: unemployed poverty income_pc num_crimes
## Number of data points: 77
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -34.524 -16.047 -1.211 14.762 47.773
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 43.9442052 10.6215451 4.137 9.41e-05 ***
## unemployed 1.6595643 0.5516481 3.008 0.00362 **
## poverty -0.8454322 0.3387132 -2.496 0.01485 *
## income_pc -0.0004019 0.0002061 -1.950 0.05503 .
## num_crimes -0.0002593 0.0004219 -0.615 0.54078
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 19.75 on 72 degrees of freedom
## Multiple R-squared: 0.2616
## Adjusted R-squared: 0.2206
## F-statistic: 6.378 on 4 and 72 DF, p-value: 0.0001887
## ***Extra Diagnostic information
## Residual sum of squares: 28086.43
## Sigma(hat): 19.35163
## AIC: 684.7577
## AICc: 685.9577
## BIC: 647.8834
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: bisquare
## Adaptive bandwidth: 47 (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 3.2834e+00 2.7103e+01 7.1212e+01 7.7152e+01 82.5058
## unemployed -2.4946e+00 2.1574e-01 5.9227e-01 1.0672e+00 1.7842
## poverty -1.5299e+00 -1.0936e+00 -7.4628e-01 1.5307e-01 0.6302
## income_pc -8.7518e-04 -6.6329e-04 -4.4458e-04 -1.7976e-04 0.0001
## num_crimes -7.3418e-04 -1.9225e-04 7.4912e-05 4.9944e-04 0.0010
## ************************Diagnostic information*************************
## Number of data points: 77
## Effective number of parameters (2trace(S) - trace(S'S)): 18.62269
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 58.37731
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 632.778
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 607.0131
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 579.7205
## Residual sum of squares: 9859.168
## R-square value: 0.7408074
## Adjusted R-square value: 0.6566824
##
## ***********************************************************************
## Program stops at: 2024-05-01 10:26:56.367214library(cowplot)##
## Attaching package: 'cowplot'## The following object is masked from 'package:ggthemes':
##
## theme_map#Map it (variable = poverty)
#Add the coefficient for poverty to the model data
model$gwr_exp_coeff <- lm.gwr$SDF$poverty
#There's a built in function to plot this, but first transform the sf object into an sp object instead
model1 <- as(model, Class = "Spatial")
#Now plot it
p1 <- spplot(model1, "gwr_exp_coeff")
plot_grid(p1,
labels = c('A) Poverty'),
hjust = -0.3, vjust = 4, label_size = 10)
library(cowplot)
#Map another variable (variable = unemployment)
model$gwr_exp_coeff <- lm.gwr$SDF$unemployed
model1 <- as(model, Class = "Spatial")
#Plot it
p2 <- spplot(model1, "gwr_exp_coeff")
plot_grid(p2,
labels = c('B) Unemployed'),
hjust = -0.3, vjust = 4, label_size = 10)
#Plot both maps side by side to compare
library(cowplot)
cowplot::plot_grid(p1, p2,
labels = c('A) Poverty', 'B) Unemployed'),
hjust = -0.3, vjust = 4, label_size = 10)
#What patterns do you notice here?:
#Looking at both maps, the poverty rate in the northern areas of Chicago are on the lower side, however there seems to be a very high percentage of people that are unemployed in these same areas. In map (A) we can see that as we travel further south, the poverty rate increases and then slightly decreases as we reach the southern tip of Chicago. As for map (B), we see that the unemployment rate in the southern area of Chicago is more or less 0. As we travel north the percentage fluctuates ending up at 2.5% at the most northern region of Chicago.