Modelo de regresión logística.

Alcance y objetivo del modelo.

El presente ejercicio se plantea ajustar un modelo que tenga el mejor performance con la muestra completa, con la finalidad de ser evaluado ante la muestra de Test correspondiente al examen final del curso.

Ajuste del modelo
  • Visualizaciones del modelo de regresion linea simple

Gráfica no. 1 (Logit simple glm(Loan.Status ~. ))

¿Qué se puede concluir del anterior GRID?

Ajuste del modelo no. 1

En este caso nos referimos al modelo no.1 solo interamente para este rmarkdown.

modelo_2 = glm(data = Train_intern_1, 
              formula = Loan.Status~., family = "binomial")
summary(modelo_2)
## 
## Call:
## glm(formula = Loan.Status ~ ., family = "binomial", data = Train_intern_1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.7204   0.4615   0.6004   0.7036   1.2926  
## 
## Coefficients:
##                                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   -4.104e+00  1.040e+00  -3.948 7.89e-05 ***
## TermShort Term                 1.235e-01  6.919e-02   1.784 0.074373 .  
## Years.in.current.job1 year     2.187e-01  1.472e-01   1.486 0.137395    
## Years.in.current.job10+ years  5.944e-02  1.079e-01   0.551 0.581821    
## Years.in.current.job2 years   -7.238e-02  1.267e-01  -0.571 0.567855    
## Years.in.current.job3 years    5.016e-02  1.347e-01   0.372 0.709541    
## Years.in.current.job4 years    1.631e-01  1.510e-01   1.080 0.280111    
## Years.in.current.job5 years    1.533e-01  1.463e-01   1.047 0.294951    
## Years.in.current.job6 years   -8.367e-02  1.502e-01  -0.557 0.577417    
## Years.in.current.job7 years   -1.018e-01  1.471e-01  -0.692 0.489201    
## Years.in.current.job8 years   -1.833e-01  1.536e-01  -1.193 0.232975    
## Years.in.current.job9 years    1.999e-01  1.682e-01   1.188 0.234703    
## Years.in.current.jobn/a        3.956e-01  1.914e-01   2.067 0.038761 *  
## Home.OwnershipHome Mortgage    8.726e-01  7.341e-01   1.189 0.234583    
## Home.OwnershipOwn Home         8.793e-01  7.375e-01   1.192 0.233133    
## Home.OwnershipRent             5.619e-01  7.321e-01   0.768 0.442762    
## Purposebuy a car               1.176e+00  3.980e-01   2.954 0.003137 ** 
## Purposebuy house               4.665e-02  3.121e-01   0.149 0.881175    
## Purposedebt consolidation      3.590e-01  1.885e-01   1.905 0.056757 .  
## Purposeeducational expenses   -9.569e-01  7.061e-01  -1.355 0.175344    
## Purposehome improvements       1.564e-01  2.117e-01   0.739 0.460106    
## Purposemajor_purchase         -3.545e-01  4.314e-01  -0.822 0.411275    
## Purposemedical bills           3.592e-01  2.743e-01   1.309 0.190434    
## Purposemoving                 -5.668e-01  7.242e-01  -0.783 0.433810    
## Purposeother                   4.502e-01  2.004e-01   2.246 0.024673 *  
## Purposerenewable_energy        1.094e+01  1.875e+02   0.058 0.953449    
## Purposesmall_business         -8.151e-01  4.677e-01  -1.743 0.081408 .  
## Purposetake a trip             2.100e-01  4.292e-01   0.489 0.624677    
## Purposevacation                5.294e-01  8.183e-01   0.647 0.517671    
## Purposewedding                -9.450e-01  8.407e-01  -1.124 0.260990    
## Current.Loan.Amount           -8.448e-07  1.977e-07  -4.272 1.93e-05 ***
## Credit.Score                   5.441e-03  1.054e-03   5.163 2.43e-07 ***
## Annual.Income                  4.766e-07  5.842e-08   8.158 3.40e-16 ***
## Monthly.Debt                  -1.371e-05  3.586e-06  -3.825 0.000131 ***
## Years.of.Credit.History        3.558e-03  4.436e-03   0.802 0.422444    
## Months.since.last.delinquent   1.160e-02  3.386e-03   3.424 0.000616 ***
## Number.of.Open.Accounts       -8.350e-03  6.297e-03  -1.326 0.184861    
## Number.of.Credit.Problems     -5.232e-02  5.226e-02  -1.001 0.316799    
## Current.Credit.Balance         1.313e-07  2.426e-07   0.541 0.588415    
## Maximum.Open.Credit            1.651e-07  1.050e-07   1.573 0.115680    
## Bankruptcies                   8.830e-02  5.044e-02   1.751 0.080024 .  
## Tax.Liens                      3.933e-02  4.955e-02   0.794 0.427360    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9260.7  on 9275  degrees of freedom
## Residual deviance: 8967.0  on 9234  degrees of freedom
## AIC: 9051
## 
## Number of Fisher Scoring iterations: 11

Salida de consola y código no.1 (Ajuste del modelo no.2 )

Update del modelo con las variables estadísticamente significativas
modelo_2 = update(object = modelo_2, formula. = .~. -Term -Home.Ownership -Years.in.current.job -Number.of.Open.Accounts -Number.of.Credit.Problems -Current.Credit.Balance -Tax.Liens -Maximum.Open.Credit -Bankruptcies)
summary(modelo_2)
## 
## Call:
## glm(formula = Loan.Status ~ Purpose + Current.Loan.Amount + Credit.Score + 
##     Annual.Income + Monthly.Debt + Years.of.Credit.History + 
##     Months.since.last.delinquent, family = "binomial", data = Train_intern_1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4011   0.4828   0.6174   0.6967   1.3146  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -4.169e+00  6.632e-01  -6.285 3.27e-10 ***
## Purposebuy a car              1.131e+00  3.955e-01   2.859 0.004246 ** 
## Purposebuy house             -2.408e-02  3.107e-01  -0.078 0.938219    
## Purposedebt consolidation     3.301e-01  1.868e-01   1.767 0.077237 .  
## Purposeeducational expenses  -9.619e-01  7.031e-01  -1.368 0.171286    
## Purposehome improvements      2.280e-01  2.093e-01   1.089 0.276023    
## Purposemajor_purchase        -3.288e-01  4.284e-01  -0.767 0.442791    
## Purposemedical bills          2.964e-01  2.729e-01   1.086 0.277442    
## Purposemoving                -6.217e-01  7.178e-01  -0.866 0.386484    
## Purposeother                  4.232e-01  1.992e-01   2.124 0.033679 *  
## Purposerenewable_energy       1.124e+01  1.875e+02   0.060 0.952199    
## Purposesmall_business        -9.071e-01  4.672e-01  -1.942 0.052190 .  
## Purposetake a trip           -7.602e-02  3.727e-01  -0.204 0.838361    
## Purposevacation               4.775e-01  8.182e-01   0.584 0.559456    
## Purposewedding               -1.060e+00  8.402e-01  -1.262 0.206998    
## Current.Loan.Amount          -7.472e-07  1.785e-07  -4.186 2.84e-05 ***
## Credit.Score                  6.624e-03  9.186e-04   7.211 5.57e-13 ***
## Annual.Income                 4.846e-07  5.763e-08   8.409  < 2e-16 ***
## Monthly.Debt                 -1.257e-05  3.227e-06  -3.894 9.85e-05 ***
## Years.of.Credit.History       1.008e-02  4.252e-03   2.370 0.017803 *  
## Months.since.last.delinquent  1.108e-02  3.341e-03   3.315 0.000915 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9260.7  on 9275  degrees of freedom
## Residual deviance: 9035.0  on 9255  degrees of freedom
## AIC: 9077
## 
## Number of Fisher Scoring iterations: 11

Salida de consola y código no.2 (Update al modelo no.2)

El modelo 1 tiene un AIC de 9’077, la version que contempla todos los predictores era de 9051

Anova del modelo
anova(modelo_2)

Salida de consola y código no.3 (ANOVA Modelo no. 2)

Selección de variables en modelo de regresion logistica por críterio de información AIC
glmulti.logistic.out <-
    glmulti(Loan.Status ~ Purpose + Current.Loan.Amount + Credit.Score + 
    Annual.Income + Monthly.Debt + Years.of.Credit.History + 
    Months.since.last.delinquent, data = Train_intern_1,
            level = 1,               # No interaction considered
            method = "h",            # Exhaustive approach
            crit = "aic",            # AIC as criteria
            confsetsize = 5,         # Keep 5 best models
            plotty = F, report = F,  # No plot or interim reports
            fitfunction = "glm",     # glm function
            family = binomial)       # binomial family for logistic regression

## Show 5 best models (Use @ instead of $ for an S4 object)
glmulti.logistic.out@formulas
## [[1]]
## Loan.Status ~ 1 + Purpose + Current.Loan.Amount + Credit.Score + 
##     Annual.Income + Monthly.Debt + Years.of.Credit.History + 
##     Months.since.last.delinquent
## <environment: 0x55db66a87b50>
## 
## [[2]]
## Loan.Status ~ 1 + Purpose + Current.Loan.Amount + Credit.Score + 
##     Annual.Income + Monthly.Debt + Months.since.last.delinquent
## <environment: 0x55db66a87b50>
## 
## [[3]]
## Loan.Status ~ 1 + Current.Loan.Amount + Credit.Score + Annual.Income + 
##     Monthly.Debt + Years.of.Credit.History + Months.since.last.delinquent
## <environment: 0x55db66a87b50>
## 
## [[4]]
## Loan.Status ~ 1 + Purpose + Current.Loan.Amount + Credit.Score + 
##     Annual.Income + Monthly.Debt + Years.of.Credit.History
## <environment: 0x55db66a87b50>
## 
## [[5]]
## Loan.Status ~ 1 + Current.Loan.Amount + Credit.Score + Annual.Income + 
##     Monthly.Debt + Months.since.last.delinquent
## <environment: 0x55db66a87b50>

Salida de consola y código no.4 (Selección de variables)-

De entrada, el mejor modelo por AIC que propone la funcion glmulti() es el mismo modelo de regresión logística que se le pasó por parametro, por lo que podemos hablar de que ese al momento es el modelo ideal, con 8 features. Sin embargo en segundo y tercer lugar proponé modelos con 7 features, por lo que se evaluará contra validación cruzada la viabilidad de irnos por un modelo más ligero.

  • Elaborando modelos alternativos a modelo_1

  • Modelo 2.1

modelo_2_1 = glm(data = Train_intern_1, 
              formula = Loan.Status~ Purpose + Current.Loan.Amount + Credit.Score + 
    Annual.Income + Monthly.Debt + Months.since.last.delinquent, family = "binomial")
summary(modelo_2_1)
## 
## Call:
## glm(formula = Loan.Status ~ Purpose + Current.Loan.Amount + Credit.Score + 
##     Annual.Income + Monthly.Debt + Months.since.last.delinquent, 
##     family = "binomial", data = Train_intern_1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4055   0.4859   0.6198   0.6943   1.3567  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -4.121e+00  6.624e-01  -6.221 4.94e-10 ***
## Purposebuy a car              1.137e+00  3.954e-01   2.876  0.00402 ** 
## Purposebuy house             -2.017e-02  3.106e-01  -0.065  0.94822    
## Purposedebt consolidation     3.410e-01  1.867e-01   1.827  0.06771 .  
## Purposeeducational expenses  -9.534e-01  7.026e-01  -1.357  0.17483    
## Purposehome improvements      2.425e-01  2.091e-01   1.160  0.24624    
## Purposemajor_purchase        -3.416e-01  4.285e-01  -0.797  0.42530    
## Purposemedical bills          3.170e-01  2.726e-01   1.163  0.24485    
## Purposemoving                -6.047e-01  7.182e-01  -0.842  0.39979    
## Purposeother                  4.291e-01  1.991e-01   2.155  0.03117 *  
## Purposerenewable_energy       1.123e+01  1.875e+02   0.060  0.95226    
## Purposesmall_business        -9.177e-01  4.670e-01  -1.965  0.04943 *  
## Purposetake a trip           -8.788e-02  3.728e-01  -0.236  0.81362    
## Purposevacation               4.339e-01  8.177e-01   0.531  0.59566    
## Purposewedding               -1.089e+00  8.404e-01  -1.296  0.19498    
## Current.Loan.Amount          -7.167e-07  1.780e-07  -4.027 5.65e-05 ***
## Credit.Score                  6.779e-03  9.156e-04   7.404 1.32e-13 ***
## Annual.Income                 4.941e-07  5.752e-08   8.591  < 2e-16 ***
## Monthly.Debt                 -1.234e-05  3.225e-06  -3.826  0.00013 ***
## Months.since.last.delinquent  1.075e-02  3.337e-03   3.221  0.00128 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9260.7  on 9275  degrees of freedom
## Residual deviance: 9040.6  on 9256  degrees of freedom
## AIC: 9080.6
## 
## Number of Fisher Scoring iterations: 11

El modelo 2.1 tiene un AIC de 9080.6, mientras que el modelo 2 tiene el AIC de 9077. Salida de consola y código no.5 (Ajustando modelo no. 2.1)-

  • Modelo 2.2
#Modelo 2.2
modelo_2_2 = glm(data = Train_intern_1, 
              formula = Loan.Status~ Current.Loan.Amount + Credit.Score + Annual.Income + 
    Monthly.Debt + Years.of.Credit.History + Months.since.last.delinquent, family = "binomial")
summary(modelo_2_2)
## 
## Call:
## glm(formula = Loan.Status ~ Current.Loan.Amount + Credit.Score + 
##     Annual.Income + Monthly.Debt + Years.of.Credit.History + 
##     Months.since.last.delinquent, family = "binomial", data = Train_intern_1)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3537   0.4982   0.6211   0.6967   1.1448  
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  -4.002e+00  6.394e-01  -6.260 3.85e-10 ***
## Current.Loan.Amount          -7.596e-07  1.717e-07  -4.423 9.72e-06 ***
## Credit.Score                  6.819e-03  8.902e-04   7.661 1.85e-14 ***
## Annual.Income                 4.748e-07  5.661e-08   8.387  < 2e-16 ***
## Monthly.Debt                 -1.163e-05  3.188e-06  -3.646 0.000266 ***
## Years.of.Credit.History       1.057e-02  4.235e-03   2.497 0.012520 *  
## Months.since.last.delinquent  1.125e-02  3.330e-03   3.379 0.000728 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 9260.7  on 9275  degrees of freedom
## Residual deviance: 9068.9  on 9269  degrees of freedom
## AIC: 9082.9
## 
## Number of Fisher Scoring iterations: 4

El modelo 2.2 tiene un AIC de 9’080.6, mientras que el modelo 2.1 tiene el AIC de 9’077. Salida de consola y código no.7 (Ajustando modelo no. 2.2)-

Validación cruzada y MSPR

De los 5 modelos previamente propuestos por selección de variables con glmulti en regresion logistica y con estadístico AIC, se evaluarán los primeros 3, cada uno por validación cruzada y por MSPR para solo quedarnos con el mejor.

  • Validación por k-fold para el modelo 2.1
MSE_kf_2_1 <- cv.glm(data = Test_intern_1,
               glmfit = modelo_2_1, 
               K = 50)
MSE_kf_2_1$delta
## [1] 0.1665053 0.1587315

Salida de consola y código no.9 (MSPR por K-Pliegues del modelo no. 2.1)-

Interpretando $delta:

El primer valor representa el error medio cuadrático promedio para los K pliegues (MSPR o Mean Squared Prediction Error), y el segundo valor representa la desviación estándar de los errores medios cuadráticos de los K pliegues.

  • Validación por k-fold para el modelo 2.2
MSE_kf_2_2 <- cv.glm(data = Test_intern_1,
               glmfit = modelo_2_2, 
               K = 50)
MSE_kf_2_2$delta
## [1] 0.1662972 0.1600346

Salida de consola y código no. 10 (MSPR por K-Pliegues del modelo no. 2.2)-

Concluimos que la alternativa 2.2 con solo 6 features presenta ligeramente un MSPR menor que el modelo base, mientras que el modelo 2.1 presenta la ventaja de tener una predictora menos, por lo que dichos modelos pasaran a la etapa de remuestreo y el estudio de la matriz de confusión con la finalidad de proponer el de mejor performance en estas 2 instancias (validación cruzada donde el modelo 2.1 fue el mejor y matriz de confusión).

Modelo no. 2.1
Matriz de confusión para el modelo base con 6 features
  • Respuesta del modelo en el conjunto de prueba
prediccion_1 = predict(object = modelo_2_1,
                    newdata = Test_intern_1, type = "response")
df_prediccion_1 <- data.frame(preds = prediccion_1)

dist_1 = ggplot(df_prediccion_1, aes(y="preds", x=prediccion_1))+
  geom_jitter(alpha=0.5, color="greenyellow")+
  geom_violin(alpha=0.5)+
  geom_boxplot(width = 0.3, fill = "white", alpha = 0.5)+
  scale_x_continuous(labels=c("40%", "60%", "80%", "100%"), name = "Respuesta")+
  scale_y_discrete(labels=c(""), name=NULL)+
  ggtitle("Predicciones para el conjunto Test en el modelo 2.1")+
  theme_bw()

dist_1

Gráfica no. 2 (Distribución de la predicción para el modelo no 2.1)-

  • punto de corte ideal para maximizar la sensibilidad (vs Test)
prediccion_1 = predict(object = modelo_2_1,
                    newdata = Test_intern_1, type = "response")

Sensibilidad_1 <- vector()
Corte_1 <- seq(0.05, 0.95, by = 0.001)

for (i in 1:length(Corte_1)) {
  Prediccion <- as.factor(ifelse(prediccion_1 >= Corte_1[i], yes = 1, no = 0))
  
  Sensibilidad_1[i] <- tryCatch(
    {
      # confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), Prediccion)$overall[1]
      confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), 
                             Prediccion)$byClass[1]
      
    },
    error = function(e) {
      0 # O cualquier valor predeterminado
    }
  )
}

# plot(x=Corte_1, y = Sensibilidad_1, type = "l", main ="Punto de corte óptimo para maximizar la sensibilidad (vs Test)", xlab="Corte", ylab="Exactitud")
df_sensibilidad = data.frame(corte= Corte_1, sensibilidad = Sensibilidad_1)
sens1 = ggplot(df_sensibilidad, aes(x=Corte_1, y=Sensibilidad_1))+
  geom_line()+
  geom_point(aes(x = Corte_1[which.max(Sensibilidad_1)], y = Sensibilidad_1[which.max(Sensibilidad_1)]), 
             shape = 21, size = 5, fill = "white", color = "blue") +
  geom_text(aes(x = Corte_1[which.max(Sensibilidad_1)], y = Sensibilidad_1[which.max(Sensibilidad_1)] - 0.07, label = glue("Max \n{Corte_1[which.max(Sensibilidad_1)]}")),
            hjust = 1) +
  scale_y_continuous(breaks=c(0.25, 0.5, 0.75, 1), labels = c("25%", "50%", "75%", "100%"),
                     name = "Sensibilidad")+
  scale_x_continuous("Punto de corte")+
  ggtitle("Punto de corte ideal para máximizar la sensibilidad (modelo 2.1)")+
  theme_bw()
sens1

Gráfica no. 3 (punto de corte ideal para maximizar la sensibilidad (vs Test))-

  • punto de corte ideal para maximizar la especificidad (vs Test)
prediccion_1 = predict(object = modelo_2_1,
                    newdata = Test_intern_1, type = "response")

Especificidad_1 <- vector()
Corte_1 <- seq(0.05, 0.95, by = 0.001)

for (i in 1:length(Corte_1)) {
  Prediccion <- as.factor(ifelse(prediccion_1 >= Corte_1[i], yes = 1, no = 0))
  
  Especificidad_1[i] <- tryCatch(
    {
      # confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), Prediccion)$overall[1]
      confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), 
                             Prediccion)$byClass[2]
    },
    error = function(e) {
      0 # O cualquier valor predeterminado
    }
  )
}

# plot(x=Corte_1, y = Especificidad_1, type = "l", main ="Punto de corte óptimo para maximizar la especificidad (vs Test)", xlab="Corte", ylab="Exactitud")
df_especificidad = data.frame(corte= Corte_1, especificidad = Especificidad_1)
espe1 = ggplot(df_especificidad, aes(x=Corte_1, y=Especificidad_1))+
  geom_line()+
  geom_point(aes(x = Corte_1[which.max(Especificidad_1)], y = Especificidad_1[which.max(Especificidad_1)]), 
             shape = 21, size = 5, fill = "white", color = "blue") +
  geom_text(aes(x = Corte_1[which.max(Especificidad_1)], y = Especificidad_1[which.max(Especificidad_1)] - 0.03, label = glue("Max \n{Corte_1[which.max(Especificidad_1)]}")),
            hjust = 1.5) +
  scale_y_continuous(breaks=c(0.25, 0.50, 0.75, 0.95, 1), labels = c("25%", "50%", "75%","95%", "100%"),
                     name = "Especificidad")+
  scale_x_continuous("Punto de corte")+
  ggtitle("Punto de corte ideal para máximizar la especificidad (modelo 2.1)")+
  theme_bw()
espe1

Gráfica no. 4 (punto de corte ideal para maximizar la especificidad (vs Test))-

  • punto de corte ideal para maximizar la exactitud (vs Test)
prediccion_1 = predict(object = modelo_2_1,
                    newdata = Test_intern_1, type = "response")

Exactitud_1 <- vector()
Corte_1 <- seq(0.05, 0.95, by = 0.001)

for (i in 1:length(Corte_1)) {
  Prediccion <- as.factor(ifelse(prediccion_1 >= Corte_1[i], yes = 1, no = 0))
  
  Exactitud_1[i] <- tryCatch(
    {
      confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), Prediccion)$overall[1]
    },
    error = function(e) {
      0 # O cualquier valor predeterminado
    }
  )
}

# plot(x=Corte_1, y = Exactitud_1, type = "l", main ="Punto de corte óptimo para maximizar el accuracy (vs Test)", xlab="Corte", ylab="Exactitud")
df_exactitud = data.frame(corte= Corte_1, exactitud = Exactitud_1)
exac1 = ggplot(df_especificidad, aes(x=Corte_1, y=Exactitud_1))+
  geom_line()+
  geom_point(aes(x = Corte_1[which.max(Exactitud_1)], y = Exactitud_1[which.max(Exactitud_1)]), 
             shape = 21, size = 5, fill = "white", color = "blue") +
  geom_text(aes(x = Corte_1[which.max(Exactitud_1)], y = Exactitud_1[which.max(Exactitud_1)] - 0.03, label = glue("Max \n{Corte_1[which.max(Exactitud_1)]}")),
            hjust = 1.3) +
  scale_y_continuous(breaks=c(0.2, 0.4, 0.5, 0.6, 0.7, 0.8), labels = c("20%", "40%", "50%", "60%", "70%", "80%"),
                     name = "Exactitud")+
  scale_x_continuous("Punto de corte")+
  ggtitle("Punto de corte ideal para máximizar la exactitud (modelo 1.1)")+
  theme_bw()
exac1

Gráfica no. 5 (punto de corte ideal para maximizar la exactitud (vs Test))-

corte_optimo_1 = which(Exactitud_1==max(Exactitud_1))
Corte_1[corte_optimo_1][1]
## [1] 0.383
Corte_modelo_2_1 = Corte_1[corte_optimo_1][1]

Salida de consola y código no. 11 (Corte ótimo para máximizar la exactitud)-

Prediccion_1<-as.factor(ifelse(prediccion_1>=Corte_1[corte_optimo_1][1],yes = 1, no = 0))


confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), 
                             Prediccion_1)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0    1  605
##          1    0 2486
##                                           
##                Accuracy : 0.8043          
##                  95% CI : (0.7899, 0.8182)
##     No Information Rate : 0.9997          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.0027          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 1.0000000       
##             Specificity : 0.8042705       
##          Pos Pred Value : 0.0016502       
##          Neg Pred Value : 1.0000000       
##              Prevalence : 0.0003234       
##          Detection Rate : 0.0003234       
##    Detection Prevalence : 0.1959897       
##       Balanced Accuracy : 0.9021352       
##                                           
##        'Positive' Class : 0               
## 
cfmx1 = confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), 
                             Prediccion_1)

Salida de consola y código no. 12 (Matriz de confusión para el modelo no 1.2)-

mosaic(cfmx1$table, shade = T, colorize = T,
       gp = gpar(fill = matrix(c("green3", "red2", "red2", "green3"), 2, 2)),
       main = "Matriz de confusión (modelo 2.1)")

Salida de consola y código no. 13 (Matriz de confusión ilustrada para el modelo no 1.2)-

Curva ROC modelo 2.1
##### Comparar modelo ya ajustado (no.1) con uno nuevo y muy simple
Prediccion_1_roc<-predict(object = modelo_2_1,
                    newdata = Test_intern_1, type = "response", 
                    se.fit = T)

predicciones_roc_1=prediction(Prediccion_1_roc$fit,Test_intern_1$Loan.Status)
Desempeno_roc_1<-performance(predicciones_roc_1,'tpr','fpr')

## Cgenerando nuevo modelo simple
modelo_simple_prueba <- glm(data = Train_intern_1, 
              formula = Loan.Status~Purpose+Bankruptcies, family = "binomial")

Prediccion_1_prueba<-predict(object = modelo_simple_prueba,
                    newdata = Test_intern_1, type = "response", 
                    se.fit = T)

predicciones1_prueba<-prediction(Prediccion_1_prueba$fit,Test_intern_1$Loan.Status)
Desempeno1_prueba<-performance(predicciones1_prueba,'tpr','fpr')

#ejecutar en conjunto
plot(Desempeno_roc_1,col="darkgreen")
plot(Desempeno1_prueba, col = "firebrick", add = T)
legend(x=.7, y=.4, legend = c("Simple","Completo"),
      lwd = 2, col = c("red","darkgreen"),
      title = "Modelos")
title("Curva ROC del modelo 2.1")

Gráfica no. 6 (Curva ROC modelo 2.1)-

Identificando los coeficientes modelo 2.1
nombres_vars = names(coef(modelo_2_1))
resumen = summary(modelo_2_1)

# Iterar por cada variable e imprimir su nombre y coeficiente
for (i in seq_along(nombres_vars)) {
  cat(paste(nombres_vars[i], ": ", resumen$coefficients[i, "Estimate"], "\n"))
}
## (Intercept) :  -4.12071498684856 
## Purposebuy a car :  1.13731763232773 
## Purposebuy house :  -0.0201690250725846 
## Purposedebt consolidation :  0.341045944606258 
## Purposeeducational expenses :  -0.953383735885413 
## Purposehome improvements :  0.242492905344161 
## Purposemajor_purchase :  -0.341616661009678 
## Purposemedical bills :  0.317032149542507 
## Purposemoving :  -0.60469119609014 
## Purposeother :  0.429061241732213 
## Purposerenewable_energy :  11.2255206940418 
## Purposesmall_business :  -0.91766985323541 
## Purposetake a trip :  -0.0878816912220709 
## Purposevacation :  0.43389064926999 
## Purposewedding :  -1.0891896737311 
## Current.Loan.Amount :  -7.16715220095716e-07 
## Credit.Score :  0.00677876856845841 
## Annual.Income :  4.94146651307591e-07 
## Monthly.Debt :  -1.23366879380709e-05 
## Months.since.last.delinquent :  0.0107501234588733

Salida de consola y código no. 14 (Coeficientes modelo no 2.1)-

Intervalos de confianza modelo no.1
suppressWarnings(
confint(object = modelo_2_1, level = 0.95))
## Waiting for profiling to be done...
##                                      2.5 %        97.5 %
## (Intercept)                  -5.415933e+00 -2.818757e+00
## Purposebuy a car              4.032216e-01  1.969489e+00
## Purposebuy house             -6.220912e-01  5.995633e-01
## Purposedebt consolidation    -3.440497e-02  6.991294e-01
## Purposeeducational expenses  -2.343939e+00  4.958366e-01
## Purposehome improvements     -1.742786e-01  6.470503e-01
## Purposemajor_purchase        -1.162913e+00  5.317129e-01
## Purposemedical bills         -2.142365e-01  8.571879e-01
## Purposemoving                -2.029354e+00  8.681899e-01
## Purposeother                  3.079703e-02  8.129298e-01
## Purposerenewable_energy      -1.322280e+01            NA
## Purposesmall_business        -1.831335e+00  1.743265e-02
## Purposetake a trip           -8.008136e-01  6.696588e-01
## Purposevacation              -1.013842e+00  2.360537e+00
## Purposewedding               -2.817091e+00  6.379996e-01
## Current.Loan.Amount          -1.065388e-06 -3.676268e-07
## Credit.Score                  4.979516e-03  8.569336e-03
## Annual.Income                 3.825349e-07  6.080412e-07
## Monthly.Debt                 -1.865794e-05 -6.015303e-06
## Months.since.last.delinquent  4.212787e-03  1.729535e-02

Salida de consola y código no. 15 (Intervalos de confianza para el modelo no 2.1)-

ODDS Ratios del modelo 2.1
  • Identificando los odds ratios \(e^{B_{k}}\)
nombres_vars = names(coef(modelo_2_1))
resumen = summary(modelo_2_1)

# Iterar por cada variable e imprimir su nombre y el num de euler elevado al coeficiente
for (i in seq_along(nombres_vars)) {
  cat(paste(nombres_vars[i], ": ", exp(resumen$coefficients[i, "Estimate"]), "\n"))
}
## (Intercept) :  0.0162329039789215 
## Purposebuy a car :  3.11839246145354 
## Purposebuy house :  0.980033009155944 
## Purposedebt consolidation :  1.40641785669244 
## Purposeeducational expenses :  0.385434605508224 
## Purposehome improvements :  1.27442220754439 
## Purposemajor_purchase :  0.710620561071819 
## Purposemedical bills :  1.37304671403958 
## Purposemoving :  0.546243082597199 
## Purposeother :  1.53581508749365 
## Purposerenewable_energy :  75020.8004981217 
## Purposesmall_business :  0.399448731675834 
## Purposetake a trip :  0.915869225682716 
## Purposevacation :  1.54325010336399 
## Purposewedding :  0.336489049178136 
## Current.Loan.Amount :  0.999999283285037 
## Credit.Score :  1.00680179642421 
## Annual.Income :  1.00000049414677 
## Monthly.Debt :  0.999987663388159 
## Months.since.last.delinquent :  1.01080811365035

Salida de consola y código no. 16 (Odds ratios para el modelo no 2.1)-

Evaluación supuestos para el modelo 2.1
  • VIF
#vif(modelo_2_1)
vif_modelo <- vif(modelo_2_1)
datos_vif <- data.frame(VIF = vif_modelo)
ggplot(datos_vif, aes(x = row.names(datos_vif), y = VIF.GVIF)) +
  geom_bar(stat = "identity", fill = "blue", width = 0.5) +
  coord_flip() +
  ggtitle("Valores de VIF para cada variable del modelo 2.1") +
  xlab("Variables del modelo") +
  ylab("Valor de VIF")+
  theme_bw()

Gráfica no. 7 (VIF modelo 2.1)-

  • La variable de interés es de respuesta binaria. Este supuesto se corrobora al inicio de este rmarkdown.
  • Las variables independientes estan relacionadas linealmente con la variable dependiente. Este supuesto se corroboró en el EDA.
Bootstraping

Empleamos como estadístico a evalua el % de aciertos (suma de verdaderos positivos y verdaderos negativos sobre el numero de observaciones en Train para este caso)

modelo.fun <- function(datos, subconjunto){
  #Generamos un vector de Predicción para la probabilidad de fully charged (1) para un nuevo vector random [llamado predicciones]
  predicciones<-predict(modelo_2_1, newdata = datos[subconjunto,])
  #Con dicho vector [predicciones], separar por clases para fully paid (1) y charged off(0) para generar el nuevo vector [Prediccion]
  Predicciones <- as.factor(ifelse(predicciones >= Corte_1[corte_optimo_1][1],
                                yes = 1, no = 0))
  # Instanciamos la matriz de confusion para la variable de respuesta       (loan.status) con el vector de probabilidades llamado predicciones
  cfmtx = confusionMatrix(data = as.factor(datos$Loan.Status), 
                             Predicciones)
  
  # Generamos un promedio de el porcentaje de aciertos general
  class_prop = (cfmtx$table[1] + cfmtx$table[4]) / length(Train_intern_1$Loan.Status)
  
  
  # Regresamos dicho promedio
  return(class_prop)
}


Bootstrap<-boot(data = Train_intern_1, 
                statistic = modelo.fun,## <- el parámetro 
                #### statistic requiere una función que reciba dos 
                #### parámetros: la muestra y los índices
                R = 10000) #10mil muestras simuladas
Bootstrap
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Train_intern_1, statistic = modelo.fun, R = 10000)
## 
## 
## Bootstrap Statistics :
##      original      bias     std. error
## t1* 0.7999138 -0.00298481 0.0008464804

Salida de consola y código no. 17 (Bootstraping en el modelo no 2.1)-

mean = mean(Bootstrap$t[,1])
bootstrap_vector_modelo_2_1 = Bootstrap$t[,1]

hist(bootstrap_vector_modelo_2_1, breaks = 100,
     xlab = "% de aciertos",
     main = "Bootstrap de la predicción para el modelo 2.1",
     ylab = "Frecuencia")
abline(v = quantile(x = Bootstrap$t[,1], c(.025,.975)),col = "red")
abline(v = mean(Bootstrap$t[,1]), col = "darkgreen", lwd = 2)

Gráfica no. 8 (Bootstraping en el modelo 1.2)-

Promedio de la estimación con 10mil muestras a partir de Train:

mean(Bootstrap$t[,1])
## [1] 0.7969289
Modelo no. 2.2
Matriz de confusión para el modelo base con 6 features
  • Respuesta del modelo en el conjunto de prueba
prediccion_2 = predict(object = modelo_2_2,
                    newdata = Test_intern_1, type = "response")
df_prediccion_2 <- data.frame(preds = prediccion_2)

dist2 = ggplot(df_prediccion_2, aes(y="preds", x=prediccion_2))+
  geom_jitter(alpha=0.5, color="greenyellow")+
  geom_violin(alpha=0.5)+
  geom_boxplot(width = 0.3, fill = "white", alpha = 0.5)+
  scale_x_continuous(breaks =c(0.6, 0.7, 0.8, 0.9),labels=c("60%", "70%", "80%", "90%"), name = "Respuesta")+
  scale_y_discrete(labels=c(""), name=NULL)+
  ggtitle("Predicciones del modelo 2.2 para el conjunto Test")+
  theme_bw()
dist2

Gráfica no. 9 (Distribución de la predicción para el modelo no 2.2)-

  • punto de corte ideal para maximizar la sensibilidad (vs Test)
prediccion_2 = predict(object = modelo_2_2,
                    newdata = Test_intern_1, type = "response")

Sensibilidad_2 <- vector()
Corte_1 <- seq(0.05, 0.95, by = 0.001)

for (i in 1:length(Corte_1)) {
  Prediccion <- as.factor(ifelse(prediccion_2 >= Corte_1[i], yes = 1, no = 0))
  
  Sensibilidad_2[i] <- tryCatch(
    {
      # confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), Prediccion)$overall[1]
      confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), 
                             Prediccion)$byClass[1]
    },
    error = function(e) {
      0 # O cualquier valor predeterminado
    }
  )
}

# plot(x=Corte_1, y = Sensibilidad_2, type = "l", main ="Punto de corte óptimo para maximizar la sensibilidad (vs Test)", xlab="Corte", ylab="Exactitud")
df_sensibilidad = data.frame(corte= Corte_1, sensibilidad = Sensibilidad_2)
sens2 = ggplot(df_sensibilidad, aes(x=Corte_1, y=Sensibilidad_2))+
  geom_line()+
  geom_point(aes(x = Corte_1[which.max(Sensibilidad_2)], y = Sensibilidad_2[which.max(Sensibilidad_2)]), 
             shape = 21, size = 5, fill = "white", color = "blue") +
  geom_text(aes(x = Corte_1[which.max(Sensibilidad_2)], y = Sensibilidad_2[which.max(Sensibilidad_2)] - 0.07, label = glue("Max \n{Corte_1[which.max(Sensibilidad_2)]}")),
            hjust = 1.3) +
  scale_y_continuous(breaks=c(0.25, 0.5, 0.75, 1), labels = c("25%", "50%", "75%", "100%"),
                     name = "Sensibilidad")+
  scale_x_continuous("Punto de corte")+
  ggtitle("Punto de corte ideal para máximizar la sensibilidad (modelo 2.2)")+
  theme_bw()
sens2

Gráfica no. 10 (punto de corte ideal para maximizar la sensibilidad (vs Test) en modelo no 2.2)-

  • punto de corte ideal para maximizar la especificidad (vs Test)
prediccion_2 = predict(object = modelo_2_2,
                    newdata = Test_intern_1, type = "response")

Especificidad_2 <- vector()
Corte_1 <- seq(0.05, 0.95, by = 0.001)

for (i in 1:length(Corte_1)) {
  Prediccion <- as.factor(ifelse(prediccion_2 >= Corte_1[i], yes = 1, no = 0))
  
  Especificidad_2[i] <- tryCatch(
    {
      # confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), Prediccion)$overall[1]
      confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), 
                             Prediccion)$byClass[2]
    },
    error = function(e) {
      0 # O cualquier valor predeterminado
    }
  )
}

# plot(x=Corte_1, y = Especificidad_2, type = "l", main ="Punto de corte óptimo para maximizar la especificidad (vs Test)", xlab="Corte", ylab="Exactitud")
df_especificidad = data.frame(corte= Corte_1, especificidad = Especificidad_2)
espe2 = ggplot(df_especificidad, aes(x=Corte_1, y=Especificidad_2))+
  geom_line()+
  geom_point(aes(x = Corte_1[which.max(Especificidad_2)], y = Especificidad_2[which.max(Especificidad_2)]), 
             shape = 21, size = 5, fill = "white", color = "blue") +
  geom_text(aes(x = Corte_1[which.max(Especificidad_2)], y = Especificidad_2[which.max(Especificidad_2)] - 0.03, label = glue("Max \n{Corte_1[which.max(Especificidad_2)]}")),
            hjust = 1.4) +
  scale_y_continuous(breaks=c(0.25, 0.50, 0.75), labels = c("25%", "50%", "75%"),
                     name = "Especificidad")+
  scale_x_continuous("Punto de corte")+
  ggtitle("Punto corte ideal para máximizar la especificidad (modelo 2.2)")+
  theme_bw()
espe2

Gráfica no. 11 (punto de corte ideal para maximizar la especificidad (vs Test) en modelo no 2.2)-

  • punto de corte ideal para maximizar la exactitud (vs Test)
prediccion_2 = predict(object = modelo_2_2,
                    newdata = Test_intern_1, type = "response")

Exactitud_2 <- vector()
Corte_1 <- seq(0.05, 0.95, by = 0.001)

for (i in 1:length(Corte_1)) {
  Prediccion <- as.factor(ifelse(prediccion_2 >= Corte_1[i], yes = 1, no = 0))
  
  Exactitud_2[i] <- tryCatch(
    {
      confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), Prediccion)$overall[1]
    },
    error = function(e) {
      0 # O cualquier valor predeterminado
    }
  )
}

# plot(x=Corte_1, y = Exactitud_2, type = "l", main ="Punto de corte óptimo para maximizar la exactitud (vs Test)", xlab="Corte", ylab="Exactitud")
df_exactitud = data.frame(corte= Corte_1, exactitud = Exactitud_2)
exac2 = ggplot(df_especificidad, aes(x=Corte_1, y=Exactitud_2))+
  geom_line()+
  geom_point(aes(x = Corte_1[which.max(Exactitud_2)], y = Exactitud_2[which.max(Exactitud_2)]), 
             shape = 21, size = 5, fill = "white", color = "blue") +
  geom_text(aes(x = Corte_1[which.max(Exactitud_2)], y = Exactitud_2[which.max(Exactitud_2)] - 0.05, label = glue("Max \n{Corte_1[which.max(Exactitud_2)]}")),
            hjust = 1.3) +
  scale_y_continuous(breaks=c(0.2, 0.4, 0.6, 0.8), labels = c("20%", "40%", "60%", "80%"),
                     name = "Exactitud")+
  scale_x_continuous("Punto de corte")+
  ggtitle("Punto de corte ideal para máximizar la exactitud (modelo 2.2)")+
  theme_bw()
exac2

Gráfica no. 12 (punto de corte ideal para maximizar la exactitud (vs Test) en modelo no 2.2)-

corte_optimo_1 = which(Exactitud_2==max(Exactitud_2))
Corte_1[corte_optimo_1][1]
## [1] 0.573
Corte_modelo_2_2 = Corte_1[corte_optimo_1][1]

Salida de consola y código no. 18 (Corte ótimo para máximizar la exactitud en modelo no 2.2)-

Filtro por EDA al vector de predicciones

filtro manual La variable Loan.Status adoptara el valor de 1 (Fully Paid) cuando: - El registro presente más de 50 unidades en la variable Number.of.Open.Accounts. - El registro presente un valor mayor a 50M en la variable Current.Loan.Amount - El registro presente un valor mayor a 100M en la variable Maximum.Open.Credit - El registro presente un valor mayor a 125K en la variable Monthly.Debt

#Obteniendo los indices de Test_intern_1 en donde se cumplan los requisitos para eplicar 0 o 1 por filto manual

indices <- which(Test_intern_1$Number.of.Open.Accounts > 50 | Test_intern_1$Current.Loan.Amount > 50000000 | Test_intern_1$Maximum.Open.Credit > 100000000 | Test_intern_1$Monthly.Debt > 125000)
indices
## integer(0)

La variable Loan.Status adoptara el valor de 0 (Charged off) cuando: - El registro presente un valor mayor a 6K en la variable Credit.Score

indices <- which(Test_intern_1$Credit.Score > 6000000)
indices
## integer(0)

Dado que ningún registro presenta dichas condiciones del filtro manual para aplicar directamente y a criterios la probabilidad de 1 o 0 respectivamente, no se sufre ninguna variación en la matriz de confusión en si se aplica dicho filtro o no.

Prediccion_2<-as.factor(ifelse(prediccion_2>=Corte_1[corte_optimo_1][1],yes = 1, no = 0))
confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), 
                             Prediccion_2)
## Confusion Matrix and Statistics
## 
##           Reference
## Prediction    0    1
##          0    2  604
##          1    0 2486
##                                           
##                Accuracy : 0.8047          
##                  95% CI : (0.7902, 0.8185)
##     No Information Rate : 0.9994          
##     P-Value [Acc > NIR] : 1               
##                                           
##                   Kappa : 0.0053          
##                                           
##  Mcnemar's Test P-Value : <2e-16          
##                                           
##             Sensitivity : 1.0000000       
##             Specificity : 0.8045307       
##          Pos Pred Value : 0.0033003       
##          Neg Pred Value : 1.0000000       
##              Prevalence : 0.0006468       
##          Detection Rate : 0.0006468       
##    Detection Prevalence : 0.1959897       
##       Balanced Accuracy : 0.9022654       
##                                           
##        'Positive' Class : 0               
## 
cfmx2 = confusionMatrix(data = as.factor(Test_intern_1$Loan.Status), 
                             Prediccion_2)

Salida de consola y código no. 19 (Matriz de confusión para el modelo no 2.2)-

mosaic(cfmx2$table, shade = T, colorize = T,
       gp = gpar(fill = matrix(c("green3", "red2", "red2", "green3"), 2, 2)),
       main = "Matriz de confusión (modelo 2.2)")

Salida de consola y código no. 20 (Matriz de confusión ilustrada para el modelo no 2.2)-

Curva ROC modelo 2.2
##### Comparar modelo ya ajustado (no.1) con uno nuevo y muy simple
Prediccion_2_roc<-predict(object = modelo_2_2,
                    newdata = Test_intern_1, type = "response", 
                    se.fit = T)

predicciones_roc_2=prediction(Prediccion_2_roc$fit,Test_intern_1$Loan.Status)
Desempeno_roc_2<-performance(predicciones_roc_2,'tpr','fpr')

#ejecutar en conjunto
plot(Desempeno_roc_2,col="darkgreen")
plot(Desempeno1_prueba, col = "firebrick", add = T)
legend(x=.7, y=.4, legend = c("Simple","Completo"),
      lwd = 2, col = c("red","darkgreen"),
      title = "Modelos")
title("Curva ROC del modelo 2.2")

Gráfica no. 13 (Curva ROC modelo 2.2)-

Identificando los coeficientes modelo 2.2
nombres_vars = names(coef(modelo_2_2))
resumen = summary(modelo_2_2)

# Iterar por cada variable e imprimir su nombre y coeficiente
for (i in seq_along(nombres_vars)) {
  cat(paste(nombres_vars[i], ": ", resumen$coefficients[i, "Estimate"], "\n"))
}
## (Intercept) :  -4.00242427128793 
## Current.Loan.Amount :  -7.59611478635483e-07 
## Credit.Score :  0.00681907623649969 
## Annual.Income :  4.74824737452753e-07 
## Monthly.Debt :  -1.16253487015612e-05 
## Years.of.Credit.History :  0.0105742451654072 
## Months.since.last.delinquent :  0.0112506811614992

Salida de consola y código no. 21 (Coeficientes modelo no 2.2)-

Intervalos de confianza modelo 2.2
suppressWarnings(
confint(object = modelo_2_2, level = 0.95)
)
## Waiting for profiling to be done...
##                                      2.5 %        97.5 %
## (Intercept)                  -5.252685e+00 -2.745875e+00
## Current.Loan.Amount          -1.096237e-06 -4.229663e-07
## Credit.Score                  5.069649e-03  8.559729e-03
## Annual.Income                 3.649823e-07  5.869423e-07
## Monthly.Debt                 -1.787591e-05 -5.376094e-06
## Years.of.Credit.History       2.308956e-03  1.891024e-02
## Months.since.last.delinquent  4.727480e-03  1.778216e-02

Salida de consola y código no. 22 (Intervalos de confianza para el modelo no 1.3)-

ODDS Ratios del modelo 2.2
  • Identificando los odds ratios \(e^{B_{k}}\)
nombres_vars = names(coef(modelo_2_2))
resumen = summary(modelo_2_2)

# Iterar por cada variable e imprimir su nombre y el num de euler elevado al coeficiente
for (i in seq_along(nombres_vars)) {
  cat(paste(nombres_vars[i], ": ", exp(resumen$coefficients[i, "Estimate"]), "\n"))
}
## (Intercept) :  0.0182712905891308 
## Current.Loan.Amount :  0.99999924038881 
## Credit.Score :  1.00684237907469 
## Annual.Income :  1.00000047482485 
## Monthly.Debt :  0.999988374718873 
## Years.of.Credit.History :  1.01063035007713 
## Months.since.last.delinquent :  1.01131420809168

Salida de consola y código no. 23 (Odds ratios para el modelo no 2.2)-

Interpretación Odds ratios en el modelo no 2.2

\(-\)

Evaluación supuestos para el modelo 2.2
  • VIF
vif(modelo_2_2)
##          Current.Loan.Amount                 Credit.Score 
##                     1.449714                     1.092051 
##                Annual.Income                 Monthly.Debt 
##                     1.594329                     1.448372 
##      Years.of.Credit.History Months.since.last.delinquent 
##                     1.038071                     1.008373
vif_modelo2 <- vif(modelo_2_2)
datos_vif <- data.frame(VIF = vif_modelo)
ggplot(datos_vif, aes(x = row.names(datos_vif), y = vif_modelo2)) +
  geom_bar(stat = "identity", fill = "blue", width = 0.5) +
  coord_flip() +
  ggtitle("Valores de VIF para cada variable del modelo 2.2") +
  xlab("Variables del modelo") +
  ylab("Valor de VIF")+
  theme_bw()

Gráfica no. 14 (VIF modelo 2.2)-

  • No existen problemas de multicolinealidad.
  • La variable de interés es de respuesta binaria.
  • Las variables independientes estan relacionadas linealmente con la variable dependiente.
Bootstraping para el modelo 2.2

Empleamos como estadístico a evalua el % de aciertos (suma de verdaderos positivos y verdaderos negativos sobre el numero de observaciones en Train para este caso)

modelo.fun <- function(datos, subconjunto){
  #Generamos un vector de Predicción para la probabilidad de fully charged (1) para un nuevo vector random [llamado predicciones]
  predicciones<-predict(modelo_2_2, newdata = datos[subconjunto,])
  #Con dicho vector [predicciones], separar por clases para fully paid (1) y charged off(0) para generar el nuevo vector [Prediccion]
  Predicciones <- as.factor(ifelse(predicciones >= Corte_1[corte_optimo_1][1],
                                yes = 1, no = 0))
  # Instanciamos la matriz de confusion para la variable de respuesta       (loan.status) con el vector de probabilidades llamado predicciones
  cfmtx = confusionMatrix(data = as.factor(datos$Loan.Status), 
                             Predicciones)
  
  # Generamos un promedio de el porcentaje de aciertos general
  class_prop = (cfmtx$table[1] + cfmtx$table[4]) / length(Train_intern_1$Loan.Status)
  
  
  # Regresamos dicho promedio
  return(class_prop)
}


Bootstrap<-boot(data = Train_intern_1, 
                statistic = modelo.fun,## <- el parámetro 
                #### statistic requiere una función que reciba dos 
                #### parámetros: la muestra y los índices
                R = 10000) #10mil muestras simuladas
Bootstrap
## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Train_intern_1, statistic = modelo.fun, R = 10000)
## 
## 
## Bootstrap Statistics :
##      original       bias    std. error
## t1* 0.7988357 -0.004603849 0.001083998

Salida de consola y código no. 24 (Bootstraping en el modelo no 2.2)-

mean = mean(Bootstrap$t[,1])
bootstrap_vector_modelo_2_2 = Bootstrap$t[,1]

hist(bootstrap_vector_modelo_2_2, breaks = 100,
     xlab = "% de aciertos",
     main = "Bootstrap de la predicción para el modelo 2.2",
     ylab = "Frecuencia")
abline(v = quantile(x = Bootstrap$t[,1], c(.025,.975)),col = "red")
abline(v = mean(Bootstrap$t[,1]), col = "darkgreen", lwd = 2)

Gráfica no. 15 (Bootstraping en el modelo 2.2)-

Promedio de la estimación con 10mil muestras a partir de Train:

mean(Bootstrap$t[,1])
## [1] 0.7942319

Comparación de remuestreos en los modelos

df <- data.frame(valor = c(bootstrap_vector_modelo_2_1, bootstrap_vector_modelo_2_2),
                 grupo = factor(rep(c("Modelo 2.1", "Modelo 2.2"), each = length(bootstrap_vector_modelo_2_1))))

medias <- tapply(df$valor, df$grupo, mean)

ggplot(df, aes(x = valor, fill = grupo)) +
  geom_histogram(alpha = 0.6, position = "identity", bins = 100) +
  scale_fill_manual(values = c("blue", "red")) +
  labs(x = "Valor", y = "Frecuencia", fill = "Grupo")+
  ggtitle("Comparación para cada modelo en el remuestreo (10k simulaciones)")+
  geom_vline(xintercept = mean(bootstrap_vector_modelo_2_1), color = "blue")+
  geom_vline(xintercept = mean(bootstrap_vector_modelo_2_2), color = "red")+
  theme_bw()

plot(Desempeno_roc_1,col="darkblue")
plot(Desempeno_roc_2,col="darkgreen", add = T)
plot(Desempeno1_prueba, col = "firebrick", add = T)
legend(x=.6, y=.5, legend = c("Simple","Modelo 2.1", "Modelo 2.2"),
      lwd = 2, col = c("red", "darkblue","darkgreen"),
      title = "Modelos")
title("Curvas ROC modelado 2")

Gráfica no. 17 (Compraración ROC del modelo 1.1 y modelo 1.2)-

Conclusiones Generales

Dado que las pruebas con MSPR y accuracy por matriz de confusión fueron tecnicamente similares, decidimos optar por quedarnos con el modelo_2_1 con 6 features porque es un modelo con menos que nos arroja una media del estadśitico de prueba ligeramente mayor.

Del mejor modelo de la comparación, podemos interpretar de la anterior salida de consola los ODDS ratios que:

  • Si la variable Current.Loan.Amount aumenta en 1,000 unidades, este evento estara asociado con una disminución del 7.5% en las probabilidades del pago al corriente (Fully Paid).
  • Si la variable Credit.Score aumenta en 100 unidades, este evento estara asociado con un aumento del 68%% en las probabilidades del pago al corriente (Fully Paid).
  • Si la variable Annual.Income aumenta en 1,000 unidades, este evento estara asociado con un aumento del 4.7% en las probabilidades del pago al corriente (Fully Paid).
  • Si la variable Monthly.Debt aumenta en 100 unidades, este evento estara asociado con una disminución del 11.6% en las probabilidades del pago al corriente (Fully Paid).
  • Si la variable Years.of.Credit.History aumenta en 1 unidades, este evento estara asociado con un aumento del 1% en las probabilidades del pago al corriente (Fully Paid).
  • Si la variable Months.since.last.delinquent aumenta en 1 unidad, este evento estara asociado con un aumento del 1.1% en las probabilidades del pago al corriente (Fully Paid).

EL resto de variables nominales se encuentran representadas en el intercepto, aunque este no es estadísticamente significativo.

De lo que podemos concluir que las variables con mayor incidencia en este modelo son: - la variable Years.of.Credit.History - la variable Months.since.last.delinquent

La variable Loan.Status adoptara el valor de 0 (Charged off) cuando: - El registro presente un valor mayor a 6K en la variable Credit.Score

Exportando modelo al directorio
cfmx_modelo2 = cfmx2
save(modelo_2_2, cfmx_modelo2, Corte_modelo_2_2, file = "../models/Modelo_2.RData")