Objetivos

  1. Regresion Lineal logística

Se utilizará el dataset \(\textbf{Default}\). Es una base de datos que contiene 10000 registros de clientes de un banco. La idea es poder inferir, en base a los registros, que clientes no podrán pagar sus cuota de tarjeta de credito y caer en default.

names(Default)
[1] "default"       "student"       "balance"       "income"        "coded_default"
attach(Default)
The following objects are masked from Default (pos = 4):

    balance, coded_default, default, income, student
head(Default)
summary(Default)
 default    student       balance           income      coded_default   
 No :9667   No :7056   Min.   :   0.0   Min.   :  772   Min.   :0.0000  
 Yes: 333   Yes:2944   1st Qu.: 481.7   1st Qu.:21340   1st Qu.:0.0000  
                       Median : 823.6   Median :34553   Median :0.0000  
                       Mean   : 835.4   Mean   :33517   Mean   :0.0333  
                       3rd Qu.:1166.3   3rd Qu.:43808   3rd Qu.:0.0000  
                       Max.   :2654.3   Max.   :73554   Max.   :1.0000  

Vamos a dicotomizar la variable Default y agregarla al dataset con el nombre de coded_default. Vemos como la regresion lineal estándar no sirve para modelar lo que sucede con los datos, obtenemos valores fuera de rango y no se ajusta a las variaciones bruzcas, ademas que el error cometido en el rango entre extremos es muy alto.

Default$coded_default <- ifelse(Default$default=="Yes",1,0)
Default %>% ggplot(aes(x=balance,y=coded_default))+
  geom_point(aes(col=Default$default))+
  geom_smooth(method = "lm")

De la recta de regresión, podemos concluir que a medida que crece que balance de la tarjeta de crédito (deuda), mayor es la probabilidad de entrar en default. Por las limitaciones de MCO en este caso, solo podemos extraer conclusiones cualitativas.

fit_lm <- lm(coded_default ~ balance, data=Default)
summary(fit_lm)

Call:
lm(formula = coded_default ~ balance, data = Default)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.23533 -0.06939 -0.02628  0.02004  0.99046 

Coefficients:
                Estimate   Std. Error t value            Pr(>|t|)    
(Intercept) -0.075191959  0.003354360  -22.42 <0.0000000000000002 ***
balance      0.000129872  0.000003475   37.37 <0.0000000000000002 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1681 on 9998 degrees of freedom
Multiple R-squared:  0.1226,    Adjusted R-squared:  0.1225 
F-statistic:  1397 on 1 and 9998 DF,  p-value: < 0.00000000000000022
\(\beta_0\) \(\beta_1\)
-0.075192 0.0001299

\(\beta_1\) es cuanto aumenta la propabilidad de caer en default cuando aumenta en una unidad el balance (deuda)

Formulación de la regresión logística

Simil a teoria de desición. Dado un balance, quiero saber la probabilidad de caer en default. Se puede definir un umbral de probabilidad a traves del cual se toma la desición.

Trabajamos con el logaritmos de Odd (chance), que se define como \(\frac{p}{1-p}\). Donde \(p = P(caer\_en\_mora)\)

\[ P(default = Yes|balance) \]

\[ P(default = Yes|balance) = \beta_0 + \beta_1 balance \] \[ P(X) = \beta_0 + \beta_1 X \]

\[ ln(\frac{p}{1-p}) = \beta_0 + \beta_1 X \]

\[ \frac{p}{1-p} = e^{\beta_0 + \beta_1 X} \]

\[ p = e^{\beta_0 + \beta_1 X} (1-p) \]

\[ p = e^{\beta_0 + \beta_1 X} - p e^{\beta_0 + \beta_1 X} \]

\[ p + p e^{\beta_0 + \beta_1 X} = e^{\beta_0 + \beta_1 X} \] \[ p (1 + e^{\beta_0 + \beta_1 X}) = e^{\beta_0 + \beta_1 X} \]

\[ p = \frac{e^{\beta_0 + \beta_1 X}}{(1 + e^{\beta_0 + \beta_1 X})} \]

Calculamos los parámetros \(\beta_0\) y \(\beta_1\) de la regresion logística

logit_reg <- glm(default ~ balance, data=Default, family = binomial)
summary(logit_reg)$coefficients
                 Estimate   Std. Error   z value
(Intercept) -10.651330614 0.3611573721 -29.49221
balance       0.005498917 0.0002203702  24.95309
                                                                                                                                                                                                           Pr(>|z|)
(Intercept) 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003623124
balance     0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001976601736462045296666174727056650651646940960033647670370745

Odds aumentan 0.0054989 por cada dolar que tengo en balance (deuda). Podemos resumir los resultados del modelo de la siguiente manera

\[ ln Odds = \beta_0 + \beta_1 \cdot balance \] donde \(\beta_0 =\) -10.6513306 y \(\beta_1 =\) 0.0054989

\[ ln Odds = -10.6513306 + 0.0054989 \cdot balance \]

\[ Odds = e^{-10.6513306 + 0.0054989 \cdot balance} \] \[ p = \frac{e^{-10.6513306 + 0.0054989 \cdot balance}}{1+e^{-10.6513306 + 0.0054989 \cdot balance}} \]

beta0 <- logit_reg$coefficients[1]
beta1 <- logit_reg$coefficients[2]
p <- exp(beta0+beta1*Default$balance)/(1+exp(beta0+beta1*Default$balance))
plot(Default$balance,p)

Al aumentar el balance en la tarjeta de credito aumenta la probabilidad \(p\) de caer en mora.

Armamos una funcion para calcular la probabilidad.

prob_default <- function(beta0, beta1, balance){
  p <- exp(beta0+beta1*balance)/(1+exp(beta0+beta1*balance))
  return(as.numeric(p))
}

prob_default(beta0, beta1, 1500)
[1] 0.08294762
prob_default(beta0, beta1, 1800)
[1] 0.320107
prob_default(beta0, beta1, 2000)
[1] 0.5857694
prob_default(beta0, beta1, 2500)
[1] 0.9567259

Interpretacion de los coeficientes

summary(logit_reg)$coefficients
                 Estimate   Std. Error   z value
(Intercept) -10.651330614 0.3611573721 -29.49221
balance       0.005498917 0.0002203702  24.95309
                                                                                                                                                                                                           Pr(>|z|)
(Intercept) 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000003623124
balance     0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001976601736462045296666174727056650651646940960033647670370745

Funciones para hallar coeficientes e intervalos de confianza para regresiones. Me interesa si en el intervalo se encuentra \(e^{0} = 1\), si no pertence entonces puede asegurar significancia estadística (test de hipótesis).

coef(logit_reg)
  (Intercept)       balance 
-10.651330614   0.005498917 
confint(logit_reg)
Waiting for profiling to be done...
                    2.5 %       97.5 %
(Intercept) -11.383288936 -9.966565064
balance       0.005078926  0.005943365

Generamos otro modelo contemplando todas las variables del dataframe para evaluar que estadisticos son relevantes para el análisis de la mora.

logit_reg2 <- glm(default ~ balance + student + income, data = default, family = binomial)
summary(logit_reg2)$coefficients
                   Estimate     Std. Error    z value
(Intercept) -10.86904519617 0.492255515606 -22.080088
balance       0.00573650526 0.000231894519  24.737563
studentYes   -0.64677580664 0.236252528745  -2.737646
income        0.00000303345 0.000008202615   0.369815
                                                                                                                                                   Pr(>|z|)
(Intercept) 0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004911279576424167001347828332502319
balance     0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004219578
studentYes  0.006188063286482599519022773648657675948925316333770751953125000000000000000000000000000000000000000000000000000000000000000000000000000000000
income      0.711520342121125359824418410426005721092224121093750000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
LS0tCnRpdGxlOiAiUmVncmVzacOzbiBMaW5lYWwgVlJTIChyZWdyZXNpw7NuIGxvZ8Otc3RpY2EpIgpvdXRwdXQ6CiAgaHRtbF9ub3RlYm9vazogZGVmYXVsdAogIHBkZl9kb2N1bWVudDogZGVmYXVsdAotLS0KCmBgYHtyIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0UsIGluY2x1ZGU9RkFMU0UsIHBhZ2VkLnByaW50PUZBTFNFfQojaW5zdGFsbC5wYWNrYWdlcygnSVNMUicpCiNpbnN0YWxsLnBhY2thZ2VzKCdNQVNTJykKI2luc3RhbGwucGFja2FnZXMoJ2NhcicpCiNpbnN0YWxsLnBhY2thZ2VzKCdxdWFudHJlZycpCiNpbnN0YWxsLnBhY2thZ2VzKCdzcGFyc2VNJykKCmxpYnJhcnkoU3BhcnNlTSkKbGlicmFyeShjYXIpCmxpYnJhcnkoTUFTUykKbGlicmFyeShJU0xSKQpsaWJyYXJ5KGxtdGVzdCkKbGlicmFyeShkcGx5cikKbGlicmFyeShnZ3Bsb3QyKQoKb3B0aW9ucyhzY2lwZW4gPSA5OTk5OSkgIyBjYW1iaW8gZm9ybWF0byBkZSBzaG93IGRlY2ltYWxzCnNldHdkKCcvaG9tZS9qZ29uemFsZXovZGV2L2RhdGFTY2llbmNlJykKYGBgCgojIyBPYmpldGl2b3MKCjEuIFJlZ3Jlc2lvbiBMaW5lYWwgbG9nw61zdGljYQoKU2UgdXRpbGl6YXLDoSBlbCBkYXRhc2V0ICRcdGV4dGJme0RlZmF1bHR9JC4gRXMgdW5hIGJhc2UgZGUgZGF0b3MgcXVlIGNvbnRpZW5lIGByIGxlbmd0aChEZWZhdWx0JGJhbGFuY2UpYCByZWdpc3Ryb3MgZGUgY2xpZW50ZXMgZGUgdW4gYmFuY28uIExhIGlkZWEgZXMgcG9kZXIgaW5mZXJpciwgZW4gYmFzZSBhIGxvcyByZWdpc3Ryb3MsIHF1ZSBjbGllbnRlcyBubyBwb2Ryw6FuIHBhZ2FyIHN1cyBjdW90YSBkZSB0YXJqZXRhIGRlIGNyZWRpdG8geSBjYWVyIGVuIGRlZmF1bHQuCgpgYGB7cn0KbmFtZXMoRGVmYXVsdCkKYXR0YWNoKERlZmF1bHQpCmBgYAoKYGBge3J9CmhlYWQoRGVmYXVsdCkKYGBgCgpgYGB7cn0Kc3VtbWFyeShEZWZhdWx0KQpgYGAKClZhbW9zIGEgZGljb3RvbWl6YXIgbGEgdmFyaWFibGUgRGVmYXVsdCB5IGFncmVnYXJsYSBhbCBkYXRhc2V0IGNvbiBlbCBub21icmUgZGUgY29kZWRfZGVmYXVsdC4gVmVtb3MgY29tbyBsYSByZWdyZXNpb24gbGluZWFsIGVzdMOhbmRhciBubyBzaXJ2ZSBwYXJhIG1vZGVsYXIgbG8gcXVlIHN1Y2VkZSBjb24gbG9zIGRhdG9zLCBvYnRlbmVtb3MgdmFsb3JlcyBmdWVyYSBkZSByYW5nbyB5IG5vIHNlIGFqdXN0YSBhIGxhcyB2YXJpYWNpb25lcyBicnV6Y2FzLCBhZGVtYXMgcXVlIGVsIGVycm9yIGNvbWV0aWRvIGVuIGVsIHJhbmdvIGVudHJlIGV4dHJlbW9zIGVzIG11eSBhbHRvLgoKYGBge3J9CkRlZmF1bHQkY29kZWRfZGVmYXVsdCA8LSBpZmVsc2UoRGVmYXVsdCRkZWZhdWx0PT0iWWVzIiwxLDApCkRlZmF1bHQgJT4lIGdncGxvdChhZXMoeD1iYWxhbmNlLHk9Y29kZWRfZGVmYXVsdCkpKwogIGdlb21fcG9pbnQoYWVzKGNvbD1EZWZhdWx0JGRlZmF1bHQpKSsKICBnZW9tX3Ntb290aChtZXRob2QgPSAibG0iKQpgYGAKCkRlIGxhIHJlY3RhIGRlIHJlZ3Jlc2nDs24sIHBvZGVtb3MgY29uY2x1aXIgcXVlIGEgbWVkaWRhIHF1ZSBjcmVjZSBxdWUgYmFsYW5jZSBkZSBsYSB0YXJqZXRhIGRlIGNyw6lkaXRvIChkZXVkYSksIG1heW9yIGVzIGxhIHByb2JhYmlsaWRhZCBkZSBlbnRyYXIgZW4gZGVmYXVsdC4gUG9yIGxhcyBsaW1pdGFjaW9uZXMgZGUgTUNPIGVuIGVzdGUgY2Fzbywgc29sbyBwb2RlbW9zIGV4dHJhZXIgY29uY2x1c2lvbmVzIGN1YWxpdGF0aXZhcy4KCmBgYHtyfQpmaXRfbG0gPC0gbG0oY29kZWRfZGVmYXVsdCB+IGJhbGFuY2UsIGRhdGE9RGVmYXVsdCkKc3VtbWFyeShmaXRfbG0pCmBgYAoKIHwgJFxiZXRhXzAkIHwgJFxiZXRhXzEkIHwgCiB8LS0tLS0tLS0tLS18LS0tLS0tLS0tLS18CiB8YHIgZml0X2xtJGNvZWZmaWNpZW50c1sxXWB8YHIgZml0X2xtJGNvZWZmaWNpZW50c1syXWB8CiAKICRcYmV0YV8xJCBlcyBjdWFudG8gYXVtZW50YSBsYSBwcm9wYWJpbGlkYWQgZGUgY2FlciBlbiBkZWZhdWx0IGN1YW5kbyBhdW1lbnRhIGVuIHVuYSB1bmlkYWQgZWwgYmFsYW5jZSAoZGV1ZGEpCgojIyBGb3JtdWxhY2nDs24gZGUgbGEgcmVncmVzacOzbiBsb2fDrXN0aWNhCgpTaW1pbCBhIHRlb3JpYSBkZSBkZXNpY2nDs24uIERhZG8gdW4gYmFsYW5jZSwgcXVpZXJvIHNhYmVyIGxhIHByb2JhYmlsaWRhZCBkZSBjYWVyIGVuIGRlZmF1bHQuIFNlIHB1ZWRlIGRlZmluaXIgdW4gdW1icmFsIGRlIHByb2JhYmlsaWRhZCBhIHRyYXZlcyBkZWwgY3VhbCBzZSB0b21hIGxhIGRlc2ljacOzbi4KClRyYWJhamFtb3MgY29uIGVsIGxvZ2FyaXRtb3MgZGUgT2RkIChjaGFuY2UpLCBxdWUgc2UgZGVmaW5lIGNvbW8gJFxmcmFje3B9ezEtcH0kLiBEb25kZSAkcCA9IFAoY2FlclxfZW5cX21vcmEpJCAKCiQkClAoZGVmYXVsdCA9IFllc3xiYWxhbmNlKQokJAoKJCQKUChkZWZhdWx0ID0gWWVzfGJhbGFuY2UpID0gXGJldGFfMCArIFxiZXRhXzEgYmFsYW5jZQokJAokJApQKFgpID0gXGJldGFfMCArIFxiZXRhXzEgWAokJAoKJCQKbG4oXGZyYWN7cH17MS1wfSkgPSBcYmV0YV8wICsgXGJldGFfMSBYCiQkCgokJApcZnJhY3twfXsxLXB9ID0gZV57XGJldGFfMCArIFxiZXRhXzEgWH0KJCQKCiQkCnAgPSBlXntcYmV0YV8wICsgXGJldGFfMSBYfSAoMS1wKQokJAoKJCQKcCA9IGVee1xiZXRhXzAgKyBcYmV0YV8xIFh9IC0gcCBlXntcYmV0YV8wICsgXGJldGFfMSBYfQokJAoKJCQKcCArIHAgZV57XGJldGFfMCArIFxiZXRhXzEgWH0gPSBlXntcYmV0YV8wICsgXGJldGFfMSBYfQokJAokJApwICgxICsgZV57XGJldGFfMCArIFxiZXRhXzEgWH0pID0gZV57XGJldGFfMCArIFxiZXRhXzEgWH0KJCQKCiQkCnAgPSBcZnJhY3tlXntcYmV0YV8wICsgXGJldGFfMSBYfX17KDEgKyBlXntcYmV0YV8wICsgXGJldGFfMSBYfSl9CiQkCgpDYWxjdWxhbW9zIGxvcyBwYXLDoW1ldHJvcyAkXGJldGFfMCQgeSAkXGJldGFfMSQgZGUgbGEgcmVncmVzaW9uIGxvZ8Otc3RpY2EKCmBgYHtyfQpsb2dpdF9yZWcgPC0gZ2xtKGRlZmF1bHQgfiBiYWxhbmNlLCBkYXRhPURlZmF1bHQsIGZhbWlseSA9IGJpbm9taWFsKQpzdW1tYXJ5KGxvZ2l0X3JlZykkY29lZmZpY2llbnRzCmBgYAoKT2RkcyBhdW1lbnRhbiBgciBsb2dpdF9yZWckY29lZmZpY2llbnRzWzJdYCBwb3IgY2FkYSBkb2xhciBxdWUgdGVuZ28gZW4gYmFsYW5jZSAoZGV1ZGEpLiBQb2RlbW9zIHJlc3VtaXIgbG9zIHJlc3VsdGFkb3MgZGVsIG1vZGVsbyBkZSBsYSBzaWd1aWVudGUgbWFuZXJhCgokJApsbiBPZGRzID0gXGJldGFfMCArIFxiZXRhXzEgXGNkb3QgYmFsYW5jZQokJApkb25kZSAkXGJldGFfMCA9JCBgciBsb2dpdF9yZWckY29lZmZpY2llbnRzWzFdYCB5ICRcYmV0YV8xID0kIGByIGxvZ2l0X3JlZyRjb2VmZmljaWVudHNbMl1gCgokJApsbiBPZGRzID0gYHIgbG9naXRfcmVnJGNvZWZmaWNpZW50c1sxXWAgKyBgciBsb2dpdF9yZWckY29lZmZpY2llbnRzWzJdYCBcY2RvdCBiYWxhbmNlCiQkCgokJApPZGRzID0gZV57YHIgbG9naXRfcmVnJGNvZWZmaWNpZW50c1sxXWAgKyBgciBsb2dpdF9yZWckY29lZmZpY2llbnRzWzJdYCBcY2RvdCBiYWxhbmNlfQokJAokJApwID0gXGZyYWN7ZV57YHIgbG9naXRfcmVnJGNvZWZmaWNpZW50c1sxXWAgKyBgciBsb2dpdF9yZWckY29lZmZpY2llbnRzWzJdYCBcY2RvdCBiYWxhbmNlfX17MStlXntgciBsb2dpdF9yZWckY29lZmZpY2llbnRzWzFdYCArIGByIGxvZ2l0X3JlZyRjb2VmZmljaWVudHNbMl1gIFxjZG90IGJhbGFuY2V9fQokJAoKCmBgYHtyfQpiZXRhMCA8LSBsb2dpdF9yZWckY29lZmZpY2llbnRzWzFdCmJldGExIDwtIGxvZ2l0X3JlZyRjb2VmZmljaWVudHNbMl0KcCA8LSBleHAoYmV0YTArYmV0YTEqRGVmYXVsdCRiYWxhbmNlKS8oMStleHAoYmV0YTArYmV0YTEqRGVmYXVsdCRiYWxhbmNlKSkKcGxvdChEZWZhdWx0JGJhbGFuY2UscCkKYGBgCgpBbCBhdW1lbnRhciBlbCBiYWxhbmNlIGVuIGxhIHRhcmpldGEgZGUgY3JlZGl0byBhdW1lbnRhIGxhIHByb2JhYmlsaWRhZCAkcCQgZGUgY2FlciBlbiBtb3JhLgoKQXJtYW1vcyB1bmEgZnVuY2lvbiBwYXJhIGNhbGN1bGFyIGxhIHByb2JhYmlsaWRhZC4KCmBgYHtyfQpwcm9iX2RlZmF1bHQgPC0gZnVuY3Rpb24oYmV0YTAsIGJldGExLCBiYWxhbmNlKXsKICBwIDwtIGV4cChiZXRhMCtiZXRhMSpiYWxhbmNlKS8oMStleHAoYmV0YTArYmV0YTEqYmFsYW5jZSkpCiAgcmV0dXJuKGFzLm51bWVyaWMocCkpCn0KCnByb2JfZGVmYXVsdChiZXRhMCwgYmV0YTEsIDE1MDApCnByb2JfZGVmYXVsdChiZXRhMCwgYmV0YTEsIDE4MDApCnByb2JfZGVmYXVsdChiZXRhMCwgYmV0YTEsIDIwMDApCnByb2JfZGVmYXVsdChiZXRhMCwgYmV0YTEsIDI1MDApCmBgYAoKIyMgSW50ZXJwcmV0YWNpb24gZGUgbG9zIGNvZWZpY2llbnRlcwoKYGBge3J9CnN1bW1hcnkobG9naXRfcmVnKSRjb2VmZmljaWVudHMKYGBgCgpGdW5jaW9uZXMgcGFyYSBoYWxsYXIgY29lZmljaWVudGVzIGUgaW50ZXJ2YWxvcyBkZSBjb25maWFuemEgcGFyYSByZWdyZXNpb25lcy4gTWUgaW50ZXJlc2Egc2kgZW4gZWwgaW50ZXJ2YWxvIHNlIGVuY3VlbnRyYSAkZV57MH0gPSAxJCwgc2kgbm8gcGVydGVuY2UgZW50b25jZXMgcHVlZGUgYXNlZ3VyYXIgc2lnbmlmaWNhbmNpYSBlc3RhZMOtc3RpY2EgKHRlc3QgZGUgaGlww7N0ZXNpcykuCgpgYGB7cn0KY29lZihsb2dpdF9yZWcpCmNvbmZpbnQobG9naXRfcmVnKQpgYGAKCkdlbmVyYW1vcyBvdHJvIG1vZGVsbyBjb250ZW1wbGFuZG8gdG9kYXMgbGFzIHZhcmlhYmxlcyBkZWwgZGF0YWZyYW1lIHBhcmEgZXZhbHVhciBxdWUgZXN0YWRpc3RpY29zIHNvbiByZWxldmFudGVzIHBhcmEgZWwgYW7DoWxpc2lzIGRlIGxhIG1vcmEuCgpgYGB7cn0KbG9naXRfcmVnMiA8LSBnbG0oZGVmYXVsdCB+IGJhbGFuY2UgKyBzdHVkZW50ICsgaW5jb21lLCBkYXRhID0gZGVmYXVsdCwgZmFtaWx5ID0gYmlub21pYWwpCnN1bW1hcnkobG9naXRfcmVnMikkY29lZmZpY2llbnRzCmBgYAoK