library(tidyverse)

library(magrittr)

library(lubridate)

library(forcats)

library(caret)

library(stringr)

setwd('C:/Users/Laurens/Dropbox/projects/ecmwf/code/ReadingBus/Data/')

claims <- read_csv('busClaims.csv')

apiKey <- readline('Enter API key:')

routes <- c('2','3','4','5','6','7','9','11','12','13','14','15','16','17',

'19a','19b','19c','21','22','23','24','25','26','27','28','33','500') %>%

factor

claims %<>%

mutate(

route = Route %>%

# Take part before slash

str_split('/') %>%

lapply(FUN = function(x) {x[[1]]}) %>%

unlist,

# Remove trailing spaces

route = route %>%

str_split(boundary('word')) %>%

lapply(FUN = function(x) {x[[1]]}) %>%

unlist,

# Remove leading zeros

route = route %>%

str_split('^[0]{1,9}') %>%

lapply(FUN = function(x) {x[[length(x)]]}) %>%

unlist %>% factor,

route = route %>%

fct_recode(

`2` = '2A',

`4` = 'X4',

`6` = '6A',

`27` = '29')) %>%

# Create proper date field

mutate(dates = `Accident Date` %>% as.Date(format = '%d/%m/%y')) %>%

# Only claims in 2015 and 2016

#filter(year(dates) > 2014 & year(dates) < 2017) %>%

# Only valid routes

filter(route %in% routes) %>%

# Make bad flag

mutate(claimed = 1) %>%

select(one_of('dates', 'route', 'claimed'))

modeldata <- data_frame(dates = seq(ymd('2001-01-01'), ymd('2016-12-31'), by = 'days')) %>%

# Add routes

cbind(setNames(lapply(routes, function(x) x=NA), routes)) %>%

# One row per route per day

gather(key = route, value = claim, -dates) %>%

# Set claims to zero

mutate(claim = 0) %>%

# Merge claims data

full_join(claims, by = c('dates', 'route')) %>%

# Get number of claims per day per route

group_by(dates, route) %>%

summarise(claims = sum(claim, na.rm = TRUE) + sum(claimed, na.rm = TRUE)) %>%

# Make factor binary numeric

mutate(claims = ifelse(claims > 0, 1, 0)) %>%

ungroup

modeldata %>%

filter(claims == 'yes') %>%

ggplot(aes(x = dates)) +

geom_histogram(binwidth = 1)

envidata_orig <- read_rds('D:/projects/ecmwf/evidata.rds') %>%

mutate(dates = timestamp %>% as.Date) %>%

unique

envidata <- envidata_orig %>%

dplyr::select(dates) %>% unique %>%

# 6 am data

left_join(

envidata_orig %>%

filter(hour(timestamp) == 6 | hour(timestamp) == 7) %>%

set_colnames(str_c('v06h_', colnames(.))) %>%

dplyr::select(-v06h_timestamp),

by = c('dates' = 'v06h_dates')) #%>%

# # 12 pm data

# left_join(

# envidata_orig %>%

# filter(hour(timestamp) == 12 | hour(timestamp) == 13) %>%

# set_colnames(str_c('v12h_', colnames(.))) %>%

# dplyr::select(-v12h_timestamp),

# by = c('dates' = 'v12h_dates')) %>%

# # 6 pm data

# left_join(

# envidata_orig %>%

# filter(hour(timestamp) == 18 | hour(timestamp) == 19) %>%

# set_colnames(str_c('v18h_', colnames(.))) %>%

# dplyr::select(-v18h_timestamp),

# by = c('dates' = 'v18h_dates'))

modeldata <- modeldata %>%

# Join to original bus data

left_join(envidata, by = c('dates' = 'dates')) %>%

# Change flag variable

mutate(claims = ifelse(claims > 0, 'yes', 'no') %>% factor,

# Change route variable

route = route %>% factor,

# Get weekday

weekday = dates %>% weekdays %>% factor,

year = dates %>% year %>% factor) %>%

# Keep only complete cases

filter(complete.cases(.))

set.seed(2345)

kfolds <- createFolds(modeldata$claims, k = 3)

data_train <- modeldata[-kfolds[[3]],] %>%

dplyr::select(-matches('timestamp'), -matches('dates'))

data_test <- modeldata[kfolds[[3]],] %>%

select(-matches('timestamp'), -matches('dates'))

mod_tune <- expand.grid(

alpha = c(0, 0.5, 1),

lambda = c(0.000001, 0.00003))

mod_eval <- trainControl(

method = "repeatedcv",

number = 5,

repeats = 1,

classProbs = TRUE,

verboseIter = TRUE,

summaryFunction = twoClassSummary)

mod_train <- caret::train(

claims ~ .,

data = data_train,

method = "glmnet",

family = "binomial",

metric = "ROC",

tuneGrid = mod_tune,

trControl = mod_eval)

write_rds(mod_train, path = 'D:/projects/ecmwf/trained_big.rds')

varImp(mod_train) %>% plot(scale = FALSE)

plot(varImp(mod_train, scale = FALSE))

mod_pred <- predict(mod_train, newdata = data_test, type = 'prob')

mod_pred_val <- bind_cols(mod_pred, data_test) %>% select(yes, claims)

ModelMetrics::auc(actual = mod_pred_val$claims, predicted = mod_pred_val$yes)

library(AUC)

roc(mod_pred_val$yes, mod_pred_val$claims) %>% plot