lrgv_23_grandparents <- get_acs(geography = "tract", county = c("Cameron"), table = "S1002", state = "TX", year = 2023, survey = "acs5", geometry = TRUE) %>% left_join(., all_acs_vars, by = c("variable"="name")) %>% separate_wider_delim(., cols = "label", delim = "!!", names_sep="", too_few = "align_start") lrgv_23_grandparents %>% filter(variable == "S1002_C01_001") %>% st_as_sf(.) %>% tm_shape(.) + tm_polygons(col = c("estimate"), popup.vars = c("Name" = "label3", "NAME", "Estimate"="estimate"), alpha = 0.5, style = "jenks") + tm_basemap("OpenStreetMap")

lrgv_23_grandparents <- get_acs(geography = "tract", county = c("Cameron"), table = "S1002", state = "TX", year = 2023, survey = "acs5", geometry = TRUE) %>% left_join(., all_acs_vars, by = c("variable"="name")) %>% separate_wider_delim(., cols = "label", delim = "!!", names_sep="", too_few = "align_start") custom_Bville<- read_sf(dsn = "C:/Users/daniel.pinon/Downloads/my_custom_brownsville_shape/map.geojson") %>% as.data.frame(.) %>% mutate(location = c("inner_Brownsville", "outer_Brownsville")) %>% st_as_sf(.) %>% st_transform(., crs = 4269) lrgv_23_grandparents %>% filter(variable == "S1002_C01_001") %>% st_as_sf(.) %>% st_intersection(., st_buffer(cam_places %>% filter(NAME == "Brownsville"), dist = 10)) %>% tm_shape(.) + tm_polygons(col = c("estimate"), popup.vars = c("Name" = "label3", "NAME", "Estimate"="estimate"), alpha = 0.5, style = "jenks") + tm_basemap("OpenStreetMap")

lrgv_23_grandparents <- get_acs(geography = "tract", county = c("Cameron"), table = "S1002", state = "TX", year = 2023, survey = "acs5", geometry = TRUE) %>% left_join(., all_acs_vars, by = c("variable"="name")) %>% separate_wider_delim(., cols = "label", delim = "!!", names_sep="", too_few = "align_start") custom_Bville<- read_sf(dsn = "C:/Users/daniel.pinon/Downloads/my_custom_brownsville_shape/map.geojson") %>% as.data.frame(.) %>% mutate(location = c("inner_Brownsville", "outer_Brownsville")) %>% st_as_sf(.) %>% st_transform(., crs = 4269) lrgv_23_grandparents %>% filter(variable == "S1002_C01_001") %>% st_as_sf(.) %>% st_intersection(., st_buffer(st_union(custom_Bville), dist = 10)) %>% tm_shape(.) + tm_polygons(col = c("estimate"), popup.vars = c("Name" = "label3", "NAME", "Estimate"="estimate"), alpha = 0.5, style = "jenks") + tm_basemap("OpenStreetMap")

This is a table detailing the incorporated and unincorporated communities throughout Cameron County using data spanning from 2023; this represents the most up-to-date information available publicly, as far as I know. A Census-Designated Place (CDP), per the U.S. Census Bureau, is "...[a] statistical equivalent of [an] incorporated place and represent[s] unincorporated communities that do not have a legally defined boundary or an active, functioning governmental structure." Source: https://www.census.gov/programs-surveys/bas/information/cdp.html ... code: library(easypackages) libraries(c("readxl", "ggmap", "ggiraph", "ggforce", "ggcorrplot", "ggthemes", "ggsignif", "ggsflabel", "ggrepel", "ggpubr", "ggsci", "glue", "gt", "janitor", "maptools", "mapview", "magrittr", "plyr", "prettyunits", "progress", "progressr", "psych", "rgeos", "rio", "rms", "Hmisc", "robustbase", "rspat", "s2", "sfheaders", "sfweight", "snakecase", "smoothr", "sp", "spatial", "spatialEco", "spatstat", "spatstat.linnet", "spatstat.model", "rpart", "spatstat.explore", "nlme", "spatstat.random", "spatstat.geom", "spatstat.data", "spdep", "sf", "spData", "abind", "summarytools", "terra", "tidycensus", "tidylog", "tidyselect", "lubridate", "forcats", "stringr", "dplyr", "purrr", "readr", "tidyr", "tibble", "ggplot2", "tidyverse", "tigris", "tmap", "vctrs", "viridis", "viridisLite", "vroom", "waldo", "wk", "stats", "graphics", "grDevices", "utils", "datasets", "methods", "base", "haven", "foreign", "survey", "srvyr", "sitrep", "questionr", "srvyr", "stringr")) ##tigris package, pulls Census data TX_places<- places(state = "TX", cb = FALSE, year = 2023) cam_pop_tracts<- get_acs("tract", table = "B01001", state = "TX", county = "Cameron", year= 2022, survey = "acs5", geometry = TRUE) # colonias data available through the Texas OAG database via filtering # https://texasoag.maps.arcgis.com/apps/webappviewer/index.html?id=1bc9c4f7b1da47dd8fc535fbd17dc060 # the file I'm using comes from the Cameron County DOT though, I'm pretty sure, circa 05/2022 colonias_sf<- read_sf(dsn = "E:/COLONIAS/COLONIAS.shp") ## Step 1: (after loading everything in) make a "cookie cutter" to filter the places data cam_pop_tracts %>% st_as_sf(.) %>% st_union(.) -> cam_union ## Step 2: Cookie Cut! st_intersection(TX_places, cam_union) -> cam_places ## Step 3: Make the table of places cam_places %>% mutate(TYPE = str_remove(NAMELSAD, paste0(NAME)) %>% str_squish(.) %>% toupper(.)) %>% st_drop_geometry(.) %>% arrange(., desc(TYPE)) %>% gt(.) %>% cols_align("center") %>% opt_stylize("4")

This is a map detailing the incorporated and unincorporated communities throughout Cameron County using data spanning from 2022 to 2023; this represents the most up-to-date information available publicly, as far as I know. A Census-Designated Place (CDP), per the U.S. Census Bureau, is "...[a] statistical equivalent of [an] incorporated place and represent[s] unincorporated communities that do not have a legally defined boundary or an active, functioning governmental structure." Source: https://www.census.gov/programs-surveys/bas/information/cdp.html ... code: library(easypackages) libraries(c("readxl", "ggmap", "ggiraph", "ggforce", "ggcorrplot", "ggthemes", "ggsignif", "ggsflabel", "ggrepel", "ggpubr", "ggsci", "glue", "gt", "janitor", "maptools", "mapview", "magrittr", "plyr", "prettyunits", "progress", "progressr", "psych", "rgeos", "rio", "rms", "Hmisc", "robustbase", "rspat", "s2", "sfheaders", "sfweight", "snakecase", "smoothr", "sp", "spatial", "spatialEco", "spatstat", "spatstat.linnet", "spatstat.model", "rpart", "spatstat.explore", "nlme", "spatstat.random", "spatstat.geom", "spatstat.data", "spdep", "sf", "spData", "abind", "summarytools", "terra", "tidycensus", "tidylog", "tidyselect", "lubridate", "forcats", "stringr", "dplyr", "purrr", "readr", "tidyr", "tibble", "ggplot2", "tidyverse", "tigris", "tmap", "vctrs", "viridis", "viridisLite", "vroom", "waldo", "wk", "stats", "graphics", "grDevices", "utils", "datasets", "methods", "base", "haven", "foreign", "survey", "srvyr", "sitrep", "questionr", "srvyr", "stringr")) ##tigris package, pulls Census data TX_places<- places(state = "TX", cb = FALSE, year = 2023) cam_pop_tracts<- get_acs("tract", table = "B01001", state = "TX", county = "Cameron", year= 2022, survey = "acs5", geometry = TRUE) # colonias data available through the Texas OAG database via filtering # https://texasoag.maps.arcgis.com/apps/webappviewer/index.html?id=1bc9c4f7b1da47dd8fc535fbd17dc060 # the file I'm using comes from the Cameron County DOT though, I'm pretty sure, circa 05/2022 colonias_sf<- read_sf(dsn = "E:/COLONIAS/COLONIAS.shp") ## Step 1: (after loading everything in) make a "cookie cutter" to filter the places data cam_pop_tracts %>% st_as_sf(.) %>% st_union(.) -> cam_union ## Step 2: Cookie Cut! st_intersection(TX_places, cam_union) -> cam_places ## Step 3: Make the map of places (unincorporated) cam_places %>% mutate(TYPE = str_remove(NAMELSAD, paste0(NAME)) %>% str_squish(.) %>% toupper(.)) %>% st_as_sf(.) %>% tm_shape(.) + tm_polygons(col = "TYPE", popup.vars = c("NAME", "TYPE"), alpha = 0.6) + tm_basemap("OpenStreetMap") -> bigboss_CDP_map ## Step 4 (optional): Show relevant colonias too based on size using "tm_bubbles" bigboss_CDP_map + tm_shape(st_make_valid(st_as_sf(colonias_sf))) + tm_bubbles(size="EST_POP", fill = "COLONIA_NM", col="purple", alpha = 0.4, popup.vars=c("COLONIA_NM", "EST_POP", "INCORPORAT"))

Note: This map IS stratified by age (i.e. early to mid to late life) as was done in the Lancet Commission reports -- all risk factors were combined in aggregate *after taking into account the age group for which the risk factor is relevant*. Moreover, "age stratified" in the sense that it's used here does NOT reflect a weighting of all risk factors per age group, as might be done but, instead, reflects inclusion of estimates ONLY belonging to the age groups identified by the Lancet Commission such that no greater that 12 risk factors -- which were included in the construction of the composite index -- are considered for any one PHR.

Note: this map is NOT stratified by age (i.e. early to mid to late life) as was done in the Lancet Commission reports -- all risk factors were combined in aggregate, *independent of age*. A separate version which *does* take age into account has also been produced.

Alternate breaks into 4 buckets!

library(easypackages) libraries(c("readxl", "gtsummary", "ggmap", "ggiraph", "ggforce", "ggcorrplot", "ggthemes", "ggsignif", "ggsflabel", "ggrepel", "ggpubr", "ggsci", "glue", "gt", "janitor", "maptools", "mapview", "magrittr", "plyr", "prettyunits", "progress", "progressr", "psych", "rgeos", "rio", "rms", "Hmisc", "robustbase", "rspat", "s2", "sfheaders", "sfweight", "snakecase", "smoothr", "sp", "spatial", "spatialEco", "spatstat", "spatstat.linnet", "spatstat.model", "rpart", "spatstat.explore", "nlme", "spatstat.random", "spatstat.geom", "spatstat.data", "spdep", "sf", "spData", "abind", "summarytools", "terra", "tidycensus", "tidylog", "tidyselect", "lubridate", "forcats", "stringr", "dplyr", "purrr", "readr", "tidyr", "tibble", "ggplot2", "tidyverse", "tigris", "tmap", "vctrs", "viridis", "viridisLite", "vroom", "waldo", "wk", "stats", "graphics", "grDevices", "utils", "datasets", "methods", "base", "haven", "foreign", "survey", "srvyr", "sitrep", "questionr", "srvyr", "stringr", "readxl", "gtsummary")) acs_2022_vars <- load_variables(2022, "acs5", cache = TRUE) score_test <- read_excel("2024_socialisolation_vals_maybe_seniorhealthrankigs_mapranks.xlsx") shr_2024 <- read_excel("2024_socialisolation_vals_maybe_seniorhealthrankigs.xlsx") states_w_abbr <- read_excel("states_w_abbr.xlsx",col_names = FALSE) ##do USA variables first #############NOTE: this code is meant to confirm the findings of the 2024 Senior Report from the American Health Rankings score_test <- read_excel("2024_socialisolation_vals_maybe_seniorhealthrankigs_mapranks.xlsx") shr_2024 <- read_excel("2024_socialisolation_vals_maybe_seniorhealthrankigs.xlsx") states_w_abbr <- read_excel("states_w_abbr.xlsx",col_names = FALSE) get_acs("state", table = "B18101", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% .[,-8] %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> disability_USAt #Poverty get_acs("state", table = "B17001", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> poverty_USAt #Marital Status (i.e. Never Married, Separated, Widowed, Divorced) get_acs("state", table = "B12002", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> maritalstatus_USAt #Living Alone 65+ (why they didn't check for multicollinearity between this and marital status or whatever I don't know but w/e) get_acs("state", table = "B11007", year= 2022, survey = "acs5", summary_var = "B11007_002", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> living_alone_65_USAt #Disability get_acs("state", table = "B18101", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% .[,-8] %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> disability_USAt #Independent Living get_acs("state", table = "B18107", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> USA_independentlivingt #... #USA Disability scores (1) #confirmed #standard output disability_USAt %>% #single out only the elderly (65+) filter(label4 == "65 to 74 years:" | label4 == "75 years and over:") %>% #remove all the "grand total" rows filter(!is.na(label5)) %>% #calculate the "grand total" that we need to use by summing the disabled + not disabled across men and women transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% #calculate the percentage of those that are disable and sum them up, and then divide by total elders to get the total proportion of disabled filter(label5 == "With a disability") %>% mutate(prop = 100 * (estimate/total_elders)) %>% ##this feels REALLY stupid, but it's what they did, I swear to God transform(total_disabled = ave(.$prop, .$NAME, FUN=sum)) %>% #calculate the z_score for the total proportion of elders disabled by state as compared to national average mutate(z_score = (total_disabled - mean(total_disabled)/sd(total_disabled))) %>% dplyr::select(NAME, total_disabled, z_score, geometry) %>%left_join(., states_w_abbr %>% dplyr::rename(., state=1, abbr=2), by=c("NAME" = "state")) %>% unique(.) %>% left_join(., shr_2024 %>% filter(Measure == "Disabled"), by = c("abbr" = "State")) %>% filter(abbr != "DC" & abbr != "PR") %>% unique(.) %>% arrange(desc(-total_disabled)) %>% mutate(total_disabled = round(total_disabled, 1), rank = rank(total_disabled, ties.method = "min"), total_prop = total_disabled) %>% dplyr::select(NAME, abbr, Measure, total_prop, rank, z_score, Value, Rank, geometry) %>% st_as_sf(.) -> disability_USA_Zt #USA Poverty scores (2) #confirmed poverty_USAt %>% filter(label5 == "65 to 74 years" | label5 == "75 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(grepl("below", label3)) %>% mutate(prop = 100 * (estimate/total_elders)) %>% transform(total_poverty = ave(.$prop, .$NAME, FUN=sum)) %>% mutate(z_score = (total_poverty - mean(total_poverty)/sd(total_poverty))) %>% left_join(., states_w_abbr %>% dplyr::rename(., state=1, abbr=2), by=c("NAME" = "state")) %>% unique(.) %>% left_join(., shr_2024 %>% filter(Measure == "Living in Poverty"), by = c("abbr" = "State")) %>% filter(abbr != "DC" & abbr != "PR") %>% arrange(desc(-total_poverty)) %>% mutate(total_poverty = round(total_poverty, 1), rank = rank(total_poverty, ties.method = "min"), total_prop = total_poverty) %>% dplyr::select(NAME, abbr, Measure, total_prop, rank, z_score, Value, Rank, geometry) %>% st_as_sf(.) ->poverty_USA_Zt #USA Living Alone scores (3) #confirmed!! living_alone_65_USAt %>% filter(variable == "B11007_003") %>% mutate(prop = 100 * (estimate/summary_est)) %>% mutate(age = "65+") %>% mutate(z_score = (prop - mean(prop))/sd(prop)) %>% left_join(., states_w_abbr %>% dplyr::rename(., state=1, abbr=2), by=c("NAME" = "state")) %>% unique(.) %>% arrange(desc(-prop)) %>% filter(abbr != "DC" & abbr != "PR") %>% mutate(total_prop = round(prop, 1), rank = rank(prop, ties.method = "min")) %>% filter(abbr != "DC" & abbr != "PR") %>% left_join(., shr_2024 %>% filter(Measure == "Living Alone"), by = c("abbr" = "State")) %>% dplyr::select(NAME, abbr, Measure, total_prop, rank, z_score, Value, Rank, geometry) %>% st_as_sf(.) -> living_alone_65_USA_Zt #USA Not Married scores (4) #confirmed! maritalstatus_USAt %>% filter(label5 == "65 to 74 years" | label5 == "75 to 84 years" | label5 == "85 years and over" | label6 == "65 to 74 years" | label6 == "75 to 84 years" | label6 == "85 years and over" | label7 == "65 to 74 years" | label7 == "75 to 84 years" | label7 == "85 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label4 == "Never married:") %>% transform(not_married_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (not_married_sum/total_elders)) %>% mutate(z_score = (prop - mean(prop))/sd(prop)) %>% dplyr::select(NAME, not_married_sum, prop, z_score, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% left_join(., states_w_abbr %>% dplyr::rename(., state=1, abbr=2), by=c("NAME" = "state")) %>% unique(.) %>% filter(abbr != "DC" & abbr != "PR") %>% arrange(desc(-not_married_sum)) %>% mutate(total_prop = round(total_prop, 1), rank = rank(total_prop, ties.method = "min")) %>% left_join(., shr_2024 %>% filter(Measure == "Never Married"), by = c("abbr" = "State")) %>% dplyr::select(NAME, abbr, Measure, total_prop, rank, z_score, Value, Rank, geometry) %>% st_as_sf(.) -> USA_not_married_Zt #USA Widowed, Separated, or Divorced scores (5) #confirmed! ##minor note: this is technically an inaccuracy; it should be including "B12002_155" through "B12002_157" because "Other" for that category is functionally equivalent to "Separated" for our purposes -- their being apart owing to marital discord or not wasn't the point, 'not being together physically' was, but also who cares; source: https://www2.census.gov/programs-surveys/acs/tech_docs/subject_definitions/2023_ACSSubjectDefinitions.pdf maritalstatus_USAt %>% filter(label5 == "65 to 74 years" | label5 == "75 to 84 years" | label5 == "85 years and over" | label6 == "65 to 74 years" | label6 == "75 to 84 years" | label6 == "85 years and over" | label7 == "65 to 74 years" | label7 == "75 to 84 years" | label7 == "85 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label4 == "Widowed:" | label4 == "Divorced:" | label6 == "Separated:") %>% transform(widow_separate_divorced_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% transform(div_sep_wid_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (div_sep_wid_sum/total_elders)) %>% mutate(z_score = (prop - mean(prop))/sd(prop)) %>% dplyr::select(NAME, div_sep_wid_sum, prop, z_score, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% left_join(., states_w_abbr %>% dplyr::rename(., state=1, abbr=2), by=c("NAME" = "state")) %>% unique(.) %>% filter(abbr != "DC" & abbr != "PR") %>% arrange(desc(-div_sep_wid_sum)) %>% mutate(total_prop = round(total_prop, 1), rank = rank(total_prop, ties.method = "min")) %>% left_join(., shr_2024 %>% filter(Measure == "Divorced, Separated or Widowed"), by = c("abbr" = "State")) %>% dplyr::select(NAME, abbr, Measure, total_prop, rank, z_score, Value, Rank, geometry) %>% st_as_sf(.) -> USA_widow_separate_divorced_Zt #USA Independent Living scores (6) #confirmed! USA_independentlivingt %>% filter(label4 == "65 to 74 years:" | label4 == "75 years and over:") %>% filter(!is.na(label5)) %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label5 == "With an independent living difficulty") %>% transform(diff_indep_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% transform(diff_indep_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (diff_indep_sum/total_elders)) %>% mutate(z_score = (prop - mean(prop))/sd(prop)) %>% dplyr::select(NAME, diff_indep_sum, prop, z_score, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% left_join(., states_w_abbr %>% dplyr::rename(., state=1, abbr=2), by=c("NAME" = "state")) %>% unique(.) %>% filter(abbr != "DC" & abbr != "PR") %>% arrange(desc(-diff_indep_sum)) %>% mutate(total_prop = round(total_prop, 1), rank = rank(total_prop, ties.method = "min")) %>% left_join(., shr_2024 %>% filter(Measure == "Independent Living Difficulty"), by = c("abbr" = "State")) %>% dplyr::select(NAME, abbr, Measure, total_prop, rank, z_score, Value, Rank, geometry) %>% st_as_sf(.) -> USA_independentliving_Zt #rbind(poverty_USA_Z, living_alone_65_USA_Z, USA_widow_separate_divorced_Z, USA_not_married_Z, USA_widow_separate_divorced_Z, disability_USA_Z, USA_independentliving_Z) %>% mutate(z_score = replace(z_score, z_score >2, 2)) %>% transform(composite_score = ave(.$z_score, .$NAME, FUN=sum)) %>% unique(.) %>% mutate(z_score_percentile = percent_rank(z_score) * 100) %>% mutate(FLAG=all.equal(as.numeric(total_prop), as.numeric(Value))) -> equality_check ##old one, don't use rbind(poverty_USA_Zt, living_alone_65_USA_Zt, USA_widow_separate_divorced_Zt, USA_not_married_Zt, USA_widow_separate_divorced_Zt, disability_USA_Zt, USA_independentliving_Zt) %>% unique(.) %>% mutate(FLAG=all.equal(as.numeric(total_prop), as.numeric(Value))) %>% group_by(Measure) %>% mutate(my_z_score = (total_prop - mean(total_prop))/sd(total_prop), their_z_score = (Value - mean(Value))/sd(Value), z_FLAG=all.equal(as.numeric(my_z_score), as.numeric(their_z_score))) %>% ungroup(.) %>% transform(my_state_mean_z = ave(.$my_z_score, .$abbr, FUN=mean), their_state_mean_Z = ave(.$their_z_score, .$abbr, FUN=mean)) %>% mutate(my_possible_composite = my_state_mean_z %>% scales::rescale(., to=c(1,100)) %>% round(.), my_possible_rank = dense_rank(my_possible_composite)) %>% unique(.) %>% dplyr::select(abbr, my_possible_composite, my_possible_rank, my_state_mean_z) %>% unique(.) %>% left_join(., score_test, by=c("abbr"="State")) %>% arrange(., -desc(my_possible_composite)) %>% mutate(my_possible_rank_two = 1:50) -> composite_score_equality_check composite_score_equality_check %>% filter(abbr != "HI" & abbr != "AK") %>% st_bbox(.) ->mainland_bbox tmap_mode("view") composite_score_equality_check %>% tm_shape(., bbox = mainland_bbox) + tm_polygons(col = "my_possible_composite", palette = "blues", breaks = c(0,37, 48, 57, 70, 100), popup.vars = c("abbr", "my_possible_composite", "Value", "Rank", "Measure")) -> test_map_w_bounds #compare with https://www.americashealthrankings.org/explore/measures/isolationrisk_sr_b

############ This map is a "modern" iteration of the Senior Health Report from "America's Health Rankings" by the United Health Foundation using 2023-2018 U.S Census data to build its values. As of 02/17/2025, the Senior Health Report using this data has not been released. Some methodology data can be found here: https://www.americashealthrankings.org/about/methodology/rankings The maps themselves can be found here: https://www.americashealthrankings.org/explore/measures/isolationrisk_sr_b/CA#measure-trend-summary #################################### original code below: library(easypackages) libraries(c("readxl", "gtsummary", "ggmap", "ggiraph", "ggforce", "ggcorrplot", "ggthemes", "ggsignif", "ggsflabel", "ggrepel", "ggpubr", "ggsci", "glue", "gt", "janitor", "maptools", "mapview", "magrittr", "plyr", "prettyunits", "progress", "progressr", "psych", "rgeos", "rio", "rms", "Hmisc", "robustbase", "rspat", "s2", "sfheaders", "sfweight", "snakecase", "smoothr", "sp", "spatial", "spatialEco", "spatstat", "spatstat.linnet", "spatstat.model", "rpart", "spatstat.explore", "nlme", "spatstat.random", "spatstat.geom", "spatstat.data", "spdep", "sf", "spData", "abind", "summarytools", "terra", "tidycensus", "tidylog", "tidyselect", "lubridate", "forcats", "stringr", "dplyr", "purrr", "readr", "tidyr", "tibble", "ggplot2", "tidyverse", "tigris", "tmap", "vctrs", "viridis", "viridisLite", "vroom", "waldo", "wk", "stats", "graphics", "grDevices", "utils", "datasets", "methods", "base", "haven", "foreign", "survey", "srvyr", "sitrep", "questionr", "srvyr", "stringr", "readxl", "gtsummary")) ##note: switch to 2022 if you want to replicate the SHR map here: https://assets.americashealthrankings.org/app/uploads/rosi2024_all.pdf # acs_2022_vars <- load_variables(2022, "acs5", cache = TRUE) # ##Poverty # get_acs("county", table = "B17001", year= 2022, state="TX", survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> poverty_USA # ##Marital Status (i.e. Never Married, Separated, Widowed, Divorced) # get_acs("county", table = "B12002", state = "TX", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> maritalstatus_USA # ##Living Alone 65+ (why they didn't check for multicollinearity between this and marital status or whatever I don't know but w/e) # get_acs("county", table = "B11007", state="TX", year= 2022, survey = "acs5", summary_var = "B11007_002", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> living_alone_65_USA # ##Disability # get_acs("county", table = "B18101", state="TX", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% .[,-8] %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> disability_USA # ##Independent Living # get_acs("county", table = "B18107", state="TX", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> USA_independentliving #########uncomment the singles above to replicate acs_2023_vars <- load_variables(2023, "acs5", cache = TRUE) #Poverty get_acs("county", table = "B17001", year= 2023, state="TX", survey = "acs5", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> poverty_USA #Marital Status (i.e. Never Married, Separated, Widowed, Divorced) get_acs("county", table = "B12002", state = "TX", year= 2023, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> maritalstatus_USA #Living Alone 65+ (why they didn't check for multicollinearity between this and marital status or whatever I don't know but w/e) get_acs("county", table = "B11007", state="TX", year= 2023, survey = "acs5", summary_var = "B11007_002", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> living_alone_65_USA #Disability get_acs("county", table = "B18101", state="TX", year= 2023, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% .[,-8] %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> disability_USA #Independent Living get_acs("county", table = "B18107", state="TX", year= 2023, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> USA_independentliving ## current state z-composite draft (2/15/24, 1:22pm) #USA Disability scores (1) disability_USA %>% #single out only the elderly (65+) filter(label4 == "65 to 74 years:" | label4 == "75 years and over:") %>% #remove all the "grand total" rows filter(!is.na(label5)) %>% #calculate the "grand total" that we need to use by summing the disabled + not disabled across men and women transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% #calculate the percentage of those that are disable and sum them up, and then divide by total elders to get the total proportion of disabled filter(label5 == "With a disability") %>% mutate(prop = 100 * (estimate/total_elders)) %>% ##this feels REALLY [adjective], but it's what they did, I swear to God transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Disabled") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% unique(.) %>% arrange(desc(-total_prop)) %>% st_as_sf(.) ->disability_USA_Z #USA Poverty scores (2) poverty_USA %>% filter(label5 == "65 to 74 years" | label5 == "75 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(grepl("below", label3)) %>% mutate(prop = 100 * (estimate/total_elders)) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% arrange(desc(-total_prop)) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Poverty") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) ->poverty_USA_Z #USA Living Alone scores (3) living_alone_65_USA %>% filter(variable == "B11007_003") %>% mutate(prop = 100 * (estimate/summary_est)) %>% mutate(age = "65+") %>% unique(.) %>% arrange(desc(-prop)) %>% mutate(total_prop = prop, rank = rank(prop, ties.method = "min"), measure = "Living Alone") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) -> living_alone_65_USA_Z #USA Not Married scores (4) maritalstatus_USA %>% filter(label5 == "65 to 74 years" | label5 == "75 to 84 years" | label5 == "85 years and over" | label6 == "65 to 74 years" | label6 == "75 to 84 years" | label6 == "85 years and over" | label7 == "65 to 74 years" | label7 == "75 to 84 years" | label7 == "85 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label4 == "Never married:") %>% transform(not_married_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (not_married_sum/total_elders)) %>% dplyr::select(NAME, not_married_sum, prop, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% unique(.) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Not Married") %>% arrange(desc(-total_prop)) %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) -> USA_not_married_Z #USA Widowed, Separated, or Divorced scores (5) maritalstatus_USA %>% filter(label5 == "65 to 74 years" | label5 == "75 to 84 years" | label5 == "85 years and over" | label6 == "65 to 74 years" | label6 == "75 to 84 years" | label6 == "85 years and over" | label7 == "65 to 74 years" | label7 == "75 to 84 years" | label7 == "85 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label4 == "Widowed:" | label4 == "Divorced:" | label6 == "Separated:") %>% transform(widow_separate_divorced_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% transform(div_sep_wid_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (div_sep_wid_sum/total_elders)) %>% dplyr::select(NAME, div_sep_wid_sum, prop, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% unique(.) %>% arrange(desc(-div_sep_wid_sum)) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Widowed, Divorced, or Separated") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) -> USA_widow_separate_divorced_Z #USA Independent Living scores (6) USA_independentliving %>% filter(label4 == "65 to 74 years:" | label4 == "75 years and over:") %>% filter(!is.na(label5)) %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label5 == "With an independent living difficulty") %>% transform(diff_indep_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% transform(diff_indep_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (diff_indep_sum/total_elders)) %>% dplyr::select(NAME, diff_indep_sum, prop, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% unique(.) %>% arrange(desc(-diff_indep_sum)) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Difficulty Living Independently") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) -> USA_independentliving_Z rbind(poverty_USA_Z, living_alone_65_USA_Z, USA_widow_separate_divorced_Z, USA_not_married_Z, USA_widow_separate_divorced_Z, disability_USA_Z, USA_independentliving_Z) %>% unique(.) %>% group_by(measure) %>% mutate(my_z_score = (total_prop - mean(total_prop))/sd(total_prop)) %>% ungroup(.) %>% transform(my_state_mean_z = ave(.$my_z_score, .$NAME, FUN=mean)) %>% mutate(my_possible_composite = my_state_mean_z %>% scales::rescale(., to=c(1,100)) %>% round(.), my_possible_rank = dense_rank(my_possible_composite)) %>% unique(.) %>% dplyr::select(NAME, my_possible_composite, my_possible_rank, my_state_mean_z) %>% unique(.) -> z_composites_county z_composites_county %>% dplyr::rename(., County = 1, `Composite Score` = 2, `County Rank` = 3 ) %>% tm_shape(.) + tm_polygons(col = "Composite Score", palette = blues9, breaks = c(1,34, 39, 45, 51, 100), popup.vars = c("County", "Composite Score", "County Rank"), as.count=TRUE) + tm_layout(title = "Risk of Social Isolation by County (5-Year Composite Score) (2023-2019)") + tm_credits("Aggregate Index using U.S. Census American Community Survey Values. Thanks to United Heath Foundation for original concept.")

library(easypackages) libraries(c("readxl", "gtsummary", "ggmap", "ggiraph", "ggforce", "ggcorrplot", "ggthemes", "ggsignif", "ggsflabel", "ggrepel", "ggpubr", "ggsci", "glue", "gt", "janitor", "maptools", "mapview", "magrittr", "plyr", "prettyunits", "progress", "progressr", "psych", "rgeos", "rio", "rms", "Hmisc", "robustbase", "rspat", "s2", "sfheaders", "sfweight", "snakecase", "smoothr", "sp", "spatial", "spatialEco", "spatstat", "spatstat.linnet", "spatstat.model", "rpart", "spatstat.explore", "nlme", "spatstat.random", "spatstat.geom", "spatstat.data", "spdep", "sf", "spData", "abind", "summarytools", "terra", "tidycensus", "tidylog", "tidyselect", "lubridate", "forcats", "stringr", "dplyr", "purrr", "readr", "tidyr", "tibble", "ggplot2", "tidyverse", "tigris", "tmap", "vctrs", "viridis", "viridisLite", "vroom", "waldo", "wk", "stats", "graphics", "grDevices", "utils", "datasets", "methods", "base", "haven", "foreign", "survey", "srvyr", "sitrep", "questionr", "srvyr", "stringr", "readxl", "gtsummary")) ##note: switch to 2022 if you want to replicate the SHR map here: https://assets.americashealthrankings.org/app/uploads/rosi2024_all.pdf acs_2022_vars <- load_variables(2022, "acs5", cache = TRUE) ##Poverty get_acs("county", table = "B17001", year= 2022, state="TX", survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> poverty_USA ##Marital Status (i.e. Never Married, Separated, Widowed, Divorced) get_acs("county", table = "B12002", state = "TX", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> maritalstatus_USA ##Living Alone 65+ (why they didn't check for multicollinearity between this and marital status or whatever I don't know but w/e) get_acs("county", table = "B11007", state="TX", year= 2022, survey = "acs5", summary_var = "B11007_002", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> living_alone_65_USA ##Disability get_acs("county", table = "B18101", state="TX", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% .[,-8] %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> disability_USA ##Independent Living get_acs("county", table = "B18107", state="TX", year= 2022, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2022_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> USA_independentliving ##########uncomment the singles above to replicate # # acs_2023_vars <- load_variables(2023, "acs5", cache = TRUE) # #Poverty # get_acs("county", table = "B17001", year= 2023, state="TX", survey = "acs5", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> poverty_USA # #Marital Status (i.e. Never Married, Separated, Widowed, Divorced) # get_acs("county", table = "B12002", state = "TX", year= 2023, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> maritalstatus_USA # #Living Alone 65+ (why they didn't check for multicollinearity between this and marital status or whatever I don't know but w/e) # get_acs("county", table = "B11007", state="TX", year= 2023, survey = "acs5", summary_var = "B11007_002", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> living_alone_65_USA # #Disability # get_acs("county", table = "B18101", state="TX", year= 2023, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% .[,-8] %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> disability_USA # #Independent Living # get_acs("county", table = "B18107", state="TX", year= 2023, survey = "acs5", geometry = TRUE) %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") -> USA_independentliving ## current state z-composite draft (2/15/24, 1:22pm) #USA Disability scores (1) disability_USA %>% #single out only the elderly (65+) filter(label4 == "65 to 74 years:" | label4 == "75 years and over:") %>% #remove all the "grand total" rows filter(!is.na(label5)) %>% #calculate the "grand total" that we need to use by summing the disabled + not disabled across men and women transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% #calculate the percentage of those that are disable and sum them up, and then divide by total elders to get the total proportion of disabled filter(label5 == "With a disability") %>% mutate(prop = 100 * (estimate/total_elders)) %>% ##this feels REALLY [adjective], but it's what they did, I swear to God transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Disabled") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% unique(.) %>% arrange(desc(-total_prop)) %>% st_as_sf(.) ->disability_USA_Z #USA Poverty scores (2) poverty_USA %>% filter(label5 == "65 to 74 years" | label5 == "75 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(grepl("below", label3)) %>% mutate(prop = 100 * (estimate/total_elders)) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% arrange(desc(-total_prop)) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Poverty") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) ->poverty_USA_Z #USA Living Alone scores (3) living_alone_65_USA %>% filter(variable == "B11007_003") %>% mutate(prop = 100 * (estimate/summary_est)) %>% mutate(age = "65+") %>% unique(.) %>% arrange(desc(-prop)) %>% mutate(total_prop = prop, rank = rank(prop, ties.method = "min"), measure = "Living Alone") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) -> living_alone_65_USA_Z #USA Not Married scores (4) maritalstatus_USA %>% filter(label5 == "65 to 74 years" | label5 == "75 to 84 years" | label5 == "85 years and over" | label6 == "65 to 74 years" | label6 == "75 to 84 years" | label6 == "85 years and over" | label7 == "65 to 74 years" | label7 == "75 to 84 years" | label7 == "85 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label4 == "Never married:") %>% transform(not_married_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (not_married_sum/total_elders)) %>% dplyr::select(NAME, not_married_sum, prop, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% unique(.) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Not Married") %>% arrange(desc(-total_prop)) %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) -> USA_not_married_Z #USA Widowed, Separated, or Divorced scores (5) maritalstatus_USA %>% filter(label5 == "65 to 74 years" | label5 == "75 to 84 years" | label5 == "85 years and over" | label6 == "65 to 74 years" | label6 == "75 to 84 years" | label6 == "85 years and over" | label7 == "65 to 74 years" | label7 == "75 to 84 years" | label7 == "85 years and over") %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label4 == "Widowed:" | label4 == "Divorced:" | label6 == "Separated:") %>% transform(widow_separate_divorced_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% transform(div_sep_wid_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (div_sep_wid_sum/total_elders)) %>% dplyr::select(NAME, div_sep_wid_sum, prop, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% unique(.) %>% arrange(desc(-div_sep_wid_sum)) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Widowed, Divorced, or Separated") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) -> USA_widow_separate_divorced_Z #USA Independent Living scores (6) USA_independentliving %>% filter(label4 == "65 to 74 years:" | label4 == "75 years and over:") %>% filter(!is.na(label5)) %>% transform(total_elders = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(age = "65+") %>% filter(label5 == "With an independent living difficulty") %>% transform(diff_indep_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% transform(diff_indep_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(prop = 100 * (diff_indep_sum/total_elders)) %>% dplyr::select(NAME, diff_indep_sum, prop, geometry) %>% unique(.) %>% transform(total_prop = ave(.$prop, .$NAME, FUN=sum)) %>% unique(.) %>% arrange(desc(-diff_indep_sum)) %>% mutate(rank = rank(total_prop, ties.method = "min"), measure = "Difficulty Living Independently") %>% dplyr::select(NAME, measure, total_prop, rank, geometry) %>% st_as_sf(.) -> USA_independentliving_Z rbind(poverty_USA_Z, living_alone_65_USA_Z, USA_widow_separate_divorced_Z, USA_not_married_Z, USA_widow_separate_divorced_Z, disability_USA_Z, USA_independentliving_Z) %>% unique(.) %>% group_by(measure) %>% mutate(my_z_score = (total_prop - mean(total_prop))/sd(total_prop)) %>% ungroup(.) %>% transform(my_state_mean_z = ave(.$my_z_score, .$NAME, FUN=mean)) %>% mutate(my_possible_composite = my_state_mean_z %>% scales::rescale(., to=c(1,100)) %>% round(.), my_possible_rank = dense_rank(my_possible_composite)) %>% unique(.) %>% dplyr::select(NAME, my_possible_composite, my_possible_rank, my_state_mean_z) %>% unique(.) -> z_composites_county ########## note: use this code to build the map for States instead of counties ## ########## when computing for counties, remove the filter for Puerto Rico and DC # rbind(poverty_USA_Z, living_alone_65_USA_Z, USA_widow_separate_divorced_Z, USA_not_married_Z, USA_widow_separate_divorced_Z, disability_USA_Z, USA_independentliving_Z) %>% filter(NAME != "Puerto Rico" & NAME != "District of Columbia") %>% unique(.) %>% group_by(measure) %>% mutate(my_z_score = (total_prop - mean(total_prop))/sd(total_prop)) %>% ungroup(.) %>% dplyr::select(NAME, my_z_score, geometry) %>% unique(.) %>% transform(my_county_mean_z = ave(.$my_z_score, .$NAME, FUN=mean)) %>% dplyr::select(NAME, my_county_mean_z, geometry) %>% unique(.) %>% mutate(my_possible_composite = my_county_mean_z %>% scales::rescale(., to=c(1,100)) %>% round(.), my_possible_rank = dense_rank(my_possible_composite)) %>% unique(.) -> z_composites_state ## z_composites_state %>% dplyr::rename(., State = 1, `Composite Score` = 4, `State Rank` = 5) %>% tm_shape(.) + tm_polygons(col = "Composite Score", palette = blues9, breaks = c(1, 37, 48, 57, 70, 100), popup.vars = c("State", "Composite Score", "State Rank")) + tm_layout(title = "Risk of Social Isolation by State (5-Year Composite Score) (2023-2019)") ################ #z_composites_county %>% dplyr::rename(., County = 1, `Composite Score` = 2, `County Rank` = 3 ) %>% tm_shape(.) + tm_polygons(col = "Composite Score", palette = blues9, breaks = c(1,35, 40, 46, 52, 100), popup.vars = c("County", "Composite Score", "County Rank")) + tm_layout(title = "Risk of Social Isolation by County (5-Year Composite Score) (2023-2019)") + tm_credits("Aggregate Index using U.S. Census American Community Survey Values. Thanks to United Heath Foundation for original concept.") z_composites_county %>% dplyr::rename(., County = 1, `Composite Score` = 2, `County Rank` = 3 ) %>% tm_shape(.) + tm_polygons(col = "Composite Score", palette = blues9, breaks = c(1,34, 39, 45, 51, 100), popup.vars = c("County", "Composite Score", "County Rank"), as.count=TRUE) + tm_layout(title = "Risk of Social Isolation by County (5-Year Composite Score) (2022-2018)") + tm_credits("Aggregate Index using U.S. Census American Community Survey Values. Thanks to United Heath Foundation for original concept.")

This map should describe the proportion -- as a percentage -- of those 65 and older across the Rio Grande Valley who report living alone. The denominator for this proportion should be interpreted as the *total number* of those 65 and older in each tract. The data is grouped via Jenks Natural Breaks classification which, per the CDC, "...clusters data into groups that minimize the within-group variance and maximize the between-group variance." When hovering over any given tract, one should be able to see: a) the Census Tract b) the total population estimate for those aged 65 and older c) the margin-of-error (MOE) for the total population estimate mentioned above, d) the sum of those living alone in each tract (which was originally delineated by sex in variables "B09020_015" and "B09020_018", and e) the proportion of individuals living alone in each tract as a percentage, where bullet 'd' above is divided by bullet 'b', the total population estimate >=65y/o for that area. Thank you! .. Relevant code: get_acs("tract", table = "B09020", state = "TX", county = c("Cameron","Hidalgo", "Willacy", "Starr"), year= 2023, survey = "acs5", summary_var = "B09020_001", geometry = TRUE) -> living_alone_65_LRGV_tracts_alt living_alone_65_LRGV_tracts_alt %>% left_join(., acs_2023_vars, by = c("variable"="name")) %>% separate_wider_delim(cols="label", delim = "!!", names_sep = "", too_few = "align_start") %>% filter(label7 == "Living alone") %>% transform(tract_alone_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(proportion_of_population = round(100 * (tract_alone_sum/summary_est),2), tract_alone_sum = round(tract_alone_sum,2)) %>% dplyr::rename(., `Total 65+`=6, `Total 65+ MOE`=7, `Sum Living Alone`=18, `Proportion Living Alone`=19) %>% st_as_sf(.) %>% filter(!st_is_empty(.)) %>% tm_shape(.) + tm_polygons(col = "Proportion Living Alone", style = "jenks", fill_alpha = 0.30, palette = "reds", title = "Density of Isolation", popup.vars = c("NAME", "Sum Living Alone", "Total 65+", "Total 65+ MOE", "Proportion Living Alone")) + tm_layout(title = "Proportion of those aged 65+ who are Living Alone by Census Tract across the Rio Grande Valley") + tm_credits("Source: U.S. Census Bureau. 'Relationship by Household Type (Including Living Alone) for the Population 65 Years and Over.' American Community Survey, ACS 5-Year Estimates Detailed Tables, Table B09020, 2019-2023, https://data.census.gov/table/ACSDT1Y2022.B09020. Accessed on February 3, 2025.") + tm_basemap("OpenStreetMap")

U.S Census Table B18104, "Sex by Age by Cognitive Difficulty", for all values where variables indicate "With a cognitive difficulty" across census tracts in Cameron County. (American Community Survey, 5-year Average Estimates, 2023-2019)

U.S Census Table B18104, "Sex by Age by Cognitive Difficulty", for all values where variables indicate "With a cognitive difficulty" across census tracts. (American Community Survey, 5-year Average Estimates, 2023-2019)

A collection of four maps; from left to right, row-wise: 1) Cameron County ZCTAs, cropped to include -- and be bounded by -- the U.S. Census Tracts for Cameron County 2) Adjusted Population Estimates for fractionated U.S. Census Tracts arbitrarily divided by ZCTAs (adjusted population estimates sum to 421,110; original 2022 population estimates equal 421,854 -- 99.8% of county reflected) 3) The estimated population estimate by ZCTA based on the adjusted population estimates generated for Map 2. 4) The original population estimates as presented by the ACS-5 2022 population estimates (table B01001) .. **Be sure to adjust the map's usability by playing with the layers! ** ** the gray icon of stacked squares! ** To use: use the cursor to click and select the fragments in each map. Use the hovering red circle in each other map section to compare and contrast as desired. Thanks! .. library("easypackages") libraries(c("tidyverse", "tmap", "sf", "tidycensus", "tigris", "readxl")) cam_pop_tracts<- get_acs("tract", table = "B01001", state = "TX", county = "Cameron", year= 2022, survey = "acs5", geometry = TRUE) national_zctas<- zctas() crosswalk_2022 <- read_excel("C:/Users/daniel.pinon/Downloads/zip_tract_122022.xlsx") // note: you'll need to download the crosswalk file here: // https://www.huduser.gov/portal/datasets/usps_crosswalk.html cam_ZCTAs_22<- left_join(crosswalk_2022 %>% filter(USPS_ZIP_PREF_STATE == "TX") %>% filter(grepl("48061", .$TRACT)), national_zctas, by=c("ZIP"="ZCTA5CE20")) unique_intersected_zcta_tract_areas<- st_intersection(st_as_sf(cam_pop_tracts %>% filter(variable == "B01001_001")), st_as_sf(left_join(crosswalk_2022 %>% filter(USPS_ZIP_PREF_STATE == "TX") %>% filter(grepl("48061", .$TRACT)), national_zctas, by=c("ZIP"="ZCTA5CE20")) , crs = st_crs(cam_pop_tracts %>% filter(variable == "B01001_001")))) %>% dplyr::select(ZIP, GEOID, NAME, estimate, geometry) %>% cbind(.,st_area(.)) %>% unique(.) cam_pop_ests_across_ZCTA<- st_as_sf(cam_pop_tracts %>% filter(variable == "B01001_001")) %>% dplyr::select(GEOID, NAME, estimate, geometry) %>% cbind(.,st_area(.)) %>% st_drop_geometry(.) %>% left_join(., as.data.frame(unique_intersected_zcta_tract_areas), by=c("GEOID"="GEOID")) %>% unique(.) %>% dplyr::rename(., orig_area=4, intersected_area=8) %>% drop_na(.) %>% mutate(diff = intersected_area/orig_area, new_est = as.numeric(diff*estimate.x), new_est = round(new_est)) %>% group_by(ZIP) %>% mutate(zip_pop = sum(new_est)) %>% st_as_sf(.) tmap_mode("view") cam_zips_maps<- tm_shape(cam_pop_ests_across_ZCTA %>% dplyr::rename(., "New Population Estimate"=11, "New Zip Population Estimate"=12, "Original Population Estimate"=3) %>% .[,-1]) + tm_polygons(col=c("ZIP","New Population Estimate", "New Zip Population Estimate", "Original Population Estimate"), alpha=0.7) + tm_facets(nrow=2, sync=TRUE)

This is a set of maps that show which parts of Cameron County have the highest proportion of people aged 45 and older; these would be most relevant to secondary prevention and, to some extent, tertiary prevention. The top map, "Percentage (%) of Population age 45 and over", illustrates the proportion of people who are 65 and older as divided by the population in that tract, multiplied by 100. The bottom map, "Raw # of Population age 45 and over", shows the numerator for that proportion calculation; essentially, how many people 65 and older are *actually* living in that Census tract. These can be considered a map of the 'high risk' people with respect to ADRD in Cameron County when considering age alone. .. (you'll need 'tidyverse', 'tm', 'tidycensus', 'sf', and 'cwi' from DataHaven for ACS labeling) Code:: get_acs("tract", table = "B01001", state = "TX", county = "Cameron", year= 2022, survey = "acs5", geometry = TRUE) %>% label_acs(.) %>% separate_wider_delim(., cols = "label", delim = "!!", names_sep = "", too_few = "align_start") %>% separate_wider_delim(., cols = "variable", delim = "_", names_sep = "", too_few = "align_start") %>% mutate(variable2 = as.numeric(variable2)) %>% filter(variable2 != "1" & variable2 != "2" & variable2!= "26") %>% filter(variable2 > 6 & variable2 < 26 | variable2 > 30 & variable2 <=49) %>% transform(., tract_pop_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(tract_pop_prop = estimate/tract_pop_sum, test = if_else(variable2 > 14 & variable2 < 26 | variable2>38 & variable2<=49,"atrisk","norm")) %>% group_by(test, NAME) %>% mutate(atrisk_sum = sum(estimate)) %>% transform(., tract_pop_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(tract_pop_prop = estimate/tract_pop_sum, test = if_else(variable2 > 14 & variable2 < 26 | variable2>38 & variable2<=49,"atrisk","norm")) %>% group_by(test, NAME) %>% mutate(atrisk_sum = sum(estimate), prop_atrisk = 100 * atrisk_sum/tract_pop_sum) -> cam_pop_tracts_votingage cam_pop_tracts_votingage %>% filter(test == "atrisk") %>% st_as_sf(.) %>% dplyr::select(prop_atrisk, test, atrisk_sum, geometry, NAME) %>% drop_na(.) %>% unique(.) %>% mutate(prop_atrisk = round(prop_atrisk, 1)) %>% dplyr::rename(., "Percentage (%) of Population age 45 and over"=1, "Raw # of Population age 45 and over"=3) %>% tm_shape(.) + tm_polygons(col=c("Percentage (%) of Population age 45 and over", "Raw # of Population age 45 and over"), alpha=0.6, style="jenks", palette="Reds") + tm_layout(legend.outside = "TRUE") + tm_facets(nrow = 2, sync = TRUE)

This is a set of maps that show which parts of Cameron County have the highest proportion of people aged 65 and older; these would be most relevant to tertiary prevention of ADRD, though also tie to secondary prevention to some extent. The top map, "Percentage (%) of Population age 65 and over", illustrates the proportion of people who are 65 and older as divided by the population in that tract, multiplied by 100. The bottom map, "Raw # of Population age 65 and over", shows the numerator for that proportion calculation; essentially, how many people 65 and older are *actually* living in that Census tract. .. (you'll need 'tidyverse', 'tm', 'tidycensus', 'sf', and 'cwi' from DataHaven for ACS labeling) Code:: get_acs("tract", table = "B01001", state = "TX", county = "Cameron", year= 2022, survey = "acs5", geometry = TRUE) %>% label_acs(.) %>% separate_wider_delim(., cols = "label", delim = "!!", names_sep = "", too_few = "align_start") %>% separate_wider_delim(., cols = "variable", delim = "_", names_sep = "", too_few = "align_start") %>% mutate(variable2 = as.numeric(variable2)) %>% filter(variable2 != "1" & variable2 != "2" & variable2!= "26") %>% filter(variable2 > 6 & variable2 < 26 | variable2 > 30 & variable2 <=49) %>% transform(., tract_pop_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(tract_pop_prop = estimate/tract_pop_sum, test = if_else(variable2 > 14 & variable2 < 26 | variable2>38 & variable2<=49,"atrisk","norm")) %>% group_by(test, NAME) %>% mutate(atrisk_sum = sum(estimate)) %>% transform(., tract_pop_sum = ave(.$estimate, .$NAME, FUN=sum)) %>% mutate(tract_pop_prop = estimate/tract_pop_sum, test = if_else(variable2 > 19 & variable2 < 26 | variable2>43 & variable2<=49,"65plus","norm")) %>% group_by(test, NAME) %>% mutate(atrisk_sum = sum(estimate), prop_atrisk = 100 * atrisk_sum/tract_pop_sum) -> cam_pop_tracts_65plus cam_pop_tracts_65plus %>% filter(test == "65plus") %>% st_as_sf(.) %>% dplyr::select(prop_atrisk, test, atrisk_sum, geometry, NAME) %>% drop_na(.) %>% unique(.) %>% mutate(prop_atrisk = round(prop_atrisk, 1)) %>% dplyr::rename(., "Percentage (%) of Population age 65 and over"=1, "Raw # of Population age 65 and over"=3) %>% tm_shape(.) + tm_polygons(col=c("Percentage (%) of Population age 65 and over", "Raw # of Population age 65 and over"), alpha=0.6, style="jenks", palette="Reds") + tm_layout(legend.outside = "TRUE") + tm_facets(nrow = 2, sync = TRUE)

Includes a three-year rolling average!

ugly variable names (and ugly in general) but functional!

Ugly but functional!

base_CAM_fullDIAG_freqtable[,c(1,7:11)] %>% .[order(.$prop_17, decreasing = TRUE),] %>% gt(.) %>% data_color() where base_CAM_fullDIAG_freqtable is, more or less, <- = full_join(base_CAM_17_fullDIAG_freq, base_CAM_18_fullDIAG_freq, by=c("CCSR_desc"="CCSR_desc")) %>% dplyr::rename(., freq_17 = 2, prop_total_diag_17=3, freq_18 = 4, prop_total_diag_18=5) %>% right_join(., base_CAM_19_fullDIAG_freq, by=c("CCSR_desc"="CCSR_desc")) %>% right_join(., base_CAM_20_fullDIAG_freq, by=c("CCSR_desc"="CCSR_desc")) %>% right_join(., base_CAM_21_fullDIAG_freq, by=c("CCSR_desc"="CCSR_desc")) %>% dplyr::rename(., freq_19 = 6, prop_total_diag_19=7,freq_20 = 8, prop_total_diag_20=9, freq_21 = 10, prop_total_diag_21=11) %>% as.data.frame(.) %>% mutate(diff_18 = prop_total_diag_18-prop_total_diag_17,diff_19 = prop_total_diag_19-prop_total_diag_18, diff_20 = prop_total_diag_20-prop_total_diag_19, diff_21 = prop_total_diag_21-prop_total_diag_20, trend_18 = sign(prop_total_diag_18-prop_total_diag_17),trend_19 = sign(prop_total_diag_19-prop_total_diag_18), trend_20 = sign(prop_total_diag_20-prop_total_diag_19), trend_21 = sign(prop_total_diag_21-prop_total_diag_20)) %>% .[,c(1,2,4,6,8,10)] %>% mutate(prop_17 = round(100 * (freq_17/{freq_17 %>% na.omit(.) %>% sum(.)}),2), prop_18 = round(100 * (freq_18/{freq_18 %>% na.omit(.) %>% sum(.)}),2), prop_19 = round(100 * (freq_19/{freq_19 %>% na.omit(.) %>% sum(.)}),2), prop_20 = round(100 * (freq_20/{freq_20 %>% na.omit(.) %>% sum(.)}),2), prop_21 = round(100 * (freq_21/{freq_21 %>% na.omit(.) %>% sum(.)}),2), diff_18 = prop_18-prop_17, diff_19 = prop_19-prop_18, diff_20 = prop_19-prop_18, diff_21 = prop_20-prop_19, trend18 = sign(diff_18), trend19 = sign(diff_19), trend20 = sign(diff_20), trend21 = sign(diff_21))

HIPAA compliant! (big zips/zctas only)

base_CAM_fullDIAG_freqtable[,c(1,7:11)] %>% .[order(.$prop_17, decreasing = TRUE),] %>% gt(.) %>% data_color() where base_CAM_fullDIAG_freqtable is, more or less, <- = full_join(base_CAM_17_fullDIAG_freq, base_CAM_18_fullDIAG_freq, by=c("CCSR_desc"="CCSR_desc")) %>% dplyr::rename(., freq_17 = 2, prop_total_diag_17=3, freq_18 = 4, prop_total_diag_18=5) %>% right_join(., base_CAM_19_fullDIAG_freq, by=c("CCSR_desc"="CCSR_desc")) %>% right_join(., base_CAM_20_fullDIAG_freq, by=c("CCSR_desc"="CCSR_desc")) %>% right_join(., base_CAM_21_fullDIAG_freq, by=c("CCSR_desc"="CCSR_desc")) %>% dplyr::rename(., freq_19 = 6, prop_total_diag_19=7,freq_20 = 8, prop_total_diag_20=9, freq_21 = 10, prop_total_diag_21=11) %>% as.data.frame(.) %>% mutate(diff_18 = prop_total_diag_18-prop_total_diag_17,diff_19 = prop_total_diag_19-prop_total_diag_18, diff_20 = prop_total_diag_20-prop_total_diag_19, diff_21 = prop_total_diag_21-prop_total_diag_20, trend_18 = sign(prop_total_diag_18-prop_total_diag_17),trend_19 = sign(prop_total_diag_19-prop_total_diag_18), trend_20 = sign(prop_total_diag_20-prop_total_diag_19), trend_21 = sign(prop_total_diag_21-prop_total_diag_20)) %>% .[,c(1,2,4,6,8,10)] %>% mutate(prop_17 = round(100 * (freq_17/{freq_17 %>% na.omit(.) %>% sum(.)}),2), prop_18 = round(100 * (freq_18/{freq_18 %>% na.omit(.) %>% sum(.)}),2), prop_19 = round(100 * (freq_19/{freq_19 %>% na.omit(.) %>% sum(.)}),2), prop_20 = round(100 * (freq_20/{freq_20 %>% na.omit(.) %>% sum(.)}),2), prop_21 = round(100 * (freq_21/{freq_21 %>% na.omit(.) %>% sum(.)}),2), diff_18 = prop_18-prop_17, diff_19 = prop_19-prop_18, diff_20 = prop_19-prop_18, diff_21 = prop_20-prop_19, trend18 = sign(diff_18), trend19 = sign(diff_19), trend20 = sign(diff_20), trend21 = sign(diff_21))

I don't know if the viewers can actually see this? If you can, hello! This is State Data, but it's agglomerated, so there's no HIPAA violations here to my knowledge -- nor is there any health information attributable to any one particular individual or group of individuals (and the margin of error associated with the Census values essentially guarantees there's a 10% bigger window than was reported that encompasses the true population value, which, in some cases, falls BELOW zero -- so, yeah)

lrgv_blocks_CAM_localmoran[lrgv_blocks_CAM_localmoran$Pr.z....E.Ii...Sim <=0.05,] %>% tm_shape(.) + tm_fill(col = "mean", alpha = 0.2) + tm_borders() + tm_layout(main.title = "Local Moran's I Population Estimate Clusters and Outliers for p-values <= 0.05") + tm_minimap()

**NOTE: Cameron County has a transparency applied to increase see-through visibility -- please click on the block you're interested in for a view of the percentage of households. ** You are strongly encouraged to share the stacked squares on the leftmost aspect of the map/screen to enable view of cities and places. ... This is a map of households in the Lower Rio Grande Valley (Hidalgo County, Willacy County, and Cameron County) that possess at least one computer, including: desktops, laptops, smartphones, tablets, and/or any other portable computer, wireless or otherwise. The color palette has been inverted to highlight areas with the *least* access, such that the darkest red regions are the ones with the lowest proportion of computer/technology access. There are limitations to census data, including undercount or omission of populations that don't respond to U.S Census surveys (such as undocumented people and/or people with low confidence in government.) ... R Code used to generate the map: lrgv_super[[24]] %>% filter(!st_is_empty(.) & variable == "B28010_002") %>% mutate(pct = 100 * (estimate/summary_est)) %>% tm_shape(.) + tm_fill(col = "pct", style = "jenks", palette = "-Reds", title = "Percent of Households", alpha = 0.2) + tm_borders() + tm_layout(main.title = "LRGV Distribution of Households w/ at least ONE type of Computing Device", main.title.fontface = "bold.italic") + tm_minimap()