04_LC02_model_train_eval.Rmd

---
title: "LC02 chl-a model"
author: "Sam Sillen"
date: "2023-12-08"
output: html_document
---

# Packages
```{r setup, include=FALSE}
library(tidyverse)
library(feather)
library(viridis)
library(sf)
library(maps)
library(magrittr)
library(purrr)
library(data.table)
library(ggthemes)
library(dplyr)
library(ggplot2)
library(mapview)
library(fs)
library(httr)
library(leaflet)
library(nhdplusTools)
library(foreign)
library(CAST)
library(caret)
library(sp)
library(xgboost)
library(Metrics)
library(parallelly)
library(doParallel)
library(ggpmisc)

knitr::opts_chunk$set(echo = TRUE)
```

# Holdout function for test / train 
```{r functions, echo=FALSE}

# This function randomly samples match-ups across different locations, times, and concentrations for splitting training/validation data

holdout <- function(x) {

  x <- x %>%
  group_by(long_group, time_group) %>% # Split up into spatial and temporal groups
  dplyr::mutate(mag = cut(value, quantile(
  x = value,
  c(0, 0.2, 0.4, 0.6, 0.8, 0.9, 1),
  include.lowest = T
  )),
  mag = factor(
  mag,
  labels = c( 0.2, 0.4, 0.6, 0.8, 0.9, 1)
  )) %>%
  ungroup()
  
  set.seed(22)
  
  train <- x %>%
  group_by(time_group, long_group, mag) %>%
  sample_frac(.8) %>% # 80% of data will be used for training
  ungroup() %>%
  dplyr::mutate(.partitions = 1)
  
  validate <- x %>%
   anti_join(train) %>%
   dplyr::mutate(.partitions = 2)

  out <- train %>%
  bind_rows(validate) 
    
  return(out)
}
```

# Filter matchup data to roi 
```{r chla}
# SF for HUC2 basins from USDA NHD WBD - includes Mid Atlantic (HUC2_02), Great Lakes (HUC2_04), Ohio (HUC2_05), Tennessee (HUC2_06) and Upper Mississippi (HUC2_07)
huc2_basins <- read_sf(dsn="C:/Users/samsi/OneDrive - University of Pittsburgh/nhd_data/HUC2_Merge", layer="HUC2_Merge") 
st_transform(huc2_basins, crs ="+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

huc2_basins <- huc2_basins %>% filter(name != 'Mid Atlantic Region')

# Load in LC02 matchup database
matchup <- read_feather("C:/Users/samsi/OneDrive - University of Pittsburgh/OhioRiver_SR/Matchups/LC02_matchup_algal_mask.feather")

matchup <- matchup %>% 
filter(harmonized_parameter == 'chl.a')  %>% 
mutate(uniqueID = row_number())

pnts <- matchup %>% # Create new df with matchup ID, lat and long columns
  select(uniqueID, long, lat) 

sp_pnts = st_as_sf(pnts, coords=c('long', 'lat')) 
st_crs(sp_pnts)=4326
sp_pnts = st_transform(sp_pnts, crs=st_crs(huc2_basins)) # Make points spatial and set coordinate system to same as the HUC polygon

# Perform intersection between HUC2 shapefile and matchup points
Int=st_intersects(sp_pnts, huc2_basins)

# Create a new column in the sp_pnts data frame and write intersection results to it
sp_pnts$Intersect=lengths(Int)>0 #Finds only countries where intersection is TRUE

# Combine intersection results with original matchup dataset
matchup <- cbind(matchup, sp_pnts$Intersect) 

# Rename newly added column
matchup <- matchup %>% rename(Intersect = "sp_pnts$Intersect")
```

# Applying filters to improve model performance 
```{r}
# Finish cleaning the data by filtering out matchups that are outside of HUC boundaries and other filters
matchup_filter <- matchup %>%
  filter(algal_mask == 1) %>% 
  filter(value < 200) %>% 
  filter(hillShadow == 1 | is.na(hillShadow)) %>%
  filter(value > 0.1) %>% # Minimum accepted chlorophyll value
  filter(Intersect == TRUE) %>% 
  filter(pCount_algal_mask > 30) %>% 
  filter(characteristicName != 'Chlorophyll a, corrected for pheophytin') %>%
  filter(Surface_temp_kelvin > 270) %>%
  select(date:gn.gn, lat, long, uniqueID, hue, dw, saturation, type, SiteID) # Select only the important vars, otherwise calling the holdout function will produce an error


test <- matchup_filter %>%
  filter(Green <0.045, value > 30) 


test1 <- matchup_filter %>% 
  filter(Green > 0.035 & value < 5)

sr_tests <- rbind(test, test1)

matchup_filter_final <- matchup_filter %>% 
  filter(!uniqueID %in% test2$uniqueID)

```

# Setting up base model (no ffs, hypertuning)
```{r}
# Make splits for training/validation

df <- matchup_filter_final %>%
  mutate(lat_group = cut_number(lat, 2, right= F),
         long_group = cut_number(long, 2, right=F),
         date = lubridate::ymd(date),
         julian = as.numeric(julian.Date(date)),
         space_group = paste0(lat_group,long_group),
         time_group = cut_number(julian, 3, right=F),
        value = log(value)) %>%
         holdout() %>% 
         ungroup() %>%
         as.data.frame() %>%
  filter_all(all_vars(!is.infinite(.))) %>%
  filter_all(all_vars(!is.nan(.)))

df <- df %>% distinct(date, Blue, Red, Green, value, .keep_all = TRUE) # double check for duplicates


# Training data #80% of data to training 
train <- df %>%
  filter(.partitions ==1) %>% 
  ungroup() %>%
  as.data.frame() %>%
  filter_all(all_vars(!is.infinite(.))) %>%
  filter_all(all_vars(!is.nan(.))) 


# Validation data 10% of data to validation
validate <- df %>%
  filter(.partitions ==2) %>%
  ungroup() %>%
  as.data.frame()%>%
  filter_all(all_vars(!is.infinite(.))) %>%
  filter_all(all_vars(!is.nan(.))) 


# Make validation row index so you can rejoin later
val.cols <- df %>% 
  filter(.partitions ==2) %>%
  ungroup() %>%
  filter_all(all_vars(!is.infinite(.))) %>%
  filter_all(all_vars(!is.nan(.))) 

# Select spectral indices/bands to use as predictors in model.
features_1 <- df %>%
  select(NR:gn.gn, GCI, IRG, SABI, KIVU, GB, GNDVI, EVI, KAB, KRL, Surface_temp_kelvin, dw, hue, saturation) %>%
  names(.)

# Create cross validation folds for spatial-temporal cross validation
folds <- CreateSpacetimeFolds(train, spacevar = "long_group", timevar = "time_group", k = 2)

# Set training parameters
train_control_final <- caret::trainControl(
  method = "cv",
  savePredictions = T,
  returnResamp = 'final',
  index = folds$index,
  indexOut = folds$indexOut,
  verboseIter = T,
  allowParallel = TRUE,
  p = 0.8 #80% of the data is used to predict the other 10%
  )
  
grid_final <-expand.grid(
  nrounds = 50,
  alpha = 0,
  lambda =1,
  eta = 0.3
)

# Make a model
model <- caret::train(
  x = train[,features_1],
  y = train$value,
  trControl = train_control_final,
  tuneLength = 1, 
  method = "xgbLinear",
   importance = T,
  verbose = TRUE
)

# Use model to make predictions over validation data
pred<- predict(model, validate[,features_1])
actual <- (validate$value)
uniqueID <- val.cols$uniqueID

output <- tibble(Predicted = pred, Actual = actual, uniqueID = uniqueID) %>%
  mutate(Actual = exp(Actual), Predicted =  exp(Predicted)) %>%
  left_join(df, by="uniqueID") %>%
  mutate(residual = Actual - Predicted,
         year = year(date),
         month = month(date),
         obs = ifelse(abs(residual) > 15, "bad", "good"))

# Calcualate error metrics 
evals <- output %>%
  mutate(Actual = (Actual), 
         Predicted = (Predicted)) %>%
  summarise(rmse = rmse(Actual, Predicted),
            mae = mae(Actual, Predicted),
            mape = mape(Actual, Predicted),
            bias = bias(Actual, Predicted),
            p.bias = percent_bias(Actual, Predicted),
            smape = smape(Actual, Predicted)) 

print(evals)
```

# Validation plot
```{r}

# Plot in log space

ggplot(output , aes(x = log(Actual), y = log(Predicted))) + 
  geom_point(colour = 'black', alpha = 0.7) +
  geom_abline(slope=1, intercept = 0, color = 'black')+
  xlab("Measured  (ug/L)") +
  ylab("Predicted (ug/L)")+
  labs(title = 'Base Model')+
  stat_poly_eq(colour = 'black') +
  theme_few() +
    theme(axis.text = element_text(colour = 'black', size = 20), axis.title = element_text(size = 20, colour = 'black'),plot.title = element_text(size = 30, colour = 'black'),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), 
            panel.border = element_blank(),
            legend.background = element_rect(fill = '#292b30'), 
            legend.title = element_text(colour = 'white'))

#Light theme ; non-log space

ggplot(output , aes(x = (Actual), y = (Predicted))) + 
  geom_point(colour = 'black', alpha = 0.7) +
  xlim(0, 200) +
  ylim(0, 200) +
  geom_abline(slope=1, intercept = 0, color = 'black')+
  xlab("Measured  (ug/L)") +
  ylab("Predicted (ug/L)")+
  #scale_x_continuous(trans = "log10") +
  #scale_y_continuous(trans = "log10") +
  #labs(title = 'Base Model')+
  stat_poly_eq(colour = 'black') +
  theme_few() +
    theme(axis.text = element_text(colour = 'black', size = 20), axis.title = element_text(size = 20, colour = 'black'),plot.title = element_text(size = 30, colour = 'black'),
           # panel.grid.major = element_blank(),
          #  panel.grid.minor = element_blank(), 
            panel.border = element_blank(),
            #legend.background = element_rect(fill = '#ffe098ff'), 
            legend.title = element_text(colour = 'white'))

#ggsave("/Users/samsillen/OneDrive - University of Pittsburgh/OhioRiver_SR/Plots/ASLO_model_results.png", width=8, height =6, units="in", dpi=330)
```

# FFS for selecting best combination of predictor variables 
```{r}
set.seed(10)

folds <- CreateSpacetimeFolds(train,
  spacevar = "long_group",
  timevar = "time_group" )
  
control <- trainControl(
  method = "cv",
  savePredictions = 'none',
  returnResamp = 'final',
  index = folds$index,
  indexOut = folds$indexOut,
  p = 0.8)
  
# Do initial feature selection with conservative hyperparameters
tuneGrid1 <- expand.grid(
  nrounds = 300,
  eta = .1,
  lambda = 0,
  alpha = 0)

# Set it up to run in parallel. This can take 1-2 days.
cl <- makePSOCKcluster(availableCores() - 4)
registerDoParallel(cl)

ffs <- ffs(df[,features_1], df$value, method = 'xgbLinear', metric = 'RMSE', tuneGrid = tuneGrid1, Control = control, verbose = T)

on.exit(stopCluster(cl))

ffsResults <- ffs$perf_all

# Save the results
write_feather(ffsResults, "C:/Users/samsi/OneDrive - University of Pittsburgh/OhioRiver_SR/Models/final_ffs_results.feather")
#
ffsResults %>%
  group_by(nvar) %>%
  summarise(RMSE = median(RMSE),
            SE = median(SE)) %>%
  ggplot(.) + geom_line(aes(x = nvar, y = RMSE)) +
  geom_errorbar(aes(x = nvar, ymin = RMSE - SE, ymax = RMSE + SE), color = 'red')

#ggsave(paste0('figs/rfeRMSE_', iter, '.png'), device = 'png', width = 6, height = 4, units = 'in')
```

# Hyperparamter tuning 
```{r}
#Sometimes the following function may be necessary due to a lingering parallel operation

#unregister_dopar <- function() {
#    env <- foreach:::.foreachGlobals
#  rm(list=ls(name=env), pos=env)
#}

#unregister_dopar()

ffsResults <- read_feather("C:/Users/samsi/OneDrive - University of Pittsburgh/OhioRiver_SR/Models/final_ffs_results.feather")

ffs_features <- ffsResults[ffsResults$RMSE == min(ffsResults$RMSE),]  %>%
  dplyr::select(-c(nvar, RMSE, SE)) %>%
  slice_head(n = 1) %>%
  select(-var13) %>% 
  paste(.)

grid_base <- expand.grid(
  nrounds = seq(100,500,100),
  alpha = c(0.01, 0.1, 0.5, 1),
  lambda = c(0.01, 0.1, 0.5, 1),
  eta = c(0.05, 0.1, 0.3)
)

set.seed(10)

folds <- CreateSpacetimeFolds(train, spacevar = "long_group", timevar = "time_group" , k=5)

train_control <- caret::trainControl(
  method = "cv",
  savePredictions = T,
  returnResamp = 'final',
  index = folds$index,
  indexOut = folds$indexOut,
  allowParallel = TRUE,
  p = 0.9,
  )
  
base_model <- caret::train(
  x = train[,ffs_features],
  y = train$value,
  trControl = train_control,
  tuneGrid = grid_base,
  method = "xgbLinear",
  verbose = TRUE,
 # preProcess = c('center', 'scale'),
  importance = F
)

base_model$bestTune
 
train_control_final <- caret::trainControl(
  method = "cv",
  savePredictions = T,
  returnResamp = 'final',
  index = folds$index,
  indexOut = folds$indexOut,
  allowParallel = TRUE,
  p = 0.9,
  )
  
grid_final <- expand.grid(
  nrounds = base_model$bestTune$nrounds,
  alpha = base_model$bestTune$alpha,
  lambda = base_model$bestTune$lambda,
  eta = base_model$bestTune$eta
)

final_model <- caret::train(
  x = train[,ffs_features],
  y = train$value,
  trControl = train_control_final,
  tuneGrid = grid_final,
  method = "xgbLinear",
 # preProcess = c('center', 'scale'),
  importance = F
)

final_model$bestTune

saveRDS(final_model, file = "C:/Users/samsi/OneDrive - University of Pittsburgh/OhioRiver_SR/Models/CHOIR_final_model.rDS")

```

# Evaluate final model
```{r eval}
# load final model 

final_model <- readRDS("C:/Users/samsi/OneDrive - University of Pittsburgh/OhioRiver_SR/Models/CHOIR_final_model.rDS")

# evaulate the model
pred<- predict(final_model, validate[,ffs_features])
actual <- (validate$value)
uniqueID <- val.cols$uniqueID

output <- tibble(Predicted = pred, Actual = actual,uniqueID = uniqueID) %>%
  mutate(Actual = exp(Actual), Predicted = exp(Predicted)) %>%
  left_join(df, by="uniqueID") %>%
  mutate(residual = Actual - Predicted,
         year = year(date),
         month = month(date),
         obs = ifelse(abs(residual) > quantile(abs(residual), .7, na.rm=T), "bad", "good"))

# calculate error metrics 
evals <- output %>%
  mutate(Actual = (Actual), 
         Predicted = (Predicted)) %>%
  summarise(rmse = rmse(Actual, Predicted),
            mae = mae(Actual, Predicted),
            mape = mape(Actual, Predicted),
            bias = bias(Actual, Predicted),
            p.bias = percent_bias(Actual, Predicted),
            smape = smape(Actual, Predicted)) 

print(evals)
```

# Final model plots
```{r}

# Plot in log space

ggplot(output , aes(x = log(Actual), y = log(Predicted))) + 
  geom_point(colour = 'black', alpha = 0.7) +
  geom_abline(slope=1, intercept = 0, color = 'black')+
  ylim(-2, 6) + 
  xlim(-2, 6) +
  xlab("Measured  (ug/L)") +
  ylab("Predicted (ug/L)")+
 # labs(title = 'Base Model')+
  stat_poly_eq(colour = 'black') +
  theme_bw() +
    theme(axis.text = element_text(colour = 'black', size = 20), axis.title = element_text(size = 20, colour = 'black'),plot.title = element_text(size = 30, colour = 'black'),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), 
            legend.background = element_rect(fill = '#292b30'), 
            legend.title = element_text(colour = 'white'))

ggsave("C:/Users/samsi/OneDrive - University of Pittsburgh/OhioRiver_SR/Plots/final_model_plot.jpeg", width=8, height =6, units="in", dpi=330)

#Light theme ; non-log space

ggplot(output , aes(x = (Actual), y = (Predicted))) + 
  geom_point(colour = 'black', alpha = 0.7) +
  xlim(0, 200) +
  ylim(0, 200) +
  geom_abline(slope=1, intercept = 0, color = 'black')+
  xlab("Measured  (ug/L)") +
  ylab("Predicted (ug/L)")+
  #scale_x_continuous(trans = "log10") +
  #scale_y_continuous(trans = "log10") +
  #labs(title = 'Base Model')+
  stat_poly_eq(colour = 'black') +
  theme_few() +
    theme(axis.text = element_text(colour = 'black', size = 20), axis.title = element_text(size = 20, colour = 'black'),plot.title = element_text(size = 30, colour = 'black'),
           # panel.grid.major = element_blank(),
          #  panel.grid.minor = element_blank(), 
            panel.border = element_blank(),
            #legend.background = element_rect(fill = '#ffe098ff'), 
            legend.title = element_text(colour = 'white'))

```

# Explore variable importance
```{r}
library(vip)

imp <- varImp(final_model, scale = FALSE)

imp <- as.data.frame(imp$importance) %>% rownames_to_column() %>% rename(Feature = 'rowname')

#For cleaner predictor variables names in plot (note: this is outdated , need to update new feature names based on this model iteration)
#names_new <- tibble(Feature = c("R.BS", "fai", "GCI", "GB", "Surface_temp_kelvin", "IRG", "SN", 
#"BR_G", "N.RB", "RS", "RN"), namesNew = c("Red / (Blue + Swir1)" , "Floating Algal Index", "Green Chlorophyll Index", "Green / Blue", "Surface Temp Kelvin", "Red / Green", "Swir 1 / Nir", 
#"Blue - Red / Green", "Nir / (Red + Blue)", "Red / Swir1", "Red / Nir"))

#imp <- left_join(imp, names_new) %>% select(-Feature) %>% rename(Feature = "namesNew") %>% drop_na()

imp$Feature <- factor(imp$Feature, levels = imp$Feature)

ggplot(imp, aes(x = Overall)) + 
geom_bar( aes(y = Feature), stat = 'identity', fill = "#fab62f") + 
guides(fill="none")+ 
theme_bw() +
labs(title = "")+
    scale_y_discrete(limits=rev)+ 
    xlab("Gain") + 
ylab("Predictor") + 
  theme(axis.text = element_text(colour = 'white', size = 20), axis.title = element_text(size = 20, colour = 'white'),plot.title = element_text(size = 30, colour = 'white'),
            plot.background = element_rect(fill = '#292b30', color = 'NA'),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(), 
            panel.border = element_blank(),
            panel.background = element_rect(fill = '#292b30'),
            legend.background = element_rect(fill = '#292b30'), 
            legend.title = element_text(colour = 'white'), legend.text = element_text(colour = 'white'))

#ggsave("C:/Users/samsi/OneDrive - University of Pittsburgh/Ohio_River_SR/Models/final_model_Varimp.jpg", width=8, height = 6, units="in", dpi=330)
```

# Use final model to predict chlorophyll-a on RiverSR dataset 
```{r}
# read in riverSR and predict  

riverSR <- read_feather("C:/Users/samsi/OneDrive - University of Pittsburgh/OhioRiver_SR/Data/algal_mask/IL_algal_mask_pred.feather")

riverSR$pred <- predict(final_model, riverSR[,ffs_features])

# predicted value will be logged so you need to call exp 
riverSR <- riverSR %>% 
mutate(pred = exp(pred))

write_feather(riverSR, "C:/Users/samsi/OneDrive - University of Pittsburgh/OhioRiver_SR/Data/RiverSR_IL_pred.feather")
```