# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .r
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: R
#     language: R
#     name: ir
# ---

# # Determining Feature Importance With Boruta

# ## Importing Data
# Also removing some columns from the data that are unnecessary for this step 

library(here)
library(ggplot2)
library(caret)
library(Boruta)
data = read.csv(here("data","data.csv"))
data[,c("X","income","education","population","unemployment","popDensity","laborforce","poverty","diversity","x_sp","y_sp","zip_city","cncldist","borocode")] = NULL

# ## Boruta Plot Function
# This function will let us visualize the results of our Boruta models

boruta_plot = function(model, main_title="Boruta Feature Importance Plot",save_plot=F) {
  if(save_plot==T){png(here("plots", paste(main_title,".png", sep="")), width = 465, height = 225, units='mm', res = 300)}
  plot(model, xlab = "", xaxt = "n", main = main_title)
  
  lz<-lapply(1:ncol(model$ImpHistory),function(i)
    
    model$ImpHistory[is.finite(model$ImpHistory[,i]),i])
  
  names(lz) <- colnames(model$ImpHistory)
  
  Labels <- sort(sapply(lz,median))
  
  axis(side = 1,las=2,labels = names(Labels),
       
       at = 1:ncol(model$ImpHistory), cex.axis = 0.7)
  if(save_plot==T){dev.off()}
  return()
}

# ## Determining Initial Correlation

# +
library(reshape2)
data_num = data[,unlist(lapply(data, is.numeric))]
data_cor = cor(data_num, use="complete.obs")

get_upper_tri = function(df){
  df[lower.tri(df)]<- NA
  return(df)
}

upper = get_upper_tri(data_cor)

melted = melt(upper, na.rm = F)
ggplot(data = melted, aes(Var2, Var1, fill = value))+
  geom_tile(color = "white")+
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-1,1), space = "Lab", 
                       name="Pearson\nCorrelation") +
  theme_minimal()+ 
  theme(axis.text.x = element_text(angle = 90))+
  coord_fixed()
# -

# ### Preparing Data for Boruta

# +
testing = data[complete.cases(data),]
testing = testing[testing$year == 2015,]
testing[,c("ntaname_full","nta_name","nta","genus","species")] = NULL
testing$healthstatus[testing$healthstatus == 'Alive'] = 'Good'
testing$healthstatus = droplevels(testing$healthstatus)

partition <- createDataPartition(testing$healthstatus, p = .5, list = FALSE, times = 1)
train = testing[partition,]
test = testing[-partition,]
# -

# ## Primary Data Section
# ### Performing Boruta on Healthstatus in Primary Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/NYCTrees%20Feature%20Importance.png)

bor = Boruta(healthstatus ~., train, doTrace = 2)

saveRDS(bor,here("models","boruta_treehealth_variable_importance.rds"))
boruta_plot(bor, main_title = "Primary Data Healthstatus Feature Importance", save_plot = T)

# ### Performing Boruta on tree_dbh in Primary Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Boruta%20tree_dbh%20Feature%20Importance.png)

bor = Boruta(tree_dbh ~., train, doTrace = 2)
fixed_bor = TentativeRoughFix(bor)
boruta_plot(bor, main_title = "Boruta tree_dbh Feature Importance", save_plot = T)
boruta_plot(fixed_bor, main_title = "Boruta tree_dbh Feature Importance", save_plot = T)

# ### Performing Boruta on healthstatus in Primary Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Primary%20Data%20Healthstatus%20Feature%20Importance.png)

health_bor = Boruta(healthstatus ~., train, doTrace = 2)

# ## Aggregate Data Section
# ### Performing Boruta on Num_Trees in Aggregate Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Aggregated%20Data%20Num_Trees%20Feature%20Importance.png)

# +
agg_data = read.csv(here("data","Joined_aggregate_data.csv"))
agg_data[c("Borough","X.1","X.x","X.y","X")] = NULL
agg_train = agg_data[complete.cases(agg_data),]

agg_bor = Boruta(Num_Trees ~., agg_train, doTrace = 2)
agg_bor = TentativeRoughFix(agg_bor)

boruta_plot(agg_bor, main_title = "Aggregated Data Num_Trees Feature Importance", save_plot = T)
# -

# ### Performing Boruta on Income in Aggregate Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Aggregated%20Data%20Income%20Feature%20Importance.png)

agg_income = Boruta(Income ~., agg_train, doTrace = 2)
agg_income = TentativeRoughFix(agg_income)
boruta_plot(agg_income, main_title = "Aggregated Data Income Feature Importance", save_plot = T)

# ### Performing Boruta on Professional in Aggregate Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Aggregated%20Data%20Professional%20Feature%20Importance.png)

agg_prof = Boruta(Professional ~., agg_train, doTrace=2)
agg_prof = TentativeRoughFix(agg_prof)
boruta_plot(agg_prof, main_title = "Aggregated Data Professional Feature Importance", save_plot = T)

# ## Water Aggregate Data Section
# ### Performing Boruta on Trees_x_sq_mi in Water Aggregate Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Aggregated%20Water%20Data%20Trees_x_sq_mi%20Feature%20Importance.png)

# +
wtr_data = read.csv(here("data","Joined_aggregate_water_data.csv"))
wtr_data[,c("ID","X...Coordinate.y","Y...Coordinate.y","geometry","Borough","X.1","X","intersection")] = NULL
wtr_data = wtr_data[complete.cases(wtr_data),]

trees_x_bor = Boruta(Trees_x_sq_mi ~., wtr_data, doTrace = 2)
trees_x_bor = TentativeRoughFix(trees_x_bor)
boruta_plot(trees_x_bor, main_title = "Water Data Trees_x_sq_mi Feature Importance", save_plot = T)
# -

# ### Performing Boruta on Fluoride in Water Aggregate Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Aggregated%20Water%20Data%20Fluoride%20Feature%20Importance.png)

fluor_bor = Boruta(Fluoride..mg.L. ~., wtr_data, doTrace = 2)
fluor_bor = TentativeRoughFix(fluor_bor)
boruta_plot(fluor_bor, main_title = "Aggregated Water Data Fluoride Feature Importance", save_plot = T)

# ### Performing Boruta on Coliform in Water Aggregate Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Aggregated Water Data Coliform Feature Importance.png)

coli_bor = Boruta(Coliform..Quanti.Tray...MPN..100mL. ~., wtr_data, doTrace=2)
coli_bor = TentativeRoughFix(coli_bor)
boruta_plot(coli_bor, main_title = "Aggregated Water Data Coliform Feature Importance", save_plot = T)

# ### Performing Boruta on E. Coli in Water Aggregate Data
# [Output found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Aggregated%20Water%20Data%20E.Coli%20Feature%20Importance.png)

Ecoli_bor = Boruta(E.coli.Quanti.Tray...MPN.100mL. ~., wtr_data, doTrace=2)
Ecoli_bor = TentativeRoughFix(Ecoli_bor)
boruta_plot(Ecoli_bor, main_title = "Aggregated Water Data E.Coli Feature Importance", save_plot = T)

# ## PCA Section
# ### PCA on Primary Data
# First, several columns need to be translated to numeric. Only numeric columns are kept for PCA.
# Afterwards, [output can be found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Raw%20Data%20PCA.png)

# +
pca_data = data
pca_data$inf_guard = ifelse(pca_data$inf_guard == "Yes", 1, 0)
pca_data$inf_shoes = ifelse(pca_data$inf_shoes == "Yes", 1, 0)
pca_data$sidw_crack = ifelse(pca_data$sidw_crack == "Yes", 1, 0)
pca_data$healthstatus = ifelse(pca_data$healthstatus == "Good" | pca_data$healthstatus == "Alive", 1, 0)
pca_data = pca_data[,unlist(lapply(pca_data, is.numeric))]
pca_data = pca_data[complete.cases(pca_data),]

pca_data_mod = prcomp(pca_data, scale = TRUE)

eig = (pca_data_mod$sdev)^2
variance = eig*100/sum(eig)
cumvar = cumsum(variance)

eig.dat = data.frame(eig = eig, variance = variance, cumvariance = cumvar)
eig.dat

png(here("plots", "Raw Data PCA.png"), width = 465, height = 225, units='mm', res = 300)
fviz_pca_var(pca_data_mod, col.var="contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE,      # Avoid text overlapping
             title = "Raw Data PCA"
)
dev.off()
# -

# ### PCA on Aggregate Data
# Several columns were dropped and others were converted from factors to numeric for PCA.
# [Output can be found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Aggregated%20Data%20PCA.png)

# +
pca_agg = agg_train
pca_agg[,c("NTA_large","NTA_small")] = NULL
pca_agg$Hosp_5.17_Avg_yr_Num=as.numeric(pca_agg$Hosp_5.17_Avg_yr_Num)
pca_agg$Hosp_5.17_Avg_Yr_Rate_per_10K = as.numeric(pca_agg$Hosp_5.17_Avg_Yr_Rate_per_10K)

pca_agg_mod = prcomp(pca_agg, scale = TRUE)

eig = (pca_data_mod$sdev)^2
variance = eig*100/sum(eig)
cumvar = cumsum(variance)

agg_eig.dat = data.frame(eig = eig, variance = variance, cumvariance = cumvar)
agg_eig.dat

png(here("plots", "Aggregated Data PCA.png"), width = 465, height = 225, units='mm', res = 300)
fviz_pca_var(pca_agg_mod, col.var="contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE,      # Avoid text overlapping,
             title = "Aggregated Data PCA"
)
dev.off()
# -

# ### PCA on Water Aggregate Data
# Several columns were dropped and others were converted from factors to numeric for PCA. 
# Afterwards, [output can be found here](https://github.com/kbfoerster/nyctrees/blob/master/plots/Water%20Data%20PCA.png)

# +
pca_wtr_data = wtr_data
pca_wtr_data[,c("Location","Sample.Date","Sample.Site","area","NTA_large","Sample.Time")] = NULL
pca_wtr_data$Residual.Free.Chlorine..mg.L. = as.numeric(pca_wtr_data$Residual.Free.Chlorine..mg.L.)
pca_wtr_data$Turbidity..NTU. = as.numeric(pca_wtr_data$Turbidity..NTU.)
pca_wtr_data$Fluoride..mg.L. = as.numeric(pca_wtr_data$Fluoride..mg.L.)
pca_wtr_data$Coliform..Quanti.Tray...MPN..100mL. = as.numeric(pca_wtr_data$Coliform..Quanti.Tray...MPN..100mL.)
pca_wtr_data$E.coli.Quanti.Tray...MPN.100mL. = as.numeric(pca_wtr_data$E.coli.Quanti.Tray...MPN.100mL.)
pca_wtr_data$Hosp_5.17_Avg_yr_Num = as.numeric(pca_wtr_data$Hosp_5.17_Avg_yr_Num)
pca_wtr_data$Hosp_5.17_Avg_Yr_Rate_per_10K = as.numeric(pca_wtr_data$Hops_Adults_Avg_Yr_Rate_per_10K_adj)

pca_wtr_data_mod = prcomp(pca_wtr_data, scale = TRUE)

eig = (pca_data_mod$sdev)^2
variance = eig*100/sum(eig)
cumvar = cumsum(variance)

wtr_eig.dat = data.frame(eig = eig, variance = variance, cumvariance = cumvar)
wtr_eig.dat

png(here("plots", "Water Data PCA.png"), width = 465, height = 225, units='mm', res = 300)
fviz_pca_var(pca_wtr_data_mod, col.var="contrib",
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), 
             repel = TRUE, # Avoid text overlapping
             title = "Water Data PCA"
)
dev.off()

StarcoderdataJupyterScriptsDedupFiltered