################################################################################ ################################################################################ ################################################################################ #Project: Manufactured Housing ################################################################################ ################################################################################ ################################################################################ ################################################################################ #set working directory ################################################################################ #INSERT DIRECTORY ################################################################################ #load libraries ################################################################################ library(dplyr) library(tidyr) library(scales) library(fixest) library(broom) library(purrr) library(sf) library(ggplot2) library(viridis) library(survey) library(reldist) library(data.table) library(kableExtra) library(gtsummary) library(patchwork) library(cowplot) library(ggpubr) library(png) options(scipen=999) ################################################################################ #buildings data ################################################################################ files <- c( "Manufactured_Homes_Data_Buildings_Final_part1.RData", "Manufactured_Homes_Data_Buildings_Final_part2.RData", "Manufactured_Homes_Data_Buildings_Final_part3.RData" ) # Helper to identify tabular objects is_tabular <- function(x) inherits(x, c("data.frame", "data.table", "tbl_df")) # Load a tabular chunk from a .RData file (chooses the largest tabular object if multiple) load_tabular_from_rdata <- function(path) { e <- new.env(parent = emptyenv()) nms <- load(path, envir = e) if (!length(nms)) stop(sprintf("No objects found in '%s'.", path)) tabs <- nms[vapply(nms, function(nm) is_tabular(get(nm, envir = e)), logical(1))] if (!length(tabs)) { stop(sprintf("No tabular (data.frame/data.table/tibble) objects found in '%s'.", path)) } # Pick the largest by size (fallback to row count if size ties) sizes <- sapply(tabs, function(nm) object.size(get(nm, envir = e))) target <- tabs[which.max(as.numeric(sizes))] get(target, envir = e) } # Load all three parts parts <- lapply(files, load_tabular_from_rdata) # Align columns (union of all columns) before row-binding all_cols <- Reduce(union, lapply(parts, names)) parts_aligned <- lapply(parts, function(df) { missing <- setdiff(all_cols, names(df)) if (length(missing)) { for (m in missing) df[[m]] <- NA } # Reorder to the union's order df[all_cols] }) # Combine d <- do.call(rbind, parts_aligned) #summarize buildings #select 1 row per county to sum FOOTPRINT_COUNT, MH_COUNT, and NON_MH_COUNT d_summary_buildings <- d %>% group_by(GISJOIN) %>% slice(1) #124,838,836 buildings sum(d_summary_buildings$FOOTPRINT_COUNT, na.rm = TRUE) #121,502,451 (97%) other buildings sum(d_summary_buildings$NON_MH_COUNT, na.rm = TRUE) sum(d_summary_buildings$NON_MH_COUNT, na.rm = TRUE) / sum(d_summary_buildings$FOOTPRINT_COUNT, na.rm = TRUE) #3,336,385 (3%) manufactured housing unit sum(d_summary_buildings$MH_COUNT, na.rm = TRUE) sum(d_summary_buildings$MH_COUNT, na.rm = TRUE) / sum(d_summary_buildings$FOOTPRINT_COUNT, na.rm = TRUE) #summarize buildings by census region d_summary_buildings_region <- d_summary_buildings %>% group_by(CENSUS_REGION) %>% summarise(FOOTPRINT_COUNT = sum(FOOTPRINT_COUNT), MH_COUNT = sum(MH_COUNT), NON_MH_COUNT = sum(NON_MH_COUNT), MH_PCT = sum(MH_COUNT) / sum(FOOTPRINT_COUNT), NON_MH_PCT = sum(NON_MH_COUNT) / sum(FOOTPRINT_COUNT)) #summarize MHs #351,782 (28%) MHs that are located in a MH community/park sum(ifelse(d$BUILDING_TYPE == "MH - Park", 1, 0)) sum(ifelse(d$BUILDING_TYPE == "MH - Park", 1, 0)) / sum(d$MH) #924,022 (72%) MHs that are not located in a MH community/park sum(ifelse(d$BUILDING_TYPE == "MH - Not Park", 1, 0)) sum(ifelse(d$BUILDING_TYPE == "MH - Not Park", 1, 0)) / sum(d$MH) #summarize MH by census region d_summary_mh <- d %>% group_by(CENSUS_REGION) %>% summarise(MH_COUNT = sum(MH), MH_PARK_COUNT = sum(ifelse(BUILDING_TYPE == "MH - Park", 1, 0)), MH_NOT_PARK_COUNT = sum(ifelse(BUILDING_TYPE == "MH - Not Park", 1, 0)), MH_PARK_PCT = MH_PARK_COUNT / MH_COUNT, MH_NOT_PARK_PCT = MH_NOT_PARK_COUNT / MH_COUNT) ################################################################################ #summary statistics ################################################################################ #function to calculate weighted mean weighted_mean <- function(x, w) { weighted.mean(x, w, na.rm = TRUE) } #function to calculate weighted standard deviation weighted_sd <- function(x, w) { mean_w <- weighted_mean(x, w) sqrt(sum(w * (x - mean_w)^2, na.rm = TRUE) / sum(w, na.rm = TRUE)) } ################################################################################ #code creates summary statistics table #variables reported in summary statistics vars_summary <- c("POLLUTION_PM25", "EPA_DISTANCE", "FLOOD_100_500_YR", "INDUSTRY_JOBS_500", "LAND_COVER_GREEN_500", "LAND_COVER_TREES_500", "LAND_COVER_WATER_500", "SCHOOLS_DISTANCE", "SNAP_DISTANCE", "STREET_INTERSECTIONS_500", "TRANSIT_STOP_500", "POPULATION_DENSITY", "RACE_WHITE_PCT", "IMMIGRANT_PCT", "POVERTY_PCT", "UNEMPLOYMENT_PCT", "SINGLE_FAMILY_HOUSING_PCT", "HOMEOWNERSHIP_PCT", "VACANCY_PCT", "HOUSEHOLD_INCOME_MEDIAN", "HOUSING_VALUE_MEDIAN", "YEAR_STRUCTURE_MEDIAN") #function calculates summary statistics for full sample fun_calculate_weighted_stats_full <- function(data, var_names) { #loop over each variable to create individual dataframes stats <- lapply(var_names, function(v) { mean_weighted <- round(weighted_mean(data[[v]], data[["SELECTED_WEIGHT"]]), 4) sd_weighted <- round(weighted_sd(data[[v]], data[["SELECTED_WEIGHT"]]), 4) return(data.frame(variable = v, weighted_mean = mean_weighted, weighted_sd = sd_weighted)) }) #combine variable dataframes into single data frame stats_df <- do.call(rbind, stats) #convert wide dataframe to long stats_df_long <- pivot_longer(stats_df, cols = c("weighted_mean", "weighted_sd"), names_to = "stat_type", values_to = "value") #return dataframe return(stats_df_long) } #function calculates summary statistics for individual building types fun_calculate_weighted_stats_group <- function(data, group_name, var_names) { #loop over each variable to create individual dataframes #calculate statistics within building type stats <- lapply(var_names, function(v) { mean_weighted <- round(weighted_mean(data[data$BUILDING_TYPE == group_name, v], data[data$BUILDING_TYPE == group_name, "SELECTED_WEIGHT"]), 4) sd_weighted <- round(weighted_sd(data[data$BUILDING_TYPE == group_name, v], data[data$BUILDING_TYPE == group_name, "SELECTED_WEIGHT"]), 4) return(data.frame(variable = v, weighted_mean = mean_weighted, weighted_sd = sd_weighted)) }) #combine variable dataframes into single data frame stats_df <- do.call(rbind, stats) #convert wide dataframe to long stats_df_long <- pivot_longer(stats_df, cols = c("weighted_mean", "weighted_sd"), names_to = "stat_type", values_to = "value") #return dataframe return(stats_df_long) } #calculate summary statistics for full sample stats_full <- fun_calculate_weighted_stats_full(d, vars_summary) %>% dplyr::rename(All = value) #calculate summary statistics for individual building types stats_other_buildings <- fun_calculate_weighted_stats_group(d, "Other", vars_summary) %>% dplyr::rename("Other" = value) stats_mh_park <- fun_calculate_weighted_stats_group(d, "MH - Park", vars_summary) %>% dplyr::rename("MH (Community)" = value) stats_mh_not_park <- fun_calculate_weighted_stats_group(d, "MH - Not Park", vars_summary) %>% dplyr::rename("MH (Not Community)" = value) #combine summary statistics into table stats_summary <- cbind(dplyr::select(stats_full, - stat_type), dplyr::select(stats_other_buildings, "Other"), dplyr::select(stats_mh_park, "MH (Community)"), dplyr::select(stats_mh_not_park, "MH (Not Community)")) #hard code observations for summary statistics table stats_summary <- rbind(stats_summary, c("Observations", nrow(d), nrow(filter(d, BUILDING_TYPE == "Other")), nrow(filter(d, BUILDING_TYPE == "MH - Park")), nrow(filter(d, BUILDING_TYPE == "MH - Not Park"))), c("(%)", "", nrow(filter(d, BUILDING_TYPE == "Other")) / nrow(d), nrow(filter(d, BUILDING_TYPE == "MH - Park")) / nrow(d), nrow(filter(d, BUILDING_TYPE == "MH - Not Park")) / nrow(d))) #latex code that creates summary statistics table kable(stats_summary, format = "latex", booktabs = FALSE, linesep = "") ################################################################################ #regressions ################################################################################ #code estimates relative effects #function to estimate relative effects for plots fun_reg_relative_effects <- function(data, dep_vars, predictors, selected_predictor, weights) { #create empty list to store the results for each dependent variable results_list <- list() #loop over each dependent variable for (dep_var in dep_vars) { #create formula for regression with county fixed effects formula <- as.formula(paste(dep_var, " ~ ", paste(predictors, collapse = " + "), " | GISJOIN")) #create subset of data that is used in regression (i.e., remove NAs) relevant_cols <- c(dep_var, predictors, weights, "GISJOIN") data_clean <- na.omit(data[relevant_cols]) #estimate regression model <- feols(formula, data = data_clean, weights = data_clean[[weights]]) #keep coefficients and standard errors from regression coef_vals <- coef(model)[predictors] se_vals <- se(model)[predictors] #extract values for selected predictor (e.g., MH (Community) and MH (Not Community)) coef_val <- coef_vals[selected_predictor] se_val <- se_vals[selected_predictor] #calculate 95% CI for coefficient lower_ci_coef <- coef_val - (1.96 * se_val) upper_ci_coef <- coef_val + (1.96 * se_val) #calculate weighted mean of dependent variable weighted_mean_dep_var <- weighted.mean(data_clean[[dep_var]], w = data_clean[[weights]]) #calculate relative effect of dependent variable relative_effect <- (coef_val / weighted_mean_dep_var) * 100 relative_effect_se <- (se_val / weighted_mean_dep_var) * 100 #calculate 95% CI for relative effect lower_ci_relative <- relative_effect - (1.96 * relative_effect_se) upper_ci_relative <- relative_effect + (1.96 * relative_effect_se) #store results for selected predictor (e.g., MH (Community) and MH (Not Community)) dep_results <- data.table(DEPENDENT_VARIABLE = dep_var, MH_TYPE = selected_predictor, COEFFICIENT = coef_val, COEFFICIENT_SE = se_val, COEFFICIENT_CI_LOWER = lower_ci_coef, COEFFICIENT_CI_UPPER = upper_ci_coef, WEIGHTED_MEAN = weighted_mean_dep_var, RELATIVE_EFFECT = relative_effect, RELATIVE_EFFECT_SE = relative_effect_se, RELATIVE_EFFECT_CI_LOWER = lower_ci_relative, RELATIVE_EFFECT_CI_UPPER = upper_ci_relative) #append results results_list[[dep_var]] <- dep_results } #combine results final_results <- rbindlist(results_list) return(final_results) } #variables for above function vars_dep <- c("POLLUTION_PM25", "EPA_DISTANCE", "FLOOD_100_500_YR", "INDUSTRY_JOBS_500", "LAND_COVER_GREEN_500", "LAND_COVER_TREES_500", "LAND_COVER_WATER_500", "SCHOOLS_DISTANCE", "SNAP_DISTANCE", "STREET_INTERSECTIONS_500", "TRANSIT_STOP_500") vars_controls <- c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK", "POPULATION_DENSITY", "RACE_WHITE_PCT", "IMMIGRANT_PCT", "POVERTY_PCT", "UNEMPLOYMENT_PCT", "SINGLE_FAMILY_HOUSING_PCT", "HOMEOWNERSHIP_PCT", "VACANCY_PCT", "HOUSEHOLD_INCOME_MEDIAN", "HOUSING_VALUE_MEDIAN", "YEAR_STRUCTURE_MEDIAN") vars_mh <- c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK") #function to rename and order variables fun_rename_dependent_variables <- function(data) { data %>% mutate(DEPENDENT_VARIABLE = case_when(DEPENDENT_VARIABLE == "POLLUTION_PM25" ~ "Air Pollution", DEPENDENT_VARIABLE == "EPA_DISTANCE" ~ "Facilities Regulated by EPA", DEPENDENT_VARIABLE == "FLOOD_100_500_YR" ~ "Floodplains", DEPENDENT_VARIABLE == "INDUSTRY_JOBS_500" ~ "Industry Jobs", DEPENDENT_VARIABLE == "LAND_COVER_GREEN_500" ~ "Land Cover - Green Space", DEPENDENT_VARIABLE == "LAND_COVER_TREES_500" ~ "Land Cover - Trees", DEPENDENT_VARIABLE == "LAND_COVER_WATER_500" ~ "Land Cover - Water", DEPENDENT_VARIABLE == "SCHOOLS_DISTANCE" ~ "Schools", DEPENDENT_VARIABLE == "SNAP_DISTANCE" ~ "SNAP Retailers", DEPENDENT_VARIABLE == "STREET_INTERSECTIONS_500" ~ "Street Intersections", DEPENDENT_VARIABLE == "TRANSIT_STOP_500" ~ "Transit Stops", TRUE ~ DEPENDENT_VARIABLE), DEPENDENT_VARIABLE = factor(DEPENDENT_VARIABLE, levels = rev(c("Air Pollution", "Facilities Regulated by EPA", "Floodplains", "Industry Jobs", "Land Cover - Green Space", "Land Cover - Trees", "Land Cover - Water", "Schools", "SNAP Retailers", "Street Intersections", "Transit Stops")))) } #regression results #full sample relative_effects <- fun_reg_relative_effects(d, vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() #census region relative_effects_midwest <- fun_reg_relative_effects(filter(d, CENSUS_REGION == "Midwest"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() relative_effects_northeast <- fun_reg_relative_effects(filter(d, CENSUS_REGION == "Northeast"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() relative_effects_south <- fun_reg_relative_effects(filter(d, CENSUS_REGION == "South"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() relative_effects_west <- fun_reg_relative_effects(filter(d, CENSUS_REGION == "West"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() #city governance / population relative_effects_1K <- fun_reg_relative_effects(filter(d, CITY_POPULATION_BIN == "<1K"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() relative_effects_1K_10K <- fun_reg_relative_effects(filter(d, CITY_POPULATION_BIN == "1K-10K"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() relative_effects_10K_50K <- fun_reg_relative_effects(filter(d, CITY_POPULATION_BIN == "10K-50K"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() relative_effects_50K_100K <- fun_reg_relative_effects(filter(d, CITY_POPULATION_BIN == "50K-100K"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() relative_effects_100K <- fun_reg_relative_effects(filter(d, CITY_POPULATION_BIN == ">100K"), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() relative_effects_unincorporated <- fun_reg_relative_effects(filter(d, is.na(CITY_POPULATION_BIN)), vars_dep, vars_controls, vars_mh, "SELECTED_WEIGHT") %>% fun_rename_dependent_variables() ################################################################################ #code plots relative effects in individual figures #function to make plots fun_plot_relative_effects_individual <- function(data, title, y_breaks, y_limits) { ggplot(data, aes(x = DEPENDENT_VARIABLE, y = RELATIVE_EFFECT, group = MH_TYPE, color = MH_TYPE, fill = MH_TYPE, shape = MH_TYPE)) + geom_errorbar(aes(ymin = RELATIVE_EFFECT_CI_LOWER, ymax = RELATIVE_EFFECT_CI_UPPER), width = 0, size = 1) + geom_point(size = 4) + coord_flip() + xlab("Variable") + ylab("Relative Effect") + ggtitle(title) + geom_hline(yintercept = 0, linetype = 2, size = 0.5) + scale_y_continuous(breaks = y_breaks, limits = y_limits) + scale_shape_manual(name = "Building Type:", labels = c("MH (Community)", "MH (Not Community)"), values = c(16, 17)) + scale_fill_manual(name = "Building Type:", labels = c("MH (Community)", "MH (Not Community)"), values = c("blue", "red")) + scale_color_manual(name = "Building Type:", labels = c("MH (Community)", "MH (Not Community)"), values = c("blue", "red")) + theme_classic() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14), axis.title.x = element_text(size = 16, margin = margin(t = 0, r = 0, b = 20, l = 0), vjust = -5), axis.title.y = element_text(size = 16, vjust = 2), legend.text = element_text(size = 12), legend.title = element_text(size = 14), legend.position = "bottom", plot.title = element_text(size = 20, hjust = 0.5)) } #national plot f_estimates_relative <- fun_plot_relative_effects_individual( data = relative_effects, title = "Exposure to Amenities/Risks", y_breaks = seq(-50, 50, by = 25), y_limits = c(-50, 50) ) #census region plots f_estimates_relative_midwest <- fun_plot_relative_effects_individual( data = relative_effects_midwest, title = "Exposure to Amenities/Risks (Midwest)", y_breaks = seq(-50, 125, by = 25), y_limits = c(-58, 131) ) f_estimates_relative_northeast <- fun_plot_relative_effects_individual( data = relative_effects_northeast, title = "Exposure to Amenities/Risks (Northeast)", y_breaks = seq(-50, 100, by = 25), y_limits = c(-64, 103) ) f_estimates_relative_south <- fun_plot_relative_effects_individual( data = relative_effects_south, title = "Exposure to Amenities/Risks (South)", y_breaks = seq(-50, 50, by = 25), y_limits = c(-50, 50) ) f_estimates_relative_west <- fun_plot_relative_effects_individual( data = relative_effects_west, title = "Exposure to Amenities/Risks (West)", y_breaks = seq(-100, 50, by = 25), y_limits = c(-100, 50) ) #city governance/population plots f_estimates_relative_1K <- fun_plot_relative_effects_individual( data = relative_effects_1K, title = "Exposure to Amenities/Risks (<1K)", y_breaks = seq(-50, 50, by = 25), y_limits = c(-50, 50) ) f_estimates_relative_1K_10K <- fun_plot_relative_effects_individual( data = relative_effects_1K_10K, title = "Exposure to Amenities/Risks (1K-10K)", y_breaks = seq(-50, 50, by = 25), y_limits = c(-50, 50) ) f_estimates_relative_10K_50K <- fun_plot_relative_effects_individual( data = relative_effects_10K_50K, title = "Exposure to Amenities/Risks (10K-50K)", y_breaks = seq(-50, 50, by = 25), y_limits = c(-50, 50) ) f_estimates_relative_50K_100K <- fun_plot_relative_effects_individual( data = relative_effects_50K_100K, title = "Exposure to Amenities/Risks (50K-100K)", y_breaks = seq(-25, 75, by = 25), y_limits = c(-25, 75) ) f_estimates_relative_100K <- fun_plot_relative_effects_individual( data = relative_effects_100K, title = "Exposure to Amenities/Risks (>100K)", y_breaks = seq(-25, 100, by = 25), y_limits = c(-25, 100) ) f_estimates_relative_unincorporated <- fun_plot_relative_effects_individual( data = relative_effects_unincorporated, title = "Exposure to Amenities/Risks (Unincorporated)", y_breaks = seq(-50, 50, by = 25), y_limits = c(-50, 50) ) #plots #full sample f_estimates_relative #census region f_estimates_relative_midwest f_estimates_relative_northeast f_estimates_relative_south f_estimates_relative_west #city governance/population f_estimates_relative_1K f_estimates_relative_1K_10K f_estimates_relative_10K_50K f_estimates_relative_50K_100K f_estimates_relative_100K f_estimates_relative_unincorporated ################################################################################ #code plots relative effects and combines into a single figure #function to make plots fun_plot_relative_effect_grid <- function(data, title, y_breaks, y_limits, blank_y_axis = FALSE, large_labels = FALSE, show_legend = FALSE) { p <- ggplot(data, aes(x = DEPENDENT_VARIABLE, y = RELATIVE_EFFECT, group = MH_TYPE, color = MH_TYPE, fill = MH_TYPE, shape = MH_TYPE)) + geom_errorbar(aes(ymin = RELATIVE_EFFECT_CI_LOWER, ymax = RELATIVE_EFFECT_CI_UPPER), width = 0, size = ifelse(show_legend, 1, 0.6)) + geom_point(size = ifelse(show_legend, 4, 2.3)) + coord_flip() + xlab("") + ylab("") + ggtitle(title) + geom_hline(yintercept = 0, linetype = 2, size = 0.5) + scale_y_continuous(breaks = y_breaks, limits = y_limits) + scale_shape_manual(name = "Building Type:", labels = c("MH (Community)", "MH (Not Community)"), values = c(16, 17)) + scale_fill_manual(name = "Building Type:", labels = c("MH (Community)", "MH (Not Community)"), values = c("blue", "red")) + scale_color_manual(name = "Building Type:", labels = c("MH (Community)", "MH (Not Community)"), values = c("blue", "red")) + theme_classic() + theme(axis.text.x = element_text(size = ifelse(large_labels, 14, 10)), axis.text.y = if (blank_y_axis) element_blank() else element_text(size = ifelse(large_labels, 14, 10)), axis.ticks.y = if (blank_y_axis) element_blank() else element_line(), axis.title.y = element_text(size = ifelse(large_labels, 16, 10), vjust = 2), axis.title.x = element_blank(), legend.position = if (show_legend) "bottom" else "none", legend.text = element_text(size = 12), legend.title = element_text(size = 14), plot.title = element_text(size = ifelse(large_labels, 20, 16), hjust = 0.5)) return(p) } #full sample f_estimates_relative_grid <- fun_plot_relative_effect_grid( relative_effects, title = "United States", y_breaks = seq(-50, 50, by = 25), y_limits = c(-50, 50), large_labels = TRUE, show_legend = FALSE ) #census region region_inputs <- list( "Midwest" = list(data = relative_effects_midwest, title = "Midwest", y_breaks = seq(-50, 125, by = 25), y_limits = c(-58, 131), blank_y = FALSE), "Northeast" = list(data = relative_effects_northeast, title = "Northeast", y_breaks = seq(-50, 100, by = 25), y_limits = c(-64, 103), blank_y = TRUE), "South" = list(data = relative_effects_south, title = "South", y_breaks = seq(-50, 50, by = 25), y_limits = c(-50, 50), blank_y = FALSE), "West" = list(data = relative_effects_west, title = "West", y_breaks = seq(-100, 50, by = 25), y_limits = c(-100, 50), blank_y = TRUE) ) region_plots <- imap(region_inputs, function(region, name) { fun_plot_relative_effect_grid(region$data, region$title, region$y_breaks, region$y_limits, blank_y_axis = region$blank_y) }) f_estimates_relative_midwest_grid <- region_plots$Midwest f_estimates_relative_northeast_grid <- region_plots$Northeast f_estimates_relative_south_grid <- region_plots$South f_estimates_relative_west_grid <- region_plots$West #city governance/population city_inputs <- list( "1K" = list(data = relative_effects_1K, title = "City Population: <1,000", y_breaks = seq(-50, 100, by = 25), y_limits = c(-35, 106), blank_y = FALSE), "1K-10K" = list(data = relative_effects_1K_10K, title = "City Population: 1,000-10,000", y_breaks = seq(-50, 75, by = 25), y_limits = c(-50, 75), blank_y = TRUE), "10K-50K" = list(data = relative_effects_10K_50K, title = "City Population: 10,000-50,000", y_breaks = seq(-50, 75, by = 25), y_limits = c(-50, 75), blank_y = FALSE), "50K-100K" = list(data = relative_effects_50K_100K, title = "City Population: 50,000-100,000", y_breaks = seq(-25, 75, by = 25), y_limits = c(-25, 75), blank_y = TRUE), "100K" = list(data = relative_effects_100K, title = "City Population: >100,000", y_breaks = seq(-25, 100, by = 25), y_limits = c(-25, 100), blank_y = FALSE), "Unincorporated" = list(data = relative_effects_unincorporated, title = "Unincorporated Area", y_breaks = seq(-50, 75, by = 25), y_limits = c(-50, 55), blank_y = TRUE) ) city_plots <- imap(city_inputs, function(city, name) { fun_plot_relative_effect_grid(city$data, city$title, city$y_breaks, city$y_limits, blank_y_axis = city$blank_y) }) f_estimates_relative_1K_grid <- city_plots$"1K" f_estimates_relative_1K_10K_grid <- city_plots$"1K-10K" f_estimates_relative_10K_50K_grid <- city_plots$"10K-50K" f_estimates_relative_50K_100K_grid <- city_plots$"50K-100K" f_estimates_relative_100K_grid <- city_plots$"100K" f_estimates_relative_unincorporated_grid <- city_plots$"Unincorporated" #plots #full sample f_estimates_relative_grid #census region f_estimates_relative_midwest_grid f_estimates_relative_northeast_grid f_estimates_relative_south_grid f_estimates_relative_west_grid #city governance/population f_estimates_relative_1K_grid f_estimates_relative_1K_10K_grid f_estimates_relative_10K_50K_grid f_estimates_relative_50K_100K_grid f_estimates_relative_100K_grid f_estimates_relative_unincorporated_grid #combine individual plots into a single figure empty_plot <- ggplot() + theme_void() legend <- ggpubr::get_legend(f_estimates_relative_grid + theme(legend.position = "bottom")) body_plot_relative <- (f_estimates_relative_grid) / (empty_plot) / (f_estimates_relative_midwest_grid | f_estimates_relative_northeast_grid) / (empty_plot) / (f_estimates_relative_south_grid | f_estimates_relative_west_grid) / (empty_plot) / (f_estimates_relative_1K_grid | f_estimates_relative_1K_10K_grid) / (empty_plot) / (f_estimates_relative_10K_50K_grid | f_estimates_relative_50K_100K_grid) / (empty_plot) / (f_estimates_relative_100K_grid | f_estimates_relative_unincorporated_grid) + plot_layout(heights = c(2, .2, 1, .2, 1, .2, 1, .2, 1, .2, 1)) f_estimates_relative_grid <- cowplot::plot_grid(body_plot_relative, legend, ncol = 1, rel_heights = c(1, 0.08)) f_estimates_relative_grid # ggsave("Estimates_Relative_Effects.pdf", plot = f_estimates_relative_grid, width = 12, height = 18, units = "in") ################################################################################ #function to estimate relative effects for tables fun_reg_fe <- function(data, dep_vars, predictors, selected_predictor, weights) { #create formula for regression with county fixed effects formula <- as.formula(paste(dep_vars, " ~ ", paste(predictors, collapse = " + "), " | GISJOIN")) #create subset of data that is used in regression (i.e., remove NAs) relevant_cols <- c(dep_vars, predictors, weights, "GISJOIN") data_clean <- na.omit(data[relevant_cols]) #estimate regression model <- feols(formula, data = data_clean, weights = data_clean[[weights]]) } #variable: air pollution reg_pollution <- fun_reg_fe(d, "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_pollution, reg_pollution_midwest, reg_pollution_northeast, reg_pollution_south, reg_pollution_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_pollution_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_pollution_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "POLLUTION_PM25", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_pollution_1K, reg_pollution_1K_10K, reg_pollution_10K_50K, reg_pollution_50K_100K, reg_pollution_100K, reg_pollution_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: facility/site regulated by EPA reg_epa <- fun_reg_fe(d, "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_epa, reg_epa_midwest, reg_epa_northeast, reg_epa_south, reg_epa_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_epa_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_epa_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "EPA_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_epa_1K, reg_epa_1K_10K, reg_epa_10K_50K, reg_epa_50K_100K, reg_epa_100K, reg_epa_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: floodplains reg_flood <- fun_reg_fe(d, "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_flood, reg_flood_midwest, reg_flood_northeast, reg_flood_south, reg_flood_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_flood_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_flood_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "FLOOD_100_500_YR", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_flood_1K, reg_flood_1K_10K, reg_flood_10K_50K, reg_flood_50K_100K, reg_flood_100K, reg_flood_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: industry jobs reg_industry_jobs <- fun_reg_fe(d, "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_industry_jobs, reg_industry_jobs_midwest, reg_industry_jobs_northeast, reg_industry_jobs_south, reg_industry_jobs_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_industry_jobs_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_industry_jobs_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "INDUSTRY_JOBS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_industry_jobs_1K, reg_industry_jobs_1K_10K, reg_industry_jobs_10K_50K, reg_industry_jobs_50K_100K, reg_industry_jobs_100K, reg_industry_jobs_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: land cover - green space reg_land_cover_green <- fun_reg_fe(d, "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_land_cover_green, reg_land_cover_green_midwest, reg_land_cover_green_northeast, reg_land_cover_green_south, reg_land_cover_green_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_land_cover_green_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_green_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "LAND_COVER_GREEN_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_land_cover_green_1K, reg_land_cover_green_1K_10K, reg_land_cover_green_10K_50K, reg_land_cover_green_50K_100K, reg_land_cover_green_100K, reg_land_cover_green_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: land cover - trees reg_land_cover_trees <- fun_reg_fe(d, "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_land_cover_trees, reg_land_cover_trees_midwest, reg_land_cover_trees_northeast, reg_land_cover_trees_south, reg_land_cover_trees_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_land_cover_trees_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_trees_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "LAND_COVER_TREES_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_land_cover_trees_1K, reg_land_cover_trees_1K_10K, reg_land_cover_trees_10K_50K, reg_land_cover_trees_50K_100K, reg_land_cover_trees_100K, reg_land_cover_trees_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: land cover - water reg_land_cover_water <- fun_reg_fe(d, "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_land_cover_water, reg_land_cover_water_midwest, reg_land_cover_water_northeast, reg_land_cover_water_south, reg_land_cover_water_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_land_cover_water_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_land_cover_water_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "LAND_COVER_WATER_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_land_cover_water_1K, reg_land_cover_water_1K_10K, reg_land_cover_water_10K_50K, reg_land_cover_water_50K_100K, reg_land_cover_water_100K, reg_land_cover_water_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: schools reg_schools <- fun_reg_fe(d, "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_schools, reg_schools_midwest, reg_schools_northeast, reg_schools_south, reg_schools_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_schools_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_schools_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "SCHOOLS_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_schools_1K, reg_schools_1K_10K, reg_schools_10K_50K, reg_schools_50K_100K, reg_schools_100K, reg_schools_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: SNAP retailers reg_snap <- fun_reg_fe(d, "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_snap, reg_snap_midwest, reg_snap_northeast, reg_snap_south, reg_snap_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_snap_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_snap_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "SNAP_DISTANCE", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_snap_1K, reg_snap_1K_10K, reg_snap_10K_50K, reg_snap_50K_100K, reg_snap_100K, reg_snap_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: street intersections reg_intersections <- fun_reg_fe(d, "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_intersections, reg_intersections_midwest, reg_intersections_northeast, reg_intersections_south, reg_intersections_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_intersections_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_intersections_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "STREET_INTERSECTIONS_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_intersections_1K, reg_intersections_1K_10K, reg_intersections_10K_50K, reg_intersections_50K_100K, reg_intersections_100K, reg_intersections_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) #variable: transit stops reg_transit <- fun_reg_fe(d, "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_midwest <- fun_reg_fe(filter(d, CENSUS_REGION == "Midwest"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_northeast <- fun_reg_fe(filter(d, CENSUS_REGION == "Northeast"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_south <- fun_reg_fe(filter(d, CENSUS_REGION == "South"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_west <- fun_reg_fe(filter(d, CENSUS_REGION == "West"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_transit, reg_transit_midwest, reg_transit_northeast, reg_transit_south, reg_transit_west, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) reg_transit_1K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "<1K"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_1K_10K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "1K-10K"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_10K_50K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "10K-50K"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_50K_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == "50K-100K"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_100K <- fun_reg_fe(filter(d, CITY_POPULATION_BIN == ">100K"), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") reg_transit_NA <- fun_reg_fe(filter(d, is.na(CITY_POPULATION_BIN)), "TRANSIT_STOP_500", vars_controls, vars_mh, "SELECTED_WEIGHT") etable(reg_transit_1K, reg_transit_1K_10K, reg_transit_10K_50K, reg_transit_50K_100K, reg_transit_100K, reg_transit_NA, keep = c("MH_TYPE_PARK", "MH_TYPE_NOT_PARK"), fitstat = c("my", "n"), tex = TRUE) ################################################################################ #weighted CDF ################################################################################ #code plots weighted CDFs #run Weighted_CDF_Function.R that creates stat_ewcdf() function source("Weighted_CDF_Function.R") #function plots weighted CDFs (individual) fun_f_wcdf <- function(data, var, weight, title, xlab_text, breaks, limits, filter_na = TRUE, gini_x, gini_y = 0) { #get variable values (optionally filtered for NA) vals <- data[[var]] if (filter_na) { vals <- vals[!is.na(vals)] wts <- data[[weight]][!is.na(data[[var]])] } else { wts <- data[[weight]] } #calculate Gini coeffient gini_val <- round(gini(vals, weights = wts), 2) #plot ggplot(data, aes(x = .data[[var]], weights = .data[[weight]], group = BUILDING_TYPE, color = BUILDING_TYPE, fill = BUILDING_TYPE)) + stat_ewcdf() + xlab(xlab_text) + ylab("") + ggtitle(title) + scale_x_continuous(breaks = breaks, limits = limits, labels = comma) + scale_fill_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c("black", "blue", "red")) + scale_color_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c("black", "blue", "red")) + annotate("text", x = gini_x, y = gini_y, size = 5, hjust = 0.71, vjust = 0, label = paste("Gini:", gini_val)) + theme_classic() + theme(axis.text.x = element_text(size = 14, color = "black"), axis.text.y = element_text(size = 14, color = "black"), axis.ticks.x = element_line(color = "black"), axis.ticks.y = element_line(color = "black"), axis.title.x = element_text(size = 16, margin = margin(t = 0, b = 20), vjust = -5), axis.title.y = element_text(size = 16, vjust = 2), legend.text = element_text(size = 12), legend.title = element_text(size = 14), legend.position = "bottom", plot.title = element_text(size = 20, hjust = 0.5)) } #air pollution f_wcdf_pollution <- fun_f_wcdf( data = d, var = "POLLUTION_PM25", weight = "SELECTED_WEIGHT", title = "Air Pollution", xlab_text = expression(PM[2.5]~(µg/m^3)), breaks = seq(2, 12, by = 2), limits = c(2, 12), filter_na = TRUE, gini_x = 12 ) #facility/site regulated by EPA f_wcdf_epa <- fun_f_wcdf( data = d, var = "EPA_DISTANCE", weight = "SELECTED_WEIGHT", title = "Facilities Regulated by EPA", xlab_text = "Distance (km)", breaks = seq(0, 8, by = 2), limits = c(0, 8), filter_na = TRUE, gini_x = 8 ) #land cover - green space f_wcdf_land_cover_green <- fun_f_wcdf( data = d, var = "LAND_COVER_GREEN_500", weight = "SELECTED_WEIGHT", title = "Land Cover - Green Space", xlab_text = "Share (%)", breaks = seq(0, 1, by = 0.25), limits = c(0, 1), filter_na = TRUE, gini_x = 1 ) #land cover - trees f_wcdf_land_cover_trees <- fun_f_wcdf( data = d, var = "LAND_COVER_TREES_500", weight = "SELECTED_WEIGHT", title = "Land Cover - Trees", xlab_text = "Share (%)", breaks = seq(0, 1, by = 0.25), limits = c(0, 1), filter_na = TRUE, gini_x = 1 ) #land cover - water f_wcdf_land_cover_water <- fun_f_wcdf( data = d, var = "LAND_COVER_WATER_500", weight = "SELECTED_WEIGHT", title = "Land Cover - Water", xlab_text = "Share (%)", breaks = seq(0, 1, by = 0.25), limits = c(0, 1), filter_na = TRUE, gini_x = 1 ) #schools f_wcdf_schools <- fun_f_wcdf( data = d, var = "SCHOOLS_DISTANCE", weight = "SELECTED_WEIGHT", title = "Schools", xlab_text = "Distance (km)", breaks = seq(0, 20, by = 5), limits = c(0, 20), filter_na = TRUE, gini_x = 20 ) #SNAP retailers f_wcdf_snap <- fun_f_wcdf( data = d, var = "SNAP_DISTANCE", weight = "SELECTED_WEIGHT", title = "SNAP Retailers", xlab_text = "Distance (km)", breaks = seq(0, 20, by = 5), limits = c(0, 20), filter_na = TRUE, gini_x = 20 ) #street intersections f_wcdf_intersections <- fun_f_wcdf( data = d, var = "STREET_INTERSECTIONS_500", weight = "SELECTED_WEIGHT", title = "Street Intersections", xlab_text = "# per sq km", breaks = seq(0, 100, by = 25), limits = c(0, 100), filter_na = TRUE, gini_x = 100 ) # ggsave("Weighted_CDF_Air_Pollution.pdf", plot = f_wcdf_pollution, width = 9, height = 6, units = "in") # ggsave("Weighted_CDF_EPA.pdf", plot = f_wcdf_epa, width = 9, height = 6, units = "in") # ggsave("Weighted_CDF_Land_Cover_Green_Space.pdf", plot = f_wcdf_land_cover_green, width = 9, height = 6, units = "in") # ggsave("Weighted_CDF_Land_Cover_Trees.pdf", plot = f_wcdf_land_cover_trees, width = 9, height = 6, units = "in") # ggsave("Weighted_CDF_Land_Cover_Water.pdf", plot = f_wcdf_land_cover_water, width = 9, height = 6, units = "in") # ggsave("Weighted_CDF_Schools.pdf", plot = f_wcdf_schools, width = 9, height = 6, units = "in") # ggsave("Weighted_CDF_SNAP.pdf", plot = f_wcdf_snap, width = 9, height = 6, units = "in") # ggsave("Weighted_CDF_Street_Intersections.pdf", plot = f_wcdf_intersections, width = 9, height = 6, units = "in") #function to plot weighted CDFs (grid) fun_f_wcdf_grid <- function(data, var, weight, title, xlab_text, breaks, limits, filter_na = TRUE, gini_x, gini_y = 0, gini_size = 4, gini_hjust = 1) { #get variable values (optionally filtered for NA) vals <- data[[var]] if (filter_na) { vals <- vals[!is.na(vals)] wts <- data[[weight]][!is.na(data[[var]])] } else { wts <- data[[weight]] } #calculate Gini coeffient gini_val <- round(gini(vals, weights = wts), 2) ggplot(data, aes(x = .data[[var]], weights = .data[[weight]], group = BUILDING_TYPE, color = BUILDING_TYPE, fill = BUILDING_TYPE)) + stat_ewcdf() + xlab(xlab_text) + ylab("") + ggtitle(title) + scale_x_continuous(breaks = breaks, limits = limits, labels = comma) + scale_fill_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c("black", "blue", "red")) + scale_color_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c("black", "blue", "red")) + annotate("text", x = gini_x, y = gini_y, size = gini_size, hjust = gini_hjust, vjust = 0, label = paste("Gini:", gini_val)) + theme_classic() + theme(axis.text.x = element_text(size = 14, color = "black"), axis.text.y = element_text(size = 14, color = "black"), axis.ticks.x = element_line(color = "black"), axis.ticks.y = element_line(color = "black"), axis.title.y = element_text(size = 16, vjust = 2), legend.text = element_text(size = 12), legend.title = element_text(size = 14), legend.position = "bottom", plot.title = element_text(size = 20, hjust = 0.5)) } #air pollution f_wcdf_pollution_grid <- fun_f_wcdf_grid( data = d, var = "POLLUTION_PM25", weight = "SELECTED_WEIGHT", title = "Air Pollution", xlab_text = expression(PM[2.5]~(µg/m^3)), breaks = seq(2, 12, by = 2), limits = c(2, 12), gini_x = 12 ) #facility/site regulated by EPA f_wcdf_epa_grid <- fun_f_wcdf_grid( data = d, var = "EPA_DISTANCE", weight = "SELECTED_WEIGHT", title = "Facilities Regulated by EPA", xlab_text = "Distance (km)", breaks = seq(0, 8, by = 2), limits = c(0, 8), gini_x = 8 ) #land cover - green space f_wcdf_land_cover_green_grid <- fun_f_wcdf_grid( data = d, var = "LAND_COVER_GREEN_500", weight = "SELECTED_WEIGHT", title = "Land Cover - Green Space", xlab_text = "Share (%)", breaks = seq(0, 1, by = 0.25), limits = c(0, 1), filter_na = TRUE, gini_x = 1 ) #land cover - trees f_wcdf_land_cover_trees_grid <- fun_f_wcdf_grid( data = d, var = "LAND_COVER_TREES_500", weight = "SELECTED_WEIGHT", title = "Land Cover - Trees", xlab_text = "Share (%)", breaks = seq(0, 1, by = 0.25), limits = c(0, 1), filter_na = TRUE, gini_x = 1 ) #land cover - water f_wcdf_land_cover_water_grid <- fun_f_wcdf_grid( data = d, var = "LAND_COVER_WATER_500", weight = "SELECTED_WEIGHT", title = "Land Cover - Water", xlab_text = "Share (%)", breaks = seq(0, 1, by = 0.25), limits = c(0, 1), filter_na = TRUE, gini_x = 1 ) #schools f_wcdf_school_grid <- fun_f_wcdf_grid( data = d, var = "SCHOOLS_DISTANCE", weight = "SELECTED_WEIGHT", title = "Schools", xlab_text = "Distance (km)", breaks = seq(0, 20, by = 5), limits = c(0, 20), filter_na = TRUE, gini_x = 20 ) #SNAP retailers f_wcdf_snap_grid <- fun_f_wcdf_grid( data = d, var = "SNAP_DISTANCE", weight = "SELECTED_WEIGHT", title = "SNAP Retailers", xlab_text = "Distance (km)", breaks = seq(0, 20, by = 5), limits = c(0, 20), filter_na = TRUE, gini_x = 20 ) #street intersections f_wcdf_intersections_grid <- fun_f_wcdf_grid( data = d, var = "STREET_INTERSECTIONS_500", weight = "SELECTED_WEIGHT", title = "Street Intersections", xlab_text = "# per sq km", breaks = seq(0, 100, by = 25), limits = c(0, 100), filter_na = TRUE, gini_x = 100 ) empty_plot <- ggplot() + theme_void() body_plot_wcdf <- (f_wcdf_pollution_grid | f_wcdf_epa_grid | f_wcdf_land_cover_green_grid) / (f_wcdf_land_cover_trees_grid | f_wcdf_land_cover_water_grid | f_wcdf_school_grid) / (f_wcdf_snap_grid | f_wcdf_intersections_grid | empty_plot) f_wcdf <- body_plot_wcdf + plot_layout(guides = "collect") & theme(legend.position = "bottom", axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), axis.title.x = element_text(size = 12, margin = margin(t = 10)), axis.title.y = element_text(size = 12), plot.title = element_text(size = 14), legend.text = element_text(size = 14), legend.title = element_text(size = 14)) # ggsave("Weighted_CDF.pdf", plot = f_wcdf, width = 8.5, height = 11, units = "in") ################################################################################ #income decile plots ################################################################################ load(file = "Manufactured_Homes_Data_Census_Demographics.Rdata") #calculate income decile by census block group deciles_income <- quantile(d_census_demographics$HOUSEHOLD_INCOME_MEDIAN, probs = seq(0, 1, by = 0.1), na.rm = TRUE) deciles_income #create variable: weighted decile d <- d %>% mutate(WEIGHTED_DECILE_INCOME = cut(HOUSEHOLD_INCOME_MEDIAN, breaks = deciles_income, include.lowest = TRUE, labels = 1:10)) #function to calculate weighted average for each group (i.e., "Other", "MH (Community)", "MH (Not Community)") fun_weighted_mean_group <- function(data, group_name, var_names) { mean_weighted <- sapply(var_names, function(v) { weighted.mean(data[data$BUILDING_TYPE == group_name, v], data[data$BUILDING_TYPE == group_name, "SELECTED_WEIGHT"], na.rm = TRUE) }) return(mean_weighted) } #function to calculated weighted average for each group within decile fun_weighted_mean_decile <- function(decile, data, variable_names) { building_types <- c("Other", "MH - Park", "MH - Not Park") mean_list <- map(building_types, function(bt) { fun_weighted_mean_group(filter(data, WEIGHTED_DECILE_INCOME == decile), bt, variable_names) } ) #keep results df <- data.frame(Variable = variable_names, Other = mean_list[[1]], MH_PARK = mean_list[[2]], MH_NOT_PARK = mean_list[[3]], WEIGHTED_DECILE_INCOME = decile) return(df) } #apply above function over all deciles deciles_income_results <- map_dfr(1:10, fun_weighted_mean_decile, d, vars_dep) #convert wide data frame to long deciles_income_results <- deciles_income_results %>% pivot_longer(cols = Other:MH_NOT_PARK, names_to = "Group", values_to = "Value") %>% mutate(Group = factor(Group, levels = c("Other", "MH_PARK", "MH_NOT_PARK"))) #function to plot weighted average within deciles (individual) fun_plot_deciles_income <- function(data, variable, ylab_text, title_text, y_breaks, y_limits, y_lab_expr = NULL) { ggplot(filter(data, Variable == variable), aes(x = WEIGHTED_DECILE_INCOME, y = Value, group = Group, color = Group, fill = Group, shape = Group)) + geom_point(size = 4) + xlab("Income Decile") + { if (is.null(y_lab_expr)) ylab(ylab_text) else ylab(y_lab_expr) } + ggtitle(title_text) + scale_x_continuous(breaks = seq(1, 10, by = 1), limits = c(1, 10)) + scale_y_continuous(breaks = y_breaks, limits = y_limits) + scale_shape_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c(15, 16, 17)) + scale_fill_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c("black", "blue", "red")) + scale_color_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c("black", "blue", "red")) + theme_classic() + theme(axis.text.x = element_text(size = 14, color = "black"), axis.text.y = element_text(size = 14, color = "black"), axis.ticks.x = element_line(color = "black"), axis.ticks.y = element_line(color = "black"), axis.title.x = element_text(size = 16, margin = margin(b = 20), vjust = -5), axis.title.y = element_text(size = 16, vjust = 2), legend.text = element_text(size = 12), legend.title = element_text(size = 14), legend.position = "bottom", plot.title = element_text(size = 20, hjust = 0.5)) } #air pollution f_deciles_income_pollution <- fun_plot_deciles_income( data = deciles_income_results, variable = "POLLUTION_PM25", ylab_text = NULL, title_text = "Air Pollution", y_breaks = seq(9.5, 11, by = 0.5), y_limits = c(9.5, 11), y_lab_expr = expression(PM2.5~(µg/m^{3})) ) #facility/site regulated by EPA f_deciles_income_epa <- fun_plot_deciles_income( data = deciles_income_results, variable = "EPA_DISTANCE", ylab_text = "Distance (km)", title_text = "Facilities Regulated by EPA", y_breaks = seq(0, 2, by = .5), y_limits = c(.5, 2.25) ) #floodplains f_deciles_income_flood <- fun_plot_deciles_income( data = deciles_income_results, variable = "FLOOD_100_500_YR", ylab_text = "Share (%)", title_text = "Floodplains", y_breaks = seq(0.06, 0.14, by = 0.02), y_limits = c(0.055, 0.14) ) #industry jobs f_deciles_income_industry_jobs <- fun_plot_deciles_income( data = deciles_income_results, variable = "INDUSTRY_JOBS_500", ylab_text = "Share (%)", title_text = "Industry Jobs", y_breaks = seq(0, 0.8, by = 0.2), y_limits = c(0, 0.8) ) #land cover - green space f_deciles_income_land_cover_green <- fun_plot_deciles_income( data = deciles_income_results, variable = "LAND_COVER_GREEN_500", ylab_text = "Share (%)", title_text = "Land Cover - Green Space", y_breaks = seq(0.2, 0.8, by = 0.2), y_limits = c(0.2, 0.8) ) #land cover - trees f_deciles_income_land_cover_trees <- fun_plot_deciles_income( data = deciles_income_results, variable = "LAND_COVER_TREES_500", ylab_text = "Share (%)", title_text = "Land Cover - Trees", y_breaks = seq(0.1, 0.3, by = 0.05), y_limits = c(0.1, 0.3) ) #land cover - water f_deciles_income_land_cover_water <- fun_plot_deciles_income( data = deciles_income_results, variable = "LAND_COVER_WATER_500", ylab_text = "Share (%)", title_text = "Land Cover - Water", y_breaks = seq(0.005, 0.02, by = 0.005), y_limits = c(0.005, 0.021) ) #schools f_deciles_income_schools <- fun_plot_deciles_income( data = deciles_income_results, variable = "SCHOOLS_DISTANCE", ylab_text = "Distance (km)", title_text = "Schools", y_breaks = seq(1, 6, by = 1), y_limits = c(1, 6) ) #SNAP retailers f_deciles_income_snap <- fun_plot_deciles_income( data = deciles_income_results, variable = "SNAP_DISTANCE", ylab_text = "Distance (km)", title_text = "SNAP Retailers", y_breaks = seq(1, 5, by = 1), y_limits = c(0.9, 5) ) #street intersections f_deciles_income_intersections <- fun_plot_deciles_income( data = deciles_income_results, variable = "STREET_INTERSECTIONS_500", ylab_text = "# per sq km", title_text = "Street Intersections", y_breaks = seq(10, 50, by = 10), y_limits = c(10, 50) ) #transit stops f_deciles_income_transit <- fun_plot_deciles_income( data = deciles_income_results, variable = "TRANSIT_STOP_500", ylab_text = "Share (%)", title_text = "Transit Stops", y_breaks = seq(0, 0.3, by = 0.1), y_limits = c(0, 0.3) ) f_deciles_income_pollution f_deciles_income_epa f_deciles_income_flood f_deciles_income_industry_jobs f_deciles_income_land_cover_green f_deciles_income_land_cover_trees f_deciles_income_land_cover_water f_deciles_income_schools f_deciles_income_snap f_deciles_income_intersections f_deciles_income_transit # ggsave("Deciles_Income_Pollution.pdf", plot = f_deciles_income_pollution, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_EPA.pdf", plot = f_deciles_income_epa, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_Flood.pdf", plot = f_deciles_income_flood, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_Industry_Jobs.pdf", plot = f_deciles_income_industry_jobs, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_Land_Cover_Green.pdf", plot = f_deciles_income_land_cover_green, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_Land_Cover_Trees.pdf", plot = f_deciles_income_land_cover_trees, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_Land_Cover_Water.pdf", plot = f_deciles_income_land_cover_water, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_Schools.pdf", plot = f_deciles_income_schools, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_SNAP.pdf", plot = f_deciles_income_snap, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_Street_Intersections.pdf", plot = f_decile_income_intersections, width = 9, height = 6, units = "in") # ggsave("Deciles_Income_Transit_Stop.pdf", plot = f_decile_income_transit, width = 9, height = 6, units = "in") #function to plot weighted average within deciles (grid) fun_plot_deciles_income_grid <- function(data, variable, ylab_text, title_text, y_transform = identity, y_breaks = NULL, y_limits = NULL) { ggplot(filter(data, Variable == variable), aes(x = WEIGHTED_DECILE_INCOME, y = y_transform(Value), group = Group, color = Group, fill = Group, shape = Group)) + geom_point(size = 2) + xlab("") + ylab(ylab_text) + ggtitle(title_text) + scale_x_continuous(breaks = seq(1, 10, by = 1), limits = c(1, 10)) + scale_y_continuous(breaks = y_breaks, limits = y_limits) + scale_shape_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c(15, 16, 17)) + scale_fill_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c("black", "blue", "red")) + scale_color_manual(name = "Building Type:", labels = c("Other", "MH (Community)", "MH (Not Community)"), values = c("black", "blue", "red")) + theme_classic() + theme(axis.text.x = element_text(size = 14, color = "black"), axis.text.y = element_text(size = 14, color = "black"), axis.ticks.x = element_line(color = "black"), axis.ticks.y = element_line(color = "black"), axis.title.x = element_text(size = 16, margin = margin(b = 20), vjust = -5), axis.title.y = element_text(size = 16, vjust = 2), legend.text = element_text(size = 12), legend.title = element_text(size = 14), legend.position = "bottom", plot.title = element_text(size = 20, hjust = 0.5)) } #air pollution f_deciles_income_pollution_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "POLLUTION_PM25", ylab_text = expression(PM2.5~(µg/m^{3})), title_text = "Air Pollution", y_breaks = seq(9.5, 11, by = 0.5), y_limits = c(9.5, 11) ) #facility/site regulated by EPA f_deciles_income_epa_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "EPA_DISTANCE", ylab_text = "Distance (km)", title_text = "Facilities Regulated by EPA", y_breaks = seq(0, 2, by = .5), y_limits = c(.5, 2.25) ) #floodplains f_deciles_income_flood_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "FLOOD_100_500_YR", ylab_text = "Share (%)", title_text = "Floodplains", y_breaks = seq(0.06, 0.14, by = 0.02), y_limits = c(0.055, 0.14) ) #industry jobs f_deciles_income_industry_jobs_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "INDUSTRY_JOBS_500", ylab_text = "Share (%)", title_text = "Industry Jobs", y_breaks = seq(0, 0.8, by = 0.2), y_limits = c(0, 0.8) ) #land cover - green space f_deciles_income_land_cover_green_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "LAND_COVER_GREEN_500", ylab_text = "Share (%)", title_text = "Land Cover - Green Space", y_breaks = seq(0.2, 0.8, by = 0.2), y_limits = c(0.2, 0.8) ) #land cover - trees f_deciles_income_land_cover_trees_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "LAND_COVER_TREES_500", ylab_text = "Share (%)", title_text = "Land Cover - Trees", y_breaks = seq(0.1, 0.3, by = 0.05), y_limits = c(0.1, 0.3) ) #land cover - water f_deciles_income_land_cover_water_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "LAND_COVER_WATER_500", ylab_text = "Share (%)", title_text = "Land Cover - Water", y_breaks = seq(0.005, 0.02, by = 0.005), y_limits = c(0.005, 0.021) ) #schools f_deciles_income_schools_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "SCHOOLS_DISTANCE", ylab_text = "Distance (km)", title_text = "Schools", y_breaks = seq(1, 6, by = 1), y_limits = c(1, 6) ) #SNAP retailers f_deciles_income_snap_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "SNAP_DISTANCE", ylab_text = "Distance (km)", title_text = "SNAP Retailers", y_breaks = seq(1, 5, by = 1), y_limits = c(0.9, 5) ) #street intersections f_deciles_income_intersections_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "STREET_INTERSECTIONS_500", ylab_text = "# per sq km", title_text = "Street Intersections", y_breaks = seq(10, 50, by = 10), y_limits = c(10, 50) ) #transit stops f_deciles_income_transit_grid <- fun_plot_deciles_income_grid( data = deciles_income_results, variable = "TRANSIT_STOP_500", ylab_text = "Share (%)", title_text = "Transit Stops", y_breaks = seq(0, 0.3, by = 0.1), y_limits = c(0, 0.3) ) f_deciles_income_pollution_grid f_deciles_income_epa_grid f_deciles_income_flood_grid f_deciles_income_industry_jobs_grid f_deciles_income_land_cover_green_grid f_deciles_income_land_cover_trees_grid f_deciles_income_land_cover_water_grid f_deciles_income_schools_grid f_deciles_income_snap_grid f_deciles_income_intersections_grid f_deciles_income_transit_grid empty_plot <- ggplot() + theme_void() body_plot_deciles <- (f_deciles_income_pollution_grid | f_deciles_income_epa_grid | f_deciles_income_flood_grid) / (f_deciles_income_industry_jobs_grid | f_deciles_income_land_cover_green_grid | f_deciles_income_land_cover_trees_grid) / (f_deciles_income_land_cover_water_grid | f_deciles_income_schools_grid | f_deciles_income_snap_grid) / (f_deciles_income_intersections_grid | f_deciles_income_transit_grid | empty_plot) f_deciles_income <- body_plot_deciles + plot_layout(guides = "collect") & theme(legend.position = "bottom", axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10), axis.title.x = element_text(size = 12), axis.title.y = element_text(size =12, margin = margin(r = 10)), plot.margin = margin(t = 5, r = 5, b = 5, l = 5), plot.title = element_text(size = 14), legend.text = element_text(size = 14), legend.title = element_text(size = 14)) f_deciles_income # ggsave("Deciles_Income.pdf", plot = f_deciles_income, width = 8.5, height = 11, units = "in") ################################################################################ #model performance ################################################################################ #have to load package after deciles plot (because of map_dfr()) library(maps) #import data: counties load(file = "Manufactured_Homes_Data_Counties.Rdata") d_summary_county_mh <- d %>% group_by(GISJOIN) %>% summarise(MH_PCT = sum(MH_COUNT) / sum(FOOTPRINT_COUNT), MH_TYPE_PARK_PCT = sum(MH_TYPE_PARK) / sum(MH)) d_counties <- left_join(d_counties, d_summary_county_mh, by = c("GISJOIN"), na_matches = "never") #import data: states d_states <- maps::map("state", fill=TRUE, plot =FALSE) %>% st_as_sf() %>% st_transform(st_crs(d_counties)) ################################################################################ #map: pct of MH (denominator = all buildings) f_map_mh_pct <- ggplot() + geom_sf(data = d_counties, aes(fill = MH_PCT), color = NA) + ggtitle("Share of Buildings that are \n Manufactured Housing Units") + scale_fill_viridis_c(name = "Share (%)", direction = -1, breaks = seq(0, 0.3, by = 0.1), limits = c(0, 0.3)) + geom_sf(data = d_states, fill = NA, color = "black", size = 0.2) + theme_classic() + theme(line = element_blank(), legend.title = element_text(size = 16, vjust = 5), legend.text = element_text(size = 12), plot.title = element_text(size = 20, hjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks.x = element_blank(), axis.ticks.y = element_blank()) # ggsave("Map_Share_MH.pdf", plot = f_map_mh_pct, width = 9, height = 6, units = "in") f_map_mh_pct_grid <- ggplot() + geom_sf(data = d_counties, aes(fill = MH_PCT), color = NA) + ggtitle("Share of Buildings that are \n Manufactured Housing Units") + scale_fill_viridis_c(name = "Share (%)", direction = -1, breaks = seq(0, 0.3, by = 0.1), limits = c(0, 0.3)) + geom_sf(data = d_states, fill = NA, color = "black", size = 0.2) + theme_classic() + theme(line = element_blank(), legend.title = element_text(size = 10, vjust = 5), legend.text = element_text(size = 8), plot.title = element_text(size = 12, hjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks.x = element_blank(), axis.ticks.y = element_blank()) ################################################################################ #map: pct of MH types (denominator = MH) f_map_mh_park_pct <- ggplot() + geom_sf(data = d_counties, aes(fill = MH_TYPE_PARK_PCT), color = NA) + ggtitle("Share of Manufactured Housing Units that are located within MHCs") + scale_fill_viridis_c(name = "Share (%)", direction = -1, breaks = seq(0, 1, by = 0.25), limits = c(0, 1)) + geom_sf(data = d_states, fill = NA, color = "black", size = 0.2) + theme_classic() + theme(line = element_blank(), legend.title = element_text(size = 16, vjust = 5), legend.text = element_text(size = 12), plot.title = element_text(size = 20, hjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks.x = element_blank(), axis.ticks.y = element_blank()) # ggsave("Map_Share_MH_Community.pdf", plot = f_map_mh_park_pct, width = 9, height = 6, units = "in") f_map_mh_park_pct_grid <- ggplot() + geom_sf(data = d_counties, aes(fill = MH_TYPE_PARK_PCT), color = NA) + ggtitle("Share of Manufactured Housing Units \n that are located within MHCs") + scale_fill_viridis_c(name = "Share (%)", direction = -1, breaks = seq(0, 1, by = 0.2), limits = c(0, 1)) + geom_sf(data = d_states, fill = NA, color = "black", size = 0.2) + theme_classic() + theme(line = element_blank(), legend.title = element_text(size = 10, vjust = 5), legend.text = element_text(size = 8), plot.title = element_text(size = 12, hjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks.x = element_blank(), axis.ticks.y = element_blank()) ################################################################################ #plot: correlation between building footprints and ACS data #import data: building footprints and ACS load(file = "Manufactured_Homes_Data_Comparison.Rdata") #calculate correlation between Pct_MH_foot and Pct_MH_ACS correlation <- cor(d_mh_comparison$PCT_MH_FOOTPRINT, d_mh_comparison$PCT_MH_ACS, use = "complete.obs") # Plot Pct_MH_foot against Pct_MH_ACS f_correlation <- ggplot(d_mh_comparison, aes(x = PCT_MH_FOOTPRINT, y = PCT_MH_ACS)) + geom_point(size = .5) + geom_smooth(method = "lm", se = FALSE, color = "blue") + xlab("MH - Building Footprints (%)") + ylab("MH - ACS (%)") + ggtitle("Correlation between Building Footprints and ACS") + scale_x_continuous(breaks = seq(0, .25, by = .05), limits = c(0, .25)) + scale_y_continuous(breaks = seq(0, .80, by = .20), limits = c(0, .80)) + annotate(geom = "text", x = .25, y = 0, hjust = 1, label = paste("r =", round(correlation, 2)), size = 5) + theme_classic() + theme(axis.text.x = element_text(size = 14), axis.text.y = element_text(size = 14), axis.title.x = element_text(size = 16, margin = margin(t = 0, r = 0, b = 20, l = 0), vjust = -5), axis.title.y = element_text(size = 16, vjust = 2), legend.text = element_text(size = 12), legend.title = element_text(size = 14), legend.position = "bottom", plot.title = element_text(size = 20, hjust = 0.5)) f_correlation f_correlation_grid <- ggplot(d_mh_comparison, aes(x = PCT_MH_FOOTPRINT, y = PCT_MH_ACS)) + geom_point(size = .5) + geom_smooth(method = "lm", se = FALSE, color = "blue") + xlab("MH - Building Footprints (%)") + ylab("MH - ACS (%)") + ggtitle("Correlation between Building Footprints \n and American Community Survey") + scale_x_continuous(breaks = seq(0, .25, by = .05), limits = c(0, .25)) + scale_y_continuous(breaks = seq(0, .80, by = .20), limits = c(0, .80)) + annotate(geom = "text", x = .25, y = 0, hjust = 1, label = paste("r =", round(correlation, 2)), size = 2.8) + theme_classic() + theme(axis.text.x = element_text(size = 8), axis.text.y = element_text(size = 8), axis.title.x = element_text(size = 10), axis.title.y = element_text(size = 10, vjust = 2), plot.title = element_text(size = 12, hjust = 0.5)) f_correlation_grid ################################################################################ #map: difference between building footprints and ACS data f_map_mh_difference <- ggplot() + geom_sf(data = d_mh_comparison, aes(fill = PCT_MH_DIFFERENCE), color = NA) + ggtitle("Difference between Building Footprints /n and American Community Survey") + scale_fill_viridis_c(name = "Percentage", direction = -1) + geom_sf(data = d_states, fill = NA, color = "black", size = 0.2) + theme_classic() + theme(line = element_blank(), legend.title = element_text(size = 16, vjust = 5), legend.text = element_text(size = 12), plot.title = element_text(size = 20, hjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks.x = element_blank(), axis.ticks.y = element_blank()) # ggsave("MH_Correlation_Map.pdf", plot = f_map_mh_difference, width = 9, height = 6, units = "in") f_map_mh_difference_grid <- ggplot() + geom_sf(data = d_mh_comparison, aes(fill = PCT_MH_DIFFERENCE), color = NA) + ggtitle("Difference between Building Footprints \n and American Community Survey") + scale_fill_viridis_c(name = "Share (%)", direction = -1) + geom_sf(data = d_states, fill = NA, color = "black", size = 0.2) + theme_classic() + theme(line = element_blank(), legend.title = element_text(size = 10, vjust = 5), legend.text = element_text(size = 8), plot.title = element_text(size = 12, hjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks.x = element_blank(), axis.ticks.y = element_blank()) ################################################################################ #combine above figures into a single plot with (a)-(e) labels #upload latex summary table latex_table_image <- readPNG("Model_Performance.png") %>% rasterGrob(interpolate = TRUE) %>% ggdraw() # Create a plot layout grid.newpage() pushViewport(viewport(layout = grid.layout(3, 2))) latex_table_image_label <- ggdraw(latex_table_image) + draw_plot_label("A", x = 0.01, y = 0.99, vjust = 1, fontface = "bold") f_map_mh_pct_label <- ggdraw(f_map_mh_pct_grid) + draw_plot_label("B", x = 0.01, y = 0.99, vjust = 1, fontface = "bold") f_correlation_label <- ggdraw() + draw_plot(f_correlation_grid, y = -0.07) + draw_plot_label("C", x = 0.01, y = 0.99, vjust = 1, fontface = "bold") f_map_mh_difference_label <- ggdraw(f_map_mh_difference_grid) + draw_plot_label("D", x = 0.01, y = 0.99, vjust = 1, fontface = "bold") f_map_mh_park_pct_label <- ggdraw(f_map_mh_park_pct_grid) + draw_plot_label("E", x = 0.01, y = 0.99, vjust = 1, fontface = "bold") empty_plot <- ggplot() + theme_void() body_plot_model <- (latex_table_image_label) / (f_map_mh_pct_label | f_correlation_label) / (f_map_mh_park_pct_label | f_map_mh_difference_label) f_model_performance <- body_plot_model + plot_layout(guides = "collect", heights = c(.75, 1, 1), tag_level = "new") + theme(plot.margin = margin(t = 5, r = 5, b = 5, l = 5)) # ggsave("Map_Model_Performance.pdf", plot = f_model_performance, width = 10, height = 10, units = "in") ################################################################################ #maps: access/exposure ################################################################################ #function to map relative difference in access/exposure #between MH (Community) and other buildings by county fun_map_exposure <- function(d, counties, states, var, title, legend_title, breaks, limits, colors = c("low" = "midnightblue", "mid" = "white", "high" = "red4"), midpoint = 0) { #create data frame to summarize by county (GISJOIN) and building type d_summary <- d %>% filter(MH_TYPE %in% c("Park", "Other")) %>% group_by(GISJOIN, MH_TYPE_PARK) %>% summarise(MEAN = mean(.data[[var]], na.rm = TRUE), .groups = "drop") %>% pivot_wider(names_from = MH_TYPE_PARK, values_from = MEAN) %>% mutate(DIFFERENCE = `1` - `0`) #join to d_counties counties <- left_join(counties, dplyr::select(d_summary, GISJOIN, DIFFERENCE), by = "GISJOIN", na_matches = "never") ggplot() + geom_sf(data = counties, aes(fill = DIFFERENCE), color = NA) + ggtitle(title) + scale_fill_gradient2(name = legend_title, low = colors["low"], mid = colors["mid"], high = colors["high"], midpoint = midpoint, breaks = breaks, limits = limits) + geom_sf(data = states, fill = NA, color = "black", size = 0.2) + theme_classic() + theme(line = element_blank(), legend.title = element_text(size = 16, vjust = 5), legend.text = element_text(size = 12), plot.title = element_text(size = 20, hjust = 0.5), axis.text.x = element_blank(), axis.text.y = element_blank(), axis.ticks.x = element_blank(), axis.ticks.y = element_blank()) } #air pollution f_map_exposure_pollution <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "POLLUTION_PM25", title = "Air Pollution", legend_title = expression(PM[2.5]~(µg/m^3)), breaks = seq(-4, 4, by = 2), limits = c(-4, 4), colors = c("low" = "midnightblue", "mid" = "white", "high" = "red4") ) #facility/site regulated by EPA f_map_exposure_epa <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "EPA_DISTANCE", title = "Facilities Regulated by EPA", legend_title = "Distance (km)", breaks = seq(-8, 4, by = 4), limits = c(-8, 4), colors = c("low" = "red4", "mid" = "white", "high" = "midnightblue") ) #floodplains f_map_exposure_flood <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "FLOOD_100_500_YR", title = "Floodplains", legend_title = "Share (%)", breaks = seq(-0.5, 1, by = 0.5), limits = c(-0.5, 1), colors = c("low" = "midnightblue", "mid" = "white", "high" = "red4") ) #industry jobs f_map_exposure_industry_jobs <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "INDUSTRY_JOBS_500", title = "Industry Jobs", legend_title = "Share (%)", breaks = seq(-1, 1, by = 0.5), limits = c(-1, 1), colors = c("low" = "midnightblue", "mid" = "white", "high" = "red4") ) #land cover - green space f_map_exposure_land_cover_green <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "LAND_COVER_GREEN_500", title = "Land Cover - Green Space", legend_title = "Share (%)", breaks = seq(-0.8, 0.4, by = 0.4), limits = c(-0.8, 0.4), colors = c("low" = "red4", "mid" = "white", "high" = "midnightblue") ) #land cover - trees f_map_exposure_land_cover_trees <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "LAND_COVER_TREES_500", title = "Land Cover - Trees", legend_title = "Share (%)", breaks = seq(-0.5, 0.5, by = 0.25), limits = c(-0.5, 0.5), colors = c("low" = "red4", "mid" = "white", "high" = "midnightblue") ) #land cover - water f_map_exposure_land_cover_water <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "LAND_COVER_WATER_500", title = "Land Cover - Water", legend_title = "Share (%)", breaks = seq(-0.3, 0.6, by = 0.3), limits = c(-0.3, 0.6), colors = c("low" = "red4", "mid" = "white", "high" = "midnightblue") ) #schools f_map_exposure_schools <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "SCHOOLS_DISTANCE", title = "Schools", legend_title = "Distance (km)", breaks = seq(-15, 30, by = 15), limits = c(-15, 30), colors = c("low" = "midnightblue", "mid" = "white", "high" = "red4") ) #SNAP retailers f_map_exposure_snap <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "SNAP_DISTANCE", title = "SNAP Retailers", legend_title = "Distance (km)", breaks = seq(-15, 30, by = 15), limits = c(-15, 30), colors = c("low" = "midnightblue", "mid" = "white", "high" = "red4") ) #street intersections f_map_exposure_intersections <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "STREET_INTERSECTIONS_500", title = "Street Intersections", legend_title = "# per sq km", breaks = seq(-80, 80, by = 40), limits = c(-80, 80), colors = c("low" = "red4", "mid" = "white", "high" = "midnightblue") ) #transit stops f_map_exposure_transit <- fun_map_exposure( d = d, counties = d_counties, states = d_states, var = "TRANSIT_STOP_500", title = "Transit Stops", legend_title = "Share (%)", breaks = seq(-1, 1, by = 0.5), limits = c(-1, 1), colors = c("low" = "red4", "mid" = "white", "high" = "midnightblue") ) # ggsave("Map_Exposure_Air_Pollution.pdf", plot = f_map_exposure_pollution, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_EPA.pdf", plot = f_map_exposure_epa, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_Floodplains.pdf", plot = f_map_exposure_flood, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_Industry_Jobs.pdf", plot = f_map_exposure_industry_jobs, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_Land_Cover_Green_Space.pdf", plot = f_map_exposure_land_cover_green, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_Land_Cover_Trees.pdf", plot = f_map_exposure_land_cover_trees, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_Land_Cover_Water.pdf", plot = f_map_exposure_land_cover_water, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_Schools.pdf", plot = f_map_exposure_schools, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_SNAP.pdf", plot = f_map_exposure_snap, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_Street_Intersections.pdf", plot = f_map_exposure_intersections, width = 9, height = 6, units = "in") # ggsave("Map_Exposure_Transit_Stops.pdf", plot = f_map_exposure_transit, width = 9, height = 6, units = "in") empty_plot <- ggplot() + theme_void() body_plot_map_exposure <- (f_map_exposure_pollution | f_map_exposure_epa | f_map_exposure_flood) / (f_map_exposure_industry_jobs | f_map_exposure_land_cover_green | f_map_exposure_land_cover_trees) / (f_map_exposure_land_cover_water | f_map_exposure_schools | f_map_exposure_snap) / (f_map_exposure_intersections | f_map_exposure_transit | empty_plot) f_map_exposure <- body_plot_map_exposure & theme(plot.title = element_text(size = 14), legend.position = "right", legend.key.height = unit(0.4, "cm"), legend.key.width = unit(0.2, "cm"), legend.text = element_text(size = 6), legend.title = element_text(size = 8)) # ggsave("Map_Exposure.pdf", plot = f_map_exposure, width = 10, height = 8, units = "in")