R Notebook

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.

library(sf)     
library(dplyr)   
library(spData) 
library(ggplot2)
library(ggthemes)
#Let's load in some packages and data
setwd("~/Desktop/Spatial R 370/chicago_airbnb")

Warning: The working directory was changed to /Users/kenlu/Desktop/Spatial R 370/chicago_airbnb inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

chicago <- st_read("airbnb_Chicago 2015.shp")

Reading layer `airbnb_Chicago 2015' from data source 
  `/Users/kenlu/Desktop/Spatial R 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=price_pp)) +
    scale_fill_steps(
    name = "Price per person",
    low = "lightgreen",
    high = "red",
    n.breaks = 4,
    show.limits = T)+
  theme_void()

library(spdep)

#This code creates the weights
chicago <- na.omit(chicago)
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$price_pp, ., zero.policy = TRUE)

    Moran I test under randomisation

data:  chicago$price_pp  
weights: .    

Moran I statistic standard deviate = 5.456, p-value = 2.435e-08
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.433421138      -0.015384615       0.006766499 

moran.plot(chicago$price_pp, 
           chicago_list, 
           zero.policy = TRUE, 
           labels = F,
           xlab = 'Price per Person',
           ylab = 'Lagged Price of person (of Neighbors)',
                   pch=20)

lisaRslt <- localmoran(chicago$price_pp, 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)
significanceLevel <- 0.05; # 95% confidence
meanVal <- mean(chicago$price_pp);

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$price_pp >= meanVal ~ "HH",
  `Pr(z > 0)` <= 0.05 & Ii >= 0 & chicago$price_pp < meanVal ~ "LL",

  `Pr(z > 0)` <= 0.05 & Ii < 0 & chicago$price_pp < 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 = 'lightgrey') +
  scale_fill_manual(values = c('red','NA','cyan'), name='Clusters & \nOutliers') +
  labs(title = "Price per person by neighborhood")+
  theme_minimal()

---

---

---

---

---

---

---

---

---

---

---

---

ASSIGNMENT 4

ggplot()+
  geom_sf(data = chicago, aes(fill=price_pp)) +
    scale_fill_steps(
    name = "Price per person",
    low = "lightblue",
    high = "purple",
    n.breaks = 4,
    show.limits = T)+
  theme_void()

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

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

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

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

    Global Moran I for regression residuals

data:  
model: lm(formula = price_pp ~ poverty + crowded + unemployed, data = chicago)
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 

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

summary(lm_error)

Call:errorsarlm(formula = price_pp ~ poverty + crowded + unemployed, 
    data = chicago, 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,
              listw = chicago_list,
              zero.policy = TRUE, 
              na.action = na.omit)

summary(lm_lag)

Call:lagsarlm(formula = price_pp ~ poverty + crowded + unemployed, 
    data = chicago, 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

# We can see which one is better using a Lagrange Multiplier Test
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)
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)
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)
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)
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)
test weights: listw

SARMA = 16.141, df = 2, p-value = 0.0003126

# Finally, Looking at the results of our lag model we can see that 
# Crowdedness is the most significant since the p value is < 0.05, 
# From its coefficient we can tell that as crowdedness increases, price perperson will decrease
# The other two variables are insignificant since the p values are > 0.05

Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Cmd+Option+I.

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Cmd+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.