<jupyter_start><jupyter_text>Regresión logística<jupyter_code>library(tidyverse)
library(broom)
library(modelr)
library(purrr)
library(caret)
library(glmnet)
library("pROC")
set.seed(55)
path_proc =  "..\\data\\processed\\"
path_results = "..\\data\\results\\"<jupyter_output><empty_output><jupyter_text>Datasets<jupyter_code># Datos
dev = read_csv(paste0(path_proc,"data_bal_dummy_dev.csv"))
test = read_csv(paste0(path_proc,"data_bal_dummy_test.csv"))

train1 = read_csv(paste0(path_proc,"data_bal_dummy_train1.csv"))
train2 = read_csv(paste0(path_proc,"data_bal_dummy_train2.csv"))
train3 = read_csv(paste0(path_proc,"data_bal_dummy_train3.csv"))
train4 = read_csv(paste0(path_proc,"data_bal_dummy_train4.csv"))
train5 = read_csv(paste0(path_proc,"data_bal_dummy_train5.csv"))

val1 = read_csv(paste0(path_proc,"data_bal_dummy_val1.csv"))
val2 = read_csv(paste0(path_proc,"data_bal_dummy_val2.csv"))
val3 = read_csv(paste0(path_proc,"data_bal_dummy_val3.csv"))
val4 = read_csv(paste0(path_proc,"data_bal_dummy_val4.csv"))
val5 = read_csv(paste0(path_proc,"data_bal_dummy_val5.csv"))


# Se quita BMI y id
dev$BMI = NULL
test$BMI = NULL
train1$BMI = NULL
train2$BMI = NULL
train3$BMI = NULL
train4$BMI = NULL
train5$BMI = NULL
val1$BMI = NULL
val2$BMI = NULL
val3$BMI = NULL
val4$BMI = NULL
val5$BMI = NULL
# Vector con el target
dev_target = dev$target
# Matriz con los regresores
dev_mtx = model.matrix(target~., data = dev[,c(2:57)])<jupyter_output><empty_output><jupyter_text>Búsqueda de lambda para regularización Lasso<jupyter_code># Búsquda de lambda con CV
lasso_cv=cv.glmnet(x=dev_mtx, y=dev_target, alpha=1, standardize = T,  family = "binomial", type.measure="auc")
# Plot
plot(lasso_cv)
# Selección lambda óptimo
lasso_lambda_opt = lasso_cv$lambda.min

# Entrenamiento modelo óptimo
lasso_opt = glmnet(x=dev_mtx, # Matriz de regresores
                 y=dev_target, #Vector de la variable a predecir
                 family = "binomial",
                 alpha=1, # Indicador del tipo de regularizacion
                 standardize = TRUE,  # Estandarizamos
                 lambda = lasso_lambda_opt)

# Salida estandar
lasso_opt
# Variables que permanencen
tmp_coeffs = coef(lasso_opt)
data.frame(name = tmp_coeffs@Dimnames[[1]][tmp_coeffs@i + 1], coefficient = tmp_coeffs@x)
# Selección de variables
dev = dev %>% select(AGE_P,CIGSDAY,MOD,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
test = test %>% select(AGE_P,CIGSDAY,MOD,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)

train1 = train1 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
train2 = train2 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
train3 = train3 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
train4 = train4 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
train5 = train5 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)

val1 = val1 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
val2 = val2 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
val3 = val3 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
val4 = val4 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)
val5 = val5 %>% select(AGE_P,CIGSDAY,ALC,ASISLEEP,HYPMED2,HRTEV,STREV,AASSTILL,ULCYR,AHAYFYR,SINYR,PAINECK,PAINLB,PAINFACE,FLA1AR,APOX,SEX_2,RACERPI2_4,RACERPI2_6,R_MARITL_2,target)

# Vector con el target
train1_target = train1$target
train2_target = train2$target
train3_target = train3$target
train4_target = train4$target
train5_target = train5$target

val1_target = val1$target
val2_target = val2$target
val3_target = val3$target
val4_target = val4$target
val5_target = val5$target

test_target = test$target
# Fórmula de los modelos
logit_formula = formulas(.response = ~target,
                         model_todos = ~ .,
                         model_sujeto = ~ AGE_P + SEX_2 + RACERPI2_4 + RACERPI2_6 + R_MARITL_2,
                         model_habitos = ~ CIGSDAY + ALC + ASISLEEP,
                         model_dolor = ~ PAINECK + PAINLB + PAINFACE + FLA1AR,
                         model_med = ~ HYPMED2,
                         model_enf_ant = ~ HRTEV + STREV + AASSTILL + ULCYR + AHAYFYR + SINYR + APOX,
                         model_sujeto_habitos = ~ AGE_P + SEX_2 + RACERPI2_4 + RACERPI2_6 + R_MARITL_2 + CIGSDAY + ALC + ASISLEEP,
                         model_sujeto_dolor = ~ AGE_P + SEX_2 + RACERPI2_4 + RACERPI2_6 + R_MARITL_2 + PAINECK + PAINLB + PAINFACE + FLA1AR,
                         model_sujeto_med = ~ AGE_P + SEX_2 + RACERPI2_4 + RACERPI2_6 + R_MARITL_2 + HYPMED2,
                         model_sujeto_enf_ant = ~ AGE_P + SEX_2 + RACERPI2_4 + RACERPI2_6 + R_MARITL_2 + HRTEV + STREV + AASSTILL + ULCYR + AHAYFYR + SINYR + APOX   
                          )
# Generación de modelos
modelos_train = tibble(logit_formula) %>%
          mutate(model = names(logit_formula),
          expression = paste(logit_formula),
          mod = map(logit_formula, ~glm(.,family = 'binomial', data = dev))) 
# Las medidas de evaluación para cada modelo
modelos_train = modelos_train %>% mutate(glance = map(mod, broom::glance))
# Se filtran medidas de evaluacion de interes
modelos_train %>% unnest(glance, .drop = TRUE) %>%
            mutate(perc_explained_dev = 1-deviance/null.deviance) %>% 
            select(expression, null.deviance, deviance, perc_explained_dev) %>% 
            arrange(deviance)
# Regresión logística del mejor modelo para cada train
model_train1 = glm(target ~ ., family = 'binomial', data = train1)
model_train2 = glm(target ~ ., family = 'binomial', data = train2)
model_train3 = glm(target ~ ., family = 'binomial', data = train3)
model_train4 = glm(target ~ ., family = 'binomial', data = train4)
model_train5 = glm(target ~ ., family = 'binomial', data = train5)
# Predicciones sobre el dataset de validacion
y_pred1 = predict.glm(model_train1, newdata = val1, type="response")
# Dataframe que aloja las predicciones
val_values1 = data.frame(y_pred1, as.factor(val1_target))
colnames(val_values1) = c("y_pred","y_real")
# Vector de puntos de corte desde el mínimo al máximo de nuestras probabilidades con intervalo de 0.001
cutoffs1 = seq(
  min(y_pred1),
  max(y_pred1),
  0.001
)

# Predicciones sobre el dataset de validacion
y_pred2 = predict.glm(model_train2, newdata = val2, type="response")
# Dataframe que aloja las predicciones
val_values2 = data.frame(y_pred2, as.factor(val2_target))
colnames(val_values2) = c("y_pred","y_real")
# Vector de puntos de corte desde el mínimo al máximo de nuestras probabilidades con intervalo de 0.001
cutoffs2 = seq(
  min(y_pred2),
  max(y_pred2),
  0.001
)

# Predicciones sobre el dataset de validacion
y_pred3 = predict.glm(model_train3, newdata = val3, type="response")
# Dataframe que aloja las predicciones
val_values3 = data.frame(y_pred1, as.factor(val3_target))
colnames(val_values3) = c("y_pred","y_real")
# Vector de puntos de corte desde el mínimo al máximo de nuestras probabilidades con intervalo de 0.001
cutoffs3 = seq(
  min(y_pred3),
  max(y_pred3),
  0.001
)

# Predicciones sobre el dataset de validacion
y_pred4 = predict.glm(model_train4, newdata = val4, type="response")
# Dataframe que aloja las predicciones
val_values4 = data.frame(y_pred4, as.factor(val4_target))
colnames(val_values4) = c("y_pred","y_real")
# Vector de puntos de corte desde el mínimo al máximo de nuestras probabilidades con intervalo de 0.001
cutoffs4 = seq(
  min(y_pred4),
  max(y_pred4),
  0.001
)

# Predicciones sobre el dataset de validacion
y_pred5 = predict.glm(model_train5, newdata = val5, type="response")
# Dataframe que aloja las predicciones
val_values5 = data.frame(y_pred5, as.factor(val5_target))
colnames(val_values5) = c("y_pred","y_real")
# Vector de puntos de corte desde el mínimo al máximo de nuestras probabilidades con intervalo de 0.001
cutoffs5 = seq(
  min(y_pred5),
  max(y_pred5),
  0.001
)
# Función que clasifica las observaciones para diferentes puntos de corte y devuelve métricas de performance asociadas
prediction_metrics = function(cutoff, predictions = val_values1) {
      table_res = predictions %>%
      mutate(predicted_class = if_else(y_pred > cutoff, 1, 0) %>% as.factor(),
              real_class = factor(y_real)
              )
    
    confusionMatrix(table_res$predicted_class, table_res$real_class, positive = "1") %>%
      tidy() %>%
      select(term, estimate) %>%
      filter(term %in% c('accuracy', 'f1', 'specificity', 'precision','sensitivity')) %>%
      mutate(cutoff = cutoff)
  }

# Se aplica la función al vector
logit_pred1 = map_dfr(cutoffs1, prediction_metrics) %>% mutate(term = as.factor(term))
# Función que clasifica las observaciones para diferentes puntos de corte y devuelve métricas de performance asociadas
prediction_metrics = function(cutoff, predictions = val_values2) {
      table_res = predictions %>%
      mutate(predicted_class = if_else(y_pred > cutoff, 1, 0) %>% as.factor(),
              real_class = factor(y_real)
              )
    
    confusionMatrix(table_res$predicted_class, table_res$real_class, positive = "1") %>%
      tidy() %>%
      select(term, estimate) %>%
      filter(term %in% c('accuracy', 'f1', 'specificity', 'precision','sensitivity')) %>%
      mutate(cutoff = cutoff)
  }

# Se aplica la función al vector
logit_pred2 = map_dfr(cutoffs2, prediction_metrics) %>% mutate(term = as.factor(term))
# Función que clasifica las observaciones para diferentes puntos de corte y devuelve métricas de performance asociadas
prediction_metrics = function(cutoff, predictions = val_values3) {
      table_res = predictions %>%
      mutate(predicted_class = if_else(y_pred > cutoff, 1, 0) %>% as.factor(),
              real_class = factor(y_real)
              )
    
    confusionMatrix(table_res$predicted_class, table_res$real_class, positive = "1") %>%
      tidy() %>%
      select(term, estimate) %>%
      filter(term %in% c('accuracy', 'f1', 'specificity', 'precision','sensitivity')) %>%
      mutate(cutoff = cutoff)
  }

# Se aplica la función al vector
logit_pred3 = map_dfr(cutoffs3, prediction_metrics) %>% mutate(term = as.factor(term))
# Función que clasifica las observaciones para diferentes puntos de corte y devuelve métricas de performance asociadas
prediction_metrics = function(cutoff, predictions = val_values4) {
      table_res = predictions %>%
      mutate(predicted_class = if_else(y_pred > cutoff, 1, 0) %>% as.factor(),
              real_class = factor(y_real)
              )
    
    confusionMatrix(table_res$predicted_class, table_res$real_class, positive = "1") %>%
      tidy() %>%
      select(term, estimate) %>%
      filter(term %in% c('accuracy', 'f1', 'specificity', 'precision','sensitivity')) %>%
      mutate(cutoff = cutoff)
  }

# Se aplica la función al vector
logit_pred4 = map_dfr(cutoffs4, prediction_metrics) %>% mutate(term = as.factor(term))
# Función que clasifica las observaciones para diferentes puntos de corte y devuelve métricas de performance asociadas
prediction_metrics = function(cutoff, predictions = val_values5) {
      table_res = predictions %>%
      mutate(predicted_class = if_else(y_pred > cutoff, 1, 0) %>% as.factor(),
              real_class = factor(y_real)
              )
    
    confusionMatrix(table_res$predicted_class, table_res$real_class, positive = "1") %>%
      tidy() %>%
      select(term, estimate) %>%
      filter(term %in% c('accuracy', 'f1', 'specificity', 'precision','sensitivity')) %>%
      mutate(cutoff = cutoff)
  }

# Se aplica la función al vector
logit_pred5 = map_dfr(cutoffs5, prediction_metrics) %>% mutate(term = as.factor(term))
logit_pred1 %>% 
  tidyr::spread(., term, estimate) %>% # se pasan las métricas a columnas
  arrange(-accuracy) %>% # se ordena por accuracy máximo a mínimo
  filter(abs(sensitivity - specificity)==min(abs(sensitivity - specificity)) ) %>% # se busca la mínima diferencia entre sensibilidad y especificidad
  select(-c("precision","f1"))
logit_pred2 %>% 
  tidyr::spread(., term, estimate) %>% # se pasan las métricas a columnas
  arrange(-accuracy) %>% # se ordena por accuracy máximo a mínimo
  filter(abs(sensitivity - specificity)==min(abs(sensitivity - specificity)) ) %>% # se busca la mínima diferencia entre sensibilidad y especificidad
  select(-c("precision","f1"))
logit_pred3 %>% 
  tidyr::spread(., term, estimate) %>% # se pasan las métricas a columnas
  arrange(-accuracy) %>% # se ordena por accuracy máximo a mínimo
  filter(abs(sensitivity - specificity)==min(abs(sensitivity - specificity)) ) %>% # se busca la mínima diferencia entre sensibilidad y especificidad
  select(-c("precision","f1"))
logit_pred4 %>% 
  tidyr::spread(., term, estimate) %>% # se pasan las métricas a columnas
  arrange(-accuracy) %>% # se ordena por accuracy máximo a mínimo
  filter(abs(sensitivity - specificity)==min(abs(sensitivity - specificity)) ) %>% # se busca la mínima diferencia entre sensibilidad y especificidad
  select(-c("precision","f1"))
logit_pred5 %>% 
  tidyr::spread(., term, estimate) %>% # se pasan las métricas a columnas
  arrange(-accuracy) %>% # se ordena por accuracy máximo a mínimo
  filter(abs(sensitivity - specificity)==min(abs(sensitivity - specificity)) ) %>% # se busca la mínima diferencia entre sensibilidad y especificidad
  select(-c("precision","f1"))
# Cutoff
sel_cutoff1 = 0.4785447
# Clasifico utilizando el punto de corte
val_values1$predicted_class = if_else(val_values1$y_pred >= sel_cutoff1, 1, 0) %>% as.factor()
# Creo la matriz de confusión
confusionMatrix(val_values1$predicted_class, val_values1$y_real, positive = "1")
# ROC-AUC
roc_obj = roc(as.numeric(val_values1$y_real), as.numeric(val_values1$predicted_class))
auc(roc_obj)
# Cutoff
sel_cutoff2 = 0.5229288
# Clasifico utilizando el punto de corte
val_values2$predicted_class = if_else(val_values2$y_pred >= sel_cutoff2, 1, 0) %>% as.factor()
# Creo la matriz de confusión
confusionMatrix(val_values2$predicted_class, val_values2$y_real, positive = "1")
# ROC-AUC
roc_obj = roc(as.numeric(val_values2$y_real), as.numeric(val_values2$predicted_class))
auc(roc_obj)
# Cutoff
sel_cutoff3 = 0.4701566 
# Clasifico utilizando el punto de corte
val_values3$predicted_class = if_else(val_values3$y_pred >= sel_cutoff3, 1, 0) %>% as.factor()
# Creo la matriz de confusión
confusionMatrix(val_values3$predicted_class, val_values3$y_real, positive = "1")
# ROC-AUC
roc_obj = roc(as.numeric(val_values3$y_real), as.numeric(val_values3$predicted_class))
auc(roc_obj)
# Cutoff
sel_cutoff4 = 0.5030193 
# Clasifico utilizando el punto de corte
val_values4$predicted_class = if_else(val_values4$y_pred >= sel_cutoff4, 1, 0) %>% as.factor()
# Creo la matriz de confusión
confusionMatrix(val_values4$predicted_class, val_values4$y_real, positive = "1")
# ROC-AUC
roc_obj = roc(as.numeric(val_values4$y_real), as.numeric(val_values4$predicted_class))
auc(roc_obj)
# Cutoff
sel_cutoff5 = 0.4858599 
# Clasifico utilizando el punto de corte
val_values5$predicted_class = if_else(val_values5$y_pred >= sel_cutoff5, 1, 0) %>% as.factor()
# Creo la matriz de confusión
confusionMatrix(val_values5$predicted_class, val_values5$y_real, positive = "1")
# ROC-AUC
roc_obj = roc(as.numeric(val_values5$y_real), as.numeric(val_values5$predicted_class))
auc(roc_obj)
# Modelo entrenado con todo el dataset dev
model_dev = glm(target ~ ., family = 'binomial', data = dev)
# Se extraen los coeficientes
tidy(model_dev)
# Predicciones sobre el dataset de test
y_pred = predict.glm(model_dev, newdata = test, type="response")

# Dataframe que aloja las predicciones
test_values = data.frame(y_pred, as.factor(test_target))
colnames(test_values) = c("y_pred","y_real")

# Clasifico utilizando el punto de corte
test_values$predicted_class = if_else(test_values$y_pred >= sel_cutoff4, 1, 0) %>% as.factor()

# Creo la matriz de confusión
confusionMatrix(test_values$predicted_class, test_values$y_real, positive = "1")

# ROC-AUC
roc_obj = roc(as.numeric(test_values$y_real), as.numeric(test_values$predicted_class))
auc(roc_obj)
ggplot(test_values, aes(x = y_real, y = y_pred, group = y_real, fill = factor(y_real), color = factor(y_real))) + 
        theme_grey() + theme(legend.position = "none") +
        geom_violin() + 
        labs(title = 'Violin plot para el modelo', y = 'Probabilidad predicha', x = "Clase")
write_csv(test_values, paste0(path_results,"test_bal_dummy_pred.csv"))
write_csv(val_values1, paste0(path_results,"val_bal_dummy_pred1.csv"))
write_csv(val_values2, paste0(path_results,"val_bal_dummy_pred2.csv"))
write_csv(val_values3, paste0(path_results,"val_bal_dummy_pred3.csv"))
write_csv(val_values4, paste0(path_results,"val_bal_dummy_pred4.csv"))
write_csv(val_values5, paste0(path_results,"val_bal_dummy_pred5.csv"))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12058397414465619, 'fraction_numerical': 0.02685835283628664, 'mean_word_length': 5.308260105448155, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 103, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '545896', 'license': 'MIT', 'n_tokens_mistral': 7234, 'n_tokens_neox': 6927, 'n_words': 1265, 'path': 'Trabajo-especializacion/notebooks/04_RL_parametros-lasso_R-glmnet.ipynb', 'repo_name': 'ponybiam/DMKD-UBA', 'size': '17946'}