R_simple_ml_example.Rmd

---
title: "Untitled"
output: html_document
date: '2022-07-16'
---

```{r setup, include=FALSE}
library(readxl)
library(dplyr)
library(tidyverse)
library(tidymodels)      # for the recipes package, along with the rest of tidymodels
library(vip)
library(skimr)           # for variable summaries
library(rpart.plot)
```

```{r}

# file location to read
file <- "C:\\Users\\kelvi\\Desktop\\HMRC_Senior_Officials_Travel_July_to_September_2022.csv" 
# load file
file <- read_csv(file)

# drop not needed columns
file <- dplyr::select(file,-`Purpose of trip`,-`Total Cost of Use of Official Secure Car`,-`Other (including hospitality given)`)

# remove any "Nil Return" from dates
file <- file %>%
  filter(`Start date of trip` !="Nil return")

# for Class of travel replace Nil return with Standard
file <- file %>%
  mutate(`Class of travel`=ifelse(`Class of travel`=="Nil return","Standard",`Class of travel`))

# for accommodation/meals() replace Nil return with 0 and make numeric
file <- file %>%
  mutate(`Accommodation/Meals ()`=ifelse(`Accommodation/Meals ()`=="Nil return","0",`Accommodation/Meals ()`))

file$`Accommodation/Meals ()` <- as.numeric(file$`Accommodation/Meals ()`)

file$`Duration of Visit (Days)` <- as.numeric(file$`Duration of Visit (Days)`)

# make total cost numeric and rename column
file <- rename(file,c('total cost'='Total cost, including all visas, accommodation, travel, meals etc. (£)'))
file$`total cost` <- as.numeric(file$`total cost`)

# make new column that has 'transport cost'
file$`transport cost` <- file$`total cost` - file$`Accommodation/Meals ()`

# change any NAs to 0
file$`transport cost` <- replace_na(file$`transport cost`,0)

file$`transport cost` <- as.numeric(file$`transport cost`)

# filter to cost > 0
file <- file %>%
  filter(`transport cost` > 0) %>%
  #filter(`Mode of transport` != "Nil Return") %>%
  filter (`Mode of transport` == "Rail" | `Mode of transport` == "Car") #%>%
  #filter(`Class of travel` == 'Standard')

# make start date a date
file$`Start date of trip` <- as.Date(file$`Start date of trip`,format = "%d/%m/%Y") # UK date format

# to make things easier for machine learning, let's remove Mode of transport categories with multiple entries

# first we need to find instances with commas
file$`find_text` <- grepl(",",file$`Mode of transport`) # returns TRUE or FALSE

# filter to only FALSE and drop the columns not needed for machine learning
file <- file %>%
  filter(`find_text`==FALSE)

file <- select(file,-`find_text`, -`Accommodation/Meals ()`, -`total cost`, -`Senior Officials Name`)

# now set the data up for machine learning

# Fix the random numbers by setting the seed 
# This enables the analysis to be reproducible when random numbers are used 
set.seed(222)
# Put 3/4 of the data into the training set 
data_split <- initial_split(file, prop = 0.7)

# Create data frames for the two sets:
train_data <- training(data_split)
train_data <- na.omit(train_data) # remove any NAs
test_data  <- testing(data_split)

# use a tidymodels recipe to predict the transport cost
file_rec <- 
  recipe(`transport cost` ~ ., data = train_data) %>% 
  step_date(`Start date of trip`) %>%
  step_rm(`Start date of trip`) %>%
  step_other(Destination, threshold = 0.01) %>% # remove destination values that occur <1% of the time
  step_dummy(all_nominal_predictors()) %>% # change categorical values to numeric
  #step_dummy(all_nominal(), -all_outcomes()) %>% 
  #step_impute_median(all_predictors()) %>%  
  step_zv(all_predictors())%>% # remove columns with a single value
  #step_nzv(all_numeric_predictors()) %>% # remove columns with a single value
  step_normalize(all_predictors()) # centre and scale

###### decision tree model #####
tree_spec <- decision_tree() %>%
  set_engine("rpart") %>%
  set_mode("regression")

tree_wflow <- 
  workflow() %>% 
  add_model(tree_spec) %>%
  add_recipe(file_rec)

#show hyperparameters
extract_parameter_set_dials(file_rec)

#val set - for candidate models
val_set <- validation_split(file, prop = 0.50, strata = `transport cost`)

#train_folds = vfold_cv(train_data, strata = `transport cost`)

#create 25 candidate models
set.seed(22)
tree_res <- 
  tune_grid(tree_wflow,
            val_set,
            #resamples = train_folds,
            grid = 25,
            #control = control_resamples(save_pred = TRUE),
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(rmse))

# show top 5 of 25 models
tree_res %>% 
  show_best(metric = "rmse")

#select the best model
tree_best <- 
  tree_res %>% 
  select_best(metric = "rmse")
tree_best

tree_res %>% 
  collect_predictions()

tree_auc1 <- 
  tree_res %>% 
  collect_predictions(parameters = tree_best)

# finalise the model for fitting
final_mlp_wflow <- tree_wflow %>%
  finalize_workflow(tree_best)

final_mlp_fit <- final_mlp_wflow %>%
  fit(train_data) # fit on train_data

feeder_vip <-
  extract_fit_parsnip(final_mlp_fit) %>% 
  vip()
feeder_vip

### or

#feeder_vip <-
#  extract_fit_engine(final_mlp_fit) %>%
#  vi()

# view the predicted data
predict(final_mlp_fit, train_data) # predict on new, test or train data

flights_aug <- 
  augment(final_mlp_fit, train_data)

# graph of test predictions vs actual data: 
#create line graph
ggplot(flights_aug,aes(x=`.pred`,y=`transport cost`))+
geom_point(color="#CD2456")+
#scale_x_date(date_labels="%d/%m/%Y",date_breaks="1 week")+ #can be 1 day, 1 week, 1 month
#geom_text(data=filter(df,`yvalue`==last(`yvalue`)), aes(label=round(`yvalue`,0)),vjust=-0.5)+
#scale_color_manual(values=alpha(c("#CD2456","14022E")),name="xyz")+
#geom_point(aes(shape=`column`))+
#geom_vline(aes(xintercept=as.numeric(Date[c(90)])))+
#geom_vline(aes(xintercept=as.integer(as.POSIXct("2020-05-05"))),linetype=4)+
theme_classic()+
#ylim(0,max((`yaxis`)+5))+
#facet_grid(~col)+
#coord_flip()+
xlim(0,max(flights_aug$.pred,flights_aug$`transport cost`))+
ylim(0,max(flights_aug$.pred,flights_aug$`transport cost`))+
geom_abline(slope = 1, intercept = 0, color="#CD2456", linetype="dashed")+
theme(legend.position="none")+
theme(axis.text.x=element_text(angle=0,size=10))+
labs(x="Train",y="Test",title="Real vs predicted data on Train dataset")

#### repeat for random forest model ############

rfr_model <- 
  rand_forest(mtry = tune(), min_n = tune(), trees = 10000) %>% 
  set_engine("ranger",  importance = "impurity") %>% # for importance plot this line is needed
  #set_engine("ranger",num.threads=8) %>% 
  #set_engine("randomForest") %>%
  #set_mode("unknown")
  #set_mode("classification")
  set_mode("regression")

# use a tidymodels recipe to predict the transport cost
# random forest models do not need dummy values or normalised values
file_rec1 <- 
  recipe(`transport cost` ~ ., data = train_data) %>% 
  step_date(`Start date of trip`) %>%
  step_rm(`Start date of trip`) %>%
  step_other(Destination, threshold = 0.01) %>% # remove destination values that occur <1% of the time
  step_dummy(all_nominal_predictors()) %>% # change categorical values to numeric
  #step_dummy(all_nominal(), -all_outcomes()) %>% 
  step_impute_median(all_predictors()) %>%  # impute missing values with median values
  step_zv(all_predictors())%>% # remove columns with a single value
  #step_nzv(all_numeric_predictors()) %>% # remove columns with a single value
  step_normalize(all_predictors()) # centre and scale
  
rfr_wflow <- 
  workflow() %>% 
  add_model(rfr_model) %>%
  add_recipe(file_rec1)

#show hyperparameters
extract_parameter_set_dials(file_rec1)

#val set - for candidate models
#val_set <- validation_split(file, prop = 0.80, strata = `transport cost`)

extract_parameter_set_dials(rfr_model)

#create 25 candidate models
set.seed(22)
rfr_res <- 
  rfr_wflow %>% 
  tune_grid(val_set,
            grid = 25,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(rmse))

# show top 5 of 25 models
rfr_res %>% 
  show_best(metric = "rmse")

#select the best model
rfr_best <- 
  rfr_res %>% 
  select_best(metric = "rmse")
rfr_best

rfr_res %>% 
  collect_predictions()

rfr_auc1 <- 
  rfr_res %>% 
  collect_predictions(parameters = rfr_best)

# get metrics
#rfr_met <- rfr_best %>%
#  collect_metrics()

# finalise the model for fitting
final_mlp_wflow <- rfr_wflow %>%
  finalize_workflow(rfr_best)

final_mlp_fit <- final_mlp_wflow %>%
  fit(train_data) # fit on train_data

feeder_vip <-
  extract_fit_parsnip(final_mlp_fit) %>% 
  vip()
feeder_vip

### or as a table

feeder_vip <-
  extract_fit_engine(final_mlp_fit) %>%
  vi()

# graph this
feeder_vip[1:10,] %>%
  ggplot(aes(Importance, fct_reorder(Variable, Importance))) + 
  geom_col(fill="#CD2456") +
  labs(y = NULL, title = "Top 10 feature importances") +
  theme(legend.position = "none")+
  theme_classic()

feeder_vip

# view the predicted data
predict(final_mlp_fit, train_data) # predict on new, test or train data

flights_aug <- 
  augment(final_mlp_fit, train_data)

# graph of test predictions vs actual data: 
#create line graph
ggplot(flights_aug,aes(x=`.pred`,y=`transport cost`))+
geom_point(color="#CD2456")+
#scale_x_date(date_labels="%d/%m/%Y",date_breaks="1 week")+ #can be 1 day, 1 week, 1 month
#geom_text(data=filter(df,`yvalue`==last(`yvalue`)), aes(label=round(`yvalue`,0)),vjust=-0.5)+
#scale_color_manual(values=alpha(c("#CD2456","14022E")),name="xyz")+
#geom_point(aes(shape=`column`))+
#geom_vline(aes(xintercept=as.numeric(Date[c(90)])))+
#geom_vline(aes(xintercept=as.integer(as.POSIXct("2020-05-05"))),linetype=4)+
theme_classic()+
#ylim(0,max((`yaxis`)+5))+
#facet_grid(~col)+
#coord_flip()+
xlim(0,max(flights_aug$.pred,flights_aug$`transport cost`))+
ylim(0,max(flights_aug$.pred,flights_aug$`transport cost`))+
geom_abline(slope = 1, intercept = 0, color="#CD2456", linetype="dashed")+
theme(legend.position="none")+
theme(axis.text.x=element_text(angle=0,size=10))+
labs(x="Train",y="Test",title="Real vs predicted data on Train dataset")

```