green_space_demographic_churn.Rmd

---
title: 'Green Space & Residential Demographic Patterns: 1990 to 2010'
author: "Francine Stephens"
date: "Last Updated: `r format(Sys.time(), '%B %d, %Y')`"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
    theme: flatly
    highlight: espresso
---

```{r setup, include=F, warning=F, message=F}
knitr::opts_chunk$set(
    echo = TRUE, 
  message = FALSE,
  warning = FALSE
)

## LIBRARIES
packages <- c(
  "readr",
  "tidyverse",
  "sf",
  "ggplot2", 
  "ggbeeswarm",
  "ggstream",
  "scales",
  "units",
  "ggcorrplot",
  "modelsummary",
  "stargazer",
  "panelr",
  "tigris",
  "censusapi", 
  "tidycensus",
  "RColorBrewer",
  "hrbrthemes",
  "viridis", 
  "wesanderson",
  "leaflet", 
  "tmap"
)
lapply(packages, library, character.only = T)

## PATHS
setwd("~/Projects/public_space")
wd <- getwd()
green_space_path <- "/green_space/gee/"
imperv_space_path <- "/impervious_space/"
census_data_path <- "/census"
shp_repo <- "C:/Users/Franc/Documents/Shapefile_Repository"
shp2010_path <- "/2010USA_CensusGeog_Shp"
cbsa_path <- "/tl_2010_us_cbsa10/"
metdiv_path <- "/tl_2010_us_metdiv10/"
cbg10_path <- "/us_blck_grp_2010/"
places10_path <- "/tl2010_us_place_2010/"

## APIs
census_api_key("99ccb52a629609683f17f804ca875115e3f0804c",  overwrite = T)
Sys.setenv(CENSUS_KEY="99ccb52a629609683f17f804ca875115e3f0804c")

## SHAPEFILE IMPORTS
cbg10 <- st_read(paste0(shp_repo,
                        shp2010_path,
                        cbg10_path,
                        "US_blck_grp_2010.shp"), 
                 quiet = F)

```

## Data Preparation

* EVI scores come from Landsat satellite image mosaics from 1990, 2000, and 2010, 2019.
* Decennial U.S. Census data for each decade and 2019 ACS. 
  * Race
  * Housing density
  * Housing tenure - % renter-occupied housing units
  * Housing occupancy - % vacant housing units
* Other Environmental Data:
  * Precipitation - mm in previous 6 months
  * Eco-regions - 12 classes. Classes are based on the regional ecosystem - soil, landform, major vegetation types, climate.
  * Impervious space - % land cover impervious

```{r data imports, warning=F, message=F}
## DATA IMPORTS-----------------------------------------------------------------
full_cbg_data_l <- readRDS("all_vars_on_cbg10_LONG.rds")
full_cbg_data_chg_w <- readRDS("chg_vars_on_cbg10_WIDE.rds")


## PREP EVI DATA
full_cbg_data_l_evi_focus <- full_cbg_data_l %>% 
  filter(evi >= 0)

```


## Metro-Micro Area Comparison

Greenness (EVI/NDVI) may take on different meanings/levels in metropolitan and micropolitan areas. Previous work has suggested that they shouldn't be studied together. If we look at the distributions of EVI at the place level across metro and micropolitan areas, there do not appear to be striking differences. 

* The middle of the distribution tends to hover around 0.45 and increases over time. 
* Noticeably, the distribution for metro areas tends to spread more over time, especially as of 2020. Although the variation in greenness in micropolitan areas is appreciable in 2020 compared to previous decades, there is more clumping of values across the distribution in micropolitan areas, whereas the the tails of the metro area distribution are not as dense.

```{r violin evi metro and micro, message=F, warning=F}
place_evi_by_year <- full_cbg_data_l %>%
  ungroup() %>%
  group_by(PLACEFP_DERIVED, PLACE_DERIVED_NM, wave, in_metro) %>% 
  filter(!is.na(in_metro)) %>%
  summarize(evi_mean = mean(evi, na.rm = T),
            evi_min = min(evi, na.rm = T),
            evi_max = max(evi, na.rm = T),
            n = n())

metro_micro_place_evi_comparison <- place_evi_by_year %>%
  ungroup() %>%
  filter(!is.na(evi_mean)) %>%
  group_by(in_metro, wave) %>%
  summarize_at(vars(evi_mean), list(mean = mean, sd = sd, median = median, min = min, max = max))

violin_evi_by_metro_status <- place_evi_by_year %>%
  ggplot(aes(fill = as.factor(wave), y = evi_mean, x = in_metro)) + 
    geom_violin(position = "dodge", alpha = 0.5, outlier.colour = "transparent") +
    scale_fill_viridis(discrete = T, name = "") +
    theme_minimal()  +
    theme(legend.position="right") +
    labs(title = "Place-level EVI by Metro & Micro Area",
         y = "Average EVI",
         x = "")
 
violin_evi_by_metro_status
  
```

The summary statistics indicate noticeable increases from 2010 to 2020 for metro and micro areas at both the center and tails of the distributions. 

```{r evi metro and micro sum stats, message=F, warning=F}
metro_micro_place_evi_comparison <- place_evi_by_year %>%
  ungroup() %>%
  filter(!is.na(evi_mean)) %>%
  group_by(in_metro, wave) %>%
  summarize_at(vars(evi_mean), list(mean = mean, sd = sd, median = median, min = min, max = max))

knitr::kable(metro_micro_place_evi_comparison,
             digits = 3,
             caption = "Summary of Place-Level EVI by Metro-Micro Area",
             col.names = c("Region Type", "Year", "Mean", "SD", "Median", "Min", "Max")
             )
```


## Neighborhood Racial Composition & Greenness

The figure below illustrates the relationship between the average neighborhood racial composition and the EVI level for each decade within metropolitan areas.

The average neighborhood with an EVI of 0.5 was 83% white in 1990, 74% white in 2000, and 65% in 2010. The average share of Hispanics in a neighborhood of that same greenness level was 3% in 1990 and then 12% in 2010.  

If we look at a neighborhood with a low average EVI such as 0.1, the average share of whites in 1990 was 48% and then fell to 38% in 2010. 

Lastly, at the high end of the greenness scale around 0.8, the average share of whites in a neighborhood was 96% in 1990 and 93% in 2010. The average percent of Hispanics in a highly green neighborhood was less than 1% in 1990 and 3% in 2010, while the share of Black residents hovered around 2% across this time period. 

Overall, as the EVI increases, the average percent of whites also increased. However, as of 2010, minorities - particularly, Hispanics and Blacks, make up larger average shares of greener neighborhoods. This pattern appears across each of the years.


```{r evi and racial comp at cbg level, message=F, warning=F}
## CBG EVI & AVG RACIAL COMPOSITION METRO AREAS
average_cbg_composition_by_evi <- full_cbg_data_l %>% 
  ungroup() %>%
  filter(evi >= 0 & evi <= 1,
         !is.na(evi),
         in_metro == "Metropolitan Area") %>%
  mutate(wave = as.character(wave),
         evi = round(evi, 2)) %>%
  select(wave, evi, WHITE_NH_PCT:HISPANIC_PCT) %>%
  group_by(wave, evi) %>%
  summarise(across(where(is.numeric), ~ mean(.x, na.rm = TRUE))) %>%
  pivot_longer(
   cols = ends_with("_PCT"),
   names_to = "race",
   values_to = "percent"
 ) %>%
  mutate(race = str_remove(race, "_PCT"))


ggplot(average_cbg_composition_by_evi, aes(x=evi, y=percent, fill=race)) + 
   geom_area(alpha=0.6 , size=.5, colour="white") +
    facet_wrap(~ wave, nrow = 4) +
    scale_fill_brewer(palette = "Set2") +
    theme_minimal() + 
    labs(title = "Average Neighborhood Composition by Level of Greenness",
         y = "Average Neighborhood Composition (%)",
         x = "EVI",
         fill = "Race")

```

How do metro areas compare to micropolitan areas with respect to racial/ethnic composition and greenness? The multiple area graphs below show the breakdown by year and type of statistical region. Compared to the metropolitan areas, the neighborhoods in micropolitan areas did not experience the same increase in the average share of minority residents in neighborhoods with the highest levels of greenness. 

```{r test metro v micro areas evi and race comp, message=F, warning=F}
## CBG EVI & AVG RACIAL COMPOSITION METRO V. MICRO AREAS
average_cbg_composition_by_evi_all_cbsa <- full_cbg_data_l_evi_focus %>% 
  ungroup() %>%
  filter(evi <= 1,
         !is.na(evi),
         !is.na(in_metro)) %>%
  mutate(wave = as.character(wave),
         evi = round(evi, 2)) %>%
  select(wave, in_metro, evi, WHITE_NH_PCT:HISPANIC_PCT) %>%
  group_by(wave, evi, in_metro) %>%
  summarise(across(where(is.numeric), ~ mean(.x, na.rm = TRUE))) %>%
  pivot_longer(
   cols = ends_with("_PCT"),
   names_to = "race",
   values_to = "percent",
 ) %>%
  mutate(race = str_remove(race, "_PCT"))


ggplot(average_cbg_composition_by_evi_all_cbsa, aes(x=evi, y=percent, fill=race)) + 
   geom_area(alpha=0.6 , size=.5, colour="white") +
    facet_wrap(vars(wave, in_metro), nrow = 4) +
    scale_fill_brewer(palette = "Set2") +
    theme_minimal() + 
    labs(title = "Average Neighborhood Composition by Level of Greenness",
         y = "Average Neighborhood Composition (%)",
         x = "EVI",
         fill = "Race")


```


**1990 Anomaly - High EVI and high average % Black residents**

In the 1990 area chart above, the racial breakdown for high EVI values (i.e., above 0.8) changes drastically. According to the chart, the average neighborhood with an EVI of 0.84 is 95% Black. A deeper investigation into this data point reveals that only one census block group fell between 0.84 and 0.85 EVI score in 1990, and this particular census block had a racial composition that was 95% Black. 

This census block group is in Summerfield, Maryland, which is a suburb of Prince George's County. It is located in the Washington DC metro area. The census block group experienced demographic changes and decreasing greenness after 1990. Notably, the share of Black residents declined, while white and Hispanic residents increased. Despite these changes, Black residents were still the majority racial group in the neighborhood. EVI dropped to around 0.4/0.5 levels, which mirrors its neighboring census block group in Summerfield. 

```{r identify place cbgs with high evi and black 1990, message=F, warning=F}
## 1990 ANOMALY 
knitr::kable(full_cbg_data_l %>%
  filter(wave == 1990,
         evi > 0.84, 
         evi < 0.85) %>%
    ungroup() %>%
  select(PLACE_NM, cbsatitle, evi, WHITE_NH_PCT, BLACK_NH_PCT),
  digits = 2,
  caption = "1990 Anomaly",
  col.names = c("Place", "Metro", "EVI", " % White", "% Black")
)


knitr::kable(full_cbg_data_l_evi_focus %>%
  filter(PLACEFP_DERIVED == "2475810") %>%
  arrange(wave, evi) %>%
  mutate(wave = as.factor(wave)) %>%
  ungroup() %>%
  relocate(wave, .after = GEOID10) %>%
  select(GEOID10, wave, evi, BLACK_NH_PCT, POP_DENS),
  digits = 2,
  format.args = list(big.mark = ','),
  caption = "Summerfield, MD Census Block Group EVI and Demographics",
  col.names = c("CBG", "Year", "EVI", "% Black", "People/Sq. Mi.")
)

```


## Within Place Comparisons

```{r maps spatial visuals}
# place_evi_comp <- place_evi_by_year %>% 
#   mutate(evi_range = evi_max - evi_min) 
# 
# place_evi_comp %>%
#   ungroup() %>% 
#   filter(PLACEFP_DERIVED %in% top100_places$PLACEFP) %>%
#   arrange(evi_range) %>%
#   slice(1:50)

```

```{r northeast connected dotplots}
## BEESWARM PLOTS OF WITHIN PLACE DISTRIBUTIONS 
top100_places_within_time_comp <- full_cbg_data_l %>% 
  ungroup() %>% 
  filter(PLACEFP_DERIVED %in% top100_places$PLACEFP) %>% 
  mutate(state_fips = str_sub(PLACEFP_DERIVED, end = 2)) %>% 
  left_join(., state_id_regions, by = "state_fips") %>%
  mutate(PLACE_NM = str_remove(PLACE_DERIVED_NM, " city"),
         PLACE_ST_NM = str_c(PLACE_NM, state_ab, sep = ", "),
         year = as.character(wave),
         year_abb = str_sub(year, start = 3),
         year_abb = paste0("'", year_abb)
         )

# NORTHEAST
top100_places_within_time_comp %>% 
  filter(region == "Northeast") %>% 
  ggplot(., aes(x = year, y = evi, color = year)) + 
  geom_quasirandom() +
  theme_minimal() +
  scale_color_manual(values = wes_palette("Cavalcanti1")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace in Northeastern Cities",
    x = "Year",
    y = "EVI",
    caption = "Note: Each point represents a block group in the city") + 
  theme(legend.position = "none") +
  facet_wrap(~ PLACE_ST_NM) 

```

```{r imprev space analysis}
top100_places_within_time_comp %>% 
  filter(region == "Northeast") %>% 
  ggplot(., aes(x = year, y = perc_imperv, color = year)) + 
  geom_quasirandom() +
  theme_minimal() +
  scale_color_manual(values = wes_palette("Cavalcanti1")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace in Northeastern Cities",
    x = "Year",
    y = "EVI",
    caption = "Note: Each point represents a block group in the city") + 
  theme(legend.position = "none") +
  facet_wrap(~ PLACE_ST_NM) 


```


```{r midatlantic}
# MidAtlantic
top100_places_within_time_comp %>% 
  filter(region == "Mid-Atlantic") %>% 
  ggplot(., aes(x = year, y = evi, color = year)) + 
  geom_quasirandom() +
  theme_minimal() +
  scale_color_manual(values = wes_palette("Cavalcanti1")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace in Mid-Atlantic Cities",
    x = "Year",
    y = "EVI",
    caption = "Note: Each point represents a block group in the city") + 
  theme(legend.position = "none") +
  facet_wrap(~ PLACE_ST_NM)

```


```{r midwest}
# MIDWEST
top100_places_within_time_comp %>% 
  filter(region == "Midwest") %>% 
  mutate(PLACE_ST_NM = gsub("\\s*\\([^\\)]+\\)", "", PLACE_ST_NM)) %>%
  ggplot(., aes(x = year, y = evi, color = year)) + 
  geom_quasirandom() +
  theme_minimal() +
  scale_color_manual(values = wes_palette("Cavalcanti1")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace in Midwestern Cities",
    x = "Year",
    y = "EVI",
    caption = "Note: Each point represents a block group in the city") + 
  theme(legend.position = "none") +
  facet_wrap(~ PLACE_ST_NM)

```


```{r south}
# SOUTH
top100_places_within_time_comp %>% 
  filter(region == "South") %>% 
  mutate(PLACE_ST_NM = str_remove(PLACE_ST_NM, " urban county"), 
         PLACE_ST_NM = str_remove(PLACE_ST_NM, " government"),
         PLACE_ST_NM = str_remove(PLACE_ST_NM, " metropolitan"),
         PLACE_ST_NM = str_remove(PLACE_ST_NM, " metro"),
         PLACE_ST_NM = gsub("\\s*\\([^\\)]+\\)", "", PLACE_ST_NM)
         ) %>%
  ggplot(., aes(x = year, y = evi, color = year)) + 
  geom_quasirandom() +
  theme_minimal() +
  scale_color_manual(values = wes_palette("Cavalcanti1")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace in Southern Cities",
    x = "Year",
    y = "EVI",
    caption = "Note: Each point represents a block group in the city") + 
  theme(legend.position = "none") +
  facet_wrap(~ PLACE_ST_NM)

```

```{r Southwest}
# SOUTHWEST
top100_places_within_time_comp %>% 
  filter(region == "Southwest") %>% 
  ggplot(., aes(x = year, y = evi, color = year)) + 
  geom_quasirandom() +
  theme_minimal() +
  scale_color_manual(values = wes_palette("Cavalcanti1")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace in Southwestern Cities",
    x = "Year",
    y = "EVI",
    caption = "Note: Each point represents a block group in the city") + 
  theme(legend.position = "none") +
  facet_wrap(~ PLACE_ST_NM)

```


```{r mountainwest}
# MOUNTAINWEST
top100_places_within_time_comp %>% 
  filter(region == "West") %>% 
  ggplot(., aes(x = year, y = evi, color = year)) + 
  geom_quasirandom() +
  theme_minimal() +
  scale_color_manual(values = wes_palette("Cavalcanti1")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace in Mountain-West Cities",
    x = "Year",
    y = "EVI",
    caption = "Note: Each point represents a block group in the city") + 
  theme(legend.position = "none") +
  facet_wrap(~ PLACE_ST_NM)

```


```{r pacific west}
# PACIFIC WEST
top100_places_within_time_comp %>% 
  filter(region == "Pacific") %>% 
  mutate(PLACE_ST_NM = str_remove(PLACE_ST_NM, " municipality"),
         PLACE_ST_NM = str_remove(PLACE_ST_NM, " Urban")) %>%
  ggplot(., aes(x = year, y = evi, color = year)) + 
  geom_quasirandom() +
  theme_minimal() +
  scale_color_manual(values = wes_palette("Cavalcanti1")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace in Pacific-West Cities",
    x = "Year",
    y = "EVI",
    caption = "Note: Each point represents a block group in the city") + 
  theme(legend.position = "none") +
  facet_wrap(~ PLACE_ST_NM)

```

```{r greenspace evi area}
## PREP GREENSPACE DATA
evi_area_full_cbg <- evi_area_decades_cbg %>% 
  rename(green_area = "area") %>%
  left_join(., cbg_all_geoids, by = "GISJOIN") %>%
  select(GISJOIN:cbsatitle, top100_place, in_metro) %>%
  left_join(., cbg_area, by = "GEOID10") %>%
  mutate(perc_green = green_area/area, 
         perc_green = if_else(perc_green > 1, 
                              1,
                              perc_green)
  ) %>%
  left_join(., full_cbg_data, by = c("GISJOIN", "wave"))

top100_places_area_within_time_comp <- evi_area_full_cbg %>%
  filter(PLACEFP_DERIVED.x %in% top100_places$PLACEFP) %>% 
  mutate(state_fips = str_sub(PLACEFP_DERIVED.x, end = 2)) %>% 
  left_join(., state_id_regions, by = "state_fips") %>%
  mutate(PLACE_NM = str_remove(PLACE_DERIVED_NM.x, " city"),
         PLACE_ST_NM = str_c(PLACE_NM.x, state_ab, sep = ", "),
         year = as.character(wave),
         majority_minority = if_else(WHITE_NH_PCT <50,
                                     "Majority Nonwhite",
                                     "Majority white")
  )

```


```{r ne area boxplot}
# NORTHEAST
top100_places_area_within_time_comp %>% 
  filter(region == "Northeast",
         !is.na(majority_minority)) %>% 
ggplot(., aes(x=year, y=perc_green.x, fill=majority_minority)) + 
    geom_boxplot() +
    scale_fill_manual(values = wes_palette("Chevalier1")) + 
  guides(fill=guide_legend(title="Racial Composition")) +
  theme_minimal() + 
  labs(
    title = "Distribution of Greenspace Area in Northeastern Cities",
    x = "Year",
    y = "% Green Landcover") +
    facet_wrap(~PLACE_ST_NM)

```

```{r midatl area box}
# MIDATLANTIC
top100_places_area_within_time_comp %>% 
  filter(region == "Mid-Atlantic",
         !is.na(majority_minority)) %>% 
ggplot(., aes(x=year, y=perc_green, fill=majority_minority)) + 
    geom_boxplot() +
    scale_fill_manual(values = wes_palette("Chevalier1")) + 
  guides(fill=guide_legend(title="Racial Composition")) +
  theme_minimal() + 
  theme(legend.position = "none") +
  labs(
    title = "Distribution of Greenspace Area in Mid-Atlantic Cities",
    x = "Year",
    y = "% Green Landcover") +
    facet_wrap(~PLACE_ST_NM)

```

```{r midwest area box}
# MIDWEST
top100_places_area_within_time_comp %>% 
  filter(region == "Midwest",
         !is.na(majority_minority)) %>% 
  mutate(PLACE_ST_NM = gsub("\\s*\\([^\\)]+\\)", "", PLACE_ST_NM)) %>%
ggplot(., aes(x=year, y=perc_green, fill=majority_minority)) + 
    geom_boxplot() +
    scale_fill_manual(values = wes_palette("Chevalier1")) + 
  guides(fill=guide_legend(title="Racial Composition")) +
  theme_minimal() + 
  theme(legend.position = "none") +
  labs(
    title = "Distribution of Greenspace Area in Midwestern Cities",
    x = "Year",
    y = "% Green Landcover") +
    facet_wrap(~PLACE_ST_NM)

```

```{r south area boxes}
# SOUTH
top100_places_area_within_time_comp %>% 
  filter(region == "South",
         !is.na(majority_minority)) %>% 
  mutate(PLACE_ST_NM = str_remove(PLACE_ST_NM, " urban county"), 
         PLACE_ST_NM = str_remove(PLACE_ST_NM, " government"),
         PLACE_ST_NM = str_remove(PLACE_ST_NM, " metropolitan"),
         PLACE_ST_NM = str_remove(PLACE_ST_NM, " metro"),
         PLACE_ST_NM = str_remove(PLACE_ST_NM, " County"),
         PLACE_ST_NM = gsub("\\s*\\([^\\)]+\\)", "", PLACE_ST_NM)
         ) %>%
ggplot(., aes(x=year, y=perc_green, fill=majority_minority)) + 
    geom_boxplot() +
    scale_fill_manual(values = wes_palette("Chevalier1")) + 
  guides(fill=guide_legend(title="Racial Composition")) +
  theme_minimal() + 
  theme(legend.position = "none") + 
  labs(
    title = "Distribution of Greenspace Area in Southern Cities",
    x = "Year",
    y = "% Green Landcover") +
    facet_wrap(~PLACE_ST_NM)

```

```{r sw area box}
# SOUTHWESTERN
top100_places_area_within_time_comp %>% 
  filter(region == "Southwest",
         !is.na(majority_minority)) %>% 
ggplot(., aes(x=year, y=perc_green, fill=majority_minority)) + 
    geom_boxplot() +
    scale_fill_manual(values = wes_palette("Chevalier1")) + 
  guides(fill=guide_legend(title="Racial Composition")) +
  theme_minimal() + 
  theme(legend.position = "none") + 
  labs(
    title = "Distribution of Greenspace Area in Southwestern Cities",
    x = "Year",
    y = "% Green Landcover") +
    facet_wrap(~PLACE_ST_NM)

```

```{r mw area box}
# MOUNTAIN WEST
top100_places_area_within_time_comp %>% 
  filter(region == "West",
         !is.na(majority_minority)) %>% 
ggplot(., aes(x=year, y=perc_green, fill=majority_minority)) + 
    geom_boxplot() +
    scale_fill_manual(values = wes_palette("Chevalier1")) + 
  guides(fill=guide_legend(title="Racial Composition")) +
  theme_minimal() + 
  theme(legend.position = "none") + 
  labs(
    title = "Distribution of Greenspace Area in Mountain-west Cities",
    x = "Year",
    y = "% Green Landcover") +
    facet_wrap(~PLACE_ST_NM)

```

```{r pacific area box}
# PACIFIC WEST
top100_places_area_within_time_comp %>% 
  filter(region == "Pacific",
         !is.na(majority_minority)) %>% 
  mutate(PLACE_ST_NM = str_remove(PLACE_ST_NM, " municipality"),
         PLACE_ST_NM = str_remove(PLACE_ST_NM, "Urban ")) %>%
ggplot(., aes(x=year, y=perc_green, fill=majority_minority)) + 
    geom_boxplot() +
    scale_fill_manual(values = wes_palette("Chevalier1")) + 
  guides(fill=guide_legend(title="Racial Composition")) +
  theme_minimal() + 
  theme(legend.position = "none") + 
  labs(
    title = "Distribution of Greenspace Area in Pacific-west Cities",
    x = "Year",
    y = "% Green Landcover") +
    facet_wrap(~PLACE_ST_NM)

```




## Line Graphs

### Impervious space v. green-space
```{r imp v green space lines}
## Southwest imperv and green-space
full_cbg_data %>%
ungroup() %>% 
  filter(PLACEFP_DERIVED %in% top100_places$PLACEFP) %>% 
  mutate(state_fips = str_sub(PLACEFP_DERIVED, end = 2)) %>% 
  left_join(., state_id_regions, by = "state_fips") %>% 
  mutate(PLACE_NM = str_remove(PLACE_DERIVED_NM, " city"),
         PLACE_ST_NM = str_c(PLACE_NM, state_ab, sep = ", ")
  ) %>%  
  group_by(PLACE_ST_NM, region, wave) %>%
  summarize(impervious = sum(imperv_area),
            green = sum(green_area),
            area = sum(area)
            ) %>% 
  ungroup() %>% 
  mutate(perc_imperv = (impervious/area) * 100,
         perc_green = (green/area) * 100) %>%
  pivot_longer(
               cols = starts_with("perc_"),
               names_to = "landcover",
               values_to = "percent"
               ) %>%
  filter(region == "Mid-Atlantic") %>% 
  ggplot(aes(wave, percent, group = landcover)) +
  geom_line(aes(color=landcover)) + 
  geom_point(aes(color=landcover)) +
  scale_color_manual(values = wes_palette("Chevalier1")
                     ) +
  theme_minimal() + 
  facet_wrap(~ PLACE_ST_NM)

```

### Impervious space by race comp
```{r line imperv and green area}
## Southwest percent imperv
full_cbg_data %>%
ungroup() %>% 
  filter(PLACEFP_DERIVED %in% top100_places$PLACEFP) %>% 
  mutate(state_fips = str_sub(PLACEFP_DERIVED, end = 2)) %>% 
  left_join(., state_id_regions, by = "state_fips") %>% 
  mutate(PLACE_NM = str_remove(PLACE_DERIVED_NM, " city"),
         PLACE_ST_NM = str_c(PLACE_NM, state_ab, sep = ", "),
         majority_minority = if_else(WHITE_NH_PCT < 50,
                                     "Majority Nonwhite",
                                     "Majority white"),
         race_comp = case_when(
           majority_minority == "Majority white" ~ "Majority white",
           BLACK_NH_PCT > 50  ~ "Majority Black",
           ASIAN_PACIFIC_NH_PCT > 50 ~ "Majority Asian",
           HISPANIC_PCT > 50 ~ "Majority Latinx",
           TRUE  ~ "Other"
         )
  ) %>%  
  group_by(PLACE_ST_NM, region, race_comp, wave) %>%
  summarize(impervious = sum(imperv_area),
            green = sum(green_area),
            area = sum(area)
            ) %>% 
  ungroup() %>% 
  mutate(perc_imperv = (impervious/area) * 100,
         perc_green = (green/area) * 100) %>%
  pivot_longer(
               cols = starts_with("perc_"),
               names_to = "landcover",
               values_to = "percent"
               ) %>%
  filter(region == "Pacific" & landcover == "perc_imperv") %>% 
  ggplot(aes(wave, percent, group = race_comp)) +
  geom_line(aes(color=race_comp)) + 
  geom_point(aes(color=race_comp)) +
  scale_color_manual(values = wes_palette("Darjeeling1")) +
  theme_minimal() + 
  facet_wrap(~ PLACE_ST_NM)

```

```{r map green space per pixel}
bexar_all_data <- full_cbg_data %>% 
  ungroup() %>% 
  #filter(PLACEFP_DERIVED %in% top100_places$PLACEFP) %>% 
  mutate(state_fips = str_sub(PLACEFP_DERIVED, end = 2)) %>% 
  left_join(., state_id_regions, by = "state_fips") %>% 
  mutate(PLACE_NM = str_remove(PLACE_DERIVED_NM, " city"),
         PLACE_ST_NM = str_c(PLACE_NM, state_ab, sep = ", "),
         STATECO = str_sub(GEOID10, end = 5)
  ) %>%
  filter(STATECO == "48029")

bexar_shp <- cbg10 %>%
  select(GEOID10, STATEFP10, COUNTYFP10) %>% 
  right_join(bexar_all_data, by ="GEOID10") 

# %>% 
#     ggplot() + 
#     geom_sf(aes(fill=perc_green)) + 
#   facet_wrap(~ wave)


tm_shape(bexar_shp) +
  tm_polygons("perc_green", 
              style="quantile", 
              palette = "Greens",
              title="Bexar County") + 
      tm_facets(by = "wave")

tmap_mode("view")
tmap_last()

```


```{r bexar gray space}
## Impervious space
imperv_bexar <- tm_shape(bexar_shp) +
  tm_polygons("perc_imperv", 
              style="quantile", 
              palette="Greys",
              title="Bexar County Impervious Space",
              border.alpha = 0.2) + 
      tm_facets(by = "wave")
tmap_save(imperv_bexar, "Bexar_impervious.png", width=1920, height=1080, asp=0)

```



```{r map green space per pixel comal co}
comal_all_data <- full_cbg_data %>% 
  ungroup() %>% 
  #filter(PLACEFP_DERIVED %in% top100_places$PLACEFP) %>% 
  mutate(state_fips = str_sub(PLACEFP_DERIVED, end = 2)) %>% 
  left_join(., state_id_regions, by = "state_fips") %>% 
  mutate(PLACE_NM = str_remove(PLACE_DERIVED_NM, " city"),
         PLACE_ST_NM = str_c(PLACE_NM, state_ab, sep = ", "),
         STATECO = str_sub(GEOID10, end = 5)
  ) %>%
  filter(STATECO == "48091")

comal_shp <- cbg10 %>%
  select(GEOID10, STATEFP10, COUNTYFP10) %>% 
  right_join(comal_all_data, by ="GEOID10") 

tm_shape(comal_shp) +
  tm_polygons("perc_imperv", 
              style="quantile", 
              palette = "Greys",
              title="Comal County",
              border.alpha = 0.2) + 
      tm_facets(by = "wave")

tmap_mode("view")
tmap_last()

# SATX-METRO
satx_metro <- full_cbg_data %>% 
  ungroup() %>% 
  mutate(state_fips = str_sub(PLACEFP_DERIVED, end = 2)) %>% 
  left_join(., state_id_regions, by = "state_fips") %>%  
  mutate(PLACE_NM = str_remove(PLACE_DERIVED_NM, " city"),
         PLACE_ST_NM = str_c(PLACE_NM, state_ab, sep = ", "),
         STATECO = str_sub(GEOID10, end = 5)
  ) %>% 
  left_join(., cbsa_co_cw_reformat, by = c("STATECO" = "state_co")) %>%
  filter(cbsacode == "41700")


satx_metro_shp <- cbg10 %>%
  select(GEOID10, STATEFP10, COUNTYFP10) %>% 
  right_join(satx_metro, by ="GEOID10") 

tm_shape(satx_metro_shp) +
  tm_polygons("perc_imperv", 
              style="quantile", 
              palette = "Greys",
              title="San Antonio, TX Metro",
              border.alpha = 0.2) + 
      tm_facets(by = "wave")

```



```{r philly imperv}

philly_metro <- full_cbg_data %>% 
  ungroup() %>% 
  mutate(state_fips = str_sub(PLACEFP_DERIVED, end = 2)) %>% 
  left_join(., state_id_regions, by = "state_fips") %>%  
  mutate(PLACE_NM = str_remove(PLACE_DERIVED_NM, " city"),
         PLACE_ST_NM = str_c(PLACE_NM, state_ab, sep = ", "),
         STATECO = str_sub(GEOID10, end = 5)
  ) %>% 
  left_join(., cbsa_co_cw_reformat, by = c("STATECO" = "state_co")) %>%
  filter(metrodivisioncode == "37964")

philly_metro_shp <- cbg10 %>%
  select(GEOID10, STATEFP10, COUNTYFP10) %>% 
  right_join(philly_metro, by ="GEOID10") 

tm_shape(philly_metro_shp) +
  tm_polygons("perc_imperv", 
              style="quantile", 
              palette = "Greys",
              title="Philadelphia Metro",
              border.alpha = 0.2) + 
      tm_facets(by = "wave")

```

## Multivariate Analysis

* The goal is to get at changes in evi and race composition within neighborhoods and places, so run time fixed effects models.
  * Non-spatial model and spatial error model. The spatial error model would account for residual spatial autocorrelation.
* Focal predictor is % race group.
* Adjust for environmental variables and demographic factors listed above.