Telemedicine_ML_Workflow
Full text
# Optimizing Telehealth-Based Cancer Pain Management through Machine Learning ‘ suggestion: run DATA CREATION and PREPROCESSING together with TRAINING PART; then save as RData the environnement and proceed to seperately run the SENSIBILLITY ANALYSIS parts (“adding weights” without the TRAINING part and “removing cohort” + TRAINING), and separately save as RData. In fact, the “adding weights” consists of their training routine. Then you can implement for each of them the SELECTION + COMPARISON parts. ‘ # R packages #### require(readxl) require(xlsx) # for data manipulation require(janitor) library(dplyr) library(tidyr) require(tibble) require(stringr) require(stringi) # to create partition require(caret) # ML algorithms require(randomForest) require(gbm) require(class) require(e1071) require(neuralnet) # ROC and ROC statistics require(purrr) require(pROC) # table outputs library(gtsummary) library(flextable) # plot outputs require(ggplot2) require(cowplot) #### DATA CREATION #### setwd("***") # your work directory dir.create( "./tables" ) doouts <- 1 # 1 = print outputs dato <- list() dato[[1]] <- clean_names(read_excel("./data/repository_RC_29 dic.xlsx"))
dato[[2]] <- clean_names(read_excel("./data/dataset telemedicina aggiornato.xlsx")) dato_0 <- dato data_temp <- clean_names(read_excel("./data/dataset telemedicina aggiornato.xlsx", sheet = "ID")) data_temp <- cbind(dato[[2]], id = data_temp$id) data_temp$cancer_type <- tolower(data_temp$cancer_type) data_temp$pain_type <- ifelse(data_temp$pain_type == "NOCICEPTIVE", 0, 1) data_temp$metastasis <- ifelse(data_temp$metastasis == "PRESENCE", 1, 0) data_temp$b_tc_p <- ifelse(data_temp$b_tc_p == "PRESENCE", 1, 0) data_temp$bone_metastasis <- ifelse(data_temp$bone_metastasis == "PRESENCE", 1, 0) table(data_temp$cancer_type) temp <- case_when(data_temp$cancer_type == "colorectal" | data_temp$cancer_type == "gastric" ~ "gastrointestinal", data_temp$cancer_type == "ovarian" | data_temp$cancer_type == "uterin" ~ "gynecological", data_temp$cancer_type == "pancreatic" | data_temp$cancer_type == "hepathic" ~ "endocrine", data_temp$cancer_type == "bladder" ~ "bladder", data_temp$cancer_type == "breast" ~ "breast", data_temp$cancer_type == "lung" ~ "lung", data_temp$cancer_type == "utering" ~ "gynecological", data_temp$cancer_type == "prostate" ~ "prostate", data_temp$cancer_type == "other" ~ "others" ) cbind(data_temp$cancer_type, temp) table(data_temp$cancer_type, temp) data_temp$cancer_type <- temp vars <- c("cancer_type", "pain_type", "ecog") temp <- data_temp[, c("id", vars)] for(j in data_temp$id){ for(s in vars){ temp2 <- unique(data_temp[data_temp$id == j, s])
if( (NA %in% temp2) & (length(temp2) > 1) ){ temp[(temp$id == j) & is.na(temp[, s]), s] <- na.omit(temp2)[length(na.omit(temp2))] } } } data_temp[, c("id", vars)] <- temp data_temp <- data_temp %>% select(c(id, date, gender, age, ecog, disease, metastasis, bone_metastasis, cancer_type, b_tc_p, pain_type, roo)) %>% group_by(id) %>% mutate(n_televisit = row_number(), age = first(age), ecog_st = first(ecog), ecog_end = last(ecog), neoplasia_st = first(cancer_type), neoplasia_end = last(cancer_type), metastasis_st = first(metastasis), metastasis_end = last(metastasis), bone_metastasis_st = first(bone_metastasis), bone_metastasis_end = last(bone_metastasis), btcp_st = first(b_tc_p), btcp_end = last(b_tc_p), neuro_pain_st = first(pain_type), neuro_pain_end = last(pain_type), roo_st = first(roo), roo_end = last(roo), registration = 2 ) %>% rename(neuro_pain = pain_type) %>% select(!c("cancer_type", "b_tc_p", "roo")) %>% mutate(id = tolower(id), neoplasia_end = ifelse(neoplasia_end == neoplasia_st, 0, 1) ) %>% filter(disease == "CANCER") %>% select(!disease) data_temp_wide <- data_temp %>% pivot_wider( id_cols = id, names_from = n_televisit, values_from = c(date, ecog), ) data_temp <- data_temp %>% select(!c("date", "ecog", )) %>%
group_by(id) %>% filter(n_televisit == max(n_televisit)) data_temp <- merge(x = data_temp, y = data_temp_wide, by = "id") dato[[2]] <- data_temp data_temp <- dato[[1]] names(dato[[2]])[str_detect(string = names(dato[[2]]), pattern = "date_|ecog_", negate = T)] names(data_temp) data_temp <- data_temp %>% rename_with(~ stringr::str_replace(.x, pattern = "data", replacement = "date") ) %>% rename_with(~ stringr::str_replace(.x, pattern = "_visita", replacement = "" ) ) %>% rename(gender = sesso, age = eta, bone_mtx = mtx_osso, btcp_st = btcp, ecog_st = ecog, neuro_pain_st = dolore_neurop, neoplasia_st = neoplasie2, roo_st = rmed, n_televisit = n_visite ) %>% rename_with(~ stringr::str_replace(.x, pattern = "mtx", replacement = "metastasis_st") ) %>% mutate(id = paste0("1_", x1), date_18 = NA, roo_st = ifelse(roo_st == ">60", 1, 0), ecog_st = as.numeric(ecog_st), neoplasia_st = tolower(neoplasia_st), registration = 1 ) %>% mutate(neoplasia_st = if_else( neoplasia_st == "pancreas", "endocrine", if_else(neoplasia_st == "bledder", "bladder", neoplasia_st) )) %>% mutate(across(contains("date"), ~ as.Date(., format = "%d/%m/%Y"))) %>% select(c(id, age, gender, neoplasia_st, ecog_st, metastasis_st, bone_metastasis_st, btcp_st, neuro_pain_st, roo_st, n_televisit, registration, contains("date")))
names(data_temp) names(dato[[2]]) dato[[1]] <- data_temp dato <- bind_rows(dato) row_names <- gsub(dato$id, pattern = "'", replacement = "_") cbind(row_names, dato$id) dato <- dato %>% relocate(starts_with("date"), .before = "ecog_1") %>% relocate(ends_with("_end"), .after = n_televisit) %>% relocate(registration, .after = n_televisit) %>% select(!id) %>% mutate_at(c("gender", "neoplasia_st", "neoplasia_end", "registration", "roo_st", "roo_end", "metastasis_st", "metastasis_end", "bone_metastasis_st", "bone_metastasis_end", "btcp_st", "btcp_end", "neuro_pain_st", "neuro_pain_end"), as.factor) %>% mutate(gender = factor(gender, levels = c("M", "F"), labels = c("Female", "Male")), neuro_pain_st = factor(neuro_pain_st, levels = c("0", "1"), labels = c("Nociceptive", "Neuropathic/mixed")), ecog_cod_st = as.factor(ifelse(ecog_st < 3, "1-2", "3-4")), bone_metastasis_st = case_when( metastasis_st == 0 & bone_metastasis_st == 1 ~ 0, bone_metastasis_st == 0 ~ 0, bone_metastasis_st == 1 ~ 1 ) ) %>% mutate_at(c("roo_st", "roo_end", "metastasis_st", "metastasis_end", "bone_metastasis_st", "bone_metastasis_end", "btcp_st", "btcp_end", "neuro_pain_end"), .funs = recode_factor, "0" = " No", "1" = " Yes") row.names(dato) <- row_names ##### repository ##### if(doouts == 1){ write.xlsx(x = cbind(id = row.names(dato), dato), file = "./data/repository_telemedicina_25.xlsx", sheetName = "merged data", row.names = F ) } #### UNIVARIABLE ANALYSIS ####
table(dato$neoplasia_st, dato$registration) tabella <- list() tabella[[1]] <- dato %>% mutate(n_televisits_cod = cut(n_televisit, breaks = c(0, 1, max(n_televisit, na.rm = T)), include.lowest = T, right = T, labels = c("One", "More")) ) %>% select(!matches("_\\d")) %>% tbl_summary( percent = "column", type = list(all_continuous() ~ "continuous2", ecog_st ~ "continuous2", ecog_end ~ "continuous2" ), statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})" ), digits = c(everything() ~ 1, all_categorical() ~ c(0, 1)), missing_text = "(Missing)", missing = "ifany") %>% #modify_caption("Table 1: main CP characteristics") %>% bold_labels() %>% italicize_levels() tabella[[1]] tabella[[2]] <- dato %>% mutate(n_televisits_cod = cut(n_televisit, breaks = c(0, 1, max(n_televisit, na.rm = T)), include.lowest = T, right = T, labels = c("One", "More")) ) %>% select(!matches("_\\d")) %>% tbl_summary( by = "registration", percent = "column", type = list(all_continuous() ~ "continuous2", ecog_st ~ "continuous2", ecog_end ~ "continuous2" ), statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})" ), digits = c(everything() ~ 1, all_categorical() ~ c(0, 1)), missing_text = "(Missing)",missing = "ifany") %>% add_p() %>% add_q(method = "BH") %>% modify_column_hide(p.value) %>% bold_labels() %>% italicize_levels() tabella[[2]]
tabella <- tbl_merge(tbls = tabella, tab_spanner = c("Descriptive statistics", "By dataset")) %>% as_flex_table() tabella ##### final data ##### dato %>% group_by(registration) %>% summarise(date_1 = min(date_1), date_16 = max(date_16, na.rm = T), date_18 = max(date_18, na.rm = T)) dato2 <- dato %>% rownames_to_column('id') %>% mutate(n_televisit = cut(n_televisit, breaks = c(0, 1, max(n_televisit, na.rm = T)), include.lowest = T, right = T, labels = c("One", "More")) ) %>% select(c("id", "n_televisit", "age", "gender", "registration", contains("_st"))) %>% rename_with(.cols = contains("_st"), ~ str_remove(.x, "_st")) %>% column_to_rownames('id') list_var <- list("n_televisit" = "Televisits", "age" = "Age", "gender" = "Sex", "registration" = "Registration", "neoplasia" = "Tumor site", "ecog" = "ECOG", "ecog_cod" = "ECOG (class)", "metastasis" = "Metastasis", "bone_metastasis" = "Bone metastasis", "btcp" = "BTCP status", "neuro_pain" = "Type of pain", "roo" = "ROO" ) tabella <- list() p_adj <- function(y, data){ k <- 0 p_mat <- data.frame(variable = NA, p_vals = NA) data <- data %>% as.data.frame() for(j in c(1:ncol(data))[- c(y)]){
k <- k + 1 temp <- na.omit(data[, c(y, j)]) p_mat[k, 1] <- names(data)[j] z = temp[, 1] x = temp[, 2] if(is.factor(z)){ if( is.factor(data[, j])){ cat(names(data)[j], "with Pearson's X2 \n") p_mat[k, 2] <- chisq.test(x = x, y = z)$p.value }else{ cat(names(data)[j], "with Kruskal \n") p_mat[k, 2] <- kruskal.test(z ~ x)$p.value } }else{ if( is.factor(data[, j]) ){ cat(names(data)[j], "with Kruskal \n") p_mat[k, 2] <- kruskal.test(z ~ x)$p.value }else{ cat(names(data)[j], "with Spearman \n") p_mat[k, 2] <- cor.test(x = x, y = z, method = "spearman")$p.value } } } return(p_mat) } p_vals <- p_adj(y = 4, data = dato2) p_vals$p_vals[c(4, 5)] <-NA; p_vals # not including neoplasia and ECOG (continous) among tests p_vals$p_vals <- p.adjust(p = p_vals$p_vals, method = "BH") tabella <- dato2 %>% tbl_summary( by = "registration", label = list_var, percent = "column", type = list(all_continuous() ~ "continuous2", ecog ~ "continuous2"), statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})" ), digits = c(everything() ~ 1, all_categorical() ~ c(0, 1)), missing_text = "(Missing)", missing = "ifany") %>% add_p(pvalue_fun = ~style_sigfig(., digits = 3)) %>% bold_labels() %>% modify_spanning_header(all_stat_cols() ~ "**Cohort**") %>% modify_table_styling( columns = label, rows = label == c("Tumor site", "ECOG", "Metastasis"),
footnote = "not considered for the Benjamini-Hochberg correction" ) %>% italicize_levels() tabella old_p <- tabella$table_body[c("variable", "p.value")] %>% rowwise() %>% mutate(p.value = case_when(variable %in% c("ecog", "neoplasia") ~ NA, T ~ p.value) ) %>% mutate(x = !is.na(variable) & !is.na(p.value)) k <- 0 for( j in 1:length(old_p) ){ if( isTRUE(old_p$x[j]) ){ k <- k + 1 old_p[j, "p.value"] <- as.numeric(p_vals$p_vals[k]) } } tabella$table_body$p.value <- old_p$p.value tabella if(doouts == 1){ save_as_docx(tabella %>% as_flex_table(), path = "./tables/table_s1.docx") } tabella <- list() tabella[[1]] <- dato2 %>% tbl_summary( percent = "column", label = list_var, type = list(all_continuous() ~ "continuous2", ecog ~ "continuous2" ), statistic = all_continuous() ~ c("{mean} ({sd})", "{median} ({p25}, {p75})" ), digits = c(everything() ~ 1, all_categorical() ~ c(0, 1)), missing_text = "Unknown",missing = "ifany" ) %>% bold_labels() %>% italicize_levels() tabella[[1]] p_vals <- p_adj(y = 1, data = dato2) p_vals$p_vals[c(4, 5)] <-NA # not including neoplasia and ECOG (continous) among tests p_vals$p_vals <- p.adjust(p = p_vals$p_vals, method = "BH")
for (node_size in node_size_grid){ for(r_shrink in r_shrink_grid){ cat("\n n_tree = ", n_tree, "\n r_shrink = ", r_shrink, "\n node_size = ", node_size, "\n rep #", n, "\n") conto <- conto + 1 data_matrix[conto, 1: (2 + length(mod_params))] <- c(r_shrink, n_tree, node_size, n, k) data_matrix <- rbind(data_matrix, NA) set.seed(1234 + k_cv * (n - 1) + k); gbm_fit <- gbm(data = train_data_2, formula = as.numeric(n_televisit) - 1 ~ ., distribution = "bernoulli", n.trees = n_tree, interaction.depth = 2, n.minobsinnode = node_size, shrinkage = r_shrink ) predicted <- predict(gbm_fit, newdata = test_data_2[,-1], type = type_response) if(T %in% predicted == 0.5){ set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5) } predicted <- factor(predicted > 0.5, levels = c(FALSE, TRUE), labels = levels(train_data[, y])) observed <- factor(test_data_2$n_televisit, levels = levels(train_data[, y])) conf_mat <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2]) if(nlevels(dato3[,y]) > 2){ f1_temp <- mean(conf_mat$byClass[, "F1"], na.rm = T) conf_mat_distr <- paste(levels(dato3[, y]), conf_mat$byClass[, "F1"], sep = " = ", collapse = ", ") empty_class <- (conf_mat$byClass[, "Recall"] %in% c(0, NA, NaN) ) data_matrix[conto, "empty_class"] <- paste(levels(dato3[, y]), empty_class, sep = " = ", collapse = ", ")
}else{ f1_temp <- conf_mat$byClass["F1"] conf_mat_distr <- paste(c("TN", "FN", "FP", "TP"), conf_mat$table, sep = " = ", collapse = ", ") empty_class <- ifelse(0 %in% diag(conf_mat$table), 1, 0) data_matrix[conto, "empty_class"] <- empty_class } data_matrix[conto, "f1"] <- f1_temp #paste(levels(dato3[, y]), f1_temp, sep = " = ", collapse = ", ") } } } } } data_matrix_list[[2]] <- data_matrix t_gbm <- difftime(time1 = Sys.time(), time2 = t0, units = 'mins') ##### K-nearest neighbors ##### t0 <- Sys.time() cat("**** START KNN ****") n_neigh_grid <- 2:5 mod_params <- "n_neigh_grid" data_matrix <- matrix(NA, ncol = 4 + length(mod_params)); colnames(data_matrix) <- c(mod_params, c("n_rep", "k_cv", "f1", "empty_class")) type_response <- "response" conto <- 0 for(n in 1:n_rep){ for(k in 1:k_cv){ set.seed(1234 + k_cv * (n - 1) + k); train_data_2 <- createDataPartition(y = train_data$y_and_registration, p = 1 - 1/k_cv, list = F) test_data_2 <- train_data[- train_data_2, - ncol(train_data)] train_data_2 <- train_data[train_data_2, - ncol(train_data)] for(n_neigh in n_neigh_grid){
cat("\n n_neighbors = ", n_neigh, "\n rep #", n, "\n") conto <- conto + 1 data_matrix[conto, 1: (2 + length(mod_params))] <- c(n_neigh, n, k) data_matrix <- rbind(data_matrix, NA) set.seed(1234 + k_cv * (n - 1) + k); knn_fit <- knnreg(x = train_data_2[,-1], y = as.numeric(train_data_2$n_televisit) - 1, k = n_neigh, l = 0, prob = F, use.all = T ) predicted <- predict(knn_fit, newdata = test_data_2[,-1], type = type_response) if(T %in% predicted == 0.5){ set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5) } predicted <- factor(predicted > 0.5, levels = c(FALSE, TRUE), labels = levels(train_data[, y])) observed <- factor(test_data_2$n_televisit, levels = levels(train_data[, y])) conf_mat <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2]) if(nlevels(dato3[,y]) > 2){ f1_temp <- mean(conf_mat$byClass[, "F1"], na.rm = T) conf_mat_distr <- paste(levels(dato3[, y]), conf_mat$byClass[, "F1"], sep = " = ", collapse = ", ") empty_class <- (conf_mat$byClass[, "Recall"] %in% c(0, NA, NaN) ) data_matrix[conto, "empty_class"] <- paste(levels(dato3[, y]), empty_class, sep = " = ", collapse = ", ") }else{ f1_temp <- conf_mat$byClass["F1"] conf_mat_distr <- paste(c("TN", "FN", "FP", "TP"), conf_mat$table, sep = " = ", collapse = ", ") empty_class <- ifelse(0 %in% diag(conf_mat$table), 1, 0) data_matrix[conto, "empty_class"] <- empty_class
} data_matrix[conto, "f1"] <- f1_temp #paste(levels(dato3[, y]), f1_temp, sep = " = ", collapse = ", ") } } } data_matrix_list[[3]] <- data_matrix t_knn <- difftime(time1 = Sys.time(), time2 = t0, units = 'mins') ##### Support vector machine ##### t0 <- Sys.time() cat("**** START SVM ****") mod_params <- c("g_gamma", "c_cost") data_matrix <- matrix(NA, ncol = length(mod_params) + 4); colnames(data_matrix) <- c(mod_params, c("n_rep", "k_cv", "f1", "empty_class") ) g_gamma_grid <- 2 ^ seq(-15, 3, 2) c_cost_grid <- 2^ seq(-5, 15, 2) conto <- 0 for(n in 1:n_rep){ for(k in 1:k_cv){ set.seed(1234 + k_cv * (n - 1) + k); train_data_2 <- createDataPartition(y = train_data$y_and_registration, p = 1 - 1/k_cv, list = F) test_data_2 <- train_data[- train_data_2, - ncol(train_data)] train_data_2 <- train_data[train_data_2, - ncol(train_data)] f1_temp <- NA type_response <- "response" for(c_cost in c_cost_grid){ for(g_gamma in g_gamma_grid){ cat("\n r_shrink = ", r_shrink, "\n node_size = ", g_gamma, "\n cost = ", c_cost, "\n rep #", n, "\n") conto <- conto + 1
data_matrix[conto, 1: (2 + length(mod_params))] <- c(g_gamma, c_cost, n, k) data_matrix <- rbind(data_matrix, NA) set.seed(1234 + k_cv * (n - 1) + k); svm_fit <- svm(x = train_data_2[, -1], y = as.numeric(train_data_2$n_televisit) - 1, scale = F, kernel = "radial", gamma = g_gamma, cost = c_cost ) predicted <- predict(svm_fit, newdata = test_data_2[,-1]) if(T %in% predicted == 0.5){ set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5) } predicted <- factor(predicted > 0.5, levels = c(F, T), labels = levels(train_data[, y])) observed <- factor(test_data_2$n_televisit, levels = levels(train_data[, y])) conf_mat <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2]) if(nlevels(dato3[,y]) > 2){ f1_temp <- mean(conf_mat$byClass[, "F1"], na.rm = T) conf_mat_distr <- paste(levels(dato3[, y]), conf_mat$byClass[, "F1"], sep = " = ", collapse = ", ") empty_class <- (conf_mat$byClass[, "Recall"] %in% c(0, NA, NaN) ) data_matrix[conto, "empty_class"] <- paste(levels(dato3[, y]), empty_class, sep = " = ", collapse = ", ") }else{ f1_temp <- conf_mat$byClass["F1"] conf_mat_distr <- paste(c("TN", "FN", "FP", "TP"), conf_mat$table, sep = " = ", collapse = ", ") empty_class <- ifelse(0 %in% diag(conf_mat$table), 1, 0) data_matrix[conto, "empty_class"] <- empty_class }
data_matrix[conto, "f1"] <- f1_temp #paste(levels(dato3[, y]), f1_temp, sep = " = ", collapse = ", ") } } } } t_svm <- difftime(time1 = Sys.time(), time2 = t0, units = 'mins') data_matrix_list[[4]] <- data_matrix ##### Multilayer perceptron ##### t0 <- Sys.time() cat("**** START MPL ****") mod_params <- c("hidden1", "hidden2", "learning_rate") data_matrix <- matrix(NA, ncol = length(mod_params) + 4); colnames(data_matrix) <- c(mod_params, c("n_rep", "k_cv", "f1", "empty_class") ) l_rate_grid <- c(0.001, 0.005, seq(10, 200, 20) / 1000 ) hidden_grid <- 1:5 conto <- 0 for(n in 1:n_rep){ for(k in 1:k_cv){ set.seed(1234 + k_cv * (n - 1) + k); train_data_2 <- createDataPartition(y = train_data$y_and_registration, p = 1 - 1/k_cv, list = F) test_data_2 <- train_data[- train_data_2, - ncol(train_data)] train_data_2 <- train_data[train_data_2, - ncol(train_data)] test_data_2$n_televisit <- as.numeric(test_data_2$n_televisit) - 1; test_data_2 <- apply(test_data_2, MARGIN = 2, as.numeric) %>% as.data.frame train_data_2$n_televisit <- as.numeric(train_data_2$n_televisit) - 1; train_data_2 <- apply(train_data_2, MARGIN = 2, as.numeric) %>% as.data.frame f1_temp <- NA type_response <- "response" for(l_rate in l_rate_grid){
for(hidden1 in hidden_grid){ hidden2 <- 0 while(hidden2 < hidden1){ if(hidden2 == 0){ hidden2_in <- NULL }else{ hidden2_in <- hidden2 } cat("\n hidden1 = ", hidden1, "\n hidden2 = ", hidden2, "\n learning rate = ", l_rate, "\n rep #", n, "\n") conto <- conto + 1 data_matrix[conto, 1: (2 + length(mod_params))] <- c(hidden1, hidden2, l_rate, n, k) data_matrix <- rbind(data_matrix, NA) set.seed(1234 + k_cv * (n - 1) + k); mlp_fit <- neuralnet(data = train_data_2, formula = n_televisit ~ ., hidden = c(hidden1, hidden2_in), threshold = 10^(-2), learningrate = l_rate, algorithm = "backprop", err.fct = "sse" , act.fct = "tanh", linear.output = F, stepma = 10^5 ) tryCatch({ predicted <- predict(mlp_fit, newdata = test_data_2[,-1]) if(T %in% predicted == 0.5){ set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- factor( rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5), labels = levels(train_data[, y]) ) } predicted <- factor(predicted > 0.5, levels = c(F, T), labels = levels(train_data[, y])) observed <- factor(test_data_2$n_televisit, labels = levels(train_data[, y])) conf_mat <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2])
if(nlevels(dato3[,y]) > 2){ f1_temp <- mean(conf_mat$byClass[, "F1"], na.rm = T) conf_mat_distr <- paste(levels(dato3[, y]), conf_mat$byClass[, "F1"], sep = " = ", collapse = ", ") empty_class <- (conf_mat$byClass[, "Recall"] %in% c(0, NA, NaN) ) data_matrix[conto, "empty_class"] <- paste(levels(dato3[, y]), empty_class, sep = " = ", collapse = ", ") }else{ f1_temp <- conf_mat$byClass["F1"] conf_mat_distr <- paste(c("TN", "FN", "FP", "TP"), conf_mat$table, sep = " = ", collapse = ", ") empty_class <- ifelse(0 %in% diag(conf_mat$table), 1, 0) data_matrix[conto, "empty_class"] <- empty_class } }, error = function(e) { print(message("An error occurred while building previsions: ", e$message)) f1_temp <- NA conf_mat_distr <- "not implemented" empty_class <- 1 data_matrix[conto, "empty_class"] <- empty_class } ) hidden2 <- hidden2 + 1 data_matrix[conto, "f1"] <- f1_temp #paste(levels(dato3[, y]), f1_temp, sep = " = ", collapse = ", ") } } } } } data_matrix_list[[5]] <- data_matrix t_mlp <- difftime(time1 = Sys.time(), time2 = t0, units = 'mins') ##### Logistic regressor #####
glm_opt <- glm(data = train_data[, - ncol(train_data)], formula = n_televisit ~ . , family = binomial(link = "logit")) predicted_glm <- predict(glm_opt, newdata = test_data, type = "response") predicted <- predicted_glm if(T %in% predicted_glm == 0.5){ set.seed(1234); predicted[predicted == 0.5] <- factor( rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5), labels = levels(train_data[, y]) ) } predicted <- factor(predicted > 0.5, levels = c(F, T), labels = levels(train_data[, y])) observed <- factor(test_data$n_televisit, labels = levels(train_data[, y])) conf_mat_glm <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2]) } #### MODEL SELECTION #### { predicted_rf <- NA; predicted_gbm <- NA; predicted_svm <- NA; predicted_knn <- NA; predicted_mlp <- NA model_list <- list("RF" = as.numeric(predicted_rf), "GBM" = as.numeric(predicted_gbm), "KNN" = as.numeric(predicted_knn), "SVM" = as.numeric(predicted_svm), "MLP" = as.numeric(predicted_mlp), "LR" = as.numeric(predicted_glm) ) data_matrix_list_opt <- list() for(n_alg in 1:5){ data_matrix_list_opt[[n_alg]] <- data_matrix_list[[n_alg]] %>% as_tibble %>% select(!c("empty_class") ) %>% group_by(across(-c(k_cv, f1) )) %>% summarise( across( .cols = f1, # mean F1 across K folds .fns = list( mean = ~ mean(., na.rm = TRUE)
), .names = "{.col}_k_{.fn}" ) ) %>% group_by(across(!contains("f1_")) ) %>% summarise( across( .cols = contains("f1"), # mean mean F1 across N reps .fns = list( mean = ~ mean(., na.rm = TRUE) ), .names = "{.col}_r_{.fn}" ) ) %>% as.data.frame() data_matrix_list_opt[[n_alg]] <- data_matrix_list_opt[[n_alg]][which.max( data_matrix_list_opt[[n_alg]]$f1_k_mean_r_mean ), ] names(data_matrix_list_opt)[n_alg] <- names(model_list)[n_alg] } print(data_matrix_list_opt) ##### Random forest ##### n_tree_opt <- data_matrix_list_opt[[1]]$n_tree m_try_opt <- data_matrix_list_opt[[1]]$m_try node_size_opt <- data_matrix_list_opt[[1]]$node_size set.seed(1234); rf_opt <- randomForest(data = train_data[,-ncol(train_data)], n_televisit ~ ., ntree = n_tree_opt, mtry = m_try_opt, nodesize = node_size_opt, maxnodes = NULL, importance = F ) predicted_rf <- predict(rf_opt, newdata = test_data[,- ncol(train_data)], type = "prob")[, 1] predicted <- predicted_rf if(T %in% predicted == 0.5){ set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- factor( rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5), labels = levels(train_data[, y]) ) } predicted <- factor(predicted > 0.5, labels = levels(train_data[, y])) observed <- factor(test_data$n_televisit, labels = levels(train_data[, y]))
train_data_old <- train_data test_data_old <- test_data #### SENSITIVITY ANALYSIS #### ##### adding weights ##### { sens = "weights" c_weights <- NULL counter = data.frame(sort(table(train_data[, "y_and_registration"])), decreasing = T) majority = max(counter$Freq) counter$weight = majority/counter$Freq c_weights_prop = setNames(as.list(counter$weight), nm = counter$Var1) c_weights = as.numeric(c_weights_prop[as.character(train_data[, "y_and_registration"])]); c_weights_train <- data.frame(c_weights, row.names = row.names(train_data)) table(train_data[, "y_and_registration"], c_weights_train$c_weights) #c_weights_train <- c_weights_train$c_weights c_weights = setNames(as.list(counter$weight), nm = counter$Var1) c_weights = as.numeric(c_weights[as.character(test_data[, "y_and_registration"])]); c_weights_test <- data.frame(c_weights, row.names = row.names(test_data)) table(test_data[, "y_and_registration"], c_weights_test$c_weights) #c_weights_test <- c_weights_test$c_weights { n_rep <- 7 k_cv <- 5 data_matrix_list <- list() ### Random forest ### t0 <- Sys.time() cat("**** START RF **** \n") mod_params <- c("n_tree", "m_try", "node_size") data_matrix <- matrix(NA, ncol = length(mod_params) + 4); colnames(data_matrix) <- c(mod_params, c("n_rep", "k_cv", "f1", "empty_class") ) type_response <- "prob" conto <- 0
n_tree_grid <- seq(20, 100, 10) m_try_grid <- 3:6 node_size_grid <- ceiling((1 - 1/k_cv) * dim(train_data)[1]/c(20,15,10)) for(n in 1:n_rep){ for(k in 1:k_cv){ set.seed(1234 + k_cv * (n - 1) + k); train_data_2 <- createDataPartition(y = train_data$y_and_registration, p = 1 - 1/k_cv, list = F) test_data_2 <- train_data[- train_data_2, - ncol(train_data)] train_data_2 <- train_data[train_data_2, - ncol(train_data)] c_weights_2 <- c_weights_train[row.names(train_data_2), "c_weights"] f1_temp <- NA for(n_tree in n_tree_grid){ for(m_try in m_try_grid){ for(node_size in node_size_grid){ cat("\n n_tree = ", n_tree, "\n m_try = ",m_try, "\n node_size = ", node_size, "\n rep #", n, "\n") conto <- conto + 1 data_matrix[conto, 1: (2 + length(mod_params))] <- c( n_tree, m_try, node_size, n, k) data_matrix <- rbind(data_matrix, NA) set.seed(1234 + k_cv * (n - 1) + k); rf_fit <- randomForest(data = train_data_2, n_televisit ~ ., ntree = n_tree, mtry = m_try, nodesize = node_size, maxnodes = NULL, importance = F, weights = c_weights_2 ) predicted <- predict(rf_fit, newdata = test_data_2[,-1], type = type_response)[, 1] if(T %in% predicted == 0.5){ set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5) }
predicted <- factor(predicted > 0.5, levels = c(FALSE, TRUE), labels = levels(train_data[, y])) observed <- factor(test_data_2$n_televisit, levels = levels(train_data[, y])) conf_mat <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2]) if(nlevels(dato3[,y]) > 2){ f1_temp <- mean(conf_mat$byClass[, "F1"], na.rm = T) conf_mat_distr <- paste(levels(dato3[, y]), conf_mat$byClass[, "F1"], sep = " = ", collapse = ", ") empty_class <- (conf_mat$byClass[, "Recall"] %in% c(0, NA, NaN) ) data_matrix[conto, "empty_class"] <- paste(levels(dato3[, y]), empty_class, sep = " = ", collapse = ", ") }else{ f1_temp <- conf_mat$byClass["F1"] conf_mat_distr <- paste(c("TN", "FN", "FP", "TP"), conf_mat$table, sep = " = ", collapse = ", ") empty_class <- ifelse(0 %in% diag(conf_mat$table), 1, 0) data_matrix[conto, "empty_class"] <- empty_class } data_matrix[conto, "f1"] <- f1_temp #paste(levels(dato3[, y]), f1_temp, sep = " = ", collapse = ", ") } } } } } data_matrix_list[[1]] <- data_matrix t_rf <- difftime(time1 = Sys.time(), time2 = t0, units = 'mins') ### Gradient boosting machine ### t0 <- Sys.time() cat("**** START GBM ****") mod_params <- c("r_shrink", "n_tree", "node_size")
data_matrix <- matrix(NA, ncol = length(mod_params) + 4); colnames(data_matrix) <- c(mod_params, c("n_rep", "k_cv", "f1", "empty_class") ) type_response <- "response" r_shrink_grid <- c(0.001, 0.005, seq(10, 100, 10) / 1000 ) n_tree_grid <- seq(20, 100, 10) node_size_grid <- ceiling((1 - 1/k_cv) * dim(train_data)[1]/c(20,15,10)) conto <- 0 for(n in 1:n_rep){ for(k in 1:k_cv){ set.seed(1234 + k_cv * (n - 1) + k); train_data_2 <- createDataPartition(y = train_data$y_and_registration, p = 1 - 1/k_cv, list = F) test_data_2 <- train_data[- train_data_2, - ncol(train_data)] train_data_2 <- train_data[train_data_2, - ncol(train_data)] c_weights_2 <- c_weights_train[row.names(train_data_2), "c_weights"] f1_temp <- NA for(n_tree in n_tree_grid){ for (node_size in node_size_grid){ for(r_shrink in r_shrink_grid){ cat("\n n_tree = ", n_tree, "\n r_shrink = ", r_shrink, "\n node_size = ", node_size, "\n rep #", n, "\n") conto <- conto + 1 data_matrix[conto, 1: (2 + length(mod_params))] <- c(r_shrink, n_tree, node_size, n, k) data_matrix <- rbind(data_matrix, NA) set.seed(1234 + k_cv * (n - 1) + k); gbm_fit <- gbm(data = train_data_2, formula = as.numeric(n_televisit) - 1 ~ ., distribution = "bernoulli", n.trees = n_tree, interaction.depth = 2, n.minobsinnode = node_size, shrinkage = r_shrink, weights = c_weights_2 ) predicted <- predict(gbm_fit, newdata = test_data_2[,-1], type = type_response)
if(T %in% predicted == 0.5){ set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5) } predicted <- factor(predicted > 0.5, levels = c(FALSE, TRUE), labels = levels(train_data[, y])) observed <- factor(test_data_2$n_televisit, levels = levels(train_data[, y])) conf_mat <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2]) if(nlevels(dato3[,y]) > 2){ f1_temp <- mean(conf_mat$byClass[, "F1"], na.rm = T) conf_mat_distr <- paste(levels(dato3[, y]), conf_mat$byClass[, "F1"], sep = " = ", collapse = ", ") empty_class <- (conf_mat$byClass[, "Recall"] %in% c(0, NA, NaN) ) data_matrix[conto, "empty_class"] <- paste(levels(dato3[, y]), empty_class, sep = " = ", collapse = ", ") }else{ f1_temp <- conf_mat$byClass["F1"] conf_mat_distr <- paste(c("TN", "FN", "FP", "TP"), conf_mat$table, sep = " = ", collapse = ", ") empty_class <- ifelse(0 %in% diag(conf_mat$table), 1, 0) data_matrix[conto, "empty_class"] <- empty_class } data_matrix[conto, "f1"] <- f1_temp #paste(levels(dato3[, y]), f1_temp, sep = " = ", collapse = ", ") } } } } } data_matrix_list[[2]] <- data_matrix t_gbm <- difftime(time1 = Sys.time(), time2 = t0, units = 'mins')
### K-nearest neighbors ### t0 <- Sys.time() cat("**** START KNN ****") n_neigh_grid <- 2:5 mod_params <- "n_neigh_grid" data_matrix <- matrix(NA, ncol = 4 + length(mod_params)); colnames(data_matrix) <- c(mod_params, c("n_rep", "k_cv", "f1", "empty_class")) type_response <- "response" conto <- 0 for(n in 1:n_rep){ for(k in 1:k_cv){ set.seed(1234 + k_cv * (n - 1) + k); train_data_2 <- createDataPartition(y = train_data$y_and_registration, p = 1 - 1/k_cv, list = F) test_data_2 <- train_data[- train_data_2, - ncol(train_data)] train_data_2 <- train_data[train_data_2, - ncol(train_data)] c_weights_2 <- c_weights_train[row.names(train_data_2), "c_weights"] f1_temp <- NA for(n_neigh in n_neigh_grid){ cat("\n n_neighbors = ", n_neigh, "\n rep #", n, "\n") conto <- conto + 1 data_matrix[conto, 1: (2 + length(mod_params))] <- c(n_neigh, n, k) data_matrix <- rbind(data_matrix, NA) set.seed(1234 + k_cv * (n - 1) + k); knn_fit <- knnreg(x = train_data_2[,-1], y = as.numeric(train_data_2$n_televisit) - 1, k = n_neigh, l = 0, prob = F, use.all = T, weights = c_weights_2 ) predicted <- predict(knn_fit, newdata = test_data_2[,-1], type = type_response)
if(T %in% predicted == 0.5){ set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5) } predicted <- factor(predicted > 0.5, levels = c(FALSE, TRUE), labels = levels(train_data[, y])) observed <- factor(test_data_2$n_televisit, levels = levels(train_data[, y])) conf_mat <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2]) if(nlevels(dato3[,y]) > 2){ f1_temp <- mean(conf_mat$byClass[, "F1"], na.rm = T) conf_mat_distr <- paste(levels(dato3[, y]), conf_mat$byClass[, "F1"], sep = " = ", collapse = ", ") empty_class <- (conf_mat$byClass[, "Recall"] %in% c(0, NA, NaN) ) data_matrix[conto, "empty_class"] <- paste(levels(dato3[, y]), empty_class, sep = " = ", collapse = ", ") }else{ f1_temp <- conf_mat$byClass["F1"] conf_mat_distr <- paste(c("TN", "FN", "FP", "TP"), conf_mat$table, sep = " = ", collapse = ", ") empty_class <- ifelse(0 %in% diag(conf_mat$table), 1, 0) data_matrix[conto, "empty_class"] <- empty_class } data_matrix[conto, "f1"] <- f1_temp #paste(levels(dato3[, y]), f1_temp, sep = " = ", collapse = ", ") } } } data_matrix_list[[3]] <- data_matrix t_knn <- difftime(time1 = Sys.time(), time2 = t0, units = 'mins') ### Support vector machine ### t0 <- Sys.time()
cat("**** START SVM ****") mod_params <- c("g_gamma", "c_cost") data_matrix <- matrix(NA, ncol = length(mod_params) + 4); colnames(data_matrix) <- c(mod_params, c("n_rep", "k_cv", "f1", "empty_class") ) g_gamma_grid <- 2 ^ seq(-15, 3, 2) c_cost_grid <- 2^ seq(-5, 15, 2) conto <- 0 for(n in 1:n_rep){ for(k in 1:k_cv){ set.seed(1234 + k_cv * (n - 1) + k); train_data_2 <- createDataPartition(y = train_data$y_and_registration, p = 1 - 1/k_cv, list = F) test_data_2 <- train_data[- train_data_2, - ncol(train_data)] train_data_2 <- train_data[train_data_2, - ncol(train_data)] c_weights_2 <- c_weights_train[row.names(train_data_2), "c_weights"] f1_temp <- NA type_response <- "response" for(c_cost in c_cost_grid){ for(g_gamma in g_gamma_grid){ cat("\n r_shrink = ", r_shrink, "\n node_size = ", g_gamma, "\n cost = ", c_cost, "\n rep #", n, "\n") conto <- conto + 1 data_matrix[conto, 1: (2 + length(mod_params))] <- c(g_gamma, c_cost, n, k) data_matrix <- rbind(data_matrix, NA) set.seed(1234 + k_cv * (n - 1) + k); svm_fit <- svm(x = train_data_2[, -1], y = as.numeric(train_data_2$n_televisit) - 1, scale = F, kernel = "radial", gamma = g_gamma, cost = c_cost, weights = c_weights_2 ) predicted <- predict(svm_fit, newdata = test_data_2[,-1]) if(T %in% predicted == 0.5){
set.seed(1234 + k_cv * (n - 1) + k); predicted[predicted == 0.5] <- rbernoulli(n = length(predicted[predicted == 0.5]), p = 0.5) } predicted <- factor(predicted > 0.5, levels = c(F, T), labels = levels(train_data[, y])) observed <- factor(test_data_2$n_televisit, levels = levels(train_data[, y])) conf_mat <- confusionMatrix(table(observed = observed, predicted = predicted ), mode = "everything", positive = levels(train_data[, y])[2]) if(nlevels(dato3[,y]) > 2){ f1_temp <- mean(conf_mat$byClass[, "F1"], na.rm = T) conf_mat_distr <- paste(levels(dato3[, y]), conf_mat$byClass[, "F1"], sep = " = ", collapse = ", ") empty_class <- (conf_mat$byClass[, "Recall"] %in% c(0, NA, NaN) ) data_matrix[conto, "empty_class"] <- paste(levels(dato3[, y]), empty_class, sep = " = ", collapse = ", ") }else{ f1_temp <- conf_mat$byClass["F1"] conf_mat_distr <- paste(c("TN", "FN", "FP", "TP"), conf_mat$table, sep = " = ", collapse = ", ") empty_class <- ifelse(0 %in% diag(conf_mat$table), 1, 0) data_matrix[conto, "empty_class"] <- empty_class } data_matrix[conto, "f1"] <- f1_temp #paste(levels(dato3[, y]), f1_temp, sep = " = ", collapse = ", ") } } } } t_svm <- difftime(time1 = Sys.time(), time2 = t0, units = 'mins') data_matrix_list[[4]] <- data_matrix ### Multilayer perceptron ### t0 <- Sys.time()
cat("**** START MPL ****") mod_params <- c("hidden1", "hidden2", "learning_rate") data_matrix <- matrix(NA, ncol = length(mod_params) + 4); colnames(data_matrix) <- c(mod_params, c("n_rep", "k_cv", "f1", "empty_class") ) l_rate_grid <- c(0.001, 0.005, seq(10, 200, 20) / 1000 ) hidden_grid <- 1:5 conto <- 0 for(n in 1:n_rep){ for(k in 1:k_cv){ set.seed(1234 + k_cv * (n - 1) + k); train_data_2 <- createDataPartition(y = train_data$y_and_registration, p = 1 - 1/k_cv, list = F) test_data_2 <- train_data[- train_data_2, - ncol(train_data) ] #righe <- train_data[train_data_2, ] train_data_2 <- train_data[train_data_2, ] train_data_2_temp_1 <- NULL train_data_2_temp_2 <- NULL for( w in names(c_weights_prop) ){ righe <- row.names(train_data_2[train_data_2$y_and_registration %in% w, ]) set.seed(1234 + k_cv * (n - 1) + k); righe_sampled <- sample(x = righe, size = length(righe) * c_weights_prop[[w]], replace = T) train_data_2_temp_1 <- train_data_2[righe_sampled, ] dim(train_data_2_temp_1) train_data_2_temp_2 <- rbind(train_data_2[!(train_data_2$y_and_registration %in% w), ], train_data_2_temp_1, train_data_2_temp_2) dim(train_data_2_temp_2) } train_data_2 <- train_data_2_temp_2[, - ncol(train_data)] test_data_2$n_televisit <- as.numeric(test_data_2$n_televisit) - 1; test_data_2 <- apply(test_data_2, MARGIN = 2, as.numeric) %>% as.data.frame