DATA 698-Research Project

#---User-defined function(s)---#

#Adapted correlation matrix used in EDA section:
plot_corr_matrix <- function(dataframe, significance_threshold){
  title <- paste0('Correlation Matrix for significance > ',
                  significance_threshold)
  
  df_cor <- dataframe %>% mutate_if(is.character, as.factor)
  
  df_cor <- df_cor %>% mutate_if(is.factor, as.numeric)
  #run a correlation and drop the insignificant ones
  corr <- cor(df_cor, use = 'na.or.complete')
  #prepare to drop duplicates and correlations of 1     
  corr[lower.tri(corr,diag=TRUE)] <- NA 
  #drop perfect correlations
  corr[corr == 1] <- NA 
  #turn into a 3-column table
  corr <- as.data.frame(as.table(corr))
  #remove the NA values from above 
  corr <- na.omit(corr) 
  #select significant values  
  corr <- subset(corr, abs(Freq) > significance_threshold) 
  #sort by highest correlation
  corr <- corr[order(-abs(corr$Freq)),] 
  #print table
  # print(corr)
  #turn corr back into matrix in order to plot with corrplot
  mtx_corr <- reshape2::acast(corr, Var1~Var2, value.var="Freq")
  
  #plot correlations visually
  corrplot(mtx_corr,
           title=title,
           mar=c(0,0,1,0),
           method='color', 
           tl.col="black", 
           na.label= " ",
           addCoef.col = 'black',
           number.cex = .9)
}

#Evaluation metrics used in Model Building / Selection section:
eval_metrics = function(model, df, predictions, target){
  resids = df[,target] - predictions
  resids2 = resids**2
  N = length(predictions)
  r2 = as.character(round(summary(model)$r.squared, 4))
  adj_r2 = as.character(round(summary(model)$adj.r.squared, 4))
  print(adj_r2) #Adjusted R-squared
  print(as.character(round(sqrt(sum(resids2)/N), 4))) #RMSE
}

## --- Model Selection --- ##

#Create Kable table to succinctly summarize variable research
Source <- c('Bhardwaj, Amiri, Buchwald, and Amram', 'Kaminsky, Lessler, Sharfstein, and Wallace', 'Azevedo, Hallal, Wells, and Victoria', 'An, Li, and Jiang', 'Holick and Wacker', 'Bremberg', 'Ahern, Brown, and Dukas', 'Ferreira, Jardim, Jardim, Sousa, and Rosa', 'Islam, Leer, Teo, et al.', 'Hammarstrom, Janlert, and Winefield')
Type <- c('Dependent', 'Dependent', 'Dependent', 'Independent', 'Independent', 'Independent', 'Independent', 'Independent', 'Independent', 'Independent')
Variable <- c('Longevity', 'Obesity', 'Physical Activity', 'EQI (Environmental Quality)', 'Sunlight', 'Population Density', 'Nutrition', 'Alcohol Consumption', 'Income', 'Unemployment')

LRVS <- cbind(Source, Type, Variable)

LRVS %>%
  kbl(caption = "Literature Review: Variable Summary") %>%
  kable_minimal() %>%
  kable_styling(latex_options = "hold_position")

## --- Dependent Variable Creation --- ##

#Read in life expectancy data, convert to tibble, and select pertinent columns:
longevity <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/IHME_LifeExpectancy.csv")
life_table <- as_tibble(longevity)
life_table <- life_table %>% dplyr::select(1:2,13:16) %>% na.omit()

#Read in obesity data, convert to tibble, and select pertinent columns:
obesity <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/IHME_Obesity.csv")
obesity_table <- as_tibble(obesity)
obesity_table <- obesity_table %>% dplyr::select(1:2,5:6,9:10) %>% na.omit()

#Read in physical activity data, convert to tibble, and select pertinent columns:
activity <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/IHME_PhysicalActivity.csv")
act_table <- as_tibble(activity)
act_table <- act_table %>% dplyr::select(1:2,5:6,9:10) %>% na.omit()

###LIFE EXPECTANCY DATA: exploration, normalization, and compilation

#Explore life expectancy data at a county level:
glimpse(life_table)

summary(life_table)

#Extract variables of interest
m1 <- life_table$`Male life expectancy, 2010 (years)`
f1 <- life_table$`Female life expectancy, 2010 (years)`
dm1 <- life_table$`Difference in male life expectancy, 1985-2010 (years)`
df1 <- life_table$`Difference in female life expectancy, 1985-2010 (years)`

#Normalize data scale to be from 0 to 1
n_m1 = (m1-min(m1))/(max(m1)-min(m1))
n_f1 = (f1-min(f1))/(max(f1)-min(f1))
n_dm1 = (dm1-min(dm1))/(max(dm1)-min(dm1))
n_df1 = (df1-min(df1))/(max(df1)-min(df1))

#Histogram of original vs. normalized data
##Life expectancy histograms
par(mfrow=c(2,2))
hist(m1, breaks=10, xlab="Age (years)", col="lightblue", main="Male life expectancy, 2010")
hist(n_m1, breaks=10, xlab="Normalized Age (years)", col="lightblue", main="Male life expectancy, 2010")
hist(f1, breaks=10, xlab="Age (years)", col="lightblue", main="Female life expectancy, 2010")
hist(n_f1, breaks=10, xlab="Normalized Age (years)", col="lightblue", main="Female life expectancy, 2010")

##Longevity improvement histograms
par(mfrow=c(2,2))
hist(dm1, breaks=10, xlab="Age (years)", col="lightblue", main="Male longevity improvement, 1985-2010")
hist(n_dm1, breaks=10, xlab="Normalized Age (years)", col="lightblue", main="Male longevity improvement, 1985-2010")
hist(df1, breaks=10, xlab="Age (years)", col="lightblue", main="Female longevity improvement, 1985-2010")
hist(n_df1, breaks=10, xlab="Normalized Age (years)", col="lightblue", main="Female longevity improvement, 1985-2010")

#Add normalized variables together
life <- n_m1 + n_dm1 + n_f1 + n_df1

#Normalize activity to 0-1 range
n_life = (life-min(life))/(max(life)-min(life))
#head(n_life)

#Histogram of original vs. normalized data
#par(mfrow=c(1,2))
#hist(life, breaks=10, xlab="Score", col="lightblue", main="Longevity metric")
hist(n_life, breaks=10, xlab="Normalized Score", col="lightblue", main="Longevity metric")
summary(n_life) #slight left skew

###OBESITY DATA: exploration, normalization, and compilation

#Explore obesity data at a county level:
glimpse(obesity_table)

summary(obesity_table)

#Extract variables of interest
m2 <- obesity_table$`Male obesity prevalence, 2009 (%)` 
f2 <- obesity_table$`Female obesity  prevalence, 2009 (%)` 
dm2 <- obesity_table$`Difference in male obesity prevalence, 2001-2009 (percentage points)` 
df2 <- obesity_table$`Difference in female obesity prevalence, 2001-2009 (percentage points)`

#Normalize
n_m2 = (m2-min(m2))/(max(m2)-min(m2))
n_f2 = (f2-min(f2))/(max(f2)-min(f2))
n_dm2 = (dm2-min(dm2))/(max(dm2)-min(dm2))
n_df2 = (df2-min(df2))/(max(df2)-min(df2))

#Histogram of original vs. normalized data
par(mfrow=c(2,2))

hist(m2, breaks=10, xlab="Obesity rate (%)", col="lightblue", main="Male obesity prevalence, 2009")
hist(n_m2, breaks=10, xlab="Normalized obesity rate (%)", col="lightblue", main="Male obesity prevalence, 2009")
hist(f2, breaks=10, xlab="Obesity rate (%)", col="lightblue", main="Female obesity prevalence, 2009")
hist(n_f2, breaks=10, xlab="Normalized obesity rate (%)", col="lightblue", main="Female obesity prevalence, 2009")

par(mfrow=c(2,2))
hist(dm2, breaks=10, xlab="Obesity rate (%)", col="lightblue", main="Male obesity increase, 2001-2009")
hist(n_dm2, breaks=10, xlab="Normalized obesity rate (%)", col="lightblue", main="Male obesity increase, 2001-2009")
hist(df2, breaks=10, xlab="Obesity rate (%)", col="lightblue", main="Female obesity increase, 2001-2009")
hist(n_df2, breaks=10, xlab="Normalized obesity rate (%)", col="lightblue", main="Female obesity increase, 2001-2009")

#Add normalized variables together
fat <- n_m2 + n_dm2 + n_f2 + n_df2

#Normalize activity to 0-1 range
n_fat = (fat-min(fat))/(max(fat)-min(fat))

#head(n_fat)

# Histogram of original vs. normalized data
#par(mfrow=c(1,2))
#hist(fat, breaks=10, xlab="Score", col="lightblue", main="Obesity metric")
hist(n_fat, breaks=10, xlab="Normalized Score", col="lightblue", main="Obesity metric")
summary(n_fat) #right skewed

###PHYSICAL ACTIVITY DATA: exploration, normalization, and compilation

glimpse(act_table)

summary(act_table)

#Explore and normalize male physical activity data
m3 <- act_table$`Male sufficient physical activity  prevalence, 2009 (%)`
f3 <- act_table$`Female sufficient physical activity  prevalence, 2009 (%)`
dm3 <- act_table$`Difference in male sufficient physical activity prevalence, 2001-2009 (percentage points)`  
df3 <- act_table$`Difference in female sufficient physical activity prevalence, 2001-2009 (percentage points)` 

#Normalized Data
n_m3 = (m3-min(m3))/(max(m3)-min(m3))
n_f3 = (f3-min(f3))/(max(f3)-min(f3))
n_dm3 = (dm3-min(dm3))/(max(dm3)-min(dm3))
n_df3 = (df3-min(df3))/(max(df3)-min(df3))

#Histogram of original vs. normalized data
par(mfrow=c(2,2))

hist(m3, breaks=10, xlab="Physical activity rate (%)", col="lightblue", main="Male activity prevalence, 2009")
hist(n_m3, breaks=10, xlab="Normalized physical activity rate (%)", col="lightblue", main="Male activity prevalence, 2009")
hist(f3, breaks=10, xlab="Physical activity rate (%)", col="lightblue", main="Female activity prevalence, 2009")
hist(n_f3, breaks=10, xlab="Normalized obesity rate (%)", col="lightblue", main="Female activity prevalence, 2009")

par(mfrow=c(2,2))
hist(dm3, breaks=10, xlab="Physical activity rate (%)", col="lightblue", main="Male activity difference, 2001-2009")
hist(n_dm3, breaks=10, xlab="Normalized physical activity rate (%)", col="lightblue", main="Male activity difference, 2001-2009")
hist(df3, breaks=10, xlab="Physical activity rate (%)", col="lightblue", main="Female activity difference, 2009")
hist(n_df3, breaks=10, xlab="Normalized obesity rate (%)", col="lightblue", main="Female activity difference, 2009")

#Add normalized variables together
active <- n_m3 + n_dm3 + n_f3 + n_df3

#Normalize activity to 0-1 range
n_active = (active-min(active))/(max(active)-min(active))

#head(n_active)

# Histogram of original vs. normalized data
#par(mfrow=c(1,2))
#hist(active, breaks=10, xlab="Score", col="lightblue", main="Physical activity metric")
hist(n_active, breaks=10, xlab="Normalized Score", col="lightblue", main="Physical activity metric")

summary(n_active) #slight right skew

###DEPENDENT VARIABLE CREATION: sum, normalize, and plot

#Add normalized variables together
lifestyle <- n_life - n_fat + n_active

#Normalize health to 0-1 range
normalized_lifestyle = (lifestyle-min(lifestyle))/(max(lifestyle)-min(lifestyle))

#Histogram of original vs. normalized data
#par(mfrow=c(1,2))
#hist(lifestyle, breaks=10, xlab="Score", col="lightblue", main="Health metric")
hist(normalized_lifestyle, breaks=10, xlab="Normalized Score", col="lightblue", main="Health metric")

summary(normalized_lifestyle)

#head(normalized_lifestyle)

#create new df with state | county | health score
starter_df <- life_table %>% 
  dplyr::select(1:2)

starter_df$health_score <- normalized_lifestyle
healthiest_counties <- filter(starter_df, `health_score` > 0.895) #top 10
healthiest_counties <- healthiest_counties[order(-healthiest_counties$`health_score`),] #descending order

#head(starter_df)
#nrow(healthiest_counties) #10
healthiest_counties %>%
  kbl() %>%
  kable_minimal() %>%
  kable_styling(latex_options = "hold_position")

## --- Independent Variable Preprocessing --- ##

#Read in alcohol data, convert to tibble, and drop impertinent observation:
alcohol <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/AlcoholConsumption.csv")
alcohol <- as_tibble(alcohol)
alcohol <- alcohol[-1,] #drop State = National
#dim(alcohol) #3178 x 6

alcohol <- alcohol[alcohol$Location != alcohol$State,] #drop state observations from Location column
#dim(alcohol) #3178 - 3127 = 51 observations dropped

#rename columns
alcohol <- alcohol %>% rename( 
  County = Location,
  Hvy = Hvy_2012,
  Bng = Binge_2012,
  HvyPctChg = `HvyPctChg_2005-2012`,
  BngPctChg = `BingePctChg_2005-2012`)

#drop excess verbage from County column
stopwords <- c("and", "Area", "Borough", "Census", "City", "County", "Division", "Municipality", "Parish")
alcohol$County <- removeWords(alcohol$County, stopwords)
#head(alcohol)

#Read in heart data, convert to tibble, and drop impertinent observation:
heart <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/CardiovascularDisease.csv")
heart <- as_tibble(heart)
heart <- subset(heart, select = -c(`Mortality Rate, 2010*`)) #drop 2010
heart <- heart[-1,] #drop State = National
#dim(heart) #3193 x 4

# remove State == Location
heart <- heart[heart$Location != heart$State,]
dim(heart) #3193 - 3143 = 50 observations removed

# retitle columns
heart <- heart %>% rename( 
  County = Location,
  Mortality_2005 = `Mortality Rate, 2005*`,
  Mortality_2014 = `Mortality Rate, 2014*`)

# retain value EXCLUSIVELY for Mortality Rate columns
heart$Mortality_2005 <- gsub("\\s*\\([^\\)]+\\)","",as.character(heart$Mortality_2005))
heart$Mortality_2014 <- gsub("\\s*\\([^\\)]+\\)","",as.character(heart$Mortality_2014))

#convert columns to proper type
heart$Mortality_2005 <- as.double(heart$Mortality_2005)
heart$Mortality_2014 <- as.double(heart$Mortality_2014)

#drop excess verbage from County column
heart$County <- removeWords(heart$County, stopwords)
heart$County <- gsub("(.*),.*", "\\1", heart$County) #remove everything after comma

# add Chg column
heart$MortalityChg <- heart$Mortality_2014 - heart$Mortality_2005

#finalize format of df
heart <- subset(heart, select = -c(`Mortality_2005`)) #drop 2005
heart <- heart %>% rename( Mortality = Mortality_2014)

#head(heart)

#Read in education data, convert to tibble, and drop impertinent observation:
education <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/Education.csv")
education <- as_tibble(education)
education <- education[-1,] #drop State = National

education$State <- state.name[match(education$State, state.abb)] #convert state abbreviation to name

education <- education[education$`Area name` != education$State,] #drop state observations from Area name column
#dim(education) #3281 - 3233 = 48 observations dropped

#rename columns
education <- education %>% rename(
  County = `Area name`,
  LTHighSchool = `Percent of adults with less than a high school diploma, 2015-19`,
  HighSchool = `Percent of adults with a high school diploma only, 2015-19`,
  SomeCollege = `Percent of adults completing some college or associate's degree, 2015-19`,
  College = `Percent of adults with a bachelor's degree or higher, 2015-19`)

#drop excess verbage from County column
education$County <- removeWords(education$County, stopwords)

#head(education) #verify

#Read in eqi data and convert to tibble:
eqi <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/EnvironmentalQualityIndex.csv")
eqi <- as_tibble(eqi)
eqi <- subset(eqi, select = -c(3:7)) #drop indices that makeup EQI score

eqi$State <- state.name[match(eqi$State, state.abb)] #convert state abbreviation to name

#rename columns
eqi <- eqi %>% rename(
  County = County_Name,
  EQI = environmental_quality_index)

#drop excess verbage from County column
eqi$County <- removeWords(eqi$County, stopwords)

#head(eqi) #verify
#dim(eqi) #3281 x 6

#Read in food insecurity data and convert to tibble:
food <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/FoodInsecurity.csv")
food <- as_tibble(food)
#head(food)

#drop FIPS
food <- subset(food, select=-c(FIPS))
#dim(food) #3142 x 3: no need to drop observations

#convert State to full name
food$State <- state.name[match(food$State, state.abb)] #convert state abbreviation to name

#rename columns
food <- food %>% rename(
  County = `County, State`,
  FoodInsecurity = `2018 Food Insecurity Rate`)

#remove excess verbage from County
food$County <- removeWords(food$County, stopwords)
food$County <- gsub("(.*),.*", "\\1", food$County) #remove everything after comma

#drop % from Food Insecurity and convert to double
food$FoodInsecurity = as.double(gsub("[\\%,]", "", food$FoodInsecurity))

#head(food) #verify

#Read in sun data and convert to tibble:
#sun <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/Sunlight.csv")
#sun <- as_tibble(sun)
#dim(sun) #3161 x 3



# Read in sun data and convert to tibble
sun <- readr::read_csv(
  "https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/Sunlight.csv",
  show_col_types = FALSE,
  locale = readr::locale(encoding = "UTF-8"),
  name_repair = "minimal"
)

# Remove any hidden spaces from column names
names(sun) <- trimws(names(sun))

# Check column names
print(names(sun))

# Check dimensions
dim(sun)

#drop excess verbage from County column
sun$County <- removeWords(sun$County, stopwords)
sun$County <- gsub("(.*),.*", "\\1", sun$County) #remove everything after comma

#head(sun) #verify

#Read in une ployment data and convert to tibble:
unemp <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/Unemployment.csv")
unemp <- as_tibble(unemp)
unemp <- unemp[-1,] #drop State = National
unemp <- subset(unemp, select=-c(3)) #drop 2016

unemp$State <- state.name[match(unemp$Stabr, state.abb)] #convert state abbreviation to name
unemp <- unemp[unemp$area_name != unemp$State,] #drop state observations from Area name column

unemp <- unemp %>% rename(
  County = area_name,
  Unemployment = Unemployment_rate_2019,
  UnemploymentChg = `Unemployment_chg_2016-2019`) #rename columns

unemp$County <- removeWords(unemp$County, stopwords) #drop excess verbage from County column
unemp$County <- gsub("(.*),.*", "\\1", unemp$County) #remove everything after comma

#dim(unemp) #3274 - 3224 = 50 dropped
#head(unemp) #verify


#Read in wealth data and convert to tibble:
# Read wealth data
wealth <- readr::read_csv(
  "https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/Wealth.csv",
  show_col_types = FALSE
)

wealth <- as_tibble(wealth)

# Remove national observation
wealth <- wealth[-1, ]

# Convert state abbreviations to full names
wealth$State <- state.name[match(wealth$State, state.abb)]

# Remove state-level observations
wealth <- wealth[wealth$County != wealth$State, ]

# Convert columns to numeric
wealth$PovertyRate <- as.numeric(wealth$PovertyRate)
wealth$Income <- as.numeric(gsub(",", "", wealth$Income))

# Rename PovertyRate
wealth <- wealth %>%
  dplyr::rename(Poverty = PovertyRate)

# Clean county names
wealth$County <- tm::removeWords(wealth$County, stopwords)
wealth$County <- trimws(wealth$County)

# Verify
names(wealth)

class(wealth$Income)

#Read in population data and convert to tibble:
pop <- read_csv("https://raw.githubusercontent.com/BIKASHBHOWMIK15/Data-698/main/population.csv")
pop <- as_tibble(pop)

# rename columns
pop <- pop %>% rename(
  State = STNAME,
  County = CTYNAME,
  Pop_2010 = CENSUS2010POP,
  Population = POPESTIMATE2019,
  Births = BIRTHS2019,
  Deaths = DEATHS2019,
  NetMig = NETMIG2019) #rename columns

#add population change variable
pop$PopChg <- pop$Population - pop$Pop_2010

pop <- subset(pop, select=-c(3)) #drop 2010

#dim(pop) #3193 x 7
pop <- pop[pop$County != pop$State,] #drop state observations from Area name column
#dim(pop) #3193 - 3141 = 52 observations dropped

pop$County[1802] <- "Dona Ana County" #invalid UTF-8
pop$County <- removeWords(pop$County, stopwords) #drop excess verbage from County column

#head(pop)

## --- Merge df's --- ##

##1. merge health score and alcohol df's
#Trim white space
alcohol$County <- trimws(alcohol$County)
starter_df$County <- trimws(starter_df$County)
#SQL join
df <- sqldf("SELECT *
             FROM starter_df
             LEFT JOIN alcohol ON starter_df.State = alcohol.State AND starter_df.County = alcohol.County")
#Remove extra State, County columns
df <- subset(df, select=-c(4,5))

##2. merge heart to df
#Trim white space
heart$County <- trimws(heart$County)
#SQL join
df <- sqldf("SELECT *
             FROM df
             LEFT JOIN heart ON df.State = heart.State AND df.County = heart.County")
#remove extra State, County columns
df <- subset(df, select=-c(8,9))

##3. merge education to df
#Trim white space
education$County <- trimws(education$County)
#SQL join
df <- sqldf("SELECT *
             FROM df
             LEFT JOIN education ON df.State = education.State AND df.County = education.County")
#remove extra State, County columns
df <- subset(df, select=-c(10,11))

##4. merge eqi to df
#Trim white space
eqi$County <- trimws(eqi$County)
#SQL join
df <- sqldf("SELECT *
             FROM df
             LEFT JOIN eqi ON df.State = eqi.State AND df.County = eqi.County")
#remove extra State, County columns
df <- subset(df, select=-c(14,15))

##5. merge food to df
#Trim white space
food$County <- trimws(food$County)
#SQL join
df <- sqldf("SELECT *
             FROM df
             LEFT JOIN food ON df.State = food.State AND df.County = food.County")
#remove extra State, County columns
df <- subset(df, select=-c(15,16))

##6. merge sun to df
#Trim white space
sun$County <- trimws(sun$County)
#SQL join
df <- sqldf("SELECT *
             FROM df
             LEFT JOIN sun ON df.State = sun.State AND df.County = sun.County")
#remove extra State, County columns
df <- subset(df, select=-c(16,17))

##7. merge unemp to df
#Trim white space
unemp$County <- trimws(unemp$County)
#SQL join
df <- sqldf("SELECT *
             FROM df
             LEFT JOIN unemp ON df.State = unemp.State AND df.County = unemp.County")
#remove extra State, County columns
df <- subset(df, select=-c(17,18))



##8. merge wealth to df
#Trim white space
wealth$County <- trimws(wealth$County)
#SQL join
df <- sqldf("SELECT *
             FROM df
             LEFT JOIN wealth ON df.State = wealth.State AND df.County = wealth.County")
#remove extra State, County columns
df <- subset(df, select=-c(19,20))

##9. merge pop to df
#Trim white space
pop$County <- trimws(pop$County)
#SQL join
df <- sqldf("SELECT *
             FROM df
             LEFT JOIN pop ON df.State = pop.State AND df.County = pop.County")
#remove extra State, County columns
df <- subset(df, select=-c(21,22))

##verify variables and dimensions
head(df) 

dim(df) #3154 x 25

# ==========================================
# FIX MISSING INCOME VARIABLE
# ==========================================

# Check Income in wealth
print(names(wealth))

# Ensure Income is numeric
wealth$Income <- readr::parse_number(as.character(wealth$Income))

# Clean county names
wealth$County <- trimws(wealth$County)

# Check for duplicate county keys
# Check duplicate counties
duplicate_counties <- wealth %>%
  dplyr::count(State, County) %>%
  dplyr::filter(n > 1)

print(duplicate_counties, n = Inf)

# Display the actual duplicate records
wealth %>%
  dplyr::semi_join(
    duplicate_counties,
    by = c("State", "County")
  ) %>%
  dplyr::arrange(State, County) %>%
  print(n = Inf)

# Add Income to merged dataframe
df <- df %>%
  dplyr::select(-dplyr::any_of("Income")) %>%
  dplyr::left_join(
    wealth %>%
      dplyr::select(State, County, Income),
    by = c("State", "County")
  )

# Verify
print("Income" %in% names(df))

print(class(df$Income))

# ==========================================
# END INCOME FIX
# ==========================================

## --- EDA --- ##

#baseline EDA: glimpse() and summary()
glimpse(df)

#summary(df)


#drop NAs from consideration
df <- drop_na(df)
#dim(df) #3154 - 3004 = 150 dropped observations


#Select numeric variables

# Drop rows with missing values
df <- tidyr::drop_na(df)

# Select numeric variables, including health_score
df_num <- df %>%
  dplyr::select(where(is.numeric)) %>%
  as.data.frame()

# Check that health_score is present and numeric
print(names(df_num))

print(class(df_num$health_score))

# Correlation matrix
plot_corr_matrix(df_num, 0.2)

# Compute proportion of missing data per variable
Missing_Data <- sapply(
  df_num,
  function(x) sum(is.na(x)) / nrow(df_num)
)

# Compute correlation with health_score
predictors <- setdiff(names(df_num), "health_score")

HealthScore_Corr <- sapply(
  predictors,
  function(x) cor(
    df_num$health_score,
    df_num[[x]],
    use = "complete.obs"
  )
)

# Combine results
MDHSC <- data.frame(
  Variable = predictors,
  Missing_Data = Missing_Data[predictors],
  HealthScore_Corr = HealthScore_Corr
)

MDHSC %>%
  kbl(caption = "Proportion of Missing Data vs. Correlation with Health Score") %>%
  kable_minimal() %>%
  kable_styling(latex_options = "hold_position")

#Utilize Boruta for feature ranking and selection
library(Boruta)

# Perform Boruta search
boruta_output <- Boruta(health_score ~ ., data=na.omit(df_num), doTrace=0, maxRuns = 1000)

#Get significant variables including tentatives
boruta_signif <- getSelectedAttributes(boruta_output, withTentative = TRUE)
#print(boruta_signif)

# Plot variable importance
plot(boruta_output, cex.axis=.7, las=2, xlab="", main="Variable Importance")

#Feature selection (address weak target correlation and multicollinearity)
#select_df <- df_num[, c(1,4,6,8,11:14,16,18,23)]
#head(select_df) #verify

# Select variables for regression modeling
selected_vars <- c(
  "health_score",
  "Mortality",
  "College",
  "FoodInsecurity",
  "LTHighSchool",
  "Sun",
  "EQI",
  "Income",
  "PopChg",
  "Bng",
  "UnemploymentChg"
)

# Check whether any variables are missing
setdiff(selected_vars, names(df_num))

# Select variables
select_df <- df_num %>%
  dplyr::select(dplyr::all_of(selected_vars))

# Verify
dim(select_df)

names(select_df)

#Histograms for all variables
select_df %>%
  gather() %>%
  ggplot(aes(value)) +
  facet_wrap(~ key, scales = "free", ncol=4) +
  geom_histogram(bins=20,color="darkblue", fill="lightblue")

#Ensure reproducibility
set.seed(123)

#Train-test split data
dt = sort(sample(nrow(select_df), nrow(select_df)*.75))
train <- select_df[dt,]
test <- select_df[-dt,] 
#dim(train) #2253 x 11
#dim(test) # 751 x 11

###LINEAR REGRESSION (baseline model): prior to outlier handling, normalization, feature engineering

model_1 <- lm(health_score ~., data = train)
summary(model_1) #R^2 = 0.6603, vars = 10

#Imputation / removal (address NAs): 0 missing values

#Normalization: 0 impact on R^2
norm_minmax <- function(x){(x- min(x, na.rm = TRUE)) /(max(x, na.rm = TRUE)-min(x, na.rm = TRUE))}
norm_df <- as.data.frame(lapply(train, norm_minmax))
#head(norm_df) #verify

#Outlier handling: reduced R^2
#cooksD <- cooks.distance(model_1)
#influential <- as.numeric(names(cooksD)[(cooksD > (10 * mean(cooksD, na.rm = TRUE)))])
#train[influential,] #verify outliers - 20 rows
#out_df <- train[-influential,]

#Feature engineering

#combine datasets so we don't have to make features twice
train$dataset <- 'train'
test$dataset <- 'test'
final_df <- rbind(train, test)

#Creating new features (started with 1st, 3rd quartile values then adjusted)
##strong EQI and lots of sun
final_df$env_and_sun <- as.factor(ifelse(final_df$EQI > 0.90017 & final_df$Sun > 16350, 1, 0))
##high proportion college graduates and high income (3rd quartile)
final_df$college_and_income <- as.factor(ifelse(final_df$College > 20.0 & final_df$Income > 65000, 1, 0))
##low binge drinking and high income (older, more responsible?)
final_df$lowDrink_and_HiIncome <- as.factor(ifelse(final_df$Bng < 32.5 & final_df$Income > 61629, 1, 0))
##high binge drinking and high income (expendable income, socially active)
final_df$HiDrink_and_HiIncome <- as.factor(ifelse(final_df$Bng > 35.0 & final_df$Income > 61629, 1, 0))
##high food insecurity and high mortality and less than HS
final_df$LowFood_HiDeath_LowEd <- as.factor(ifelse(final_df$FoodInsecurity > 15.7 & final_df$Mortality > 500.00 & final_df$LTHighSchool > 20.0, 1, 0))
##low income and less than HS and growth in unemployment
final_df$LoIncome_LowEd_HiUnemp <- as.factor(ifelse(final_df$Income < 30000 & final_df$UnemploymentChg < -1.7 & final_df$LTHighSchool > 16.8, 1, 0))
##big PopChg and low unemployment
final_df$HiPop_LoUnemp <- as.factor(ifelse(final_df$PopChg > 1595.5 & final_df$UnemploymentChg > -1, 1, 0))

#Transformed (model 2): after handling outliers, normalization, feature engineering
train2 <- final_df %>% filter(dataset == 'train') %>% dplyr::select(-dataset)
test2 <- final_df %>% filter(dataset == 'test') %>% dplyr::select(-dataset)
dat <- final_df %>% dplyr::select(-dataset)
model_2 <- lm(health_score ~., data = train2)
#summary(model_2) #R^2 = 0.6735, vars = 17

###LINEAR REGRESSION (transformed model): features created, AIC optimized
#stepAIC(model_2)

model_3 <- lm(formula = health_score ~ Mortality + College + FoodInsecurity + 
                LTHighSchool + Sun + EQI + Income + PopChg + Bng + UnemploymentChg + 
                env_and_sun + lowDrink_and_HiIncome + HiDrink_and_HiIncome + 
                LowFood_HiDeath_LowEd + LoIncome_LowEd_HiUnemp + HiPop_LoUnemp, 
              data = train2)

#summary(model_3) #R^2 = 0.6734, vars = 16

# Predict and evaluate raw_lm model on training data
predictions = predict(model_3, newdata = train2)
eval_metrics(model_3, train2, predictions, target = 'health_score') #0.671, 0.0764

# Predict and evaluate raw_lm model on testing data
predictions = predict(model_3, newdata = test2)
eval_metrics(model_3, test2, predictions, target = 'health_score') #0.671, 0.0815

## --- Model Building --- ##

###LINEAR REGRESSION (raw data)

#Train-test split data
dt2 = sort(sample(nrow(df_num), nrow(df_num)*.75))
train_raw <- df_num[dt2,]
test_raw <- df_num[-dt2,] 
#dim(train_raw) #2253 x 23
#dim(test_raw) # 751 x 23

#Train raw_lm model:
raw_lm <- lm(health_score ~., data = train_raw)
#summary(raw_lm) #R^2 = 0.6808, vars = 22

# Predict and evaluate raw_lm model on training data
predictions = predict(raw_lm, newdata = train_raw)
eval_metrics(raw_lm, train_raw, predictions, target = 'health_score') #0.6678, 0.0759

# Predict and evaluate raw_lm model on testing data
predictions = predict(raw_lm, newdata = test_raw)
eval_metrics(raw_lm, test_raw, predictions, target = 'health_score') #0.6678, 0.0827

###LINEAR REGRESSION (raw data, AIC optimized)
#stepAIC(raw_lm)

#Train aic_raw_lm model:
aic_raw_lm <- lm(formula = health_score ~ Hvy + HvyPctChg + Bng + Mortality + 
                   HighSchool + SomeCollege + College + EQI + FoodInsecurity + 
                   Sun + Unemployment + UnemploymentChg + Poverty + Population + 
                   Births + Deaths + NetMig, data = train_raw)

#summary(aic_raw_lm) #R^2 = 0.6802, vars = 17

# Predict and evaluate aic_raw_lm model on training data
predictions = predict(aic_raw_lm, newdata = train_raw)
eval_metrics(aic_raw_lm, train_raw, predictions, target = 'health_score') #0.6671, 0.0761

# Predict and evaluate aic_raw_lm model on testing data
predictions = predict(aic_raw_lm, newdata = test_raw)
eval_metrics(aic_raw_lm, test_raw, predictions, target = 'health_score') #0.6671, 0.0828

###RIDGE REGRESSION (raw data)

#Specify column names
cr_raw = c('Hvy', 'HvyPctChg', 'BngPctChg', 'MortalityChg', 'HighSchool', 'SomeCollege', 'College', 'Unemployment', 'Poverty', 'Population', 'Births', 'Deaths', 'NetMig', 'Mortality', 'FoodInsecurity', 'LTHighSchool', 'Sun', 'EQI', 'Income', 'PopChg', 'Bng', 'UnemploymentChg', 'health_score')

#Generate dummy variables from data (if applicable)
dummies <- dummyVars(health_score ~ ., data = df_num[,cr_raw])
train_dummies = predict(dummies, newdata = train_raw[,cr_raw]) #2253 x 22
test_dummies = predict(dummies, newdata = test_raw[,cr_raw]) #751 x 22
print(dim(train_dummies)); print(dim(test_dummies))

#Create numeric model matrices
x = as.matrix(train_dummies)
y_train = train_raw$health_score
x_test = as.matrix(test_dummies)
y_test = test_raw$health_score

lambdas <- 10^seq(2, -3, by = -.1) #specify lambda sequence

#Train model
raw_ridge_reg = glmnet(x, y_train, nlambda = 25, alpha = 0, family = 'gaussian', lambda = lambdas)
summary(raw_ridge_reg)

#Compute optimal lambda
raw_cv_ridge <- cv.glmnet(x, y_train, alpha = 0, lambda = lambdas)
ol <- raw_cv_ridge$lambda.min
ol #0.003162278

#Compute R^2 from true and predicted values
eval_results <- function(true, predicted, df) {
  SSE <- sum((predicted - true)^2)
  SST <- sum((true - mean(true))^2)
  R_square <- round((1 - SSE / SST),4)
  RMSE = round(sqrt(SSE/nrow(df)),4)
  
  # Model performance metrics
  data.frame(RMSE = RMSE, Rsquare = R_square)
}

#Predict and evaluate raw_ridge_reg model on training data
predictions_train <- predict(raw_ridge_reg, s = ol, newx = x)
eval_results(y_train, predictions_train, train_raw) #RMSE: 0.0761, RSquare: 0.6734

#Predict and evaluate raw_ridge_reg model on testing data
predictions_test <- predict(raw_ridge_reg, s = ol, newx = x_test)
eval_results(y_test, predictions_test, test_raw) #RMSE: 0.0789, RSquare: 0.659

#Ridge regression lambda plot
plot(raw_ridge_reg)

###RIDGE REGRESSION (transformed data)

#Specify column names
cols_reg = c('Mortality', 'College', 'FoodInsecurity', 'LTHighSchool', 'Sun', 'EQI', 'Income', 'PopChg','Bng','UnemploymentChg','env_and_sun','lowDrink_and_HiIncome','HiDrink_and_HiIncome','LowFood_HiDeath_LowEd','LoIncome_LowEd_HiUnemp','HiPop_LoUnemp','health_score')

#Generate dummy variables from data (if applicable)
dummies <- dummyVars(health_score ~ ., data = dat[,cols_reg])
train_dummies2 = predict(dummies, newdata = train2[,cols_reg]) #2253 x 22
test_dummies2 = predict(dummies, newdata = test2[,cols_reg]) #751 x 22
print(dim(train_dummies2)); print(dim(test_dummies2))

#Create numeric model matrices
x2 = as.matrix(train_dummies2)
y_train2 = train2$health_score
x_test2 = as.matrix(test_dummies2)
y_test2 = test2$health_score

#Train model
ridge_reg = glmnet(x2, y_train2, nlambda = 25, alpha = 0, family = 'gaussian', lambda = lambdas)
summary(ridge_reg)

#Compute optimal lambda
cv_ridge2 <- cv.glmnet(x2, y_train2, alpha = 0, lambda = lambdas)
ol2 <- cv_ridge2$lambda.min
ol2 #0.003162278

#Predict and evaluate ridge_reg model on training data
predictions_train2 <- predict(ridge_reg, s = ol2, newx = x2)
eval_results(y_train2, predictions_train2, train2) #RMSE: 0.0764, RSquare: 0.6732

#Predict and evaluate ridge_reg model on training data
predictions_test2 <- predict(ridge_reg, s = ol2, newx = x_test2)
eval_results(y_test2, predictions_test2, test2) #RMSE: 0.0814, RSquare: 0.632

###LASSO REGRESSION (raw data)

# Setting alpha = 1 implements lasso regression
lasso_reg <- cv.glmnet(x, y_train, alpha = 1, lambda = lambdas, standardize = TRUE, nfolds = 5)

#Compute optimal lambda
lambda_best <- lasso_reg$lambda.min 
lambda_best #0.001

raw_lasso <- glmnet(x, y_train, alpha = 1, lambda = lambda_best, standardize = TRUE)

predictions_train <- predict(raw_lasso, s = lambda_best, newx = x)
eval_results(y_train, predictions_train, train_raw) #RMSE: 0.0767, Rsquare: 0.6681

predictions_test <- predict(raw_lasso, s = lambda_best, newx = x_test)
eval_results(y_test, predictions_test, test_raw) #RMSE: 0.0783, Rsquare: 0.6642

#Lasso regression lambda plots
op <- par(mfrow=c(1, 2))

plot(lasso_reg$glmnet.fit, "norm",   label=TRUE)
plot(lasso_reg$glmnet.fit, "lambda", label=TRUE)

par(op)

###LASSO REGRESSION (transformed data)

# Setting alpha = 1 implements lasso regression
lasso_reg2 <- cv.glmnet(x2, y_train2, alpha = 1, lambda = lambdas, standardize = TRUE, nfolds = 5)

#Compute optimal lambda
lambda_best2 <- lasso_reg2$lambda.min 
lambda_best2 #0.001

lasso_model <- glmnet(x2, y_train2, alpha = 1, lambda = lambda_best2, standardize = TRUE)

predictions_train2 <- predict(lasso_model, s = lambda_best2, newx = x2)
eval_results(y_train2, predictions_train2, train2) #RMSE: 0.0765, Rsquare: 0.6724

predictions_test2 <- predict(lasso_model, s = lambda_best2, newx = x_test2)
eval_results(y_test2, predictions_test2, test2) #RMSE: 0.0811, Rsquare: 0.6354

## --- Model Selection --- ##

#Create Kable table to succinctly summarize model optimization results
Model <- c('1', '2', '3', '4', '5', '6', '7')
Method <- c('Linear', 'Linear', 'Linear', 'Ridge', 'Ridge', 'Lasso', 'Lasso')
Data <- c('raw', 'raw(AIC)', 'transformed', 'raw', 'transformed', 'raw', 'transformed')
Var_Num <- c(22, 17, 16, 22, 16, 22, 16)
R2_train <- c(0.6678, 0.673, 0.671, 0.6734, 0.6732, 0.6681, 0.6724)
RMSE_train <- c(0.0759, 0.0759, 0.0764, 0.0761, 0.0764, 0.0767, 0.0765)
R2_test <- c(0.6678, 0.673, 0.71, 0.659, 0.632, 0.6642, 0.6354)
RMSE_test <- c(0.0827, 0.0842, 0.0815, 0.0789, 0.0814, 0.0783, 0.0811)

output <- cbind(Model, Method, Data, Var_Num, R2_train, RMSE_train, R2_test, RMSE_test)

output %>%
  kbl(caption = "Regression Model Comparison") %>%
  kable_minimal() %>%
  kable_styling(latex_options = "hold_position")