Change in a county-level ACS variable between two years

if (!require("tidyverse")) install.packages("tidyverse")
if (!require("tidycensus")) install.packages("tidycensus")
if (!require("sf")) install.packages("sf")

library(tidycensus)
library(tidyverse)
library(sf)
options(tigris_use_cache = TRUE)

# NOTE: Obtain a Census API key from 
# https://api.census.gov/data/key_signup.html
# and replace the census_api_key code line's 
# PastYourAPIKeyBetwenTheseQuoteMarks text with your
# API key. Also, delete the "#" from the front of the 
# census_api_key code.

#census_api_key("PasteYourAPIKeyBetweenTheseQuoteMarks")

# Browse the ACS variables, if you like
# Change "2021" to your preferred year
# Change "acs1" to "acs5" to retrieve the five-year ACS variables
DetailedTables_Y1 <- load_variables(2016, "acs1", cache = TRUE)
SubjectTables_Y1 <- load_variables(2016, "acs1/subject", cache = TRUE)
ProfileTables_Y1 <- load_variables(2016, "acs1/profile", cache = TRUE)
Codebook_Y1 <- rbind(DetailedTables_Y1,SubjectTables_Y1,ProfileTables_Y1)
rm(DetailedTables_Y1,SubjectTables_Y1,ProfileTables_Y1)

DetailedTables_Y2 <- load_variables(2021, "acs1", cache = TRUE)
SubjectTables_Y2 <- load_variables(2021, "acs1/subject", cache = TRUE)
ProfileTables_Y2 <- load_variables(2021, "acs1/profile", cache = TRUE)
Codebook_Y2 <- rbind(DetailedTables_Y2,SubjectTables_Y2,ProfileTables_Y2)
rm(DetailedTables_Y2,SubjectTables_Y2,ProfileTables_Y2)


# Getting the Year 1 data
Data_Y1 <- get_acs(geography = "county",
                       state = "TN",
                       variables = c(Var_Y1 = "B25031_001"),
                       year = 2016,
                       survey = "acs1",
                       output = "wide",
                       geometry = TRUE)

# Getting the Year 2 data
Data_Y2 <- get_acs(geography = "county",
                       state = "TN",
                       variables = c(Var_Y2 = "B25031_001"),
                       year = 2021,
                       survey = "acs1",
                       output = "wide",
                       geometry = FALSE)

# Merging the Year 1 and Year 2 data, using the GEOID variable
# in each dataset as a key
Mergeddata <- merge(Data_Y1, Data_Y2, by = "GEOID")
rm(Data_Y1,Data_Y2)

# Making pretty column names for the map file
Mergeddata$County <- Mergeddata$NAME.x
Mergeddata$Change <- Mergeddata$Var_Y2E - Mergeddata$Var_Y1E

# Testing the significance of each county's change

Mergeddata$Sig <- significance(Mergeddata$Var_Y1E,
                                Mergeddata$Var_Y2E,
                                Mergeddata$Var_Y1M,
                                Mergeddata$Var_Y2M,
                                 clevel = 0.90)
Mergeddata$Significance <- ifelse(Mergeddata$"Sig" == "TRUE",
                                  "Significant",
                                  "Nonsignificant")


# Making a preview map
Mergeddata %>%
  ggplot(aes(fill = Change)) + 
  geom_sf(color = NA) + 
  scale_fill_viridis_c(option = "magma") 

#Calculating error margin ranges
Mergeddata$From = Mergeddata$Var_Y2E-Mergeddata$Var_Y2M
Mergeddata$To = Mergeddata$Var_Y2E+Mergeddata$Var_Y2M

#Sorting data by change, in descending order
Mergeddata <- Mergeddata[order(-Mergeddata$Change),]

#Pruning data frame for .kml file
Mergeddata <- Mergeddata[colSums(!is.na(Mergeddata)) > 0]
Mergeddata <- subset(Mergeddata, select = -c(NAME.x,NAME.y,Sig))

#Adding meaningful variable names (optional)
#Mergeddata$CustomVariableNameY1 <- Mergeddata$Var_Y1E
#Mergeddata$CustomVariableNameY2 <- Mergeddata$Var_Y2E

#Making a .kml
library(sf)
st_write(Mergeddata,"ChangeMap.kml", append = FALSE)

#Pruning data frame for .csv file
Mergeddata <- st_drop_geometry(Mergeddata)

#Making a .csv file
write.csv(Mergeddata,"ChangeCSV.csv")  
  
##########################################################################

if (!require("plotly")) install.packages("plotly")
library(plotly)
if (!require("stringr")) install.packages("stringr")
library(stringr)

Mergeddata$Place <- word(Mergeddata$County, 1)

plot_ly(Mergeddata,
        x = ~Place,
        y = ~Median_Rent,
        type = 'bar',
        name = 'Median Rent',
        error_y = ~list(array = Var_Y2M,
                        color = '#000000')) %>% 
  layout(xaxis = list(categoryorder = "total descending"))