CDS301_Final_Project

## Warning: package 'ggplot2' was built under R version 4.3.3

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

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

## Warning: package 'sf' was built under R version 4.3.3

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

## Warning: package 'spdep' was built under R version 4.3.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.3.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Warning: package 'tmap' was built under R version 4.3.3

## Warning: package 'purrr' was built under R version 4.3.3

## Warning: package 'lubridate' was built under R version 4.3.3

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1
## ✔ readr     2.1.5

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks plotly::filter(), stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

# the spatial weights + Global Moran’s I - Income and Rent -------------------------------

# 2021
# income(2021)


Census_Joined_Income_2021 <- Census_Joined_2021 %>%
  drop_na(median_house_income)

nb_income_2021 <- poly2nb(Census_Joined_Income_2021, queen = TRUE)
listw_income_2021 <- nb2listw(nb_income_2021, style = "W")

global_moran_income_2021 <- moran.test(Census_Joined_Income_2021$median_house_income, listw = listw_income_2021)

print(global_moran_income_2021)

x_std_2021_income <- as.vector(scale(Census_Joined_Income_2021$median_house_income))
y_lag_2021_income <- lag.listw(listw_income_2021, x_std_2021_income)




# rent(2021)

Census_Joined_Rent_2021 <- Census_Joined_2021 %>%
  drop_na(median_gross_rent)

nb_rent_2021 <- poly2nb(Census_Joined_Rent_2021, queen = TRUE)
listw_rent_2021 <- nb2listw(nb_rent_2021, style = "W")

global_moran_rent_2021 <- moran.test(Census_Joined_Rent_2021$median_gross_rent, listw = listw_rent_2021)
print(global_moran_rent_2021)

x_std_2021_rent <- as.vector(scale(Census_Joined_Rent_2021$median_gross_rent))
y_lag_2021_rent <- lag.listw(listw_rent_2021, x_std_2021_rent)




# 2020
# income(2020)


Census_Joined_Income_2020 <- Census_Joined_2020 %>%
  drop_na(median_house_income)

nb_income_2020 <- poly2nb(Census_Joined_Income_2020, queen = TRUE)
listw_income_2020 <- nb2listw(nb_income_2020, style = "W")

global_moran_income_2020 <- moran.test(Census_Joined_Income_2020$median_house_income, listw = listw_income_2020)

print(global_moran_income_2020)

x_std_2020_income <- as.vector(scale(Census_Joined_Income_2020$median_house_income))
y_lag_2020_income <- lag.listw(listw_income_2020, x_std_2020_income)




# rent(2020)


Census_Joined_Rent_2020 <- Census_Joined_2020 %>%
  drop_na(median_gross_rent)

nb_rent_2020 <- poly2nb(Census_Joined_Rent_2020, queen = TRUE)
listw_rent_2020 <- nb2listw(nb_rent_2020, style = "W")

global_moran_rent_2020 <- moran.test(Census_Joined_Rent_2020$median_gross_rent, listw = listw_rent_2020)
print(global_moran_rent_2020)

x_std_2020_rent <- as.vector(scale(Census_Joined_Rent_2020$median_gross_rent))
y_lag_2020_rent <- lag.listw(listw_rent_2020, x_std_2020_rent)




# 2019
# income(2019)


Census_Joined_Income_2019 <- Census_Joined_2019 %>%
  drop_na(median_house_income)

nb_income_2019 <- poly2nb(Census_Joined_Income_2019, queen = TRUE)
listw_income_2019 <- nb2listw(nb_income_2019, style = "W")

global_moran_income_2019 <- moran.test(Census_Joined_Income_2019$median_house_income, listw = listw_income_2019)

print(global_moran_income_2019)

x_std_2019_income <- as.vector(scale(Census_Joined_Income_2019$median_house_income))
y_lag_2019_income <- lag.listw(listw_income_2019, x_std_2019_income)




# rent(2019)


Census_Joined_Rent_2019 <- Census_Joined_2019 %>%
  drop_na(median_gross_rent)

nb_rent_2019 <- poly2nb(Census_Joined_Rent_2019, queen = TRUE)
listw_rent_2019 <- nb2listw(nb_rent_2019, style = "W")

global_moran_rent_2019 <- moran.test(Census_Joined_Rent_2019$median_gross_rent, listw = listw_rent_2019)
print(global_moran_rent_2019)

x_std_2019_rent <- as.vector(scale(Census_Joined_Rent_2019$median_gross_rent))
y_lag_2019_rent <- lag.listw(listw_rent_2019, x_std_2019_rent)


# 2018
# income(2018)

Census_Joined_Income_2018 <- Census_Joined_2018 %>%
  drop_na(median_house_income)

nb_income_2018 <- poly2nb(Census_Joined_Income_2018, queen = TRUE)
listw_income_2018 <- nb2listw(nb_income_2018, style = "W")

global_moran_income_2018 <- moran.test(Census_Joined_Income_2018$median_house_income, listw = listw_income_2018)

print(global_moran_income_2018)

x_std_2018_income <- as.vector(scale(Census_Joined_Income_2018$median_house_income))
y_lag_2018_income <- lag.listw(listw_income_2018, x_std_2018_income)




# rent(2018)


Census_Joined_Rent_2018 <- Census_Joined_2018 %>%
  drop_na(median_gross_rent)

nb_rent_2018 <- poly2nb(Census_Joined_Rent_2018, queen = TRUE)
listw_rent_2018 <- nb2listw(nb_rent_2018, style = "W")

global_moran_rent_2018 <- moran.test(Census_Joined_Rent_2018$median_gross_rent, listw = listw_rent_2018)
print(global_moran_rent_2018)

x_std_2018_rent <- as.vector(scale(Census_Joined_Rent_2018$median_gross_rent))
y_lag_2018_rent <- lag.listw(listw_rent_2018, x_std_2018_rent)




# 2017
# income(2017)


Census_Joined_Income_2017 <- Census_Joined_2017 %>%
  drop_na(median_house_income)

nb_income_2017 <- poly2nb(Census_Joined_Income_2017, queen = TRUE)
listw_income_2017 <- nb2listw(nb_income_2017, style = "W")

global_moran_income_2017 <- moran.test(Census_Joined_Income_2017$median_house_income, listw = listw_income_2017)

print(global_moran_income_2017)

x_std_2017_income <- as.vector(scale(Census_Joined_Income_2017$median_house_income))
y_lag_2017_income <- lag.listw(listw_income_2017, x_std_2017_income)




# rent(2017) -> spatial unit


Census_Joined_Rent_2017 <- Census_Joined_2017 %>%
  drop_na(median_gross_rent)

nb_rent_2017 <- poly2nb(Census_Joined_Rent_2017, queen = TRUE)
listw_rent_2017 <- nb2listw(nb_rent_2017, style = "W", zero.policy = TRUE) #add zero policy to use this code

global_moran_rent_2017 <- moran.test(Census_Joined_Rent_2017$median_gross_rent, listw = listw_rent_2017)
print(global_moran_rent_2017)

x_std_2017_rent <- as.vector(scale(Census_Joined_Rent_2017$median_gross_rent))
y_lag_2017_rent <- lag.listw(listw_rent_2017, x_std_2017_rent)




# 2016
# income(2016)


Census_Joined_Income_2016 <- Census_Joined_2016 %>%
  drop_na(median_house_income)

nb_income_2016 <- poly2nb(Census_Joined_Income_2016, queen = TRUE)
listw_income_2016 <- nb2listw(nb_income_2016, style = "W")

global_moran_income_2016 <- moran.test(Census_Joined_Income_2016$median_house_income, listw = listw_income_2016)

print(global_moran_income_2016)

x_std_2016_income <- as.vector(scale(Census_Joined_Income_2016$median_house_income))
y_lag_2016_income <- lag.listw(listw_income_2016, x_std_2016_income)




# rent(2016)


Census_Joined_Rent_2016 <- Census_Joined_2016 %>%
  drop_na(median_gross_rent)

nb_rent_2016 <- poly2nb(Census_Joined_Rent_2016, queen = TRUE)
listw_rent_2016 <- nb2listw(nb_rent_2016, style = "W")

global_moran_rent_2016 <- moran.test(Census_Joined_Rent_2016$median_gross_rent, listw = listw_rent_2016)
print(global_moran_rent_2016)

x_std_2016_rent <- as.vector(scale(Census_Joined_Rent_2016$median_gross_rent))
y_lag_2016_rent <- lag.listw(listw_rent_2016, x_std_2016_rent)




# 2015
# income(2015)


Census_Joined_Income_2015 <- Census_Joined_2015 %>%
  drop_na(median_house_income)

nb_income_2015 <- poly2nb(Census_Joined_Income_2015, queen = TRUE)
listw_income_2015 <- nb2listw(nb_income_2015, style = "W")

global_moran_income_2015 <- moran.test(Census_Joined_Income_2015$median_house_income, listw = listw_income_2015)

print(global_moran_income_2015)

x_std_2015_income <- as.vector(scale(Census_Joined_Income_2015$median_house_income))
y_lag_2015_income <- lag.listw(listw_income_2015, x_std_2015_income)




# rent(2015)

Census_Joined_Rent_2015 <- Census_Joined_2015 %>%
  drop_na(median_gross_rent)

nb_rent_2015 <- poly2nb(Census_Joined_Rent_2015, queen = TRUE)
listw_rent_2015 <- nb2listw(nb_rent_2015, style = "W")

global_moran_rent_2015 <- moran.test(Census_Joined_Rent_2015$median_gross_rent, listw = listw_rent_2015)
print(global_moran_rent_2015)

x_std_2015_rent <- as.vector(scale(Census_Joined_Rent_2015$median_gross_rent))
y_lag_2015_rent <- lag.listw(listw_rent_2015, x_std_2015_rent)




# 2014
# income(2014)


Census_Joined_Income_2014 <- Census_Joined_2014 %>%
  drop_na(median_house_income)

nb_income_2014 <- poly2nb(Census_Joined_Income_2014, queen = TRUE)
listw_income_2014 <- nb2listw(nb_income_2014, style = "W")

global_moran_income_2014 <- moran.test(Census_Joined_Income_2014$median_house_income, listw = listw_income_2014)

print(global_moran_income_2014)

x_std_2014_income <- as.vector(scale(Census_Joined_Income_2014$median_house_income))
y_lag_2014_income <- lag.listw(listw_income_2014, x_std_2014_income)




# rent(2014)


Census_Joined_Rent_2014 <- Census_Joined_2014 %>%
  drop_na(median_gross_rent)

nb_rent_2014 <- poly2nb(Census_Joined_Rent_2014, queen = TRUE)
listw_rent_2014 <- nb2listw(nb_rent_2014, style = "W")

global_moran_rent_2014 <- moran.test(Census_Joined_Rent_2014$median_gross_rent, listw = listw_rent_2014)
print(global_moran_rent_2014)

x_std_2014_rent <- as.vector(scale(Census_Joined_Rent_2014$median_gross_rent))
y_lag_2014_rent <- lag.listw(listw_rent_2014, x_std_2014_rent)




# 2013
# income(2013) -> spatial unit


Census_Joined_Income_2013 <- Census_Joined_2013 %>%
  drop_na(median_house_income)

nb_income_2013 <- poly2nb(Census_Joined_Income_2013, queen = TRUE)
listw_income_2013 <- nb2listw(nb_income_2013, style = "W", zero.policy = TRUE) # spatial unit same as 2017

global_moran_income_2013 <- moran.test(Census_Joined_Income_2013$median_house_income, listw = listw_income_2013)

print(global_moran_income_2013)

x_std_2013_income <- as.vector(scale(Census_Joined_Income_2013$median_house_income))
y_lag_2013_income <- lag.listw(listw_income_2013, x_std_2013_income)




# rent(2013)


Census_Joined_Rent_2013 <- Census_Joined_2013 %>%
  drop_na(median_gross_rent)

nb_rent_2013 <- poly2nb(Census_Joined_Rent_2013, queen = TRUE)
listw_rent_2013 <- nb2listw(nb_rent_2013, style = "W")

global_moran_rent_2013 <- moran.test(Census_Joined_Rent_2013$median_gross_rent, listw = listw_rent_2013)
print(global_moran_rent_2013)

x_std_2013_rent <- as.vector(scale(Census_Joined_Rent_2013$median_gross_rent))
y_lag_2013_rent <- lag.listw(listw_rent_2013, x_std_2013_rent)




# 2012
# income(2012)


Census_Joined_Income_2012 <- Census_Joined_2012 %>%
  drop_na(median_house_income)

nb_income_2012 <- poly2nb(Census_Joined_Income_2012, queen = TRUE)
listw_income_2012 <- nb2listw(nb_income_2012, style = "W")

global_moran_income_2012 <- moran.test(Census_Joined_Income_2012$median_house_income, listw = listw_income_2012)

print(global_moran_income_2012)

x_std_2012_income <- as.vector(scale(Census_Joined_Income_2012$median_house_income))
y_lag_2012_income <- lag.listw(listw_income_2012, x_std_2012_income)




# rent(2012) -> spatial unit


Census_Joined_Rent_2012 <- Census_Joined_2012 %>%
  drop_na(median_gross_rent)

nb_rent_2012 <- poly2nb(Census_Joined_Rent_2012, queen = TRUE)
listw_rent_2012 <- nb2listw(nb_rent_2012, style = "W", zero.policy = TRUE) #spatial unit

global_moran_rent_2012 <- moran.test(Census_Joined_Rent_2012$median_gross_rent, listw = listw_rent_2012)
print(global_moran_rent_2012)

x_std_2012_rent <- as.vector(scale(Census_Joined_Rent_2012$median_gross_rent))
y_lag_2012_rent <- lag.listw(listw_rent_2012, x_std_2012_rent)




# 2011
# income(2011)


Census_Joined_Income_2011 <- Census_Joined_2011 %>%
  drop_na(median_house_income)

nb_income_2011 <- poly2nb(Census_Joined_Income_2011, queen = TRUE)
listw_income_2011 <- nb2listw(nb_income_2011, style = "W")

global_moran_income_2011 <- moran.test(Census_Joined_Income_2011$median_house_income, listw = listw_income_2011)

print(global_moran_income_2011)

x_std_2011_income <- as.vector(scale(Census_Joined_Income_2011$median_house_income))
y_lag_2011_income <- lag.listw(listw_income_2011, x_std_2011_income)




# rent(2011)


Census_Joined_Rent_2011 <- Census_Joined_2011 %>%
  drop_na(median_gross_rent)

nb_rent_2011 <- poly2nb(Census_Joined_Rent_2011, queen = TRUE)
listw_rent_2011 <- nb2listw(nb_rent_2011, style = "W")

global_moran_rent_2011 <- moran.test(Census_Joined_Rent_2011$median_gross_rent, listw = listw_rent_2011)
print(global_moran_rent_2011)

x_std_2011_rent <- as.vector(scale(Census_Joined_Rent_2011$median_gross_rent))
y_lag_2011_rent <- lag.listw(listw_rent_2011, x_std_2011_rent)

## Warning: package 'rgeoda' was built under R version 4.3.3

## Loading required package: digest

## Warning: package 'digest' was built under R version 4.3.3

## 
## Attaching package: 'rgeoda'

## The following object is masked from 'package:spdep':
## 
##     skater

## Warning in poly2nb(Census_Joined_Income, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## ℹ tmap mode set to "plot".
## 
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## 
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'breaks', 'palette' (rename to 'values'),
## 'labels' to fill.scale = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [v3->v4] `tm_legend()`: use 'tm_legend()' inside a layer function, e.g.
## 'tm_polygons(..., fill.legend = tm_legend())'

## Warning in poly2nb(Census_Joined_Income_16, queen = TRUE): neighbour object has 6 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## ℹ tmap mode set to "plot".
## 
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## 
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'breaks', 'palette' (rename to 'values'),
## 'labels' to fill.scale = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

## Warning in poly2nb(Census_Joined_Income_11, queen = TRUE): neighbour object has 6 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## ℹ tmap mode set to "plot".
## 
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## 
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'breaks', 'palette' (rename to 'values'),
## 'labels' to fill.scale = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

## Warning in poly2nb(Census_Joined_Rent, queen = TRUE): neighbour object has 3 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## ℹ tmap mode set to "plot".
## 
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## 
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'breaks', 'palette' (rename to 'values'),
## 'labels' to fill.scale = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [tm_borders()] Argument `ltw` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

## Warning in poly2nb(Census_Joined_Rent_16, queen = TRUE): neighbour object has 8 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## ℹ tmap mode set to "plot".
## 
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## 
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'breaks', 'palette' (rename to 'values'),
## 'labels' to fill.scale = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

## Warning in poly2nb(Census_Joined_Rent_11, queen = TRUE): neighbour object has 7 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## ℹ tmap mode set to "plot".
## 
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## 
## [v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of
## the visual variable `fill` namely 'breaks', 'palette' (rename to 'values'),
## 'labels' to fill.scale = tm_scale(<HERE>).
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_borders()`: use `fill_alpha` instead of `alpha`.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

## Warning in poly2nb(Census_Joined_2021, queen = TRUE): neighbour object has 3 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2020, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2019, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2018, queen = TRUE): neighbour object has 2 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2017, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2017, queen = TRUE): neighbour object has 8 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2016, queen = TRUE): neighbour object has 8 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2015, queen = TRUE): neighbour object has 7 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2014, queen = TRUE): neighbour object has 7 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2013, queen = TRUE): neighbour object has 8 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2012, queen = TRUE): some observations have no neighbours;
## if this seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2012, queen = TRUE): neighbour object has 8 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.

## Warning in poly2nb(Census_Joined_2011, queen = TRUE): neighbour object has 7 sub-graphs;
## if this sub-graph count seems unexpected, try increasing the snap argument.