AddingControls
Cindy Lee
2026-03-03
R Markdown
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.1 ✔ stringr 1.5.2
## ✔ ggplot2 4.0.0 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── 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 errorslibrary(here)## here() starts at C:/Users/ddtmd/Desktop/2026 Spring/HotelDensity_PR_GooglePlaceAPIlibrary(skimr)
library(mapview)PR_SVL22 <- read_csv("C:/Users/ddtmd/Desktop/PR_Data/2022_CencusTrack_PR_SVL/2022_CT_svi_interactive_map.csv")## Rows: 939 Columns: 160
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (6): STATE, ST_ABBR, COUNTY, LOCATION, GeoLevel, Comparison
## dbl (154): ST, STCNTY, FIPS, AREA_SQMI, E_TOTPOP, M_TOTPOP, E_HU, M_HU, E_HH...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.PR_SVL14 <- read_csv("C:/Users/ddtmd/Desktop/PR_Data/2014_CencusTrack_PR_SVL/2014_CT_svi_interactive_map.csv")## Rows: 910 Columns: 130
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (125): AFFGEOID, TRACTCE, ST, STATE, ST_ABBR, COUNTY, FIPS, LOCATION, E_...
## dbl (5): STCNTY, AREA_SQMI, F_NOVEH, Shape.STArea(), Shape.STLength()
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.#2022 data: select a subset of columns
PR_SVL22_subset <- PR_SVL22 %>%
dplyr::select(RPL_THEME1, RPL_THEME2, RPL_THEME3, RPL_THEME4, RPL_THEMES)
view(PR_SVL22_subset)
#rename the columns
library(dplyr)
PR_SVL22_subset <- PR_SVL22 %>%
dplyr::select(RPL_THEME1, RPL_THEME2, RPL_THEME3, RPL_THEME4, RPL_THEMES) %>%
rename(
`Socioeconomic Status` = RPL_THEME1,
`Household Characteristics` = RPL_THEME2,
`Racial & Ethnic Minority Status` = RPL_THEME3,
`Housing Type & Transportation` = RPL_THEME4,
`Overall Summary Ranking` = RPL_THEMES
)
print(PR_SVL22_subset)## # A tibble: 939 × 5
## `Socioeconomic Status` `Household Characteristics` Racial & Ethnic Minority…¹
## <dbl> <dbl> <dbl>
## 1 0.888 0.244 0.457
## 2 0.750 0.109 0.522
## 3 0.774 0.145 0.522
## 4 0.887 0.552 0.522
## 5 0.186 0.424 0.0684
## 6 0.643 0.690 0.522
## 7 0.565 0.646 0.0456
## 8 0.306 0.463 0.0163
## 9 0.626 0.887 0.0087
## 10 0.260 0.351 0.0456
## # ℹ 929 more rows
## # ℹ abbreviated name: ¹`Racial & Ethnic Minority Status`
## # ℹ 2 more variables: `Housing Type & Transportation` <dbl>,
## # `Overall Summary Ranking` <dbl># 2014 Data: select a subset of columns
PR_SVL14_subset <- PR_SVL14 %>%
dplyr::select(RPL_THEME1, RPL_THEME2, RPL_THEME3, RPL_THEME4, RPL_THEMES)
view(PR_SVL14_subset)
# rename the columns
library(dplyr)
PR_SVL14_subset <- PR_SVL14 %>%
dplyr::select(RPL_THEME1, RPL_THEME2, RPL_THEME3, RPL_THEME4, RPL_THEMES) %>%
rename(
`Socioeconomic Status` = RPL_THEME1,
`Household Characteristics` = RPL_THEME2,
`Racial & Ethnic Minority Status` = RPL_THEME3,
`Housing Type & Transportation` = RPL_THEME4,
`Overall Summary Ranking` = RPL_THEMES
)
print(PR_SVL14_subset)## # A tibble: 910 × 5
## `Socioeconomic Status` `Household Characteristics` Racial & Ethnic Minority…¹
## <chr> <chr> <chr>
## 1 RPL_THEME1 RPL_THEME2 RPL_THEME3
## 2 0.9717 0.7423 0.902
## 3 0.9174 0.1927 0.8976
## 4 0.9842 0.1608 0.8645
## 5 0.9333 0.3932 0.6134
## 6 0.7115 0.9604 0.7313
## 7 0.7692 0.2126 0.8656
## 8 0.6719 0.2456 0.1993
## 9 0.6244 0.4989 0.1663
## 10 0.8778 0.5297 0.5672
## # ℹ 900 more rows
## # ℹ abbreviated name: ¹`Racial & Ethnic Minority Status`
## # ℹ 2 more variables: `Housing Type & Transportation` <chr>,
## # `Overall Summary Ranking` <chr># Make sure both datasets have the key column
names(PR_SVL14)## [1] "AFFGEOID" "TRACTCE" "ST"
## [4] "STATE" "ST_ABBR" "STCNTY"
## [7] "COUNTY" "FIPS" "LOCATION"
## [10] "AREA_SQMI" "E_TOTPOP" "M_TOTPOP"
## [13] "E_HU" "M_HU" "E_HH"
## [16] "M_HH" "E_POV" "M_POV"
## [19] "E_UNEMP" "M_UNEMP" "E_PCI"
## [22] "M_PCI" "E_NOHSDP" "M_NOHSDP"
## [25] "E_AGE65" "M_AGE65" "E_AGE17"
## [28] "M_AGE17" "E_DISABL" "M_DISABL"
## [31] "E_SNGPNT" "M_SNGPNT" "E_MINRTY"
## [34] "M_MINRTY" "E_LIMENG" "M_LIMENG"
## [37] "E_MUNIT" "M_MUNIT" "E_MOBILE"
## [40] "M_MOBILE" "E_CROWD" "M_CROWD"
## [43] "E_NOVEH" "M_NOVEH" "E_GROUPQ"
## [46] "M_GROUPQ" "EP_POV" "MP_POV"
## [49] "EP_UNEMP" "MP_UNEMP" "EP_PCI"
## [52] "MP_PCI" "EP_NOHSDP" "MP_NOHSDP"
## [55] "EP_AGE65" "MP_AGE65" "EP_AGE17"
## [58] "MP_AGE17" "EP_DISABL" "MP_DISABL"
## [61] "EP_SNGPNT" "MP_SNGPNT" "EP_MINRTY"
## [64] "MP_MINRTY" "EP_LIMENG" "MP_LIMENG"
## [67] "EP_MUNIT" "MP_MUNIT" "EP_MOBILE"
## [70] "MP_MOBILE" "EP_CROWD" "MP_CROWD"
## [73] "EP_NOVEH" "MP_NOVEH" "EP_GROUPQ"
## [76] "MP_GROUPQ" "EPL_POV" "EPL_UNEMP"
## [79] "EPL_PCI" "EPL_NOHSDP" "SPL_THEME1"
## [82] "RPL_THEME1" "EPL_AGE65" "EPL_AGE17"
## [85] "EPL_DISABL" "EPL_SNGPNT" "SPL_THEME2"
## [88] "RPL_THEME2" "EPL_MINRTY" "EPL_LIMENG"
## [91] "SPL_THEME3" "RPL_THEME3" "EPL_MUNIT"
## [94] "EPL_MOBILE" "EPL_CROWD" "EPL_NOVEH"
## [97] "EPL_GROUPQ" "SPL_THEME4" "RPL_THEME4"
## [100] "SPL_THEMES" "RPL_THEMES" "F_POV"
## [103] "F_UNEMP" "F_PCI" "F_NOHSDP"
## [106] "F_THEME1" "F_AGE65" "F_AGE17"
## [109] "F_DISABL" "F_SNGPNT" "F_THEME2"
## [112] "F_MINRTY" "F_LIMENG" "F_THEME3"
## [115] "F_MUNIT" "F_MOBILE" "F_CROWD"
## [118] "F_NOVEH" "F_GROUPQ" "F_THEME4"
## [121] "F_TOTAL" "E_UNINSUR" "M_UNINSUR"
## [124] "EP_UNINSUR" "MP_UNINSUR" "Shape"
## [127] "Shape.STArea()" "Shape.STLength()" "GeoLevel"
## [130] "Comparison"names(PR_SVL22)## [1] "ST" "STATE" "ST_ABBR" "STCNTY" "COUNTY"
## [6] "FIPS" "LOCATION" "AREA_SQMI" "E_TOTPOP" "M_TOTPOP"
## [11] "E_HU" "M_HU" "E_HH" "M_HH" "E_POV150"
## [16] "M_POV150" "E_UNEMP" "M_UNEMP" "E_HBURD" "M_HBURD"
## [21] "E_NOHSDP" "M_NOHSDP" "E_UNINSUR" "M_UNINSUR" "E_AGE65"
## [26] "M_AGE65" "E_AGE17" "M_AGE17" "E_DISABL" "M_DISABL"
## [31] "E_SNGPNT" "M_SNGPNT" "E_LIMENG" "M_LIMENG" "E_MINRTY"
## [36] "M_MINRTY" "E_MUNIT" "M_MUNIT" "E_MOBILE" "M_MOBILE"
## [41] "E_CROWD" "M_CROWD" "E_NOVEH" "M_NOVEH" "E_GROUPQ"
## [46] "M_GROUPQ" "EP_POV150" "MP_POV150" "EP_UNEMP" "MP_UNEMP"
## [51] "EP_HBURD" "MP_HBURD" "EP_NOHSDP" "MP_NOHSDP" "EP_UNINSUR"
## [56] "MP_UNINSUR" "EP_AGE65" "MP_AGE65" "EP_AGE17" "MP_AGE17"
## [61] "EP_DISABL" "MP_DISABL" "EP_SNGPNT" "MP_SNGPNT" "EP_LIMENG"
## [66] "MP_LIMENG" "EP_MINRTY" "MP_MINRTY" "EP_MUNIT" "MP_MUNIT"
## [71] "EP_MOBILE" "MP_MOBILE" "EP_CROWD" "MP_CROWD" "EP_NOVEH"
## [76] "MP_NOVEH" "EP_GROUPQ" "MP_GROUPQ" "EPL_POV150" "EPL_UNEMP"
## [81] "EPL_HBURD" "EPL_NOHSDP" "EPL_UNINSUR" "SPL_THEME1" "RPL_THEME1"
## [86] "EPL_AGE65" "EPL_AGE17" "EPL_DISABL" "EPL_SNGPNT" "EPL_LIMENG"
## [91] "SPL_THEME2" "RPL_THEME2" "EPL_MINRTY" "SPL_THEME3" "RPL_THEME3"
## [96] "EPL_MUNIT" "EPL_MOBILE" "EPL_CROWD" "EPL_NOVEH" "EPL_GROUPQ"
## [101] "SPL_THEME4" "RPL_THEME4" "SPL_THEMES" "RPL_THEMES" "F_POV150"
## [106] "F_UNEMP" "F_HBURD" "F_NOHSDP" "F_UNINSUR" "F_THEME1"
## [111] "F_AGE65" "F_AGE17" "F_DISABL" "F_SNGPNT" "F_LIMENG"
## [116] "F_THEME2" "F_MINRTY" "F_THEME3" "F_MUNIT" "F_MOBILE"
## [121] "F_CROWD" "F_NOVEH" "F_GROUPQ" "F_THEME4" "F_TOTAL"
## [126] "E_DAYPOP" "E_NOINT" "M_NOINT" "E_AFAM" "M_AFAM"
## [131] "E_HISP" "M_HISP" "E_ASIAN" "M_ASIAN" "E_AIAN"
## [136] "M_AIAN" "E_NHPI" "M_NHPI" "E_TWOMORE" "M_TWOMORE"
## [141] "E_OTHERRACE" "M_OTHERRACE" "EP_NOINT" "MP_NOINT" "EP_AFAM"
## [146] "MP_AFAM" "EP_HISP" "MP_HISP" "EP_ASIAN" "MP_ASIAN"
## [151] "EP_AIAN" "MP_AIAN" "EP_NHPI" "MP_NHPI" "EP_TWOMORE"
## [156] "MP_TWOMORE" "EP_OTHERRACE" "MP_OTHERRACE" "GeoLevel" "Comparison"# Make sure both datasets have the key column
names(PR_SVL14)## [1] "AFFGEOID" "TRACTCE" "ST"
## [4] "STATE" "ST_ABBR" "STCNTY"
## [7] "COUNTY" "FIPS" "LOCATION"
## [10] "AREA_SQMI" "E_TOTPOP" "M_TOTPOP"
## [13] "E_HU" "M_HU" "E_HH"
## [16] "M_HH" "E_POV" "M_POV"
## [19] "E_UNEMP" "M_UNEMP" "E_PCI"
## [22] "M_PCI" "E_NOHSDP" "M_NOHSDP"
## [25] "E_AGE65" "M_AGE65" "E_AGE17"
## [28] "M_AGE17" "E_DISABL" "M_DISABL"
## [31] "E_SNGPNT" "M_SNGPNT" "E_MINRTY"
## [34] "M_MINRTY" "E_LIMENG" "M_LIMENG"
## [37] "E_MUNIT" "M_MUNIT" "E_MOBILE"
## [40] "M_MOBILE" "E_CROWD" "M_CROWD"
## [43] "E_NOVEH" "M_NOVEH" "E_GROUPQ"
## [46] "M_GROUPQ" "EP_POV" "MP_POV"
## [49] "EP_UNEMP" "MP_UNEMP" "EP_PCI"
## [52] "MP_PCI" "EP_NOHSDP" "MP_NOHSDP"
## [55] "EP_AGE65" "MP_AGE65" "EP_AGE17"
## [58] "MP_AGE17" "EP_DISABL" "MP_DISABL"
## [61] "EP_SNGPNT" "MP_SNGPNT" "EP_MINRTY"
## [64] "MP_MINRTY" "EP_LIMENG" "MP_LIMENG"
## [67] "EP_MUNIT" "MP_MUNIT" "EP_MOBILE"
## [70] "MP_MOBILE" "EP_CROWD" "MP_CROWD"
## [73] "EP_NOVEH" "MP_NOVEH" "EP_GROUPQ"
## [76] "MP_GROUPQ" "EPL_POV" "EPL_UNEMP"
## [79] "EPL_PCI" "EPL_NOHSDP" "SPL_THEME1"
## [82] "RPL_THEME1" "EPL_AGE65" "EPL_AGE17"
## [85] "EPL_DISABL" "EPL_SNGPNT" "SPL_THEME2"
## [88] "RPL_THEME2" "EPL_MINRTY" "EPL_LIMENG"
## [91] "SPL_THEME3" "RPL_THEME3" "EPL_MUNIT"
## [94] "EPL_MOBILE" "EPL_CROWD" "EPL_NOVEH"
## [97] "EPL_GROUPQ" "SPL_THEME4" "RPL_THEME4"
## [100] "SPL_THEMES" "RPL_THEMES" "F_POV"
## [103] "F_UNEMP" "F_PCI" "F_NOHSDP"
## [106] "F_THEME1" "F_AGE65" "F_AGE17"
## [109] "F_DISABL" "F_SNGPNT" "F_THEME2"
## [112] "F_MINRTY" "F_LIMENG" "F_THEME3"
## [115] "F_MUNIT" "F_MOBILE" "F_CROWD"
## [118] "F_NOVEH" "F_GROUPQ" "F_THEME4"
## [121] "F_TOTAL" "E_UNINSUR" "M_UNINSUR"
## [124] "EP_UNINSUR" "MP_UNINSUR" "Shape"
## [127] "Shape.STArea()" "Shape.STLength()" "GeoLevel"
## [130] "Comparison"names(PR_SVL22)## [1] "ST" "STATE" "ST_ABBR" "STCNTY" "COUNTY"
## [6] "FIPS" "LOCATION" "AREA_SQMI" "E_TOTPOP" "M_TOTPOP"
## [11] "E_HU" "M_HU" "E_HH" "M_HH" "E_POV150"
## [16] "M_POV150" "E_UNEMP" "M_UNEMP" "E_HBURD" "M_HBURD"
## [21] "E_NOHSDP" "M_NOHSDP" "E_UNINSUR" "M_UNINSUR" "E_AGE65"
## [26] "M_AGE65" "E_AGE17" "M_AGE17" "E_DISABL" "M_DISABL"
## [31] "E_SNGPNT" "M_SNGPNT" "E_LIMENG" "M_LIMENG" "E_MINRTY"
## [36] "M_MINRTY" "E_MUNIT" "M_MUNIT" "E_MOBILE" "M_MOBILE"
## [41] "E_CROWD" "M_CROWD" "E_NOVEH" "M_NOVEH" "E_GROUPQ"
## [46] "M_GROUPQ" "EP_POV150" "MP_POV150" "EP_UNEMP" "MP_UNEMP"
## [51] "EP_HBURD" "MP_HBURD" "EP_NOHSDP" "MP_NOHSDP" "EP_UNINSUR"
## [56] "MP_UNINSUR" "EP_AGE65" "MP_AGE65" "EP_AGE17" "MP_AGE17"
## [61] "EP_DISABL" "MP_DISABL" "EP_SNGPNT" "MP_SNGPNT" "EP_LIMENG"
## [66] "MP_LIMENG" "EP_MINRTY" "MP_MINRTY" "EP_MUNIT" "MP_MUNIT"
## [71] "EP_MOBILE" "MP_MOBILE" "EP_CROWD" "MP_CROWD" "EP_NOVEH"
## [76] "MP_NOVEH" "EP_GROUPQ" "MP_GROUPQ" "EPL_POV150" "EPL_UNEMP"
## [81] "EPL_HBURD" "EPL_NOHSDP" "EPL_UNINSUR" "SPL_THEME1" "RPL_THEME1"
## [86] "EPL_AGE65" "EPL_AGE17" "EPL_DISABL" "EPL_SNGPNT" "EPL_LIMENG"
## [91] "SPL_THEME2" "RPL_THEME2" "EPL_MINRTY" "SPL_THEME3" "RPL_THEME3"
## [96] "EPL_MUNIT" "EPL_MOBILE" "EPL_CROWD" "EPL_NOVEH" "EPL_GROUPQ"
## [101] "SPL_THEME4" "RPL_THEME4" "SPL_THEMES" "RPL_THEMES" "F_POV150"
## [106] "F_UNEMP" "F_HBURD" "F_NOHSDP" "F_UNINSUR" "F_THEME1"
## [111] "F_AGE65" "F_AGE17" "F_DISABL" "F_SNGPNT" "F_LIMENG"
## [116] "F_THEME2" "F_MINRTY" "F_THEME3" "F_MUNIT" "F_MOBILE"
## [121] "F_CROWD" "F_NOVEH" "F_GROUPQ" "F_THEME4" "F_TOTAL"
## [126] "E_DAYPOP" "E_NOINT" "M_NOINT" "E_AFAM" "M_AFAM"
## [131] "E_HISP" "M_HISP" "E_ASIAN" "M_ASIAN" "E_AIAN"
## [136] "M_AIAN" "E_NHPI" "M_NHPI" "E_TWOMORE" "M_TWOMORE"
## [141] "E_OTHERRACE" "M_OTHERRACE" "EP_NOINT" "MP_NOINT" "EP_AFAM"
## [146] "MP_AFAM" "EP_HISP" "MP_HISP" "EP_ASIAN" "MP_ASIAN"
## [151] "EP_AIAN" "MP_AIAN" "EP_NHPI" "MP_NHPI" "EP_TWOMORE"
## [156] "MP_TWOMORE" "EP_OTHERRACE" "MP_OTHERRACE" "GeoLevel" "Comparison"library(dplyr)
# Step 1. Make cleaned subset for 2014
PR_SVL14_subset <- PR_SVL14 %>%
dplyr::select(
STCNTY, # county FIPS
COUNTY, # county name
LOCATION, # tract description
FIPS, # tract ID
AREA_SQMI, # area (sq mi)
E_TOTPOP, # total population
RPL_THEME1,
RPL_THEME2,
RPL_THEME3,
RPL_THEME4,
RPL_THEMES
) %>%
rename(
`Socioeconomic Status` = RPL_THEME1,
`Household Characteristics` = RPL_THEME2,
`Racial & Ethnic Minority Status` = RPL_THEME3,
`Housing Type & Transportation` = RPL_THEME4,
`Overall Summary Ranking` = RPL_THEMES
)
# Step 2. Make cleaned subset for 2022
PR_SVL22_subset <- PR_SVL22 %>%
dplyr::select(
STCNTY,
COUNTY,
LOCATION,
FIPS,
AREA_SQMI,
E_TOTPOP,
RPL_THEME1,
RPL_THEME2,
RPL_THEME3,
RPL_THEME4,
RPL_THEMES
) %>%
rename(
`Socioeconomic Status` = RPL_THEME1,
`Household Characteristics` = RPL_THEME2,
`Racial & Ethnic Minority Status` = RPL_THEME3,
`Housing Type & Transportation` = RPL_THEME4,
`Overall Summary Ranking` = RPL_THEMES
)
# Step 3. Join 2014 + 2022 using FIPS (tract ID)
PR_SVL14_subset <- PR_SVL14_subset %>%
slice(-1) # Delete first row becuase it's not numeric (2022 is okay)
str(PR_SVL14_subset$`Overall Summary Ranking`)## chr [1:909] "0.8948" "0.5814" "0.6538" "0.8281" "0.8937" "0.5011" "0.4355" ...comparison <- PR_SVL14_subset %>%
mutate(FIPS = as.character(FIPS)) %>% # FIPS 문자로
rename(
Overall_2014 = `Overall Summary Ranking`
) %>%
left_join(
PR_SVL22_subset %>%
mutate(FIPS = as.character(FIPS)) %>%
rename(
Overall_2022 = `Overall Summary Ranking`
),
by = "FIPS",
suffix = c("_2014", "_2022")
) %>%
mutate(
# 여기서 진짜 숫자로 바꿔줌 (2014는 꼭 필요, 2022는 안전용)
Overall_2014 = as.numeric(Overall_2014),
Overall_2022 = as.numeric(Overall_2022),
Change = Overall_2022 - Overall_2014,
Change_Direction = case_when(
Change > 0 ~ "Increased (worse / more vulnerable)",
Change < 0 ~ "Decreased (better / less vulnerable)",
TRUE ~ "No change"
)
)
view(comparison)
table(comparison$Change_Direction)##
## Decreased (better / less vulnerable) Increased (worse / more vulnerable)
## 387 438
## No change
## 84# ---- 1. FIPS 문자형 통일 + 헤더행 제거 ----
PR_SVL14_subset <- PR_SVL14_subset %>%
mutate(FIPS = as.character(FIPS)) %>%
filter(FIPS != "FIPS", !is.na(FIPS))
PR_SVL22_subset <- PR_SVL22_subset %>%
mutate(
FIPS = format(FIPS, scientific = FALSE, trim = TRUE),
FIPS = as.character(FIPS)
) %>%
filter(FIPS != "FIPS", !is.na(FIPS))
# ---- 2. SVI 점수 numeric 변환 ----
svi_cols <- c(
"Socioeconomic Status",
"Household Characteristics",
"Racial & Ethnic Minority Status",
"Housing Type & Transportation",
"Overall Summary Ranking"
)
PR_SVL14_subset <- PR_SVL14_subset %>%
mutate(across(all_of(svi_cols), as.numeric))
PR_SVL22_subset <- PR_SVL22_subset %>%
mutate(across(all_of(svi_cols), as.numeric))
# ---- 3. 0~1 범위 밖 이상치 NA 처리 ----
PR_SVL14_subset <- PR_SVL14_subset %>%
mutate(
`Overall Summary Ranking` =
ifelse(`Overall Summary Ranking` < 0 |
`Overall Summary Ranking` > 1,
NA, `Overall Summary Ranking`)
)
PR_SVL22_subset <- PR_SVL22_subset %>%
mutate(
`Overall Summary Ranking` =
ifelse(`Overall Summary Ranking` < 0 |
`Overall Summary Ranking` > 1,
NA, `Overall Summary Ranking`)
)
# ---- 4. comparison 최종 버전 만들기 ----
comparison <- PR_SVL14_subset %>%
rename(
Overall_2014 = `Overall Summary Ranking`,
COUNTY_2014 = COUNTY,
LOCATION_2014 = LOCATION
) %>%
left_join(
PR_SVL22_subset %>%
rename(
Overall_2022 = `Overall Summary Ranking`,
COUNTY_2022 = COUNTY,
LOCATION_2022 = LOCATION
),
by = "FIPS",
suffix = c("_2014", "_2022")
) %>%
filter(!is.na(Overall_2014), !is.na(Overall_2022)) %>%
mutate(
Change = Overall_2022 - Overall_2014,
Change_Direction = case_when(
Change > 0 ~ "Increased (worse / more vulnerable)",
Change < 0 ~ "Decreased (better / less vulnerable)",
TRUE ~ "No change"
)
)
# ---- 5. 요약 / 히스토그램 ----
library(ggplot2)
summary(comparison$Change)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.69240 -0.12255 0.01160 0.01179 0.15470 0.78160hist(comparison$Change, breaks = 40)
ggplot(comparison, aes(x = Change)) +
geom_histogram(bins = 40, fill = "steelblue", color = "white") +
labs(
title = "Distribution of Change in Social Vulnerability (2022 - 2014)",
x = "Change in Overall Summary Ranking (Δ vulnerability)",
y = "Number of Census Tracts"
) +
theme_minimal()
## Mapping
#Mapping
library(sf)## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUElibrary(ggplot2)
library(dplyr)
library(tigris)## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.options(tigris_use_cache = TRUE)
# 1. Puerto Rico Tracts 불러오기
pr_tracts <- tracts("PR", year = 2022, cb = TRUE) %>%
st_transform(32620) %>% # UTM zone20N (Puerto Rico)
mutate(FIPS = as.character(GEOID))
# 2. tract + comparison 조인
pr_change_sf <- pr_tracts %>%
dplyr::left_join(
comparison %>% dplyr::select(FIPS, Change, Change_Direction),
by = "FIPS"
)
# 3. Continuous Change Map
ggplot(pr_change_sf) +
geom_sf(aes(fill = Change), color="grey40", size=0.1) +
scale_fill_distiller(palette="RdBu", direction=-1, limits=c(-0.7,0.7)) +
labs(
title="Change in Social Vulnerability (2014→2022)",
fill="Δ SVI"
) +
theme_minimal()
## 3. d
#Scatter plot
library(sf)
library(dplyr)
library(sf)
coast_pr <- st_read(
"C:/Users/ddtmd/Desktop/2026 Spring/HotelDensity_PR_GooglePlaceAPI/data_processed/pr_coastline.gpkg"
) %>%
st_transform(32620) # UTM zone 20N (Puerto Rico)## Reading layer `pr_coastline' from data source
## `C:\Users\ddtmd\Desktop\2026 Spring\HotelDensity_PR_GooglePlaceAPI\data_processed\pr_coastline.gpkg'
## using driver `GPKG'
## Simple feature collection with 16 features and 0 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: -23817.68 ymin: 1980811 xmax: 265239.9 ymax: 2052181
## Projected CRS: WGS 84 / UTM zone 20N#Creating Centroid
tract_centroids <- pr_change_sf %>%
st_centroid()## Warning: st_centroid assumes attributes are constant over geometries# centroid → coastline 거리 행렬 계산
dist_matrix <- st_distance(tract_centroids, coast_pr)
# 각 tract마다 "가장 가까운 해안까지의 거리"만 뽑기 (미터 단위)
tract_centroids$dist_to_coast_m <- apply(dist_matrix, 1, min)
# 요약 확인
summary(tract_centroids$dist_to_coast_m)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.2 2501.0 5902.5 8164.3 11632.5 27400.2library(ggplot2)
# =========================================================
library(sf)
library(dplyr)
library(ggplot2)
# 0) FIPS 타입 통일
pr_tracts <- pr_tracts %>%
mutate(FIPS = as.character(GEOID))
comparison2 <- comparison %>%
dplyr::mutate(FIPS = as.character(FIPS)) %>%
dplyr::select(FIPS, Overall_2014, Overall_2022, Change, Change_Direction)
# 1) join
analysis_sf <- pr_tracts %>%
dplyr::left_join(comparison2, by = "FIPS")
# 2) (optional) NA 확인
summary(is.na(analysis_sf$Overall_2014))## Mode FALSE TRUE
## logical 823 116summary(is.na(analysis_sf$Overall_2022))## Mode FALSE TRUE
## logical 823 116summary(is.na(analysis_sf$Change))## Mode FALSE TRUE
## logical 823 116## 2 Distance to coast 다시붙이기
sf::sf_use_s2(FALSE)## Spherical geometry (s2) switched offcoast_pr <- st_read(
"C:/Users/ddtmd/Desktop/2026 Spring/HotelDensity_PR_GooglePlaceAPI/data_processed/pr_coastline.gpkg",
quiet = TRUE
) %>% st_transform(32620)
analysis_sf <- analysis_sf %>% st_transform(32620)
cent <- st_centroid(analysis_sf)## Warning: st_centroid assumes attributes are constant over geometriesdmat <- st_distance(cent, coast_pr)
analysis_sf$dist_to_coast_m <- as.numeric(apply(dmat, 1, min))
summary(analysis_sf$dist_to_coast_m)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.2 2501.0 5902.5 8164.3 11632.5 27400.2## 3. Level plot 2014 / 2022 / Change (3개)
df <- analysis_sf %>% st_drop_geometry() %>%
filter(!is.na(dist_to_coast_m), !is.na(Overall_2014), !is.na(Overall_2022), !is.na(Change))
# (A) Level 2014
ggplot(df, aes(x = dist_to_coast_m, y = Overall_2014)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "SVI Level (2014) vs Distance to Coastline (PR)",
x = "Distance to Coast (meters)",
y = "SVI (Overall Summary Ranking, 2014)"
) +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
# (B) Level 2022
ggplot(df, aes(x = dist_to_coast_m, y = Overall_2022)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "SVI Level (2022) vs Distance to Coastline (PR)",
x = "Distance to Coast (meters)",
y = "SVI (Overall Summary Ranking, 2022)"
) +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
# (C) Change
ggplot(df, aes(x = dist_to_coast_m, y = Change)) +
geom_point(alpha = 0.6) +
geom_smooth(method = "lm", se = TRUE) +
labs(
title = "ΔSVI (2022–2014) vs Distance to Coastline (PR)",
x = "Distance to Coast (meters)",
y = "Δ SVI (2022 - 2014)"
) +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
##4. 회귀선형 확인
m14 <- lm(Overall_2014 ~ dist_to_coast_m, data = df)
m22 <- lm(Overall_2022 ~ dist_to_coast_m, data = df)
mchg <- lm(Change ~ dist_to_coast_m, data = df)
summary(m14)##
## Call:
## lm(formula = Overall_2014 ~ dist_to_coast_m, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5448 -0.2477 -0.0100 0.2393 0.5287
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.673e-01 1.537e-02 30.410 < 2e-16 ***
## dist_to_coast_m 4.210e-06 1.387e-06 3.035 0.00248 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2852 on 821 degrees of freedom
## Multiple R-squared: 0.0111, Adjusted R-squared: 0.009891
## F-statistic: 9.211 on 1 and 821 DF, p-value: 0.002481summary(m22)##
## Call:
## lm(formula = Overall_2022 ~ dist_to_coast_m, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.53016 -0.24258 0.01179 0.24680 0.49894
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.356e-01 1.536e-02 34.880 <2e-16 ***
## dist_to_coast_m -2.476e-06 1.386e-06 -1.787 0.0744 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.285 on 821 degrees of freedom
## Multiple R-squared: 0.003873, Adjusted R-squared: 0.00266
## F-statistic: 3.192 on 1 and 821 DF, p-value: 0.07436summary(mchg)##
## Call:
## lm(formula = Change ~ dist_to_coast_m, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71314 -0.13323 0.00029 0.14906 0.75157
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.827e-02 1.247e-02 5.477 5.75e-08 ***
## dist_to_coast_m -6.686e-06 1.125e-06 -5.942 4.15e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2314 on 821 degrees of freedom
## Multiple R-squared: 0.04124, Adjusted R-squared: 0.04007
## F-statistic: 35.31 on 1 and 821 DF, p-value: 4.149e-094. import hotel data
# FINAL (Copy/Paste): PR tract-level dataset
# - SVI Change (Change) 2014→2022 : from pr_change_sf
# - Distance to coast (dist_to_coast_m)
# - Hotel density (lodging_dens_km2) + Winsorize (p99)
# - Regression: Change ~ hotel_density_w + dist_to_coast_m
# =========================================================
# 패키지 로드 + s2 끄기
library(sf)
library(dplyr)
library(readr)
library(ggplot2)
sf::sf_use_s2(FALSE)
# “pr_change_sf”가 이미 있는지 확인
# ---- 0) MUST have pr_change_sf (contains Change) ----
stopifnot(exists("pr_change_sf"))
# quick check
stopifnot("Change" %in% names(pr_change_sf))
stopifnot("FIPS" %in% names(pr_change_sf))
# ---- 1) Read hotel density ----
dir_hd <- "C:/Users/ddtmd/Desktop/2026 Spring/HotelDensity_PR_GooglePlaceAPI"
#호텔 밀도 데이터 읽고 필요한 컬럼만 정리
hotel_den_raw <- read_csv(
file.path(dir_hd, "PR_tract_lodging_density.csv"),
show_col_types = FALSE
)
hotel_den <- hotel_den_raw %>%
transmute(
FIPS = as.character(GEOID),
hotel_density = as.numeric(lodging_dens_km2), # lodging density per km2
lodging_total = as.numeric(lodging_total),
area_km2 = as.numeric(area_km2)
)
#SVI 변화 데이터 + 호텔 밀도 붙이기 (조인)
# ---- 2) Join hotel density to pr_change_sf ----
analysis_sf <- pr_change_sf %>%
mutate(FIPS = as.character(FIPS)) %>%
left_join(hotel_den, by = "FIPS")
# sanity check
summary(analysis_sf$hotel_density)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.5379 0.0000 64.9653sum(is.na(analysis_sf$hotel_density)) # should be 0## [1] 05. Import distance to coastal from centroid
#해안까지 거리 만들기 (핵심 공간 계산)
# ---- 3) Distance to PR coastline (centroid -> min distance) ----
coast_pr <- st_read(
"C:/Users/ddtmd/Desktop/2026 Spring/HotelDensity_PR_GooglePlaceAPI/data_processed/pr_coastline.gpkg",
quiet = TRUE
) %>% st_transform(32620)
analysis_sf <- analysis_sf %>% st_transform(32620)
cent <- st_centroid(analysis_sf)## Warning: st_centroid assumes attributes are constant over geometriesdmat <- st_distance(cent, coast_pr)
analysis_sf$dist_to_coast_m <- as.numeric(apply(dmat, 1, min))
summary(analysis_sf$dist_to_coast_m)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 15.2 2501.0 5902.5 8164.3 11632.5 27400.26. Ready for anlysis
#호텔밀도 극단값 자르기 (winsorize p99)
# ---- 4) Winsorize hotel density at 99th percentile (no log) ----
p99 <- quantile(analysis_sf$hotel_density, 0.99, na.rm = TRUE)
analysis_sf <- analysis_sf %>%
mutate(hotel_density_w = pmin(hotel_density, p99))
quantile(analysis_sf$hotel_density, c(.95, .99), na.rm = TRUE)## 95% 99%
## 1.137027 14.438188quantile(analysis_sf$hotel_density_w, c(.95, .99), na.rm = TRUE)## 95% 99%
## 1.137027 13.922990#회귀용 데이터프레임 만들기 (geometry 제거 + NA 제거)
# ---- 5) Regression-ready df (drop geometry + remove NA) ----
df_reg <- analysis_sf %>%
st_drop_geometry() %>%
filter(
!is.na(Change),
!is.na(hotel_density_w),
!is.na(dist_to_coast_m)
)
# ---- 6) Main regression ----
m <- lm(Change ~ hotel_density_w + dist_to_coast_m, data = df_reg)
summary(m)##
## Call:
## lm(formula = Change ~ hotel_density_w + dist_to_coast_m, data = df_reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.71171 -0.13453 -0.00019 0.14989 0.75261
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.665e-02 1.292e-02 5.157 3.14e-07 ***
## hotel_density_w 2.228e-03 4.662e-03 0.478 0.633
## dist_to_coast_m -6.585e-06 1.145e-06 -5.749 1.26e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2315 on 820 degrees of freedom
## Multiple R-squared: 0.0415, Adjusted R-squared: 0.03917
## F-statistic: 17.75 on 2 and 820 DF, p-value: 2.831e-08# ---- 7) (Optional) Quick scatter checks ----
ggplot(df_reg, aes(x = hotel_density_w, y = Change)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", se = TRUE) +
theme_minimal() +
labs(title = "ΔSVI vs Hotel density (winsorized p99)", x = "Hotel density (per km2, winsorized)", y = "ΔSVI (2022-2014)")## `geom_smooth()` using formula = 'y ~ x'
ggplot(df_reg, aes(x = dist_to_coast_m, y = Change)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", se = TRUE) +
theme_minimal() +
labs(title = "ΔSVI vs Distance to coast", x = "Distance to coast (m)", y = "ΔSVI (2022-2014)")## `geom_smooth()` using formula = 'y ~ x'
#호텔밀도 → 해안거리
m_b <- lm(dist_to_coast_m ~ hotel_density_w, data = df_reg)
summary(m_b)##
## Call:
## lm(formula = dist_to_coast_m ~ hotel_density_w, data = df_reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -8691 -5546 -2268 3621 18694
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 8706.4 250.5 34.754 < 2e-16 ***
## hotel_density_w -751.0 139.6 -5.379 9.76e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7053 on 821 degrees of freedom
## Multiple R-squared: 0.03405, Adjusted R-squared: 0.03287
## F-statistic: 28.94 on 1 and 821 DF, p-value: 9.759e-08#호텔밀도만 넣은 모델
m_c <- lm(Change ~ hotel_density_w, data = df_reg)
summary(m_c)##
## Call:
## lm(formula = Change ~ hotel_density_w, data = df_reg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.70172 -0.13736 -0.00062 0.14355 0.77228
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.009320 0.008380 1.112 0.266
## hotel_density_w 0.007173 0.004670 1.536 0.125
##
## Residual standard error: 0.2359 on 821 degrees of freedom
## Multiple R-squared: 0.002865, Adjusted R-squared: 0.001651
## F-statistic: 2.359 on 1 and 821 DF, p-value: 0.12497. Test
#mediation test
# =========================================================
# Mediation (package: mediation)
# X = hotel_density_w
# M = dist_to_coast_m
# Y = Change (ΔSVI 2022-2014)
# =========================================================
# 0) packages
install_if_missing <- function(pkgs){
to_install <- pkgs[!pkgs %in% installed.packages()[, "Package"]]
if(length(to_install) > 0) install.packages(to_install)
}
install_if_missing(c("mediation", "dplyr"))
library(mediation)## Loading required package: MASS##
## Attaching package: 'MASS'## The following object is masked from 'package:dplyr':
##
## select## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack## Loading required package: mvtnorm## Loading required package: sandwich## mediation: Causal Mediation Analysis
## Version: 4.5.1library(dplyr)
# 1) make sure df_reg exists and has needed columns
stopifnot(exists("df_reg"))
stopifnot(all(c("Change","hotel_density_w","dist_to_coast_m") %in% names(df_reg)))
df_med <- df_reg %>%
filter(
!is.na(Change),
!is.na(hotel_density_w),
!is.na(dist_to_coast_m)
)
# 2) mediator model (M ~ X)
med.fit <- lm(dist_to_coast_m ~ hotel_density_w, data = df_med)
# 3) outcome model (Y ~ X + M)
out.fit <- lm(Change ~ hotel_density_w + dist_to_coast_m, data = df_med)
# 4) run mediation
# sims=2000 is usually fine; increase to 5000+ for final reporting
set.seed(123)
med.out <- mediate(
model.m = med.fit,
model.y = out.fit,
treat = "hotel_density_w",
mediator= "dist_to_coast_m",
sims = 2000,
boot = TRUE # robust to non-normality
)## Running nonparametric bootstrapsummary(med.out)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.00494542 0.00325845 0.00701564 <2e-16 ***
## ADE 0.00222778 -0.00493922 0.00951248 0.543
## Total Effect 0.00717320 0.00028132 0.01459750 0.039 *
## Prop. Mediated 0.68943029 0.26832347 4.62386151 0.039 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 823
##
##
## Simulations: 2000# 5) optional: plot ACME/ADE
plot(med.out)
Controal Variables
#Contral Variables
# ---- packages ----
install_if_missing <- function(pkgs){
to_install <- pkgs[!pkgs %in% installed.packages()[, "Package"]]
if(length(to_install) > 0) install.packages(to_install)
}
install_if_missing(c("tidycensus","tigris","sf","dplyr","stringr","jsonlite"))
library(tidycensus)
library(tigris)
library(sf)
library(dplyr)
library(stringr)
library(jsonlite)##
## Attaching package: 'jsonlite'## The following object is masked from 'package:purrr':
##
## flattenoptions(tigris_use_cache = TRUE)
sf::sf_use_s2(FALSE)
# (권장) API key 한 번만 세팅
# census_api_key("YOUR_KEY", install = TRUE, overwrite = TRUE)
# readRenviron("~/.Renviron")
#1. DP03 프로필(산업 % 비중)만 “직접 API”로 가져오는 함수
#tidycensus가 DP profile에서 404 나는 경우가 있어서, DP03만 direct API로 안정적으로 가져옴.
get_pr_dp03_profile_pct <- function(year_end){
# 우리가 원하는 DP03 "Percent" 변수 목록 (P)
vars_dp03P <- c(
pct_construction = "DP03_0034P",
pct_manufacturing = "DP03_0035P",
pct_wholesale = "DP03_0036P",
pct_retail = "DP03_0037P",
pct_transport = "DP03_0038P",
pct_information = "DP03_0039P",
pct_finance_re = "DP03_0040P",
pct_prof_sci_mgmt = "DP03_0041P",
pct_edu_health = "DP03_0042P",
pct_tourism_proxy = "DP03_0043P", # Arts/Ent/Rec/Accom/Food
pct_other_services = "DP03_0044P",
pct_public_admin = "DP03_0045P"
)
# Census API에서는 보통 Percent estimate가 "...PE" 형태로 옴
api_vars <- paste0(unname(vars_dp03P), "E") # ex: DP03_0043PE
names(api_vars) <- names(vars_dp03P)
base <- sprintf("https://api.census.gov/data/%s/acs/acs5/profile", year_end)
get_str <- paste(c("NAME", api_vars), collapse = ",")
url <- paste0(
base,
"?get=", URLencode(get_str, reserved = TRUE),
"&for=tract:*",
"&in=state:72"
)
raw <- jsonlite::fromJSON(url)
df <- as.data.frame(raw[-1, ], stringsAsFactors = FALSE)
names(df) <- raw[1, ]
out <- df %>%
dplyr::transmute(
FIPS = paste0(state, county, tract),
dplyr::across(dplyr::all_of(api_vars), ~ suppressWarnings(as.numeric(.x)) / 100)
)
# pct_* 로 rename (자동 매핑)
rename_map <- setNames(names(vars_dp03P), api_vars)
out <- out %>% dplyr::rename(!!!rename_map)
out
}
#2. PR tract controls 한 번에 가져오기 (DP03 + basic + housing + area)
pull_pr_controls <- function(year_end){
# 1) DP03 (percent industry shares)
dp03 <- get_pr_dp03_profile_pct(year_end)
# 2) basic SES
vars_basic <- c(
total_pop = "B01003_001",
med_income = "B19013_001",
pov_total = "B17001_001",
pov_below = "B17001_002"
)
basic <- tidycensus::get_acs(
geography = "tract",
variables = vars_basic,
state = "PR",
year = year_end,
survey = "acs5",
output = "wide"
) %>%
dplyr::transmute(
FIPS = GEOID,
total_pop = total_popE,
med_income = med_incomeE,
pov_rate = dplyr::if_else(pov_totalE > 0, pov_belowE / pov_totalE, NA_real_)
)
# 3) housing pressure
vars_housing <- c(
tenure_total = "B25003_001",
renters = "B25003_003",
med_gross_rent = "B25064_001",
rb_30_34 = "B25070_008",
rb_35_39 = "B25070_009",
rb_40_49 = "B25070_010",
rb_50p = "B25070_011",
rb_total = "B25070_001"
)
housing <- tidycensus::get_acs(
geography = "tract",
variables = vars_housing,
state = "PR",
year = year_end,
survey = "acs5",
output = "wide"
) %>%
dplyr::transmute(
FIPS = GEOID,
renter_share = dplyr::if_else(tenure_totalE > 0, rentersE / tenure_totalE, NA_real_),
med_gross_rent = med_gross_rentE,
rent_burden_30p = dplyr::if_else(
rb_totalE > 0,
(rb_30_34E + rb_35_39E + rb_40_49E + rb_50pE) / rb_totalE,
NA_real_
)
)
# 4) area for pop density (고정 2022 CB 추천)
tr <- tigris::tracts("PR", year = 2022, cb = TRUE) %>%
sf::st_transform(32620) %>%
dplyr::mutate(
FIPS = GEOID,
area_km2 = as.numeric(sf::st_area(geometry)) / 1e6
) %>%
sf::st_drop_geometry() %>%
dplyr::select(FIPS, area_km2)
out <- dp03 %>%
dplyr::left_join(basic, by = "FIPS") %>%
dplyr::left_join(housing, by = "FIPS") %>%
dplyr::left_join(tr, by = "FIPS") %>%
dplyr::mutate(
pop_density_km2 = dplyr::if_else(area_km2 > 0, total_pop / area_km2, NA_real_),
year_end = year_end
)
out
}
#3.
acs14 <- pull_pr_controls(2014) %>% dplyr::rename_with(~ paste0(.x, "_2014"), -FIPS)## Getting data from the 2010-2014 5-year ACS## Warning: • You have not set a Census API key. Users without a key are limited to 500
## queries per day and may experience performance limitations.
## ℹ For best results, get a Census API key at
## http://api.census.gov/data/key_signup.html and then supply the key to the
## `census_api_key()` function to use it throughout your tidycensus session.
## This warning is displayed once per session.## Getting data from the 2010-2014 5-year ACSacs22 <- pull_pr_controls(2022) %>% dplyr::rename_with(~ paste0(.x, "_2022"), -FIPS)## Getting data from the 2018-2022 5-year ACS## Getting data from the 2018-2022 5-year ACS# 확인
names(acs14)[grepl("^pct_", names(acs14))]## character(0)names(acs14)[grepl("tourism", names(acs14))]## character(0)head(acs14)## FIPS DP03_0034PE_2014 DP03_0035PE_2014 DP03_0036PE_2014
## 1 72029100200 0.077 0.133 0.021
## 2 72029100400 0.109 0.036 0.052
## 3 72029100800 0.061 0.016 0.025
## 4 72029100700 0.168 0.125 0.020
## 5 72029100101 0.039 0.048 0.033
## 6 72029100502 0.019 0.095 0.013
## DP03_0037PE_2014 DP03_0038PE_2014 DP03_0039PE_2014 DP03_0040PE_2014
## 1 0.086 0.031 0.000 0.075
## 2 0.126 0.037 0.008 0.012
## 3 0.123 0.057 0.000 0.055
## 4 0.069 0.061 0.000 0.047
## 5 0.115 0.073 0.024 0.040
## 6 0.155 0.023 0.000 0.074
## DP03_0041PE_2014 DP03_0042PE_2014 DP03_0043PE_2014 DP03_0044PE_2014
## 1 0.093 0.245 0.115 0.037
## 2 0.077 0.263 0.127 0.035
## 3 0.098 0.370 0.073 0.016
## 4 0.075 0.224 0.046 0.046
## 5 0.117 0.208 0.137 0.035
## 6 0.160 0.170 0.145 0.019
## DP03_0045PE_2014 total_pop_2014 med_income_2014 pov_rate_2014
## 1 0.087 2309 16724 0.5184062
## 2 0.116 7138 19624 0.5057651
## 3 0.020 1539 14828 0.4938272
## 4 0.113 4515 16868 0.3936642
## 5 0.114 4250 15657 0.5590663
## 6 0.126 3720 26250 0.2814516
## renter_share_2014 med_gross_rent_2014 rent_burden_30p_2014 area_km2_2014
## 1 0.3803191 405 0.5734266 3.2793940
## 2 0.2465294 759 0.9223301 3.8475950
## 3 0.2673611 463 0.8831169 12.1343082
## 4 0.2773674 415 0.7578692 19.2464521
## 5 0.3036837 456 0.8076923 4.8014769
## 6 0.2200913 672 0.6597510 0.7108695
## pop_density_km2_2014 year_end_2014
## 1 704.0935 2014
## 2 1855.1848 2014
## 3 126.8305 2014
## 4 234.5887 2014
## 5 885.1443 2014
## 6 5233.0279 2014#지금 acs14에서 DP03 컬럼을 “pct_”로 다시 rename
#너가 이미 뽑아놓은 acs14, acs22에 바로 적용해.
acs14 <- acs14 %>%
dplyr::rename(
pct_construction_2014 = DP03_0034PE_2014,
pct_manufacturing_2014 = DP03_0035PE_2014,
pct_wholesale_2014 = DP03_0036PE_2014,
pct_retail_2014 = DP03_0037PE_2014,
pct_transport_2014 = DP03_0038PE_2014,
pct_information_2014 = DP03_0039PE_2014,
pct_finance_re_2014 = DP03_0040PE_2014,
pct_prof_sci_mgmt_2014 = DP03_0041PE_2014,
pct_edu_health_2014 = DP03_0042PE_2014,
pct_tourism_proxy_2014 = DP03_0043PE_2014,
pct_other_services_2014 = DP03_0044PE_2014,
pct_public_admin_2014 = DP03_0045PE_2014
)
acs22 <- acs22 %>%
dplyr::rename(
pct_construction_2022 = DP03_0034PE_2022,
pct_manufacturing_2022 = DP03_0035PE_2022,
pct_wholesale_2022 = DP03_0036PE_2022,
pct_retail_2022 = DP03_0037PE_2022,
pct_transport_2022 = DP03_0038PE_2022,
pct_information_2022 = DP03_0039PE_2022,
pct_finance_re_2022 = DP03_0040PE_2022,
pct_prof_sci_mgmt_2022 = DP03_0041PE_2022,
pct_edu_health_2022 = DP03_0042PE_2022,
pct_tourism_proxy_2022 = DP03_0043PE_2022,
pct_other_services_2022 = DP03_0044PE_2022,
pct_public_admin_2022 = DP03_0045PE_2022
)
# 확인
names(acs14)[grepl("^pct_", names(acs14))]## [1] "pct_construction_2014" "pct_manufacturing_2014"
## [3] "pct_wholesale_2014" "pct_retail_2014"
## [5] "pct_transport_2014" "pct_information_2014"
## [7] "pct_finance_re_2014" "pct_prof_sci_mgmt_2014"
## [9] "pct_edu_health_2014" "pct_tourism_proxy_2014"
## [11] "pct_other_services_2014" "pct_public_admin_2014"names(acs14)[grepl("tourism", names(acs14))]## [1] "pct_tourism_proxy_2014"#############################################################################
#1) Join: df_reg에 acs14 + acs22 붙이기
df_full <- df_reg %>%
dplyr::mutate(FIPS = as.character(FIPS)) %>%
dplyr::left_join(acs14, by = "FIPS") %>%
dplyr::left_join(acs22, by = "FIPS")
# join 잘 됐는지 체크
dplyr::glimpse(df_full)## Rows: 823
## Columns: 63
## $ STATEFP <chr> "72", "72", "72", "72", "72", "72", "72", "72"…
## $ COUNTYFP <chr> "005", "025", "025", "025", "031", "031", "043…
## $ TRACTCE <chr> "400900", "201000", "201500", "201300", "05090…
## $ AFFGEOID <chr> "1400000US72005400900", "1400000US72025201000"…
## $ GEOID <chr> "72005400900", "72025201000", "72025201500", "…
## $ NAME <chr> "4009", "2010", "2015", "2013", "509.02", "508…
## $ NAMELSAD <chr> "Census Tract 4009", "Census Tract 2010", "Cen…
## $ STUSPS <chr> "PR", "PR", "PR", "PR", "PR", "PR", "PR", "PR"…
## $ NAMELSADCO <chr> "Aguadilla Municipio", "Caguas Municipio", "Ca…
## $ STATE_NAME <chr> "Puerto Rico", "Puerto Rico", "Puerto Rico", "…
## $ LSAD <chr> "CT", "CT", "CT", "CT", "CT", "CT", "CT", "CT"…
## $ ALAND <dbl> 1120757, 1183020, 588614, 759722, 3693709, 106…
## $ AWATER <dbl> 0, 0, 0, 0, 111446, 182694, 0, 498516, 0, 0, 0…
## $ FIPS <chr> "72005400900", "72025201000", "72025201500", "…
## $ Change <dbl> 0.0703, 0.0201, 0.1968, 0.2645, 0.0553, 0.1219…
## $ Change_Direction <chr> "Increased (worse / more vulnerable)", "Increa…
## $ hotel_density <dbl> 0.8795381, 0.0000000, 0.0000000, 0.0000000, 0.…
## $ lodging_total <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ area_km2 <dbl> 1.1369604, 1.2036730, 0.5953259, 0.7520502, 3.…
## $ dist_to_coast_m <dbl> 593.2190, 22745.5561, 23769.3409, 22759.9076, …
## $ hotel_density_w <dbl> 0.8795381, 0.0000000, 0.0000000, 0.0000000, 0.…
## $ pct_construction_2014 <dbl> 0.138, 0.000, 0.014, 0.124, 0.075, 0.022, 0.02…
## $ pct_manufacturing_2014 <dbl> 0.140, 0.112, 0.089, 0.102, 0.022, 0.032, 0.17…
## $ pct_wholesale_2014 <dbl> 0.067, 0.030, 0.037, 0.041, 0.000, 0.041, 0.01…
## $ pct_retail_2014 <dbl> 0.190, 0.250, 0.202, 0.211, 0.096, 0.184, 0.13…
## $ pct_transport_2014 <dbl> 0.015, 0.000, 0.000, 0.000, 0.071, 0.071, 0.01…
## $ pct_information_2014 <dbl> 0.031, 0.000, 0.014, 0.000, 0.023, 0.032, 0.00…
## $ pct_finance_re_2014 <dbl> 0.000, 0.043, 0.085, 0.056, 0.083, 0.079, 0.03…
## $ pct_prof_sci_mgmt_2014 <dbl> 0.079, 0.051, 0.102, 0.076, 0.090, 0.104, 0.05…
## $ pct_edu_health_2014 <dbl> 0.090, 0.100, 0.220, 0.199, 0.217, 0.197, 0.28…
## $ pct_tourism_proxy_2014 <dbl> 0.155, 0.222, 0.072, 0.056, 0.144, 0.121, 0.09…
## $ pct_other_services_2014 <dbl> 0.000, 0.120, 0.028, 0.069, 0.015, 0.021, 0.03…
## $ pct_public_admin_2014 <dbl> 0.073, 0.072, 0.128, 0.067, 0.164, 0.097, 0.06…
## $ total_pop_2014 <dbl> 2349, 2230, 2931, 2869, 5040, 4451, 2314, 2524…
## $ med_income_2014 <dbl> 10079, 11869, 23865, 25313, 29205, 29427, 1519…
## $ pov_rate_2014 <dbl> 0.74201788, 0.66995516, 0.49451679, 0.33321678…
## $ renter_share_2014 <dbl> 0.55845630, 0.85101580, 0.40265907, 0.34714004…
## $ med_gross_rent_2014 <dbl> 199, 482, 557, 625, 537, 532, 481, 248, 1220, …
## $ rent_burden_30p_2014 <dbl> 0.7052846, 0.6618037, 0.7877358, 0.6164773, 0.…
## $ area_km2_2014 <dbl> 1.1411410, 1.2098756, 0.5938731, 0.7567413, 3.…
## $ pop_density_km2_2014 <dbl> 2058.4661, 1843.1646, 4935.3982, 3791.2558, 13…
## $ year_end_2014 <dbl> 2014, 2014, 2014, 2014, 2014, 2014, 2014, 2014…
## $ pct_construction_2022 <dbl> 0.000, 0.040, 0.044, 0.000, 0.054, 0.029, 0.11…
## $ pct_manufacturing_2022 <dbl> 0.121, 0.075, 0.031, 0.042, 0.049, 0.038, 0.10…
## $ pct_wholesale_2022 <dbl> 0.048, 0.085, 0.064, 0.025, 0.026, 0.007, 0.00…
## $ pct_retail_2022 <dbl> 0.075, 0.197, 0.260, 0.027, 0.095, 0.165, 0.02…
## $ pct_transport_2022 <dbl> 0.000, 0.000, 0.067, 0.053, 0.092, 0.060, 0.01…
## $ pct_information_2022 <dbl> 0.000, 0.000, 0.000, 0.025, 0.026, 0.030, 0.00…
## $ pct_finance_re_2022 <dbl> 0.000, 0.030, 0.014, 0.028, 0.069, 0.091, 0.04…
## $ pct_prof_sci_mgmt_2022 <dbl> 0.293, 0.125, 0.084, 0.142, 0.047, 0.053, 0.08…
## $ pct_edu_health_2022 <dbl> 0.223, 0.218, 0.192, 0.344, 0.251, 0.279, 0.14…
## $ pct_tourism_proxy_2022 <dbl> 0.203, 0.128, 0.083, 0.167, 0.145, 0.110, 0.12…
## $ pct_other_services_2022 <dbl> 0.000, 0.101, 0.072, 0.037, 0.072, 0.025, 0.10…
## $ pct_public_admin_2022 <dbl> 0.036, 0.000, 0.089, 0.089, 0.065, 0.108, 0.16…
## $ total_pop_2022 <dbl> 2522, 1639, 2109, 2640, 6101, 4497, 1822, 1784…
## $ med_income_2022 <dbl> 10263, 13447, 26016, 24779, 36454, 37754, 2072…
## $ pov_rate_2022 <dbl> 0.62410785, 0.66076876, 0.39731286, 0.40660592…
## $ renter_share_2022 <dbl> 0.5263158, 0.7765727, 0.3737673, 0.3785366, 0.…
## $ med_gross_rent_2022 <dbl> 220, 431, 452, 532, 827, 513, 501, 237, 2078, …
## $ rent_burden_30p_2022 <dbl> 0.7607143, 0.6117318, 0.7810026, 0.6134021, 0.…
## $ area_km2_2022 <dbl> 1.1411410, 1.2098756, 0.5938731, 0.7567413, 3.…
## $ pop_density_km2_2022 <dbl> 2210.0687, 1354.6847, 3551.2640, 3488.6425, 16…
## $ year_end_2022 <dbl> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022…summary(is.na(df_full$pct_tourism_proxy_2014))## Mode FALSE
## logical 823summary(is.na(df_full$pct_tourism_proxy_2022))## Mode FALSE
## logical 823## 컨트롤 포함 회귀 (baseline 2014 컨트롤)
m_full <- lm(
Change ~
hotel_density_w +
dist_to_coast_m +
pop_density_km2_2014 +
med_income_2014 +
pov_rate_2014 +
renter_share_2014 +
med_gross_rent_2014 +
rent_burden_30p_2014 +
pct_tourism_proxy_2014 +
pct_manufacturing_2014 +
pct_construction_2014 +
pct_public_admin_2014,
data = df_full
)
summary(m_full)##
## Call:
## lm(formula = Change ~ hotel_density_w + dist_to_coast_m + pop_density_km2_2014 +
## med_income_2014 + pov_rate_2014 + renter_share_2014 + med_gross_rent_2014 +
## rent_burden_30p_2014 + pct_tourism_proxy_2014 + pct_manufacturing_2014 +
## pct_construction_2014 + pct_public_admin_2014, data = df_full)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.7143 -0.1327 -0.0034 0.1326 0.6720
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.456e-01 9.637e-02 5.662 2.09e-08 ***
## hotel_density_w -8.311e-03 4.530e-03 -1.835 0.066923 .
## dist_to_coast_m -4.360e-06 1.125e-06 -3.874 0.000116 ***
## pop_density_km2_2014 9.988e-06 4.391e-06 2.275 0.023188 *
## med_income_2014 -6.530e-06 1.686e-06 -3.874 0.000116 ***
## pov_rate_2014 -9.185e-01 1.085e-01 -8.468 < 2e-16 ***
## renter_share_2014 2.678e-01 7.713e-02 3.472 0.000544 ***
## med_gross_rent_2014 8.270e-06 5.710e-05 0.145 0.884877
## rent_burden_30p_2014 -8.790e-02 7.089e-02 -1.240 0.215379
## pct_tourism_proxy_2014 1.087e-01 1.651e-01 0.659 0.510288
## pct_manufacturing_2014 1.007e-01 1.462e-01 0.689 0.491119
## pct_construction_2014 1.423e-01 1.952e-01 0.729 0.466252
## pct_public_admin_2014 -8.819e-02 1.731e-01 -0.509 0.610584
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2118 on 798 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.2014, Adjusted R-squared: 0.1894
## F-statistic: 16.77 on 12 and 798 DF, p-value: < 2.2e-16# 0) baseline SVI table 준비 (FIPS + Overall_2014 only)
baseline_svi <- comparison %>%
transmute(
FIPS = as.character(FIPS),
Overall_2014 = as.numeric(Overall_2014)
)
# 1) df_full에 "baseline 2014 Overall SVI" Added! ***********
df_full <- df_full %>%
mutate(FIPS = as.character(FIPS)) %>%
left_join(baseline_svi, by = "FIPS")
# 2) 확인
"Overall_2014" %in% names(df_full)## [1] TRUEsummary(df_full$Overall_2014)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0011 0.2574 0.5000 0.5029 0.7494 1.0000m_full2 <- lm(
Change ~
hotel_density_w + dist_to_coast_m +
Overall_2014 +
pop_density_km2_2014 + med_income_2014 + pov_rate_2014 +
renter_share_2014 + med_gross_rent_2014 + rent_burden_30p_2014 +
pct_tourism_proxy_2014 + pct_manufacturing_2014 +
pct_construction_2014 + pct_public_admin_2014,
data = df_full
)
summary(m_full2)##
## Call:
## lm(formula = Change ~ hotel_density_w + dist_to_coast_m + Overall_2014 +
## pop_density_km2_2014 + med_income_2014 + pov_rate_2014 +
## renter_share_2014 + med_gross_rent_2014 + rent_burden_30p_2014 +
## pct_tourism_proxy_2014 + pct_manufacturing_2014 + pct_construction_2014 +
## pct_public_admin_2014, data = df_full)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.64279 -0.12606 -0.00341 0.12945 0.57072
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.077e-01 8.927e-02 5.687 1.81e-08 ***
## hotel_density_w -1.135e-02 4.202e-03 -2.701 0.00705 **
## dist_to_coast_m -2.937e-06 1.049e-06 -2.799 0.00525 **
## Overall_2014 -5.209e-01 4.498e-02 -11.581 < 2e-16 ***
## pop_density_km2_2014 1.252e-05 4.071e-06 3.076 0.00217 **
## med_income_2014 -7.266e-06 1.562e-06 -4.652 3.84e-06 ***
## pov_rate_2014 -2.245e-01 1.169e-01 -1.920 0.05524 .
## renter_share_2014 3.029e-01 7.147e-02 4.238 2.52e-05 ***
## med_gross_rent_2014 -2.882e-05 5.295e-05 -0.544 0.58649
## rent_burden_30p_2014 -1.108e-01 6.566e-02 -1.687 0.09190 .
## pct_tourism_proxy_2014 1.006e-01 1.528e-01 0.658 0.51047
## pct_manufacturing_2014 -1.329e-02 1.357e-01 -0.098 0.92203
## pct_construction_2014 3.130e-01 1.813e-01 1.727 0.08464 .
## pct_public_admin_2014 -9.221e-02 1.603e-01 -0.575 0.56522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1961 on 797 degrees of freedom
## (12 observations deleted due to missingness)
## Multiple R-squared: 0.3164, Adjusted R-squared: 0.3053
## F-statistic: 28.38 on 13 and 797 DF, p-value: < 2.2e-16##m_full2 : The regression results indicate that, after controlling for baseline vulnerability and socioeconomic characteristics, hotel density is negatively associated with changes in social vulnerability (β = −0.011, p = 0.007), suggesting that areas with higher hotel density tend to experience smaller increases—or slight decreases—in SVI over time. Coastal proximity is also significantly related to vulnerability change (β = −2.94e−06, p = 0.005), indicating that tracts closer to the coastline tend to experience greater increases in vulnerability. Baseline vulnerability (SVI 2014) shows a strong negative association with change (β = −0.521, p < 0.001), suggesting convergence effects where initially more vulnerable areas tend to experience decreases in SVI over time.
##baseline 포함한 mediation (권장 코드)
library(mediation)
set.seed(123)
df_med2 <- df_full %>%
dplyr::filter(
!is.na(Change),
!is.na(hotel_density_w),
!is.na(dist_to_coast_m),
!is.na(Overall_2014)
)
med.fit2 <- lm(dist_to_coast_m ~ hotel_density_w + Overall_2014, data = df_med2)
out.fit2 <- lm(Change ~ hotel_density_w + dist_to_coast_m + Overall_2014, data = df_med2)
med.out2 <- mediate(
model.m = med.fit2,
model.y = out.fit2,
treat = "hotel_density_w",
mediator = "dist_to_coast_m",
sims = 2000,
boot = TRUE
)## Running nonparametric bootstrapsummary(med.out2)##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME 0.0038582 0.0023442 0.0056569 <2e-16 ***
## ADE -0.0016882 -0.0086455 0.0059217 0.635
## Total Effect 0.0021699 -0.0044478 0.0096945 0.562
## Prop. Mediated 1.7780096 -16.5195449 13.5621255 0.562
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 823
##
##
## Simulations: 2000#Final result: Mediation analysis further shows a significant indirect effect of hotel density on vulnerability change through coastal proximity (ACME = 0.00386, p < 0.001). However, the direct effect (ADE) and the total effect are not statistically significant, indicating that hotel density does not have a significant overall effect on vulnerability change but may operate indirectly through spatial coastal dynamics.