── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.3 ✔ readr 2.1.4
✔ forcats 1.0.0 ✔ stringr 1.5.0
✔ ggplot2 3.4.3 ✔ tibble 3.2.1
✔ lubridate 1.9.2 ✔ tidyr 1.3.0
✔ purrr 1.0.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Attaching package: 'janitor'
The following objects are masked from 'package:stats':
chisq.test, fisher.test
Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
remotes::install_github('r-tmap/tmap')
Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE
To enable caching of data, set `options(tigris_use_cache = TRUE)`
in your R script or .Rprofile.
Attaching package: 'zoo'
The following objects are masked from 'package:base':
as.Date, as.Date.numeric
Nonparametric Kernel Methods for Mixed Datatypes (version 0.60-17)
[vignette("np_faq",package="np") provides answers to frequently asked questions]
[vignette("np",package="np") an overview]
[vignette("entropy_np",package="np") an overview of entropy-based methods]
Loading required package: nlme
Attaching package: 'nlme'
The following object is masked from 'package:dplyr':
collapse
This is mgcv 1.9-0. For overview type 'help("mgcv-package")'.
Loading required package: spDataSpatial Analysis
1. Subset Data
# Subset Data
## Import COD dataset
sample <- read_dta("C:/Users/User/OneDrive/Research/CDC Race COD/finaldata_240623.dta")
## Subset
sample <- sample %>%
filter(year == 2019) %>%
mutate(pr_65plus = p_65plus/p_total)
## Prepare county sf data
county <- tigris::counties() %>%
rename(county_fips = GEOID) %>%
mutate(county_fips = as.numeric(county_fips)) %>%
filter(as.numeric(STATEFP) < 60)Retrieving data for the year 2021
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%## Join data
sample_sf <- left_join(county, sample, by = "county_fips") %>%
tigris::shift_geometry()
## Prepare state overlay
state_overlay <- states(cb = TRUE, resolution = "20m") %>%
filter(as.numeric(STATEFP) < 60) %>%
tigris::shift_geometry()Retrieving data for the year 2021
|
| | 0%
|
|============= | 19%
|
|============================ | 40%
|
|================================== | 49%
|
|=========================================== | 62%
|
|================================================ | 69%
|
|=============================================================== | 90%
|
|======================================================================| 100%# important: DV cannot have NA
sample_sf <- sample_sf[!is.na(sample_sf$di_index1), ]
names(sample_sf) [1] "STATEFP" "COUNTYFP"
[3] "COUNTYNS" "county_fips"
[5] "NAME.x" "NAMELSAD"
[7] "LSAD" "CLASSFP"
[9] "MTFCC" "CSAFP"
[11] "CBSAFP" "METDIVFP"
[13] "FUNCSTAT" "ALAND"
[15] "AWATER" "INTPTLAT"
[17] "INTPTLON" "year"
[19] "abbr" "state_name"
[21] "county_name" "ajr_total_0to14_sui"
[23] "ajr_total_14to64_sui" "ajr_total_65_sui"
[25] "ajr_total_sui" "ajr_total_0to14_drug"
[27] "ajr_total_14to64_drug" "ajr_total_65_drug"
[29] "ajr_total_drug" "ajr_total_0to14_alcohol"
[31] "ajr_total_14to64_alcohol" "ajr_total_65_alcohol"
[33] "ajr_total_alcohol" "ajr_total_0to14_heart_disease"
[35] "ajr_total_14to64_heart_disease" "ajr_total_65_heart_disease"
[37] "ajr_total_heart_disease" "ajr_total_0to14_cancer"
[39] "ajr_total_14to64_cancer" "ajr_total_65_cancer"
[41] "ajr_total_cancer" "ajr_total_0to14_stroke"
[43] "ajr_total_14to64_stroke" "ajr_total_65_stroke"
[45] "ajr_total_stroke" "ajr_total_0to14_diabetes"
[47] "ajr_total_14to64_diabetes" "ajr_total_65_diabetes"
[49] "ajr_total_diabetes" "ajr_total_14to64_alzheimers"
[51] "ajr_total_65_alzheimers" "ajr_total_alzheimers"
[53] "ajr_total_0to14_clr" "ajr_total_14to64_clr"
[55] "ajr_total_65_clr" "ajr_total_clr"
[57] "ajr_total_0to14_accidents" "ajr_total_14to64_accidents"
[59] "ajr_total_65_accidents" "ajr_total_accidents"
[61] "ajr_total_0to14_chronic_liver" "ajr_total_14to64_chronic_liver"
[63] "ajr_total_65_chronic_liver" "ajr_total_chronic_liver"
[65] "ajr_total_0to14_nephrotic" "ajr_total_14to64_nephrotic"
[67] "ajr_total_65_nephrotic" "ajr_total_nephrotic"
[69] "ajr_total_0to14_other" "ajr_total_14to64_other"
[71] "ajr_total_65_other" "ajr_total_other"
[73] "ajr_total_0to14_dod" "ajr_total_14to64_dod"
[75] "ajr_total_65_dod" "ajr_total_dod"
[77] "ajr_whi_0to14_sui" "ajr_whi_14to64_sui"
[79] "ajr_whi_65_sui" "ajr_whi_sui"
[81] "ajr_whi_0to14_drug" "ajr_whi_14to64_drug"
[83] "ajr_whi_65_drug" "ajr_whi_drug"
[85] "ajr_whi_0to14_alcohol" "ajr_whi_14to64_alcohol"
[87] "ajr_whi_65_alcohol" "ajr_whi_alcohol"
[89] "ajr_whi_0to14_heart_disease" "ajr_whi_14to64_heart_disease"
[91] "ajr_whi_65_heart_disease" "ajr_whi_heart_disease"
[93] "ajr_whi_0to14_cancer" "ajr_whi_14to64_cancer"
[95] "ajr_whi_65_cancer" "ajr_whi_cancer"
[97] "ajr_whi_0to14_stroke" "ajr_whi_14to64_stroke"
[99] "ajr_whi_65_stroke" "ajr_whi_stroke"
[101] "ajr_whi_0to14_diabetes" "ajr_whi_14to64_diabetes"
[103] "ajr_whi_65_diabetes" "ajr_whi_diabetes"
[105] "ajr_whi_14to64_alzheimers" "ajr_whi_65_alzheimers"
[107] "ajr_whi_alzheimers" "ajr_whi_0to14_clr"
[109] "ajr_whi_14to64_clr" "ajr_whi_65_clr"
[111] "ajr_whi_clr" "ajr_whi_0to14_accidents"
[113] "ajr_whi_14to64_accidents" "ajr_whi_65_accidents"
[115] "ajr_whi_accidents" "ajr_whi_0to14_chronic_liver"
[117] "ajr_whi_14to64_chronic_liver" "ajr_whi_65_chronic_liver"
[119] "ajr_whi_chronic_liver" "ajr_whi_0to14_nephrotic"
[121] "ajr_whi_14to64_nephrotic" "ajr_whi_65_nephrotic"
[123] "ajr_whi_nephrotic" "ajr_whi_0to14_other"
[125] "ajr_whi_14to64_other" "ajr_whi_65_other"
[127] "ajr_whi_other" "ajr_whi_0to14_dod"
[129] "ajr_whi_14to64_dod" "ajr_whi_65_dod"
[131] "ajr_whi_dod" "ajr_blk_0to14_sui"
[133] "ajr_blk_14to64_sui" "ajr_blk_65_sui"
[135] "ajr_blk_sui" "ajr_blk_0to14_drug"
[137] "ajr_blk_14to64_drug" "ajr_blk_65_drug"
[139] "ajr_blk_drug" "ajr_blk_0to14_alcohol"
[141] "ajr_blk_14to64_alcohol" "ajr_blk_65_alcohol"
[143] "ajr_blk_alcohol" "ajr_blk_0to14_heart_disease"
[145] "ajr_blk_14to64_heart_disease" "ajr_blk_65_heart_disease"
[147] "ajr_blk_heart_disease" "ajr_blk_0to14_cancer"
[149] "ajr_blk_14to64_cancer" "ajr_blk_65_cancer"
[151] "ajr_blk_cancer" "ajr_blk_0to14_stroke"
[153] "ajr_blk_14to64_stroke" "ajr_blk_65_stroke"
[155] "ajr_blk_stroke" "ajr_blk_0to14_diabetes"
[157] "ajr_blk_14to64_diabetes" "ajr_blk_65_diabetes"
[159] "ajr_blk_diabetes" "ajr_blk_14to64_alzheimers"
[161] "ajr_blk_65_alzheimers" "ajr_blk_alzheimers"
[163] "ajr_blk_0to14_clr" "ajr_blk_14to64_clr"
[165] "ajr_blk_65_clr" "ajr_blk_clr"
[167] "ajr_blk_0to14_accidents" "ajr_blk_14to64_accidents"
[169] "ajr_blk_65_accidents" "ajr_blk_accidents"
[171] "ajr_blk_0to14_chronic_liver" "ajr_blk_14to64_chronic_liver"
[173] "ajr_blk_65_chronic_liver" "ajr_blk_chronic_liver"
[175] "ajr_blk_0to14_nephrotic" "ajr_blk_14to64_nephrotic"
[177] "ajr_blk_65_nephrotic" "ajr_blk_nephrotic"
[179] "ajr_blk_0to14_other" "ajr_blk_14to64_other"
[181] "ajr_blk_65_other" "ajr_blk_other"
[183] "ajr_blk_0to14_dod" "ajr_blk_14to64_dod"
[185] "ajr_blk_65_dod" "ajr_blk_dod"
[187] "ajr_his_0to14_sui" "ajr_his_14to64_sui"
[189] "ajr_his_65_sui" "ajr_his_sui"
[191] "ajr_his_0to14_drug" "ajr_his_14to64_drug"
[193] "ajr_his_65_drug" "ajr_his_drug"
[195] "ajr_his_0to14_alcohol" "ajr_his_14to64_alcohol"
[197] "ajr_his_65_alcohol" "ajr_his_alcohol"
[199] "ajr_his_0to14_heart_disease" "ajr_his_14to64_heart_disease"
[201] "ajr_his_65_heart_disease" "ajr_his_heart_disease"
[203] "ajr_his_0to14_cancer" "ajr_his_14to64_cancer"
[205] "ajr_his_65_cancer" "ajr_his_cancer"
[207] "ajr_his_0to14_stroke" "ajr_his_14to64_stroke"
[209] "ajr_his_65_stroke" "ajr_his_stroke"
[211] "ajr_his_0to14_diabetes" "ajr_his_14to64_diabetes"
[213] "ajr_his_65_diabetes" "ajr_his_diabetes"
[215] "ajr_his_14to64_alzheimers" "ajr_his_65_alzheimers"
[217] "ajr_his_alzheimers" "ajr_his_0to14_clr"
[219] "ajr_his_14to64_clr" "ajr_his_65_clr"
[221] "ajr_his_clr" "ajr_his_0to14_accidents"
[223] "ajr_his_14to64_accidents" "ajr_his_65_accidents"
[225] "ajr_his_accidents" "ajr_his_0to14_chronic_liver"
[227] "ajr_his_14to64_chronic_liver" "ajr_his_65_chronic_liver"
[229] "ajr_his_chronic_liver" "ajr_his_0to14_nephrotic"
[231] "ajr_his_14to64_nephrotic" "ajr_his_65_nephrotic"
[233] "ajr_his_nephrotic" "ajr_his_0to14_other"
[235] "ajr_his_14to64_other" "ajr_his_65_other"
[237] "ajr_his_other" "ajr_his_0to14_dod"
[239] "ajr_his_14to64_dod" "ajr_his_65_dod"
[241] "ajr_his_dod" "ct_total_sui"
[243] "ct_total_alcohol" "ct_total_drug"
[245] "ct_total_heart_disease" "ct_total_cancer"
[247] "ct_total_stroke" "ct_total_diabetes"
[249] "ct_total_alzheimers" "ct_total_clr"
[251] "ct_total_accidents" "ct_total_chronic_liver"
[253] "ct_total_nephrotic" "ct_total_other"
[255] "ct_whi_sui" "ct_whi_alcohol"
[257] "ct_whi_drug" "ct_whi_heart_disease"
[259] "ct_whi_cancer" "ct_whi_stroke"
[261] "ct_whi_diabetes" "ct_whi_alzheimers"
[263] "ct_whi_clr" "ct_whi_accidents"
[265] "ct_whi_chronic_liver" "ct_whi_nephrotic"
[267] "ct_whi_other" "ct_blk_sui"
[269] "ct_blk_alcohol" "ct_blk_drug"
[271] "ct_blk_heart_disease" "ct_blk_cancer"
[273] "ct_blk_stroke" "ct_blk_diabetes"
[275] "ct_blk_alzheimers" "ct_blk_clr"
[277] "ct_blk_accidents" "ct_blk_chronic_liver"
[279] "ct_blk_nephrotic" "ct_blk_other"
[281] "ct_his_sui" "ct_his_alcohol"
[283] "ct_his_drug" "ct_his_heart_disease"
[285] "ct_his_cancer" "ct_his_stroke"
[287] "ct_his_diabetes" "ct_his_alzheimers"
[289] "ct_his_clr" "ct_his_accidents"
[291] "ct_his_chronic_liver" "ct_his_nephrotic"
[293] "ct_his_other" "n_mari_nev_over65"
[295] "n_mari_sepa_over65" "n_mari_wid_over65"
[297] "n_mari_div_over65" "n_mari_sum_over65"
[299] "n_pov_over65" "n_alone_over65"
[301] "n_disab_over65" "n_livdiff_over65"
[303] "n_heardiff_over65" "n_visidiff_over65"
[305] "n_cogdiff_over65" "n_ambdiff_over65"
[307] "n_slfcardiff_over65" "n_mari_nev_under65"
[309] "n_mari_sepa_under65" "n_mari_wid_under65"
[311] "n_mari_div_under65" "n_mari_sum_under65"
[313] "n_pov_under65" "n_alone_under65"
[315] "n_disab_under65" "n_livdiff_under65"
[317] "n_heardiff_under65" "n_visidiff_under65"
[319] "n_cogdiff_under65" "n_ambdiff_under65"
[321] "n_slfcardiff_under65" "n_mari_nev_over65_std"
[323] "n_mari_nev_under65_std" "n_mari_sum_over65_std"
[325] "n_mari_sum_under65_std" "n_alone_over65_std"
[327] "n_alone_under65_std" "n_pov_over65_std"
[329] "n_pov_under65_std" "n_disab_over65_std"
[331] "n_disab_under65_std" "n_livdiff_over65_std"
[333] "n_livdiff_under65_std" "n_heardiff_over65_std"
[335] "n_heardiff_under65_std" "n_visidiff_over65_std"
[337] "n_visidiff_under65_std" "n_cogdiff_over65_std"
[339] "n_cogdiff_under65_std" "n_ambdiff_over65_std"
[341] "n_ambdiff_under65_std" "n_slfcardiff_over65_std"
[343] "n_slfcardiff_under65_std" "n_isolation_under65"
[345] "n_isolation_over65" "dissi_asi_whi"
[347] "dissi_blk_whi" "dissi_his_whi"
[349] "dissi_nonwhi_whi" "iso_asi"
[351] "iso_blk" "iso_his"
[353] "iso_nonwhi" "hc_provider"
[355] "hc_nursing" "hc_mental"
[357] "countypop" "hc_rt_provider"
[359] "hc_rt_nursing" "hc_rt_mental"
[361] "resi_own_total" "resi_own_65over"
[363] "resi_time" "medage"
[365] "p_total" "p_m"
[367] "p_f" "p_b"
[369] "p_w" "p_h"
[371] "p_65plus" "pop_density"
[373] "foreborn" "hhincome"
[375] "percapincome" "p_educ_overba"
[377] "p_pov" "NAME.y"
[379] "pop_tot" "pop_blk"
[381] "pop_asi" "pop_whi"
[383] "pop_his" "pop_nonwhi"
[385] "pr_whi" "pr_blk"
[387] "pr_his" "pr_asi"
[389] "pr_nonwhi" "di_educ"
[391] "di_fehh" "di_hhinc"
[393] "di_pov" "di_pubassi"
[395] "di_unemp" "di_index1"
[397] "rx_mic" "rx_emp"
[399] "rx_pov" "rx_hown"
[401] "rx_unemp" "rx_educ"
[403] "rx_pinc" "rx_mean"
[405] "bdscore_pr_whi" "bdscore_pr_blk"
[407] "bdscore_pr_his" "bdscore_pr_asi"
[409] "bdscore_pr_nonwhi" "bdscore_di_index1"
[411] "bdscore_rx_mic" "bdscore_rx_unemp"
[413] "bdscore_rx_pov" "bdscore_rx_hown"
[415] "bdscore_rx_educ" "bdscore_rx_pinc"
[417] "bdscore_rx_mean" "urbancode1"
[419] "urbancode2" "ajr_total_0to14_unin_traf"
[421] "ajr_total_14to64_unin_traf" "ajr_total_65_unin_traf"
[423] "ajr_total_unin_traf" "ajr_total_0to14_unin_falls"
[425] "ajr_total_14to64_unin_falls" "ajr_total_65_unin_falls"
[427] "ajr_total_unin_falls" "ajr_total_0to14_unin_gun"
[429] "ajr_total_14to64_unin_gun" "ajr_total_65_unin_gun"
[431] "ajr_total_unin_gun" "ajr_total_0to14_unin_poisoning"
[433] "ajr_total_14to64_unin_poisoning" "ajr_total_65_unin_poisoning"
[435] "ajr_total_unin_poisoning" "ajr_total_0to14_unin_oth_inj"
[437] "ajr_total_14to64_unin_oth_inj" "ajr_total_65_unin_oth_inj"
[439] "ajr_total_unin_oth_inj" "ajr_total_0to14_sui_by_gun"
[441] "ajr_total_14to64_sui_by_gun" "ajr_total_65_sui_by_gun"
[443] "ajr_total_sui_by_gun" "ajr_total_0to14_sui_not_by_gun"
[445] "ajr_total_14to64_sui_not_by_gun" "ajr_total_65_sui_not_by_gun"
[447] "ajr_total_sui_not_by_gun" "ajr_total_0to14_homi_by_gun"
[449] "ajr_total_14to64_homi_by_gun" "ajr_total_65_homi_by_gun"
[451] "ajr_total_homi_by_gun" "ajr_total_0to14_homi_not_by_gun"
[453] "ajr_total_14to64_homi_not_by_gun" "ajr_total_65_homi_not_by_gun"
[455] "ajr_total_homi_not_by_gun" "ajr_whi_0to14_unin_traf"
[457] "ajr_whi_14to64_unin_traf" "ajr_whi_65_unin_traf"
[459] "ajr_whi_unin_traf" "ajr_whi_0to14_unin_falls"
[461] "ajr_whi_14to64_unin_falls" "ajr_whi_65_unin_falls"
[463] "ajr_whi_unin_falls" "ajr_whi_0to14_unin_gun"
[465] "ajr_whi_14to64_unin_gun" "ajr_whi_65_unin_gun"
[467] "ajr_whi_unin_gun" "ajr_whi_0to14_unin_poisoning"
[469] "ajr_whi_14to64_unin_poisoning" "ajr_whi_65_unin_poisoning"
[471] "ajr_whi_unin_poisoning" "ajr_whi_0to14_unin_oth_inj"
[473] "ajr_whi_14to64_unin_oth_inj" "ajr_whi_65_unin_oth_inj"
[475] "ajr_whi_unin_oth_inj" "ajr_whi_0to14_sui_by_gun"
[477] "ajr_whi_14to64_sui_by_gun" "ajr_whi_65_sui_by_gun"
[479] "ajr_whi_sui_by_gun" "ajr_whi_0to14_sui_not_by_gun"
[481] "ajr_whi_14to64_sui_not_by_gun" "ajr_whi_65_sui_not_by_gun"
[483] "ajr_whi_sui_not_by_gun" "ajr_whi_0to14_homi_by_gun"
[485] "ajr_whi_14to64_homi_by_gun" "ajr_whi_65_homi_by_gun"
[487] "ajr_whi_homi_by_gun" "ajr_whi_0to14_homi_not_by_gun"
[489] "ajr_whi_14to64_homi_not_by_gun" "ajr_whi_65_homi_not_by_gun"
[491] "ajr_whi_homi_not_by_gun" "ajr_blk_0to14_unin_traf"
[493] "ajr_blk_14to64_unin_traf" "ajr_blk_65_unin_traf"
[495] "ajr_blk_unin_traf" "ajr_blk_0to14_unin_falls"
[497] "ajr_blk_14to64_unin_falls" "ajr_blk_65_unin_falls"
[499] "ajr_blk_unin_falls" "ajr_blk_0to14_unin_gun"
[501] "ajr_blk_14to64_unin_gun" "ajr_blk_65_unin_gun"
[503] "ajr_blk_unin_gun" "ajr_blk_0to14_unin_poisoning"
[505] "ajr_blk_14to64_unin_poisoning" "ajr_blk_65_unin_poisoning"
[507] "ajr_blk_unin_poisoning" "ajr_blk_0to14_unin_oth_inj"
[509] "ajr_blk_14to64_unin_oth_inj" "ajr_blk_65_unin_oth_inj"
[511] "ajr_blk_unin_oth_inj" "ajr_blk_0to14_sui_by_gun"
[513] "ajr_blk_14to64_sui_by_gun" "ajr_blk_65_sui_by_gun"
[515] "ajr_blk_sui_by_gun" "ajr_blk_0to14_sui_not_by_gun"
[517] "ajr_blk_14to64_sui_not_by_gun" "ajr_blk_65_sui_not_by_gun"
[519] "ajr_blk_sui_not_by_gun" "ajr_blk_0to14_homi_by_gun"
[521] "ajr_blk_14to64_homi_by_gun" "ajr_blk_65_homi_by_gun"
[523] "ajr_blk_homi_by_gun" "ajr_blk_0to14_homi_not_by_gun"
[525] "ajr_blk_14to64_homi_not_by_gun" "ajr_blk_65_homi_not_by_gun"
[527] "ajr_blk_homi_not_by_gun" "ajr_his_0to14_unin_traf"
[529] "ajr_his_14to64_unin_traf" "ajr_his_65_unin_traf"
[531] "ajr_his_unin_traf" "ajr_his_0to14_unin_falls"
[533] "ajr_his_14to64_unin_falls" "ajr_his_65_unin_falls"
[535] "ajr_his_unin_falls" "ajr_his_0to14_unin_gun"
[537] "ajr_his_14to64_unin_gun" "ajr_his_65_unin_gun"
[539] "ajr_his_unin_gun" "ajr_his_0to14_unin_poisoning"
[541] "ajr_his_14to64_unin_poisoning" "ajr_his_65_unin_poisoning"
[543] "ajr_his_unin_poisoning" "ajr_his_0to14_unin_oth_inj"
[545] "ajr_his_14to64_unin_oth_inj" "ajr_his_65_unin_oth_inj"
[547] "ajr_his_unin_oth_inj" "ajr_his_0to14_sui_by_gun"
[549] "ajr_his_14to64_sui_by_gun" "ajr_his_65_sui_by_gun"
[551] "ajr_his_sui_by_gun" "ajr_his_0to14_sui_not_by_gun"
[553] "ajr_his_14to64_sui_not_by_gun" "ajr_his_65_sui_not_by_gun"
[555] "ajr_his_sui_not_by_gun" "ajr_his_0to14_homi_by_gun"
[557] "ajr_his_14to64_homi_by_gun" "ajr_his_65_homi_by_gun"
[559] "ajr_his_homi_by_gun" "ajr_his_0to14_homi_not_by_gun"
[561] "ajr_his_14to64_homi_not_by_gun" "ajr_his_65_homi_not_by_gun"
[563] "ajr_his_homi_not_by_gun" "ct_total_homi_by_gun"
[565] "ct_total_homi_not_by_gun" "ct_total_sui_by_gun"
[567] "ct_total_sui_not_by_gun" "ct_total_unin_falls"
[569] "ct_total_unin_gun" "ct_total_unin_other"
[571] "ct_total_unin_poisoning" "ct_total_unin_traf"
[573] "ct_whi_homi_by_gun" "ct_whi_homi_not_by_gun"
[575] "ct_whi_sui_by_gun" "ct_whi_sui_not_by_gun"
[577] "ct_whi_unin_falls" "ct_whi_unin_gun"
[579] "ct_whi_unin_other" "ct_whi_unin_poisoning"
[581] "ct_whi_unin_traf" "ct_blk_homi_by_gun"
[583] "ct_blk_homi_not_by_gun" "ct_blk_sui_by_gun"
[585] "ct_blk_sui_not_by_gun" "ct_blk_unin_falls"
[587] "ct_blk_unin_gun" "ct_blk_unin_other"
[589] "ct_blk_unin_poisoning" "ct_blk_unin_traf"
[591] "ct_his_homi_by_gun" "ct_his_homi_not_by_gun"
[593] "ct_his_sui_by_gun" "ct_his_sui_not_by_gun"
[595] "ct_his_unin_falls" "ct_his_unin_gun"
[597] "ct_his_unin_other" "ct_his_unin_poisoning"
[599] "ct_his_unin_traf" "pr_65plus"
[601] "geometry" 2. Basic Map
sample_sf %>%
tm_shape() +
tm_polygons(
col = "di_index1",
palette = "-RdYlBu",
style = "fisher",
n = 5,
border.col = NA,
border.alpha = 0
) +
tm_shape(state_overlay) +
tm_borders(col = "black", lwd = 1.0)Variable(s) "di_index1" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
3. Spatial Neighborhoods
3-1. Contiguity-based neighbors - Queen: where all polygons that share at least one vertex are considered neighbors
# Create a contiguity-based neighbor list
neighbors_queen <- poly2nb(sample_sf, queen = TRUE)
summary(neighbors_queen)Neighbour list object:
Number of regions: 3136
Number of nonzero links: 18504
Percentage nonzero weights: 0.188154
Average number of links: 5.90051
3 regions with no links:
3139 3141 3143
7 disjoint connected subgraphs
Link number distribution:
0 1 2 3 4 5 6 7 8 9 10 11 13 14
3 18 35 69 265 674 1104 675 229 50 11 1 1 1
18 least connected regions:
294 309 367 391 680 1074 1140 1152 1498 1778 2131 2136 2351 3132 3134 3137 3140 3142 with 1 link
1 most connected region:
2697 with 14 links# Visualization
sample_coords <- st_coordinates(st_centroid(sample_sf))Warning: st_centroid assumes attributes are constant over geometriesplot(st_geometry(sample_sf), col = "lightgrey", main = "Map with Contiguity-Based Neighbors")
plot(neighbors_queen,
coords = sample_coords,
add = TRUE,
col = "blue",
points = FALSE)
3-2. Proximity-based neighbors - Neighbors might be defined as those that fall within a given distance threshold (e.g. all features within 2km of a given feature)
# Define distance threshold
distance_threshold <- 50000 # 50 km
# Create a distance-based neighbor list
neighbors_dist <- dnearneigh(st_centroid(sample_sf), d1 = 0, d2 = distance_threshold)Warning: st_centroid assumes attributes are constant over geometriessummary(neighbors_dist)Neighbour list object:
Number of regions: 3136
Number of nonzero links: 12424
Percentage nonzero weights: 0.1263308
Average number of links: 3.961735
446 regions with no links:
3 18 22 30 31 32 33 35 62 80 84 86 92 94 100 104 117 118 147 148 171
176 195 196 200 201 205 206 212 214 248 252 254 257 264 265 271 275 277
281 282 301 314 316 320 346 360 364 383 385 386 387 394 406 413 423 427
433 437 440 441 445 457 466 470 476 477 479 486 497 502 506 514 517 518
523 525 527 551 574 577 592 603 626 630 687 691 696 703 706 708 740 744
754 755 761 766 776 797 801 814 818 821 822 827 832 833 842 854 863 874
878 879 880 888 890 900 910 912 914 920 929 931 937 942 954 971 995 996
997 1001 1005 1008 1013 1014 1036 1037 1042 1057 1066 1067 1068 1069
1078 1083 1111 1126 1139 1142 1146 1149 1158 1160 1169 1185 1196 1199
1204 1208 1247 1255 1261 1262 1272 1283 1290 1291 1297 1325 1328 1337
1340 1345 1349 1350 1352 1357 1358 1364 1365 1368 1373 1374 1381 1401
1409 1425 1432 1433 1434 1436 1455 1457 1463 1464 1469 1470 1476 1478
1493 1495 1502 1507 1512 1518 1520 1541 1543 1551 1552 1553 1556 1564
1567 1584 1595 1598 1606 1620 1625 1632 1648 1650 1655 1658 1677 1684
1686 1687 1697 1706 1712 1715 1717 1718 1728 1730 1736 1745 1748 1758
1759 1761 1777 1780 1798 1803 1813 1814 1846 1847 1857 1858 1860 1861
1865 1871 1888 1898 1917 1924 1931 1938 1939 1942 1953 1957 1958 1984
1986 2004 2006 2009 2018 2020 2023 2032 2034 2052 2074 2086 2090 2095
2096 2113 2117 2126 2127 2135 2149 2150 2185 2187 2192 2209 2211 2226
2228 2236 2237 2255 2257 2259 2266 2268 2269 2279 2280 2284 2290 2305
2317 2342 2347 2354 2358 2363 2374 2375 2376 2378 2390 2406 2408 2410
2414 2420 2424 2425 2435 2443 2446 2447 2448 2449 2457 2475 2476 2480
2482 2492 2497 2500 2515 2518 2527 2530 2535 2542 2545 2556 2557 2560
2588 2609 2626 2647 2648 2649 2650 2652 2657 2663 2674 2679 2684 2692
2697 2699 2703 2708 2711 2721 2726 2727 2730 2737 2748 2750 2753 2760
2762 2765 2766 2773 2777 2782 2784 2786 2806 2824 2826 2829 2844 2846
2850 2859 2892 2902 2916 2926 2936 2938 2948 2961 2964 2973 2974 2995
3005 3006 3023 3026 3030 3032 3047 3055 3058 3070 3094 3095 3102 3105
3109 3110 3111 3112 3113 3118 3119 3120 3121 3123 3126 3127 3128 3130
3131 3132 3133 3134 3139 3140 3141 3142 3143
524 disjoint connected subgraphs
Link number distribution:
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
446 258 324 369 448 429 328 233 116 68 41 27 20 20 7 2
258 least connected regions:
8 24 50 55 65 72 78 114 120 127 136 137 141 150 159 170 177 187 192 194 204 217 224 227 236 272 303 312 337 338 341 368 374 377 396 397 414 419 425 429 432 436 442 472 488 491 494 499 515 563 569 586 593 598 606 609 624 632 651 658 662 673 676 680 682 692 693 702 710 746 747 748 756 768 775 795 812 815 828 835 837 846 847 857 862 869 873 885 887 891 896 902 913 961 982 989 990 991 1025 1044 1048 1058 1070 1082 1084 1113 1115 1122 1171 1193 1246 1250 1265 1266 1267 1292 1293 1303 1304 1324 1342 1360 1362 1375 1376 1386 1395 1397 1408 1413 1414 1420 1460 1462 1477 1484 1494 1511 1531 1550 1578 1583 1593 1594 1602 1610 1627 1636 1666 1670 1680 1681 1685 1734 1742 1752 1760 1786 1808 1826 1839 1843 1850 1853 1854 1863 1864 1885 1889 1892 1896 1905 1954 1963 1967 1990 2000 2022 2051 2060 2062 2079 2108 2124 2130 2143 2146 2154 2173 2217 2218 2221 2234 2241 2247 2263 2271 2273 2283 2306 2313 2324 2356 2377 2404 2405 2411 2419 2493 2496 2498 2501 2513 2537 2539 2561 2564 2618 2621 2637 2638 2675 2677 2678 2694 2696 2723 2724 2725 2729 2754 2767 2770 2798 2817 2843 2860 2885 2904 2918 2947 2953 2971 2998 3018 3025 3028 3042 3056 3061 3075 3077 3114 3116 3117 3124 3125 3137 with 1 link
2 most connected regions:
95 1202 with 15 linksplot(st_geometry(sample_sf), col = "lightgrey", main = "Map with Distance-based Neighbors")
plot(neighbors_dist,
coords = sample_coords,
add = TRUE,
col = "red",
points = FALSE)
4. Spatial weights matrix
weights <- nb2listw(neighbors_queen, style = "W", zero.policy = TRUE) # The zero.policy = TRUE argument tells the function to handle regions with no neighbors by assigning them zero weights.
weights$weights[[1]][1] 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.16666675. Spatial lags and Moran’s I
5-1. Compute Spatial Lag
# Spatial lag of "di_index1"
sample_sf$lag_estimate <- lag.listw(weights, sample_sf$di_index1, NAOK = TRUE, zero.policy = TRUE) # NAOK = TRUE: This argument allows the function to proceed even if there are NA values in the data. It will treat NA values as 0 in the calculations.5-2. Scatter Plot of Spatial Lag
sample_sf %>%
filter(di_index1 < 2000) %>%
filter(lag_estimate < 2000) %>%
ggplot(aes(x = di_index1, y = lag_estimate)) +
geom_point(alpha = 0.3, color = "darkblue") +
geom_smooth(method = "lm", color = "red", se = FALSE) + # Added a trend line
geom_abline(intercept = 0, slope = 1, color = "gray50", linetype = "dashed") + # Added y=x line
theme_minimal() +
labs(title = "County-Level Concentrated Disadvanage vs. Spatial Lag",
x = "County-Level Concentrated Disadvanage",
y = "Spatial Lag",
caption = "Data source: CDC Multiple Causes of Death") +
coord_equal() + # Makes the plot square and ensures equal scaling
theme(plot.title = element_text(hjust = 0.5), # Center the title
text = element_text(size = 12)) # Increase text size for better readability`geom_smooth()` using formula = 'y ~ x'
5-3. Moran’s I Test
moran.test(sample_sf$di_index1, weights, zero.policy = TRUE)
Moran I test under randomisation
data: sample_sf$di_index1
weights: weights
n reduced by no-neighbour observations
Moran I statistic standard deviate = 52.535, p-value < 2.2e-16
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.5572042384 -0.0003192848 0.0001126234 # na.action : a function (default na.fail), can also be na.omit or na.exclude - in these cases the weights list will be subsetted to remove NAs in the data. It may be necessary to set zero.policy to TRUE because this subsetting may create no-neighbour observations. Note that only weights lists created without using the glist argument to nb2listw may be subsetted. na.pass is not permitted because it is meaningless in a permutation test.The small p-value suggests that we reject the null hypothesis of spatial randomness in our dataset. As the statistic is positive, it suggests that our data are spatially clustered; a negative statistic would suggest spatial uniformity.
6. Local spatial autocorrelation
6-1. Compute Local Gi* Statistic
localg_weights <- nb2listw(include.self(neighbors_queen))
# Compute the local Gi* statistic
localG_values <- localG(sample_sf$di_index1, weights)
# Add the local Gi* statistic to the data frame
sample_sf$localG <- as.numeric(localG_values)6-2. Plot Local Gi* Statistic
ggplot(sample_sf) +
geom_sf(aes(fill = localG), color = NA) +
scale_fill_distiller(palette = "RdYlBu") +
theme_void() +
labs(fill = "Local Gi* statistic")
6-3. Classify and Plot Clusters
sample_sf <- sample_sf %>%
mutate(hotspot = case_when(
localG >= 1.96 ~ "High cluster",
localG <= -1.96 ~ "Low cluster",
TRUE ~ "Not significant"
))
ggplot(sample_sf) +
geom_sf(aes(fill = hotspot), color = "grey20", size = 0.1) +
scale_fill_manual(values = c("red", "blue", "grey90")) +
theme_void()
7. LISA Analysis
7-1. Compute Local Moran’s I and LISA Clusters
set.seed(1983)
sample_sf$scaled_estimate <- as.numeric(scale(sample_sf$di_index1))
dfw_lisa <- localmoran_perm(
sample_sf$scaled_estimate,
weights,
nsim = 999L,
alternative = "two.sided"
) %>%
as_tibble() %>%
set_names(c("local_i", "exp_i", "var_i", "z_i", "p_i",
"p_i_sim", "pi_sim_folded", "skewness", "kurtosis"))
sample_sf_lisa <- sample_sf %>%
select(county_fips, scaled_estimate) %>%
mutate(lagged_estimate = lag.listw(weights, scaled_estimate)) %>%
bind_cols(dfw_lisa)
sample_sf_lisaSimple feature collection with 3136 features and 12 fields
Geometry type: GEOMETRY
Dimension: XY
Bounding box: xmin: -3115585 ymin: -1702303 xmax: 2263786 ymax: 1570639
Projected CRS: USA_Contiguous_Albers_Equal_Area_Conic
First 10 features:
county_fips scaled_estimate lagged_estimate local_i exp_i
1 31039 -1.2360599 -0.28376959 0.350868105 0.003188495
2 53069 -1.1004387 -0.32637535 0.359270622 -0.014908330
3 35011 0.3833977 0.50284081 0.192849493 -0.002868632
4 31109 -0.5451106 -0.96763937 0.527638702 -0.016439838
5 31129 -0.7217293 -0.65618153 0.473736530 0.007708748
6 46099 -0.8316973 -0.98072903 0.815929840 -0.001762842
7 48327 0.0736211 0.08098845 0.005964361 0.001729908
8 6091 -0.6833035 -0.39346203 0.268939735 -0.010277699
9 21053 1.4741253 1.34956780 1.990066634 -0.034569289
10 39063 -0.7537841 -0.67067478 0.505705238 -0.006525470
var_i z_i p_i p_i_sim pi_sim_folded skewness
1 0.229411571 0.7258913 0.467905440 0.502 0.251 -0.2929378
2 0.234314158 0.7730011 0.439521755 0.434 0.217 -0.5493536
3 0.027087740 1.1891720 0.234372002 0.244 0.122 0.4374227
4 0.038077511 2.7882218 0.005299824 0.002 0.001 -0.2705622
5 0.084461653 1.6035482 0.108813675 0.084 0.042 -0.3766060
6 0.091329844 2.7057256 0.006815534 0.004 0.002 -0.4079567
7 0.000752739 0.1543388 0.877342614 0.826 0.413 0.4313343
8 0.086756772 0.9479618 0.343148900 0.374 0.187 -0.4011883
9 0.394097989 3.2251122 0.001259233 0.006 0.003 0.3949150
10 0.086211655 1.7445468 0.081063792 0.064 0.032 -0.3743221
kurtosis geometry
1 0.07892479 MULTIPOLYGON (((-45725.24 4...
2 0.67738144 MULTIPOLYGON (((-2090804 12...
3 0.54519912 MULTIPOLYGON (((-760917.9 -...
4 0.02324642 MULTIPOLYGON (((-56821.82 3...
5 -0.03152218 MULTIPOLYGON (((-172448.4 3...
6 0.07627835 MULTIPOLYGON (((-49033.05 6...
7 0.23891262 MULTIPOLYGON (((-360251.5 -...
8 0.18276602 MULTIPOLYGON (((-2069003 49...
9 -0.07438500 MULTIPOLYGON (((970049.7 -2...
10 0.42744017 MULTIPOLYGON (((1028350 478...7-2. Classify and Plot LISA Clusters
lisa_clusters <- sample_sf_lisa %>%
mutate(lisa_cluster = case_when(
p_i >= 0.05 ~ "Not significant", # adjust p value; or you can erase this to not make "not significant" place
scaled_estimate > 0 & local_i > 0 ~ "High-high",
scaled_estimate > 0 & local_i < 0 ~ "High-low",
scaled_estimate < 0 & local_i > 0 ~ "Low-low",
scaled_estimate < 0 & local_i < 0 ~ "Low-high"
))
lisa_clusters %>% tabyl(lisa_cluster) lisa_clustergeometry n percent valid_percent
High-high 410 0.1307397959 0.130864986
High-low 33 0.0105229592 0.010533035
Low-high 30 0.0095663265 0.009575487
Low-low 329 0.1049107143 0.105011171
Not significant 2331 0.7433035714 0.744015321
<NA> 3 0.0009566327 NAcolor_values <- c(`High-high` = "red",
`High-low` = "pink",
`Low-low` = "blue",
`Low-high` = "lightblue",
`Not significant` = "white")
ggplot(lisa_clusters, aes(x = scaled_estimate,
y = lagged_estimate,
fill = lisa_cluster)) +
geom_point(color = "black", shape = 21, size = 2) +
theme_minimal() +
geom_hline(yintercept = 0, linetype = "dashed") +
geom_vline(xintercept = 0, linetype = "dashed") +
scale_fill_manual(values = color_values) +
labs(x = "concentrated disadvantage (z-score)",
y = "Spatial lag of concentrated disadvantage (z-score)",
fill = "Cluster type")
7-3. Scatter Plot of LISA Clusters
ggplot(lisa_clusters, aes(fill = lisa_cluster)) +
geom_sf(color = "grey30", size = 0.1) +
theme_void() +
scale_fill_manual(values = color_values) +
labs(fill = "Cluster type")