
# compare models using cross validation

# --------------------------------------------------------------------------------------------------------------------------------
# functions

merge_all <- function(dflist, by.var) {
    Reduce(function(...) merge(..., by = by.var, all = TRUE), dflist)
}

"%!in%" <- Negate("%in%")

# --------------------------------------------------------------------------------------------------------------------------------
# options

options(stringsAsFactors = FALSE)

# --------------------------------------------------------------------------------------------------------------------------------
# packages

# installpackages("caret", dependencies = c("Imports", "Depends", "Suggests"))
packages <- c("caret", "e1071", "subselect", "ipred", "parallel", "doParallel", "corrplot",  "rowr", "ggplot2") 
sapply(packages, require, character.only = TRUE)


# --------------------------------------------------------------------------------------------------------------------------------
# load data

project_dir <- "path-to-output-directory"
data_dir <- "path-to-input-directory"

# load pre-processed data
load(file.path(data_dir, "data_PP.Rdata"))


# --------------------------------------------------------------------------------------------------------------------------------
# tuning parameters

# https://machinelearningmastery.com/tune-machine-learning-algorithms-in-r/

modelMethods <- c("cubist", "rf", "xgbLinear", "rqnc", "gamSpline", "penalized", "BstLm", 
                                "simpls", "widekernelpls", "glmnet", "gaussprPoly", "pcr", "lm")
tuningParams <- lapply(modelMethods, modelLookup)
names(tuningParams) <- modelMethods
tuningParams

# --------------------------------------------------------------------------------------------------------------------------------
control <- trainControl(method = "repeatedcv", number = 10, repeats = 10, allowParallel = TRUE)
metric <- "RMSE"

# --------------------------------------------------------------------------------------------------------------------------------
# data for PD.s.ALL response
dat_PD.s.ALL <- cv_data_PP$PD.s.ALL_PP$dataset_PP[, colnames(cv_data_PP$PD.s.ALL_PP$dataset_PP) %!in% 
    c("fips", "region",  paste("totalOntree", c("n", "i", "e", "int"), sep = "."), responseColumns[responseColumns %!in% "PD.s.ALL"])]
   
# --------------------------------------------------------------------------------------------------------------------------------
# parallel processing
cluster <- makeCluster(detectCores() - 1) # convention to leave 1 core for OS
registerDoParallel(cluster)

# -------------------------------------------------------------------------------------------------------------------------------- 

cat("starting model fitting", "\n")
   
set.seed(7)
cv_cubist_PD.s.ALL <- train(PD.s.ALL ~ . , method = "cubist", metric = metric, trControl = control, na.action = "na.exclude", 
    tuneGrid = data.frame(committees = 100, neighbors = 9),
    data = dat_PD.s.ALL)

set.seed(7)
mtry <- floor(sqrt(ncol(dat_PD.s.ALL[, !(colnames(dat_PD.s.ALL) %in% "PD.s.ALL")])))
tunegrid <- expand.grid(.mtry=mtry) # grid or random search?
cv_rf_PD.s.ALL <- train(PD.s.ALL ~ . , method = "rf", metric = metric, trControl = control, na.action = "na.exclude", 
    tuneGrid = tunegrid,
    data = dat_PD.s.ALL)
    
set.seed(7)
cv_xgbLinear_PD.s.ALL <- train(PD.s.ALL ~ . , method = "xgbLinear", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)
    
set.seed(7)
cv_rqnc_PD.s.ALL <- train(PD.s.ALL ~ . , method = "rqnc", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)

set.seed(7)
cv_gamSpline_PD.s.ALL <- train(PD.s.ALL ~ . , method = "gamSpline", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)

set.seed(7)
cv_penalized_PD.s.ALL <- train(PD.s.ALL ~ . , method = "penalized", metric = metric, trControl = control, na.action = "na.exclude", 
    tuneGrid = data.frame(lambda1 = 0, lambda2 = 1),
    data = dat_PD.s.ALL)

set.seed(7)
cv_BstLm_PD.s.ALL <- train(PD.s.ALL ~ . , method = "BstLm", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)

set.seed(7)
cv_simpls_PD.s.ALL <- train(PD.s.ALL ~ . , method = "simpls", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)

set.seed(7)
cv_widekernelpls_PD.s.ALL <- train(PD.s.ALL ~ . , method = "widekernelpls", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)

set.seed(7)
cv_glmnet_PD.s.ALL <- train(PD.s.ALL ~ . , method = "glmnet", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)

set.seed(7)
cv_gaussprPoly_PD.s.ALL <- train(PD.s.ALL ~ . , method = "gaussprPoly", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)

set.seed(7)
cv_pcr_PD.s.ALL <- train(PD.s.ALL ~ . , method = "pcr", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)

set.seed(7)
cv_lm_PD.s.ALL <- train(PD.s.ALL ~ . , method = "lm", metric = metric, trControl = control, na.action = "na.exclude", 
    data = dat_PD.s.ALL)
                    
cat("finished model fitting", "\n")


# de-register parallel processing cluster
stopCluster(cluster)
registerDoSEQ()


# --------------------------------------------------------------------------------------------------------------------------------
# save output

cv_models_comp <- list(
    "cubist" = cv_cubist_PD.s.ALL,
    "rf" = cv_rf_PD.s.ALL,
    "xgbLinear" = cv_xgbLinear_PD.s.ALL, 
    "rqnc" = cv_rqnc_PD.s.ALL,
    "gamSpline" = cv_gamSpline_PD.s.ALL,
    "penalized" = cv_penalized_PD.s.ALL,
    "BstLm" = cv_BstLm_PD.s.ALL,
    "simpls" = cv_simpls_PD.s.ALL,
    "widekernelpls" = cv_widekernelpls_PD.s.ALL,
    "glmnet" = cv_glmnet_PD.s.ALL,
    "gaussprPoly" = cv_gaussprPoly_PD.s.ALL,
    "pcr" = cv_pcr_PD.s.ALL,
    "lm" = cv_lm_PD.s.ALL
   )

save(cv_models_comp, file = file.path(project_dir, "cv_models_comp_PD.s.ALL.Rdata"), compress = "gzip")

    
