Tarea semana 9
Juan Alberto Zapata
Jun 3, 2020
options(scipen=999)
pkges<-c("mfx","pROC","tidyverse","forecast","data.table","csv")
#install.packages(pkges)
lapply(pkges,"library",character.only=T)## Loading required package: sandwich## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: MASS## Loading required package: betareg## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var## -- Attaching packages ---------------------------------------------------------- tidyverse 1.3.0 --## v ggplot2 3.3.0 v purrr 0.3.4
## v tibble 3.0.1 v dplyr 0.8.5
## v tidyr 1.1.0 v stringr 1.4.0
## v readr 1.3.1 v forcats 0.5.0## -- Conflicts ------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()
## x dplyr::select() masks MASS::select()## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Attaching package: 'data.table'## The following objects are masked from 'package:dplyr':
##
## between, first, last## The following object is masked from 'package:purrr':
##
## transpose## Warning: package 'csv' was built under R version 4.0.2## [[1]]
## [1] "mfx" "betareg" "MASS" "lmtest" "zoo" "sandwich"
## [7] "stats" "graphics" "grDevices" "utils" "datasets" "methods"
## [13] "base"
##
## [[2]]
## [1] "pROC" "mfx" "betareg" "MASS" "lmtest" "zoo"
## [7] "sandwich" "stats" "graphics" "grDevices" "utils" "datasets"
## [13] "methods" "base"
##
## [[3]]
## [1] "forcats" "stringr" "dplyr" "purrr" "readr" "tidyr"
## [7] "tibble" "ggplot2" "tidyverse" "pROC" "mfx" "betareg"
## [13] "MASS" "lmtest" "zoo" "sandwich" "stats" "graphics"
## [19] "grDevices" "utils" "datasets" "methods" "base"
##
## [[4]]
## [1] "forecast" "forcats" "stringr" "dplyr" "purrr" "readr"
## [7] "tidyr" "tibble" "ggplot2" "tidyverse" "pROC" "mfx"
## [13] "betareg" "MASS" "lmtest" "zoo" "sandwich" "stats"
## [19] "graphics" "grDevices" "utils" "datasets" "methods" "base"
##
## [[5]]
## [1] "data.table" "forecast" "forcats" "stringr" "dplyr"
## [6] "purrr" "readr" "tidyr" "tibble" "ggplot2"
## [11] "tidyverse" "pROC" "mfx" "betareg" "MASS"
## [16] "lmtest" "zoo" "sandwich" "stats" "graphics"
## [21] "grDevices" "utils" "datasets" "methods" "base"
##
## [[6]]
## [1] "csv" "data.table" "forecast" "forcats" "stringr"
## [6] "dplyr" "purrr" "readr" "tidyr" "tibble"
## [11] "ggplot2" "tidyverse" "pROC" "mfx" "betareg"
## [16] "MASS" "lmtest" "zoo" "sandwich" "stats"
## [21] "graphics" "grDevices" "utils" "datasets" "methods"
## [26] "base"maindir <- getwd()subdir <- c("/semana9/")datos <- read.csv(paste0(maindir,subdir,"UCI_Credit_Card.csv"),header = T)datos <- datos[complete.cases(datos),]Cambio a variable categorica.
cred <- datos[,-1]
cred$SEX <- as.factor(cred$SEX)
cred$MARRIAGE <- as.factor(cred$MARRIAGE)MODELO 2
Variable dependiente:
_default.payment.next.month.- pago predeterminado (1 = sí, 0 = no)
Variable independiente:
_PAY_0.- Estado de reembolso en septiembre de 2005 (-1 = pagar debidamente, 1 = retraso en el pago de un mes, 2 = retraso en el pago de dos meses, … 8 = retraso en el pago de ocho meses, 9 = retraso en el pago de nueve meses o más)
_BILL_AMT1.- Cantidad de extracto de cuenta en septiembre de 2005 (dólar NT)
_SEX.- Género (1 = masculino, 2 = femenino)
_MARRIAGE.- Estado civil (1 = casado, 2 = soltero, 3 = otros)
_LIMIT_BAL.- Cantidad de crédito otorgado en dólares NT (incluye crédito individual y familiar / complementario)
El modelo 2 se corre con toda las variables con la familia binomial con enlace lobit.
XB <- as.formula("default.payment.next.month ~ PAY_0+BILL_AMT1+factor(SEX)+factor(MARRIAGE)+LIMIT_BAL")
modelo2 <- glm(XB,data = cred,
family = binomial(link = "logit"))
summary(modelo2)##
## Call:
## glm(formula = XB, family = binomial(link = "logit"), data = cred)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.0292 -0.6957 -0.5463 -0.3156 2.8541
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.1962571601 0.5107632550 -4.300 0.0000170835589 ***
## PAY_0 0.7078746785 0.0147991913 47.832 < 0.0000000000000002 ***
## BILL_AMT1 -0.0000016192 0.0000002410 -6.720 0.0000000000182 ***
## factor(SEX)2 -0.1346710975 0.0301793170 -4.462 0.0000081060381 ***
## factor(MARRIAGE)1 1.3154488660 0.5105523184 2.577 0.00998 **
## factor(MARRIAGE)2 1.0898552620 0.5104922796 2.135 0.03277 *
## factor(MARRIAGE)3 1.2116160912 0.5273679837 2.297 0.02159 *
## LIMIT_BAL -0.0000015193 0.0000001389 -10.934 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 31705 on 29999 degrees of freedom
## Residual deviance: 28205 on 29992 degrees of freedom
## AIC: 28221
##
## Number of Fisher Scoring iterations: 4Estimate son las betas de las variables independientes del modelo 2.
Los betas con un p value de significancia de cada covariable tiene una probabilidad menor de 0.05 tienen la menor probabilidad de cometer el error tipo 1.
Estimacion con menor AIC
Te ayuda a escoger la menor cantidad de variables para obtener la menor cantidad de aic.
model.AIC <- stepAIC(modelo2)## Start: AIC=28221
## default.payment.next.month ~ PAY_0 + BILL_AMT1 + factor(SEX) +
## factor(MARRIAGE) + LIMIT_BAL
##
## Df Deviance AIC
## <none> 28205 28221
## - factor(SEX) 1 28225 28239
## - BILL_AMT1 1 28252 28266
## - factor(MARRIAGE) 3 28268 28278
## - LIMIT_BAL 1 28330 28344
## - PAY_0 1 30794 30808yhat1 es la predicción en grafico de histograma.
yhat1<-model.AIC$fitted.values
hist(yhat1)
Sensibilidad y Espesifidad
c<-seq(0.01,0.3,by=0.01)
sens<-c()
spec<-c()
for (i in 1:length(c)){
y.pred<-ifelse(model.AIC$fitted.values > c[i], yes = 1, no = 0)
spec[i]<-prop.table(table(cred$default.payment.next.month,y.pred),1)[1]
sens[i]<-prop.table(table(cred$default.payment.next.month,y.pred),1)[4]
}
o.cut<-mean(c[which(round(spec,1)==round(sens,1))],na.rm = T)
plot(c,sens,type="l",col=2,main=c("Especificidad vs Sensibilidad"),ylab=c("Especificidad/Sensibilidad"))
lines(c,spec,col=3)
abline(v=0.2124)
print(o.cut)## [1] 0.21La especifidad y sencibilidad evaluado entre el punto de corte donde las dos lineas se interceptan es de 0.21.
Matriz de Confusion
y.pred<-ifelse(model.AIC$fitted.values > o.cut, yes = 1, no = 0)
matriz_confusion <- table(cred$default.payment.next.month, y.pred,
dnn = c("observaciones", "predicciones"))Evaluamos en la matriz de confusion, donde si tiende a 0 sensibilidad y 1 especifidad
prop.table(matriz_confusion,1)## predicciones
## observaciones 0 1
## 0 0.6528420 0.3471580
## 1 0.3319771 0.6680229Donde se observa que tiene un 0.652 en sensibilidad y 0.668 especifidad.
Area bajo la curva roc
roc(cred$default.payment.next.month,yhat1)## Setting levels: control = 0, case = 1## Setting direction: controls < cases##
## Call:
## roc.default(response = cred$default.payment.next.month, predictor = yhat1)
##
## Data: yhat1 in 23364 controls (cred$default.payment.next.month 0) < 6636 cases (cred$default.payment.next.month 1).
## Area under the curve: 0.7138plot(roc(cred$default.payment.next.month,yhat1),main=c("Curva ROC"))## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
La función roc muestra el area bajo la curva donde es 0.7138
EFECTO MARGINAL
Muestra la probabilidad de caer en morosidad de las siguientes variables.
logitmfx(formula=XB, data=cred)## Call:
## logitmfx(formula = XB, data = cred)
##
## Marginal Effects:
## dF/dx Std. Err. z
## PAY_0 0.109965239802 0.002243344483 49.0184
## BILL_AMT1 -0.000000251543 0.000000037513 -6.7055
## factor(SEX)2 -0.021101484759 0.004766780170 -4.4268
## factor(MARRIAGE)1 0.212056946435 0.085010733530 2.4945
## factor(MARRIAGE)2 0.166291536624 0.076727341754 2.1673
## factor(MARRIAGE)3 0.250875657635 0.129452729584 1.9380
## LIMIT_BAL -0.000000236021 0.000000021376 -11.0414
## P>|z|
## PAY_0 < 0.00000000000000022 ***
## BILL_AMT1 0.00000000002007 ***
## factor(SEX)2 0.00000956504177 ***
## factor(MARRIAGE)1 0.01261 *
## factor(MARRIAGE)2 0.03021 *
## factor(MARRIAGE)3 0.05263 .
## LIMIT_BAL < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## dF/dx is for discrete change for the following variables:
##
## [1] "factor(SEX)2" "factor(MARRIAGE)1" "factor(MARRIAGE)2"
## [4] "factor(MARRIAGE)3"La variable PAY_0 cambia la probabilidad de 0.1099 de caer en morosidad.
La variable BILL_AMT1 cambia la probabilidad de -0.00000025 de caer en morosidad.
La variable sex cambia la probabilidad de -0.0211 de caer en morosidad.
La variable marriage en 1 cambia la probabilidad de 0.212 de caer en morosidad.
La variable marriage en 2 cambia la probabilidad de 0.16629 de caer en morosidad.
La variable marriage en 3 cambia la probabilidad de 0.2508 de caer en morosidad.
La variable LIMIT_BAL cambia la probabilidad de -0.00000023 de caer en morosidad.
RATIO DE ODDG
PAY_0: Estado de reembolso en septiembre de 2005 (-1 = pagar debidamente,
1 = retraso en el pago de un mes,
2 = retraso en el pago de dos meses,
8 = retraso en el pago de ocho meses,
9 = retraso en el pago de nueve meses o más)
ratio_oddg1 <- exp(0.7078746785)-1
ratio_oddg1## [1] 1.029673probabilidad1 <- exp(0.7078746785)/(1+exp(0.7078746785))
probabilidad1## [1] 0.6699314La chance de cambiar en morosidad es de 1.029673 La probabilidad de caer en morosidad es de 0.6699314
BILL_AMT1: Cantidad de extracto de cuenta en septiembre de 2005 (dólar NT)
ratio_oddg2 <- exp(-0.0000016192)-1
ratio_oddg2## [1] -0.000001619199probabilidad2 <- exp(-0.0000016192)/(1+exp(-0.0000016192))
probabilidad2## [1] 0.4999996La chance de BILL_AMT1 de cambiar en morosidad es de -0.0000016 La probabilidad de BILL_AMT1 de caer en morosidad es de 0.49999
SEXO: Género (1 = masculino, 2 = femenino)
ratio_oddg3 <- exp(-0.1346710975)-1
ratio_oddg3## [1] -0.1259967probabilidad3 <- exp(-0.1346710975)/(1+exp(-0.1346710975))
probabilidad3## [1] 0.466383La chance del Género de cambiar en morosidad es de -0.12599 La probabilidad del Género de caer en morosidad es de 0.46638
Estado civil (1 = casado, 2 = soltero, 3 = otros)
ratio_oddg41 <- exp(1.3154488660)-1
ratio_oddg41## [1] 2.726423probabilidad41 <- exp(1.3154488660)/(1+exp(1.3154488660))
probabilidad41## [1] 0.7884235La chance de un casado de cambiar en morosidad es de 2.7264 La probabilidad de un casado de caer en morosidad es de 0.7884
ratio_oddg42 <- exp(1.0898552620)-1
ratio_oddg42## [1] 1.973844probabilidad42 <- exp(1.0898552620)/(1+exp(1.0898552620))
probabilidad42## [1] 0.7483545La chance de un soltero de cambiar en morosidad es de 1.9738 La probabilidad de un soltero de caer en morosidad es de 0.7483
ratio_oddg43 <- exp(1.2116160912)-1
ratio_oddg43## [1] 2.358909probabilidad43 <- exp(1.2116160912)/(1+exp(1.2116160912))
probabilidad43## [1] 0.7705848La chance de otros de cambiar en morosidad es de 2.3589 La probabilidad de otros de caer en morosidad es de 0.7705
LIMIT_BAL: Cantidad de crédito otorgado en dólares NT (incluye crédito individual y familiar / complementario
ratio_oddg5 <- exp(-0.0000015193)-1
ratio_oddg5## [1] -0.000001519299probabilidad5 <- exp(-0.0000015193)/(1+exp(-0.0000015193))
probabilidad5## [1] 0.4999996La chance de la Linea de credito de cambiar en morosidad es de -0.000001519 La probabilidad de la Linea de credito de caer en morosidad es de 0.499999
Tabla de datos
id.Modelo2 <- data.frame(
Nomb= c("Estado_2005","BILL_AMT1","SEXO","Estado_civil_1","Estado_civil_2","Estado_civil_3","Cantidad_crédito"),
Estimate = c(0.7078746785,-0.0000016192,-0.1346710975,1.3154488660,1.0898552620,1.2116160912,-0.0000015193),
dF_dx = c(0.109965239802,-0.000000251543,-0.021101484759,0.212056946435,0.166291536624,0.250875657635,-0.000000236021),
RATIO_ODDG = c(ratio_oddg1,ratio_oddg2,ratio_oddg3,ratio_oddg41,ratio_oddg42,ratio_oddg43,ratio_oddg5),
PROBABILIDAD = c(probabilidad1,probabilidad2,probabilidad3,probabilidad41,probabilidad42,probabilidad43,probabilidad5),
stringsAsFactors = FALSE
)Imprimir el marco de datos.
print(id.Modelo2) ## Nomb Estimate dF_dx RATIO_ODDG PROBABILIDAD
## 1 Estado_2005 0.7078746785 0.109965239802 1.029672963802 0.6699314
## 2 BILL_AMT1 -0.0000016192 -0.000000251543 -0.000001619199 0.4999996
## 3 SEXO -0.1346710975 -0.021101484759 -0.125996673701 0.4663830
## 4 Estado_civil_1 1.3154488660 0.212056946435 2.726423274616 0.7884235
## 5 Estado_civil_2 1.0898552620 0.166291536624 1.973843613235 0.7483545
## 6 Estado_civil_3 1.2116160912 0.250875657635 2.358908571178 0.7705848
## 7 Cantidad_crédito -0.0000015193 -0.000000236021 -0.000001519299 0.4999996