R Notebook

library(sf)     

## Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; 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)
#Let's load in some packages and data
setwd("C:/Users/konst/Desktop/DIDA 370/chicago_airbnb")
chicago <- st_read("airbnb_Chicago 2015.shp")

## Reading layer `airbnb_Chicago 2015' from data source 
##   `C:\Users\konst\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

#Let's plot the data first

ggplot()+
  geom_sf(data = chicago, aes(fill=crowded))+
  scale_fill_steps(
    name = "Percent of Crowded Housing",
    low = "lightsteelblue1",
    high = "tomato1",
    n.breaks = 4,
    show.limits = T)+
  theme_void()

library(spdep)

#This code creates the weights
chicago_list <- chicago %>% 
  #first convert polygons to a neighbor object
  #R needs to know a unique ID for each neighborhood    
  #(cartodb_id)
  #Neighbors determined using Queen weights
  poly2nb(st_geometry(chicago)) %>% 
  #zero.policy just tells R not to remove 0 values
  nb2listw(zero.policy = TRUE) 

#Now we can calculate the "Global Moran's I" for the data
#start with the weights, then give R the corresponding variable
chicago_list %>% 
  moran.test(chicago$crowded, ., zero.policy = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  chicago$crowded  
## weights: .    
## 
## Moran I statistic standard deviate = 6.52, p-value = 3.516e-11
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.452376001      -0.013157895       0.005098117

moran.plot(chicago$crowded, 
           chicago_list, 
           zero.policy = TRUE, 
           labels = F,
           xlab = 'Percent of Crowded Housing',
           ylab = 'Percent Crowded hOusing of Neighbors',
                   pch=20)

lisaRslt <- localmoran(chicago$crowded, chicago_list, 
            zero.policy = TRUE, na.action = na.omit)

# Now we can derive the cluster/outlier types (COType in ArcGIS term) for each spatial feature in the data
library(magrittr)
library(tidyr)

## 
## Attaching package: 'tidyr'

## The following object is masked from 'package:magrittr':
## 
##     extract

significanceLevel <- 0.05; # 95% confidence
meanVal <- mean(chicago$crowded);

lisaRslt %<>% as_tibble() %>%
  set_colnames(c("Ii","E.Ii","Var.Ii","Z.Ii","Pr(z > 0)")) %>%
  #Create a new categorical value
  #This variable will call insignificant units "insignificant"
  #Outliers are determined by LI <0 - units higher than the mean are HL and units lower than the mean are LH
  #Clusters are units with LI >= 0, and classified as HH or LL based on whether they are above or below the mean
  mutate(coType = case_when(
  `Pr(z > 0)` > 0.05 ~ "Insignificant",
  `Pr(z > 0)` <= 0.05 & Ii >= 0 & chicago$crowded >= meanVal ~ "HH",
  `Pr(z > 0)` <= 0.05 & Ii >= 0 & chicago$crowded < meanVal ~ "LL",
  `Pr(z > 0)` <= 0.05 & Ii < 0 & chicago$crowded >= meanVal ~ "HL",
  `Pr(z > 0)` <= 0.05 & Ii < 0 & chicago$crowded < meanVal ~ "LH"))

# Now add this coType to original sf data
chicago$coType <- lisaRslt$coType %>% replace_na("Insignificant")

#Now we'll plot it!
ggplot(chicago) +
  geom_sf(aes(fill=coType),color = 'black') +
  scale_fill_manual(values = c('red','grey','NA','blue','cyan'), name='Clusters & \nOutliers') +
  labs(title = "Percent of Crowded Housing in Chicago")+
  theme_minimal()