Práctica 3: Análisis de Supervivencia
Aguirre Armada Guillermo
Colin Betanzos Angel Iván
Estrada Garcia Naúhn Gabriel
Figueroa Torres Ivan Emiliano
Maya García Daniel
Salazar Cedillo Hugo
Sandoval Mendoza Jorge
5/12/2020
library(readr)
library(survival)
library(tidyverse)
library(latex2exp)
library(ggplot2)
library(gridExtra)
library(kableExtra)
library(knitr)
library(dplyr)
library(ggplot2)
library(km.ci)
library(ggpubr)
library(goftest)
library(univariateML)Introducción
Para esta práctica trabajaremos con datos de marketing para prestamos bancarios.
Los datos están relacionados con campañas de marketing directo de una institución bancaria portuguesa. Las campañas de marketing se basaron en llamadas telefónicas. A menudo, se requería más de un contacto con el mismo cliente, para poder acceder si el producto (depósito a plazo bancario) estaría (‘sí’) o no (‘no’) suscrito.
Se cuentan con los datos de \(45,211\) personas, agrupados en las siguientes \(17\) variables:
- age: Edad del cliente.
- job: Tipo de empleo (categorica: administrativo, obrero, emprendedor, empleada doméstica, administración, jubilado, autónomo, servicios, estudiante, técnico, desempleado, desconocido).
- marital: Estado Civil (categorica: divorciado, casado, soltero, desconocido. Divorciado engloba divorcio y viudez).
- education: Nivel educativo (categorica: primaria, secundaria, curso profesional/título universitario, desconocido).
- default: ¿Tiene crédito en mora? (categorica: No, Si, desconocido).
- balance: Balance del cliente.
- housing: ¿Tiene préstamo para vivienda? (categorica: No, Si, desconocido).
- loan: ¿Tiene préstamo personal? (categorica: No, Si, desconocido).
- contact: Tipo de comunicacion para contactar (categorica: celular, telefono).
- day: Día del mes en que se hizo la llamada.
- month: mes en el que se hace el ultimo contacto (categorica: enero,…, diciembre).
- duration: Duración en segundos del ultimo contacto.
- campaign: Número de contactos realizados durante esta campaña y para este cliente(incluye el último).
- pdays: Número de días que pasaron desde la última vez que se contactó con el cliente desde una campaña anterior (-1 significa que no fue contactado previamente)
- previous: Número de contactos realizados antes de esta campaña y para este cliente.
- poutcome: Resultado de la campaña de marketing anterior (categorica: fracaso, inexistente, éxito)
- y: ¿El cliente ha suscrito un depósito a plazo? (categorico: sí, no)
Se desea estudiar el tiempo \((T)\) de la llamada y la variable de fallo será y (si el cliente suscrito realizará el depositado o no).
Análisis Exploratorio de Datos
El objetivo inicial es conocer los datos y entender su comportamiento. Dividimos las variables en numéricas y categóricas,
Variables numéricas
Usaremos histogramas para ver dos cosas, como se distribuyen las variables y que proporción de la variable \(y\) hay dentro de cada una de estas variables.
Posteriormente mostraremos un resumen estadistico para estas variables.
d <- read_csv("bank-full.csv")
d<-data.frame(d)
p1 <- ggplot(data = d) +
geom_histogram(aes(x=age,group=y,fill=y), binwidth=1,color='black') +
labs(fill="Deposito")+
labs(title="Edad",x = "Edad", y = "Frecuencia") +
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))
p2<-ggplot(data = d) +
geom_histogram(aes(x=balance,group=y,fill=y), bins=100,color='black') +
labs(fill="Deposito")+
labs(title="Balance",x = "Balance", y = "Frecuencia") +
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))
p3<-ggplot(data = d) +
geom_histogram(aes(x=duration,group=y,fill=y), bins=100,color='black') +
labs(fill="Deposito")+
labs(title="Duracion",x = "Duracion", y = "Frecuencia") +
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))
p4<-ggplot(data = d) +
geom_histogram(aes(x=campaign,group=y,fill=y), bins=100,color='black') +
labs(fill="Deposito")+
labs(title="Numero de contactos",x = "numero", y = "Frecuencia") +
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))
p5<-ggplot(data = d) +
geom_histogram(aes(x=pdays,group=y,fill=y), bins=20,color='black') +
labs(fill="Deposito")+
labs(title="Dias entre campañas",x = "dias", y = "Frecuencia") +
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))
p6<-ggplot(data = d) +
geom_histogram(aes(x=previous,group=y,fill=y), bins=50,color='black') +
labs(fill="Deposito")+
labs(title="Contactos campaña anterior",x = "contactos", y = "Frecuencia") +
theme_bw()+
theme(plot.title = element_text(hjust = 0.5))
p1;p2;p3;p4;p5;p6
#ggarrange(p1,p2,p3,p4,p5,p6,ncol=2,nrow=3, common.legend = TRUE, legend = "bottom")Como se puede ver en las graficas, a grandes rasgos siempre hay una mayor cantidad de clientes que no realizan el depósito, adicionalmente se observa que:
Edades: Podemos observar que la mayoría de los clientes tienen entre 25 y 50 años, siendo más específicos, hay más clientes que tienen 32 años.
Balance: De la gráfica, podemos observar que la mayor parte de los clientes tienen un balance 0, por lo que podemos decir que la mayoría de los clientes van al corriente con sus deudas, también podemos ver que hay pocos clientes a los que el banco les debe, pero son menos que los clientes que le deben al banco.
Duración: Vemos que la duración en segundos del último contacto se concentra en los primeros 1000 segundos, es decir, en menos de 17 minutos el cliente toma su decision con respecto a hacer o no el depósito.
Número de contactos: Podemos ver que la mayor parte de los clientes tienen entre 1 y 5 contactos, y muy pocos o prácticamente nadie tiene más de 70 contratos, es en las primeras llamadas donde se decide si se hace el depósito o no, posterior a esto es dificil hacer que el cliente lo haga.
Días entre campañas: Podemos observar que para la mayoria de los clientes, esta es su primer campaña pues existe una gran acumulacion de observaciones al rededor del 0 y el -1.
Contactos anteriores: De igual manera, como esta es la primer campaña de muchos clientes(la gran mayoria), el numero de contactos anteriores es \(0\) pues no hay registro de haber llamado a dichos clientes antes
options(digits = 4)
kable(data.frame(c("Edad", "Balance", "Duración", "N° de contactos", "Dias entre campañas", "Contactos campaña anterior"),
c(min(d$age),min(d$balance),min(d$duration),min(d$campaign),min(d$pdays),min(d$previous)),
c(quantile(d$age, prob=0.25),quantile(d$balance, prob=0.25),quantile(d$duration, prob=0.25),quantile(d$campaign, prob=0.25),quantile(d$pdays, prob=0.25),quantile(d$previous, prob=0.25)),
c(quantile(d$age, prob=0.5),quantile(d$balance, prob=0.5),quantile(d$duration, prob=0.5),quantile(d$campaign, prob=0.5),quantile(d$pdays, prob=0.5),quantile(d$previous, prob=0.5)),
c(mean(d$age),mean(d$balance),mean(d$duration),mean(d$campaign),mean(d$pdays),mean(d$previous)),
c(quantile(d$age, prob=0.75),quantile(d$balance, prob=0.75),quantile(d$duration, prob=0.75),quantile(d$campaign, prob=0.75), quantile(d$pdays, prob=0.75),quantile(d$previous, prob=0.75)),
c(max(d$age),max(d$balance),max(d$duration),max(d$campaign),max(d$pdays),max(d$previous)),
c(var(d$age),var(d$balance),var(d$duration),var(d$campaign),var(d$pdays),var(d$previous)),
c(sqrt(var(d$age)),sqrt(var(d$balance)),sqrt(var(d$duration)),sqrt(var(d$campaign)),sqrt(var(d$pdays)),sqrt(var(d$previous)))),
caption = "Resumen Estadistico"
, align = c('c', 'c', 'c','c', 'c', 'c','c', 'c', 'c'), col.names=c("Variable","Mínimo","25%","50%","Media","75%","Maximo","Varianza","Desv. Est")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Variable | Mínimo | 25% | 50% | Media | 75% | Maximo | Varianza | Desv. Est |
|---|---|---|---|---|---|---|---|---|
| Edad | 18 | 33 | 39 | 40.9362 | 48 | 95 | 1.128e+02 | 10.619 |
| Balance | -8019 | 72 | 448 | 1362.2721 | 1428 | 102127 | 9.271e+06 | 3044.766 |
| Duración | 0 | 103 | 180 | 258.1631 | 319 | 4918 | 6.632e+04 | 257.528 |
| N° de contactos | 1 | 1 | 2 | 2.7638 | 3 | 63 | 9.598e+00 | 3.098 |
| Dias entre campañas | -1 | -1 | -1 | 40.1978 | -1 | 871 | 1.003e+04 | 100.129 |
| Contactos campaña anterior | 0 | 0 | 0 | 0.5803 | 0 | 275 | 5.306e+00 | 2.303 |
Lo anteriormente descrito queda mejor evidenciado con la tabla anterior.
Variables Categóricas
Para estas variables mostramos tablas de frecuencias para observar la composicion de cada variable de acuerdo a sus categorias.
kable(data.frame(c("administrativo", "obrero", "emprendedor", "empleada doméstica", "administración", "jubilado", "autónomo", "servicios", "estudiante", "técnico", "desempleado", "desconocido"),
c(5171, 9732, 1487, 1240, 9458, 2264, 1579, 4154, 938, 7597, 1303, 288),
c(0.11,0.22,0.03,0.03,0.21,0.05,0.03,0.09,0.02,0.17,0.03,0.01)),
caption = "Tipo de Empleo"
, align = c('c', 'c', 'c'), col.names=c("Empleo","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Empleo | Cantidad de Clientes | Proporcion |
|---|---|---|
| administrativo | 5171 | 0.11 |
| obrero | 9732 | 0.22 |
| emprendedor | 1487 | 0.03 |
| empleada doméstica | 1240 | 0.03 |
| administración | 9458 | 0.21 |
| jubilado | 2264 | 0.05 |
| autónomo | 1579 | 0.03 |
| servicios | 4154 | 0.09 |
| estudiante | 938 | 0.02 |
| técnico | 7597 | 0.17 |
| desempleado | 1303 | 0.03 |
| desconocido | 288 | 0.01 |
En esta tabla se puede observar que la mayor cantidad de clientes trabajan en los sectores obrero (22%) y lo relacionado con la administración (32%), y los que trabajan en sectores desconocidos o poco remunerados cómo estudiante (2%), emprendedor (3%), empleada doméstica (3%), autónomo (3%), desempleado (3%) y desconocido (1%).
Lo cuál tiene sentido debido a que las personas que tienen mayores ingresos tienen una capacidad más amplia de poder adquirir los productos ofrecidos y poder liquidar sus deudas.
kable(data.frame(c("Divorciado", "Casado", "Soltero"),
c(5207, 27214, 12790),
c(0.1152, 0.6019, 0.2829)),
caption = "Estado Civil"
, align = c('c', 'c', 'c'), col.names=c("Estado","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Estado | Cantidad de Clientes | Proporcion |
|---|---|---|
| Divorciado | 5207 | 0.1152 |
| Casado | 27214 | 0.6019 |
| Soltero | 12790 | 0.2829 |
Se puede observar un comportamiento bastante importante referente al estado civil, siendo los que están casados la mayor cantidad de clientes que se tienen esto podría ser debido a que perciben ingresos por parte de las 2 personas de las parejas y a su vez, las deudas tambien se comparten y aumentan, seguidos por las personas solteras que tienen menor cantidad de ingresos al ser solo una persona y por ultimo los divorciados.
kable(data.frame(c("Primaria", "Secundaria", "Universidad", "Desconocido"),
c(6851,23202,13301,1857),
c(0.15153,0.51319,0.29420,0.04107)),
caption = "Nivel educativo"
, align = c('c', 'c', 'c'), col.names=c("Nivel","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Nivel | Cantidad de Clientes | Proporcion |
|---|---|---|
| Primaria | 6851 | 0.1515 |
| Secundaria | 23202 | 0.5132 |
| Universidad | 13301 | 0.2942 |
| Desconocido | 1857 | 0.0411 |
En la tabla podemos observar que la mayoría de los clientes tienen un nivel educativo de secundaria, seguido por los de universidad, los de primaria y por ultimo los que no sé sabe su nivel educativo, esto puede estar relacionado con la educación de todo el país y que la mayoría de las personas terminan la secundaria, pero ya no estudian la universidad.
kable(data.frame(c("No", "Si"),
c(44396, 815),
c(0.98, 0.02)),
caption = "Crédito en mora"
, align = c('c', 'c', 'c'), col.names=c("Situación","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Situación | Cantidad de Clientes | Proporcion |
|---|---|---|
| No | 44396 | 0.98 |
| Si | 815 | 0.02 |
La institución tiene un nivel muy bajo en cuanto a los clientes que están en situación de mora, la mayoría de sus clientes (98%) pagan por los productos que adquieren lo cual es muy bueno para la institución.
kable(data.frame(c("No", "Si"),
c(20081, 25130 ),
c(0.44, 0.56)),
caption = "Prestamo para vivienda"
, align = c('c', 'c', 'c'), col.names=c("Situación","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Situación | Cantidad de Clientes | Proporcion |
|---|---|---|
| No | 20081 | 0.44 |
| Si | 25130 | 0.56 |
El 56% de los clientes tienen préstamos para viviendas lo que nos quiere decir que una gran cantidad de nuestros clientes buscan adquirir patrimonio propio, es decir, adquirir una casa o un apartamento.
kable(data.frame(c("No", "Si"),
c(37967, 7244 ),
c(0.84, 0.16)),
caption = "Prestamo personal"
, align = c('c', 'c', 'c'), col.names=c("Situación","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Situación | Cantidad de Clientes | Proporcion |
|---|---|---|
| No | 37967 | 0.84 |
| Si | 7244 | 0.16 |
La mayoría de los clientes (84%) no adquieren préstamos personales y sólo una mínima parte (16%) de los clientes si los adquieren, podría deberse al hecho de que les interesa más tener una casa, cómo vimos anteriormente a obtener un préstamo personal.
kable(data.frame(c("Celular", "Telefono", "Desconocido"),
c(29285,2906,13020),
c(0.65,0.06,0.29)),
caption = "Método de contacto"
, align = c('c', 'c', 'c'), col.names=c("Situación","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Situación | Cantidad de Clientes | Proporcion |
|---|---|---|
| Celular | 29285 | 0.65 |
| Telefono | 2906 | 0.06 |
| Desconocido | 13020 | 0.29 |
En esta tabla es curioso el hecho de que el (29%) de los clientes tienen un método de contacto desconocido, ya que debería haber un control para poder contactar a los clientes en caso de que no paguen y la mayor parte da como contacto su número de celular, esto debido a que ahora hay mayor cantidad de personas con celular y es mucho más fácil contactar a las personas por medio de este método.
kable(data.frame(c("Fallo","Otro", "Éxito","Desconocido"),
c(4901,1840,1511,36959),
c(0.11,0.04,0.03,0.82)),
caption = "Resultado de la campaña anterior"
, align = c('c', 'c', 'c'), col.names=c("Situación","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Situación | Cantidad de Clientes | Proporcion |
|---|---|---|
| Fallo | 4901 | 0.11 |
| Otro | 1840 | 0.04 |
| Éxito | 1511 | 0.03 |
| Desconocido | 36959 | 0.82 |
Podemos ver que las campañas tienen una tasa de éxito del 3% mientras que del 82% no se tiene conocimiento y el 11% como tasa de fallo. Se podría mejorar las campañas o hacerlas más intensivas para incrementar el nivel de éxito y conseguir que los clientes hagan sus depósitos.
kable(data.frame(c("No","Si"),
c(39922, 5289),
c(0.88,0.12)),
caption = "Suscripción al depósito"
, align = c('c', 'c', 'c'), col.names=c("Situación","Cantidad de Clientes","Proporcion")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Situación | Cantidad de Clientes | Proporcion |
|---|---|---|
| No | 39922 | 0.88 |
| Si | 5289 | 0.12 |
Podemos notar que sólo el 12% de los clientes se suscribe al depósito, mientras que el 88% de los clientes no se suscribe. Esto indica una mejoria significativa con respecto a la campaña anterior, la institución debe seguir haciendo lo mismo para que los clientes realicen sus depositos más rapido.
Función de Supervivencia
El método K-M se basa en el uso de la Log-Verosimilitud para estimar los parametros que permitan calcular la función de supervivencia.
Recordemos que para una muestra aleatoria \(\{T_i\}_{i=1}^n\) con soporte discreto en \(\{u_1,u_2,...\}\) la función de supervivencia usando K-M se estima como sigue:
\[\hat S(t) = \prod_{k:u_k\leq t}\left(1-\frac{d_k}{n_k}\right)\] Donde:
\[d_k= \text{número de tiempos de fallo iguales a } u_k \\ n_k= \text{número de individuos en riesgo al tiempo } u_k\] Recordemos que al usar el método Kaplan-Meier, el intervalo puntual al \((1-\alpha)\cdot 100\%\) de confianza para \(\hat S(t_0)\) se estiman como sigue:
\[ConfInt=\left( \hat S(t_0)-z_{1-\frac{\alpha}{2}}\sqrt{\hat{Var}(\hat S(t_0))},~~\hat S(t_0)+z_{1-\frac{\alpha}{2}}\sqrt{\hat{Var}(\hat S(t_0))}\right)\] Donde \(\alpha=0.05\) y \(z_{1-\frac{\alpha}{2}}\) es el cuantil \(1-\frac{\alpha}{2}\) de una distribución Normal estandar.
Se muestra una tabla con los 10 primeros renglones del cálculo de \(\hat S(t)\):
Donde para este problema se consideró:
\[ \text{Fallo: El cliente no realizó el depósito}\\ \text{Supervivencia: El cliente realizó el depósito} \]
d1<-Surv(d$duration, d$Status)
d2<- survfit(d1~1)
#tabla de Supervivencia
a<-tibble(Tiempo = d2$time,
dk = d2$n.event,
nk = d2$n.risk,
ck = d2$n.censor,
surv = d2$surv)
#tabla con formato, muestra los primeros valores de la tabla anterior
kable(head(a) , caption = "Función de Supervivencia"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Tiempo","dk","nk","ck","S(t)")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Tiempo | dk | nk | ck | S(t) |
|---|---|---|---|---|
| 0 | 3 | 45211 | 0 | 0.9999 |
| 1 | 2 | 45208 | 0 | 0.9999 |
| 2 | 3 | 45206 | 0 | 0.9998 |
| 3 | 4 | 45203 | 0 | 0.9997 |
| 4 | 15 | 45199 | 0 | 0.9994 |
| 5 | 35 | 45184 | 0 | 0.9986 |
Gráfica
Gráfica de supervivencia estimada por Kaplan-Meier junto a sus intervalos de confianza.
#añadimos los intervalos de confianza
a<-cbind(a,VarS = a$surv^2*cumsum(a$dk/((a$nk-a$dk)*a$nk)))
a<-cbind(a, LimInf = a$surv-qnorm(0.975)*sqrt(a$VarS),LimSup = a$surv+qnorm(0.975)*sqrt(a$VarS))
a$LimInf = if_else(is.na(a$LimInf)|a$LimInf<0, 0, a$LimInf)
a$LimSup = if_else(a$LimSup>1, 1, a$LimSup)
a$LimSup = if_else(is.na(a$LimSup), 0, a$LimSup)
#grafica de S(t) ahora incluyendo los intervalos de confianza
ggplot(data=a, aes(x=Tiempo,y=surv))+
geom_step(aes(color="S(t)"),size=1)+
geom_step(aes(y=LimInf, color="IntConf"),alpha=0.4,size=0.8)+
geom_step(aes(y=LimSup, color="IntConf"),alpha=0.4,size=0.8)+
theme_bw()+
labs(title="Funcion de Supervivencia", x="Tiempo", y=TeX("$\\hat{S}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_color_manual(breaks = c("S(t)", "IntConf"),
values = c("S(t)" = "steelblue","IntConf" = "red")) +
labs(color = "Función") +
theme(legend.position = "top")
Se observa que la gráfica muestra un decrecimiento acelerado y los intervalos de confianza se ajustan bastante a la función. Más adelante veremos mas a detalle las caracteristicas de la supervivencia
Resúmenes Estadisticos
Parámetros Poblacionales
Media
Para el estimador \(\hat S(t)\) se puede calcular la media de los tiempos de falla como: \[\hat \mu = \int_0^\infty\hat S(t) dt\] En este caso no contamos con una función explicita que defina a \(\hat S(t)\) por lo que discretizaremos la integral de la siguiente manera:
\[\hat \mu = \int_0^\infty\hat S(t) dt=\sum_{i=1}^k \hat S(t_i)\cdot(t_{i+1}-t_i)\] Para este estimador, el intervalo al \((1-\alpha)\cdot 100\%\) de confianza es:
\[\hat \mu ~~\pm~~z_{1-\frac{\alpha}{2}}\sqrt{\hat{Var}(\hat \mu)}\] Donde(para esta discretización):
\[\hat{Var}(\hat \mu)=\sum_{i=1}^k\left[\sum_{j=i}^k \hat S(t_j)\cdot(t_{j+1}-t_j)\cdot \right]^2\frac{d_i}{n_i(n_i-d_i)}\] La siguiente tabla muestra los resultados de los cálculos:
Media<-function(datos){
#calculamos la integral como la suma de los rectangulos que forman la funcion de supervivencia
aux<-(lead(datos$time)-datos$time)*datos$surv
#eliminamos posibles NA's
aux<-na.omit(aux)
#la media es la suma de los rectangulos
media<-sum(aux)
#Varianza
integral<-seq(1:length(aux))
length(integral)
for(i in 1:length(aux)){
integral[i]=sum(aux[i:length(aux)])
}
integral<-integral^2
cociente<-(datos$n.event/((datos$n.risk-datos$n.event)*datos$n.risk))
length(cociente)
length(integral)
#la varianza es el producto punto de la integral al cuadrado por el cociente
varianza<-integral%*%cociente[1:length(integral)]
varianza<-varianza[1]
#Intervalo de confianza para la media
LimInf <- media - qnorm(0.975)*sqrt(varianza)
LimSup <- media + qnorm(0.975)*sqrt(varianza)
return(c(media,varianza,LimInf[1],LimSup[1]))
}
m<-Media(d2)
kable(data.frame(m[1],m[2],m[3],m[4]) , caption = "Media de los tiempos de fallo"
, align = c('c', 'c', 'c', 'c'), col.names=c("Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|
| 332.9 | 21.52 | 323.8 | 341.9 |
Entonces la media de duración de las llamadas es de \(332.9\) segundos, casi \(6\) minutos, dentro de un intervalo de confianza del \(97.5\%\) cuya longitud es de \(18\) segundos.
Este resultado es bastante similar a la media calculada directamente sobre los tiempos sin considerar la supervivencia, la cual es: \(258.16\).
Pero es más confiable el nuevo valor pues en su calculo viene implicita la supervivencia.
Cuantiles
A partir del estimador \(\hat S(t)\) se calculan los cuantiles de orden \(p\) como sigue:
\[\hat t_p= inf \{t:\hat S(t) \leq 1-p\}\] La siguiente tabla muestra algunos cuantiles:
#funcion que calcula y devuelve los cuantiles 25%, 50% y 75%
Cuantiles<-function(datos){
porcentaje <- c("25%", "50%", "75%")
#se calcula el valor de algunos cuantiles
cuantiles <- c(min(subset(datos$time,datos$surv<=1-0.25)),
min(subset(datos$time,datos$surv<=1-0.5)),
min(subset(datos$time,datos$surv<=1-0.75)))
return(list(p=porcentaje,c=cuantiles))
}
c<-Cuantiles(d2)
kable(data.frame(c$p, c$c) , caption = "Cuantiles relevantes"
, align = c('c', 'c'), col.names=c("Porcentaje","Cuantil")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | Cuantil |
|---|---|
| 25% | 104 |
| 50% | 185 |
| 75% | 345 |
25% 104 50% 185 75% 345 Como se puede observar, la primer cuarta parte de los contactos/llamadas duran \(104\) segundos(casi \(2\) minutos), la primera mitad de las llamadas duran \(185\) segundos(\(3\) minutos) y la ultima cuarta parte logra durar más de \(345\) segundos(\(5\) minutos), pese a que la duracion media de las llamadas sea de \(332.9\) segundos(casi \(6\) minutos) menos de la mitad logra llegar a este tiempo.
Sumado a lo anterior, dado que \(\hat t_{0.5}<\hat \mu\), quiere decir que la supervivencia está sesgada hacia la derecha.
Funcion de Riesgo:
Recordemos que la función de Riesgo para el estimador K-M es:
\[\hat h(u_k)=\frac{d_k}{n_k}\] La gráfica para \(\hat h(t)\) es la siguiente:
ht<-d2$n.event/d2$n.risk
#grafica de h(t)
ggplot() +
geom_line(aes(x=d2$time,y=ht), size=1,col="steelblue")+
theme_bw() +
labs(title="Funcion de Riesgo", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))
Ignorando los picos se puede observar que la función mantiene un comportamiento creciente, disparandose pasados los \(4,000\) segundos, pues pasado tanto tiempo es mas dificil de convencer al cliente para que realice el depósito.
Función de Riesgo Acumulado:
Para obtener la funcion de riesgo acumulado se usa el estimador de Nelson-Aelen el cuál es el siguiente:
\[\hat H_2(t)= \sum_{j:u_j\leq t}\hat h_j = \sum_{j:u_j\leq t}\frac{d_j}{n_j}\] El intervalo de confianza puntual al \((1-\alpha)\cdot 100\%\) para este estimador es:
\[\hat H_2(t_0)\pm z_{1-\frac{\alpha}{2}}\frac{\sqrt{\hat{Var}(\hat H_2(t_0))}}{\hat H_2(t_0)}\] Con \[\hat{Var}(\hat H_2(t))=\sum_{j:u_j\leq t}\frac{d_j}{n_j^2}\]
Tabla de los primeros valores de \(\hat H_2(t)\) incluyendo sus intervalos de confianza:
d3 <- summary(survfit(Surv(d$duration, d$Status)~1))
#se genera tabla con informacion importante sobre H(t) y sus intervalos de confianza
c<-tibble(dk = d3$n.event,
nk = d3$n.risk,
H = d3$cumhaz,
stdH = d3$std.chaz) %>%
mutate(upper = H + qnorm(0.975)*stdH/H,
lower = H - qnorm(0.975)*stdH/H)
colnames(c)<-c("dk", "nk","H(t)","Desv. Est de H(t)","Lim Inf","Lim Sup")
kable(head(c) , caption = "Función de Riesgo acumulado"
, align = c('c', 'c', 'c', 'c','c')) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| dk | nk | H(t) | Desv. Est de H(t) | Lim Inf | Lim Sup |
|---|---|---|---|---|---|
| 3 | 45211 | 0.0001 | 0e+00 | 1.1317 | -1.1315 |
| 2 | 45208 | 0.0001 | 0e+00 | 0.8766 | -0.8764 |
| 3 | 45206 | 0.0002 | 1e-04 | 0.6931 | -0.6928 |
| 4 | 45203 | 0.0003 | 1e-04 | 0.5661 | -0.5655 |
| 15 | 45199 | 0.0006 | 1e-04 | 0.3778 | -0.3766 |
| 35 | 45184 | 0.0014 | 2e-04 | 0.2503 | -0.2475 |
Gráfica de Riesgo acumulado:
#calculamos los intervalos de confianza
inf<- d2$cumhaz-qnorm(0.975)*d2$std.chaz/d2$cumhaz
sup<-d2$cumhaz+qnorm(0.975)*d2$std.chaz/d2$cumhaz
#grafica
ggplot()+
geom_step(aes(x=d2$time,y=d2$cumhaz, color="H(t)"),size=1)+
geom_step(aes(x=d2$time,y=inf, color="IntConf"),alpha=0.4,size=0.8)+
geom_step(aes(x=d2$time,y=sup, color="IntConf"),alpha=0.4,size=0.8)+
theme_bw()+
labs(title="Funcion de Riesgo Acumulada", x="Tiempo", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_color_manual(breaks = c("H(t)", "IntConf"),
values = c("H(t)" = "steelblue","IntConf" = "red")) +
labs(color = "Función") +
theme(legend.position = "right")
Segmentación de los Datos
A continuación se realizará el análisis de supervivencia sobre cada variable categórica.
Empleo
Segmentamos las supervivencia de acuerdo al tipo de empleo.
Tipo de empleo:
- Administrativo
- Obrero
- Emprendedor
- Empleada doméstica
- Administración
- Jubilado
- Autónomo
- Servicios
- Estudiante
- Técnico
- Desempleado
- Desconocido
Se muestran las gráficas de supervivencia de cada tipo.
#separamos los datos de acuerdo a los tipos de empleo
{
j1=subset(d,job=="management", select = c(duration,Status))
j2=subset(d,job=="technician", select = c(duration,Status))
j3=subset(d,job=="entrepreneur", select = c(duration,Status))
j4=subset(d,job=="blue-collar", select = c(duration,Status))
j5=subset(d,job=="unknown", select = c(duration,Status))
j6=subset(d,job=="retired", select = c(duration,Status))
j7=subset(d,job=="admin.", select = c(duration,Status))
j8=subset(d,job=="services", select = c(duration,Status))
j9=subset(d,job=="self-employed", select = c(duration,Status))
j10=subset(d,job=="unemployed", select = c(duration,Status))
j11=subset(d,job=="housemaid", select = c(duration,Status))
j12=subset(d,job=="student", select = c(duration,Status))
}
#Supervivencias
{Surv_j1_fit <- survfit(Surv(time = j1$duration, event = j1$Status)~1)
a1 <-data.frame(Tiempo = Surv_j1_fit$time, Surv = Surv_j1_fit$surv)
Surv_j2_fit <- survfit(Surv(time = j2$duration, event = j2$Status)~1)
a2 <-data.frame(Tiempo = Surv_j2_fit$time, Surv = Surv_j2_fit$surv)
Surv_j3_fit <- survfit(Surv(time = j3$duration, event = j3$Status)~1)
a3 <-data.frame(Tiempo = Surv_j3_fit$time, Surv = Surv_j3_fit$surv)
Surv_j4_fit <- survfit(Surv(time = j4$duration, event = j4$Status)~1)
a4 <-data.frame(Tiempo = Surv_j4_fit$time, Surv = Surv_j4_fit$surv)
Surv_j5_fit <- survfit(Surv(time = j5$duration, event = j5$Status)~1)
a5 <-data.frame(Tiempo = Surv_j5_fit$time, Surv = Surv_j5_fit$surv)
Surv_j6_fit <- survfit(Surv(time = j6$duration, event = j6$Status)~1)
a6 <-data.frame(Tiempo = Surv_j6_fit$time, Surv = Surv_j6_fit$surv)
Surv_j7_fit <- survfit(Surv(time = j7$duration, event = j7$Status)~1)
a7 <-data.frame(Tiempo = Surv_j7_fit$time, Surv = Surv_j7_fit$surv)
Surv_j8_fit <- survfit(Surv(time = j8$duration, event = j8$Status)~1)
a8 <-data.frame(Tiempo = Surv_j8_fit$time, Surv = Surv_j8_fit$surv)
Surv_j9_fit <- survfit(Surv(time = j9$duration, event = j9$Status)~1)
a9 <-data.frame(Tiempo = Surv_j9_fit$time, Surv = Surv_j9_fit$surv)
Surv_j10_fit <- survfit(Surv(time = j10$duration, event = j10$Status)~1)
a10 <-data.frame(Tiempo = Surv_j10_fit$time, Surv = Surv_j10_fit$surv)
Surv_j11_fit <- survfit(Surv(time = j11$duration, event = j11$Status)~1)
a11 <-data.frame(Tiempo = Surv_j11_fit$time, Surv = Surv_j11_fit$surv)
Surv_j12_fit <- survfit(Surv(time = j12$duration, event = j12$Status)~1)
a12 <-data.frame(Tiempo = Surv_j12_fit$time, Surv = Surv_j12_fit$surv)
}
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="administracion"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="tecnico"))+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="emprendedor"))+
geom_step(mapping = aes(x=a4$Tiempo, y=a4$Surv, color="obrero"))+
geom_step(mapping = aes(x=a5$Tiempo, y=a5$Surv, color="desconocido"))+
geom_step(mapping = aes(x=a6$Tiempo, y=a6$Surv, color="retirado"))+
geom_step(mapping = aes(x=a7$Tiempo, y=a7$Surv, color="administrativo"))+
geom_step(mapping = aes(x=a8$Tiempo, y=a8$Surv, color="servicios"))+
geom_step(mapping = aes(x=a9$Tiempo, y=a9$Surv, color="independiente"))+
geom_step(mapping = aes(x=a10$Tiempo, y=a10$Surv, color="desempleado"))+
geom_step(mapping = aes(x=a11$Tiempo, y=a11$Surv, color="empleado domestico"))+
geom_step(mapping = aes(x=a12$Tiempo, y=a12$Surv, color="estudiante"))+
labs(title="Supervivencia por tipo de empleo", x="Edad", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Podemos observar que la mayoría de las gráficas parecen compartarse de manera similar, usaremos pruebas de hipótesis para mostrar si su supervivencia es similar, lo haremos por parejas para encontrar todos los posibles resultados, posteriormente se agruparán los empleos con supervivencia similar.
Prueba de Hipótesis
Utilizaremos la prueba Log-Rank para comparar dos funciones de supervivencia mediante el siguiente contraste:
\[H_0: S_1(t)=S_2(t)~~~\forall t>0\\ vs \\ H_1: S_1(t)\neq S_2(t)~~p.a.~t>0\] El estadistico de esta prueba es:
\[L^2=\frac{ (\sum_{i=1}^k(d_{1i}-e_{1i}))^2}{\sum_{i=1}^kV_{1i}}\] Donde:
\[e_{1i}=n_{1i}\frac{d_i}{n_i}\\ V_{1i}=\frac{n_{1i}n_{2i}d_i(n_i-d_i)}{n_i^2(n_i-1)}\] Se rechaza \(H_0\) cuando: \[L^2>J_{1-\alpha}\] Donde \(J_{1-\alpha}\) es el cuantil \(1-\alpha\) de una distribución \(\chi^2_{(1)}\).
Los resultados son los siguientes:
#Creamos una funcion que haga la prueba Log-Rank
#recibe las 2 funciones de supervivencia que serán testeadas
LogRank <- function(S1,S2){
n1i<- S1$n.risk
n2i<- S2$n.risk
d1i<- S1$n.event
di <- S1$n.event + S2$n.event
ni <- S1$n.risk + S2$n.risk
e1i <- n1i*di/ni
v1i <- n1i*n2i*di*(ni-di)/((ni^2)*(ni-1))
#Calculamos el estadístico
L <- round((sum(d1i-e1i))^2/(sum(v1i)),3)
#Regla de decisión
if (L > qchisq(0.95, df = 1)){
conclusion <- "Se rechaza Ho"
}else{
conclusion <- "No se rechaza Ho"
}
resultado <- c(L,conclusion)
return(resultado)
}
p1<-LogRank(Surv_j1_fit,Surv_j2_fit)
p2<-LogRank(Surv_j1_fit,Surv_j3_fit)
p3<-LogRank(Surv_j1_fit,Surv_j4_fit)
p4<-LogRank(Surv_j1_fit,Surv_j5_fit)
p5<-LogRank(Surv_j1_fit,Surv_j6_fit)
p6<-LogRank(Surv_j1_fit,Surv_j7_fit)
p7<-LogRank(Surv_j1_fit,Surv_j8_fit)
p8<-LogRank(Surv_j1_fit,Surv_j9_fit)
p9<-LogRank(Surv_j1_fit,Surv_j10_fit)
p10<-LogRank(Surv_j1_fit,Surv_j11_fit)
p11<-LogRank(Surv_j1_fit,Surv_j12_fit)
p12<-LogRank(Surv_j2_fit,Surv_j3_fit)
p13<-LogRank(Surv_j2_fit,Surv_j4_fit)
p14<-LogRank(Surv_j2_fit,Surv_j5_fit)
p15<-LogRank(Surv_j2_fit,Surv_j6_fit)
p16<-LogRank(Surv_j2_fit,Surv_j7_fit)
p17<-LogRank(Surv_j2_fit,Surv_j8_fit)
p18<-LogRank(Surv_j2_fit,Surv_j9_fit)
p19<-LogRank(Surv_j2_fit,Surv_j10_fit)
p20<-LogRank(Surv_j2_fit,Surv_j11_fit)
p21<-LogRank(Surv_j2_fit,Surv_j12_fit)
p22<-LogRank(Surv_j3_fit,Surv_j4_fit)
p23<-LogRank(Surv_j3_fit,Surv_j5_fit)
p24<-LogRank(Surv_j3_fit,Surv_j6_fit)
p25<-LogRank(Surv_j3_fit,Surv_j7_fit)
p26<-LogRank(Surv_j3_fit,Surv_j8_fit)
p27<-LogRank(Surv_j3_fit,Surv_j9_fit)
p28<-LogRank(Surv_j3_fit,Surv_j10_fit)
p29<-LogRank(Surv_j3_fit,Surv_j11_fit)
p30<-LogRank(Surv_j3_fit,Surv_j12_fit)
p31<-LogRank(Surv_j4_fit,Surv_j5_fit)
p32<-LogRank(Surv_j4_fit,Surv_j6_fit)
p33<-LogRank(Surv_j4_fit,Surv_j7_fit)
p34<-LogRank(Surv_j4_fit,Surv_j8_fit)
p35<-LogRank(Surv_j4_fit,Surv_j9_fit)
p36<-LogRank(Surv_j4_fit,Surv_j10_fit)
p37<-LogRank(Surv_j4_fit,Surv_j11_fit)
p38<-LogRank(Surv_j4_fit,Surv_j12_fit)
p39<-LogRank(Surv_j5_fit,Surv_j6_fit)
p40<-LogRank(Surv_j5_fit,Surv_j7_fit)
p41<-LogRank(Surv_j5_fit,Surv_j8_fit)
p42<-LogRank(Surv_j5_fit,Surv_j9_fit)
p43<-LogRank(Surv_j5_fit,Surv_j10_fit)
p44<-LogRank(Surv_j5_fit,Surv_j11_fit)
p45<-LogRank(Surv_j5_fit,Surv_j12_fit)
p46<-LogRank(Surv_j6_fit,Surv_j7_fit)
p47<-LogRank(Surv_j6_fit,Surv_j8_fit)
p48<-LogRank(Surv_j6_fit,Surv_j9_fit)
p49<-LogRank(Surv_j6_fit,Surv_j10_fit)
p50<-LogRank(Surv_j6_fit,Surv_j11_fit)
p51<-LogRank(Surv_j6_fit,Surv_j12_fit)
p52<-LogRank(Surv_j7_fit,Surv_j8_fit)
p53<-LogRank(Surv_j7_fit,Surv_j9_fit)
p54<-LogRank(Surv_j7_fit,Surv_j10_fit)
p55<-LogRank(Surv_j7_fit,Surv_j11_fit)
p56<-LogRank(Surv_j7_fit,Surv_j12_fit)
p57<-LogRank(Surv_j8_fit,Surv_j9_fit)
p58<-LogRank(Surv_j8_fit,Surv_j10_fit)
p59<-LogRank(Surv_j8_fit,Surv_j11_fit)
p60<-LogRank(Surv_j8_fit,Surv_j12_fit)
p61<-LogRank(Surv_j9_fit,Surv_j10_fit)
p62<-LogRank(Surv_j9_fit,Surv_j11_fit)
p63<-LogRank(Surv_j9_fit,Surv_j12_fit)
p64<-LogRank(Surv_j10_fit,Surv_j11_fit)
p65<-LogRank(Surv_j10_fit,Surv_j12_fit)
p66<-LogRank(Surv_j11_fit,Surv_j12_fit)
s1<-c("administracion","administracion","administracion","administracion","administracion","administracion","administracion","administracion","administracion","administracion","administracion",
"tecnico","tecnico","tecnico","tecnico","tecnico","tecnico","tecnico","tecnico","tecnico","tecnico",
"emprendedor","emprendedor","emprendedor","emprendedor","emprendedor","emprendedor","emprendedor","emprendedor","emprendedor",
"obrero","obrero","obrero","obrero","obrero","obrero","obrero","obrero",
"desconocido","desconocido","desconocido","desconocido","desconocido","desconocido","desconocido",
"retirado","retirado","retirado","retirado","retirado","retirado",
"administrativo","administrativo","administrativo","administrativo","administrativo",
"servicios","servicios","servicios","servicios",
"independiente","independiente","independiente",
"desempleado","desempleado",
"empleado domestico")
s2<-c("tecnico","emprendedor","obrero","desconocido","retirado","administrativo","servicios","independiente","desempleado","empleado domestico","estudiante",
"emprendedor","obrero","desconocido","retirado","administrativo","servicios","independiente","desempleado","empleado domestico","estudiante",
"obrero","desconocido","retirado","administrativo","servicios","independiente","desempleado","empleado domestico","estudiante",
"desconocido","retirado","administrativo","servicios","independiente","desempleado","empleado domestico","estudiante",
"retirado","administrativo","servicios","independiente","desempleado","empleado domestico","estudiante",
"administrativo","servicios","independiente","desempleado","empleado domestico","estudiante",
"servicios","independiente","desempleado","empleado domestico","estudiante",
"independiente","desempleado","empleado domestico","estudiante",
"desempleado","empleado domestico","estudiante",
"empleado domestico","estudiante",
"estudiante")
pruebas<-data.frame(s1,s2,
c(p1[1],p2[1],p3[1],p4[1],p5[1],p6[1],p7[1],p8[1],p9[1],p10[1],
p11[1],p12[1],p13[1],p14[1],p15[1],p16[1],p17[1],p18[1],p19[1],p20[1],
p21[1],p22[1],p23[1],p24[1],p25[1],p26[1],p27[1],p28[1],p29[1],p30[1],
p31[1],p32[1],p33[1],p34[1],p35[1],p36[1],p37[1],p38[1],p39[1],p40[1],
p41[1],p42[1],p43[1],p44[1],p45[1],p46[1],p47[1],p48[1],p49[1],p50[1],
p51[1],p52[1],p53[1],p54[1],p55[1],p56[1],p57[1],p58[1],p59[1],p60[1],
p61[1],p62[1],p63[1],p64[1],p65[1],p66[1]),
c(p1[2],p2[2],p3[2],p4[2],p5[2],p6[2],p7[2],p8[2],p9[2],p10[2],
p11[2],p12[2],p13[2],p14[2],p15[2],p16[2],p17[2],p18[2],p19[2],p20[2],
p21[2],p22[2],p23[2],p24[2],p25[2],p26[2],p27[2],p28[2],p29[2],p30[2],
p31[2],p32[2],p33[2],p34[2],p35[2],p36[2],p37[2],p38[2],p39[2],p40[2],
p41[2],p42[2],p43[2],p44[2],p45[2],p46[2],p47[2],p48[2],p49[2],p50[2],
p51[2],p52[2],p53[2],p54[2],p55[2],p56[2],p57[2],p58[2],p59[2],p60[2],
p61[2],p62[2],p63[2],p64[2],p65[2],p66[2]))
colnames(pruebas)<-c("S1","S2","Estadistico","Conclusión")
pruebasDe lo anterior se observa que podemos fusionar los tipos:
* administracion-tecnico-obrero, los agrupamos en obrero.
* servicio-independiente-estudiante-desempleado, los agrupamos en servicio.
Obs: Existen más tipos cuyas supervivencias son similares pero el fin es agrupar la mayor cantidad de tipos, por lo que nos quedamos con los grupos anteriores.
La gráfica con los nuevos grupos es la siguientes
ggplot()+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="emprendedor"))+
geom_step(mapping = aes(x=a4$Tiempo, y=a4$Surv, color="obrero"))+
geom_step(mapping = aes(x=a5$Tiempo, y=a5$Surv, color="desconocido"))+
geom_step(mapping = aes(x=a6$Tiempo, y=a6$Surv, color="retirado"))+
geom_step(mapping = aes(x=a7$Tiempo, y=a7$Surv, color="administrativo"))+
geom_step(mapping = aes(x=a8$Tiempo, y=a8$Surv, color="servicios"))+
geom_step(mapping = aes(x=a11$Tiempo, y=a11$Surv, color="empleado domestico"))+
labs(title="Supervivencia por tipo de empleo", x="Edad", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Se observa que el grupo perfectamente distinguible es el de los retirados, cuya supervivencia es mayor que para los demás hasta superar los \(2000\) segundos
Media
La siguiente tabla muestra las medias junto a su varianza y los limites de sus intervalos de confianza para cada tipo de empleo:
m1<-Media(Surv_j3_fit)
m2<-Media(Surv_j4_fit)
m3<-Media(Surv_j5_fit)
m4<-Media(Surv_j6_fit)
m5<-Media(Surv_j7_fit)
m6<-Media(Surv_j8_fit)
m7<-Media(Surv_j11_fit)
medias<-data.frame(c("Emprendedor","Obrero","Desconocido","Retirado","Admin","Servicios","Empleado doméstico"),
c(m1[1],m2[1],m3[1],m4[1],m5[1],m6[1],m7[1]),
c(m1[2],m2[2],m3[2],m4[2],m5[2],m6[2],m7[2]),
c(m1[3],m2[3],m3[3],m4[3],m5[3],m6[3],m7[3]),
c(m1[4],m2[4],m3[4],m4[4],m5[4],m6[4],m7[4]))
kable(medias , caption = "Media de la duración de la llamada para cada tipo de empleo"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Empleo","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Empleo | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| Emprendedor | 303.2 | 210.34 | 274.7 | 331.6 |
| Obrero | 308.2 | 36.18 | 296.4 | 320.0 |
| Desconocido | 269.8 | 413.09 | 230.0 | 309.7 |
| Retirado | 415.7 | 279.17 | 383.0 | 448.5 |
| Admin | 336.9 | 142.96 | 313.4 | 360.3 |
| Servicios | 317.4 | 151.31 | 293.3 | 341.6 |
| Empleado doméstico | 279.2 | 159.21 | 254.4 | 303.9 |
El grupo de los retirados es el que presenta una mayor media de entre todos los demás tipos de empleo.
Cuantiles
La siguiente tabla muetras los cuantiles \(0.25, 0.5\) y \(0.75\) para las duraciones de las llamadas de acuerdo a cada tipo de empleo.
c1<-Cuantiles(Surv_j3_fit)
c2<-Cuantiles(Surv_j4_fit)
c3<-Cuantiles(Surv_j5_fit)
c4<-Cuantiles(Surv_j6_fit)
c5<-Cuantiles(Surv_j7_fit)
c6<-Cuantiles(Surv_j8_fit)
c7<-Cuantiles(Surv_j11_fit)
cuantiles<-data.frame(c1$p,c1$c,c2$c,c3$c,c4$c,c5$c,c6$c,c7$c)
kable(cuantiles , caption = "Cuantiles de las duraciones"
, align = c('c', 'c', 'c', 'c','c','c','c','c'), col.names=c("Porcentaje","Emprendedor","Obrero","Desconocido","Retirado","Admin","Servicios","Empleado doméstico")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | Emprendedor | Obrero | Desconocido | Retirado | Admin | Servicios | Empleado doméstico |
|---|---|---|---|---|---|---|---|
| 25% | 102 | 108 | 88 | 122 | 102 | 107 | 97 |
| 50% | 181 | 189 | 166 | 224 | 179 | 189 | 166 |
| 75% | 340 | 331 | 346 | 476 | 340 | 327 | 314 |
Funciones de Riesgo
Funciones de riesgo \(\hat h(t)\) y de riesgo acumulado \(\hat H(t)\) para cada tipo de empleo
#funcion de riesgo para cada tipo
ht1<-Surv_j3_fit$n.event/Surv_j3_fit$n.risk
ht2<-Surv_j4_fit$n.event/Surv_j4_fit$n.risk
ht3<-Surv_j5_fit$n.event/Surv_j5_fit$n.risk
ht4<-Surv_j6_fit$n.event/Surv_j6_fit$n.risk
ht5<-Surv_j7_fit$n.event/Surv_j7_fit$n.risk
ht6<-Surv_j8_fit$n.event/Surv_j8_fit$n.risk
ht7<-Surv_j11_fit$n.event/Surv_j11_fit$n.risk
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=Surv_j3_fit$time,y=ht1,colour="Emprendedor"))+
geom_line(aes(x=Surv_j4_fit$time,y=ht2,colour="Obrero"))+
geom_line(aes(x=Surv_j5_fit$time,y=ht3,colour="Desconocido"))+
geom_line(aes(x=Surv_j6_fit$time,y=ht4,colour="Retirado"))+
geom_line(aes(x=Surv_j7_fit$time,y=ht5,colour="Admin"))+
geom_line(aes(x=Surv_j8_fit$time,y=ht6,colour="Servicios"))+
geom_line(aes(x=Surv_j11_fit$time,y=ht7,colour="Domestico"))+
theme_bw() +
labs(title="h(t) por tipo de empleo", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p2<-ggplot()+
geom_step(aes(x=Surv_j3_fit$time,y=Surv_j3_fit$cumhaz, colour="Emprendedor"))+
geom_step(aes(x=Surv_j4_fit$time,y=Surv_j4_fit$cumhaz, colour="Obrero"))+
geom_step(aes(x=Surv_j5_fit$time,y=Surv_j5_fit$cumhaz, colour="Desconocido"))+
geom_step(aes(x=Surv_j6_fit$time,y=Surv_j6_fit$cumhaz, colour="Retirado"))+
geom_step(aes(x=Surv_j7_fit$time,y=Surv_j7_fit$cumhaz, colour="Admin"))+
geom_step(aes(x=Surv_j8_fit$time,y=Surv_j8_fit$cumhaz, colour="Servicios"))+
geom_step(aes(x=Surv_j11_fit$time,y=Surv_j11_fit$cumhaz, colour="Domestico"))+
theme_bw()+
labs(title="H(t) por tipo de empleo", x="Duración", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
#p1;p2
ggarrange(p1,p2, ncol = 2, nrow = 1)
No es posible comparar las funciones de riesgo por el comportamiento que presentan pero en cambio para \(H(t)\) que como ya se veia en la funcion de supervivencia, los retirados son los que menos riesgo acumulan de entre todos los tipos de empleo. Son los que mejor responden a la campaña y deciden realizar depósitos.
Máxima segregación
Creamos una función que permita conocer la mayor distancia entre 2 funciones de supervivencia \(S(t)\) dada una probabilidad como cota inferior.
#la funcion recibe como parametros:
#f1: dato tipo survfit
#f2: dato tipo survfit
#etiqueta1 y etiqueta2: corresponden a f1 y f2 respectivamente, cadenas que se usan para personalizar la grafica
# en este caso pueden tomar los valores {1,2,3,4} que corresponden a las 4 marcas
#de focos
#p: probabilidad cota inferior sobre la que se calcula la distancia
Maxima_segregacion<-function(f1,f2,etiqueta1,etiqueta2,p){
#obtenemos tiempo y supervivencia de los 2 conjuntos
t1<-f1$time
s1<-f1$surv
t2<-f2$time
s2<-f2$surv
#suavizamos ambas funciones de supervivencia y las evaluamos sobre el mismo conjunto
#para poder hacer la comparacion
s1aprox<-approxfun(s1,t1)
s2aprox<-approxfun(s2,t2)
df1<-data.frame(s=s1,ta1=s1aprox(s1), ta2=s2aprox(s1))
#diferencia entre los tiempos para las funciones suavizadas
df1$distancia<-abs(df1$ta1-df1$ta2)
#Eliminamos NA's
df1 <- df1[!is.na(df1$distancia),]
#Conservamos las S(t)>p
df<-df1[df1$s>p,]
#obtenemos la mayor de las distancias entre las gráficas
maxima<-max(df$distancia)
#seleccionamos el renglon que corresponda a la distancia maxima
d<-df1[df1$distancia==maxima,]
max<-as.character(round(maxima,4))
grafica<-ggplot()+
geom_step(aes(t1,s1, color=etiqueta1))+
geom_step(aes(t2,s2, color=etiqueta2))+
theme_bw()+
labs(title="Máxima Distancia",
subtitle = paste("para S(t)>",p,sep = " "),x="Tiempo", y=TeX("$\\hat{S}(t)$"))+
theme(plot.title = element_text(hjust = 0.5),plot.subtitle = element_text(hjust = 0.5))+
geom_hline(yintercept = d$s, linetype="dashed", color = "gray50", size=0.7)+
scale_colour_discrete("")+
labs(color = "Grupo") +
geom_vline(xintercept = d$ta1, linetype="dashed", color = "gray50", size=0.7)+
geom_vline(xintercept = d$ta2, linetype="dashed", color = "gray50", size=0.7)+
annotate("text", x =min(d$ta1,d$ta2), y = d$s+0.02, label = paste("Distancia=",max,sep = " "))+
theme(legend.position = "bottom", legend.direction = "horizontal")
tabla<-kable(data.frame(d$s,d$ta1,d$ta2,round(maxima,4)) , caption = "Máxima Segregacion"
, align = c('c', 'c', 'c', 'c'), col.names=c("Probabilidad","Tiempo 1","Tiempo 2","Distancia maxima")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))
return(list(tabla=tabla,grafica=grafica))
}Podemos comparar los grupos retirado y servicio, el primero es el que más se distingue de los demás, lo comparamos contra servicios el cual, tras la prueba de hipótesis es el que más grupos “absorvió”.
#empleo
maxima<-Maxima_segregacion(Surv_j6_fit,Surv_j8_fit,"Retirado","Servicios",0)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.0375 | 2027 | 1365 | 661.6 |
Estado Civil
Tipos de estado civil:
- Divorciado(incluye a los viudos)
- Casado
- Soltero
Las funciones de supervivencia son las siguientes:
#separamos los datos de acuerdo a los tipos de estado civil
m1=subset(d,marital=="married", select = c(duration,Status))
m2=subset(d,marital=="single", select = c(duration,Status))
m3=subset(d,marital=="divorced", select = c(duration,Status))
#Supervivencias
Surv_m1_fit <- survfit(Surv(time = m1$duration, event = m1$Status)~1)
a1 <-data.frame(Tiempo = Surv_m1_fit$time, Surv = Surv_m1_fit$surv)
Surv_m2_fit <- survfit(Surv(time = m2$duration, event = m2$Status)~1)
a2 <-data.frame(Tiempo = Surv_m2_fit$time, Surv = Surv_m2_fit$surv)
Surv_m3_fit <- survfit(Surv(time = m3$duration, event = m3$Status)~1)
a3 <-data.frame(Tiempo = Surv_m3_fit$time, Surv = Surv_m3_fit$surv)
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="Casado"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Soltero"))+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="Divorciado"))+
labs(title="Supervivencia por estado civil", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Pruebas de Hipótesis
Veamos si es posible agrupar algunas supervivencias.
prueba1<-LogRank(Surv_m1_fit,Surv_m2_fit)
prueba2<-LogRank(Surv_m1_fit,Surv_m3_fit)
prueba3<-LogRank(Surv_m2_fit,Surv_m3_fit)
s1<-c("Casado","Casado","Soltero")
s2<-c("Soltero","Divorciado", "Divorciado")
pruebas<-data.frame(s1,s2,c(prueba1[1],prueba2[1],prueba3[1]),
c(prueba1[2],prueba2[2],prueba3[2]))
colnames(pruebas)<-c("S1","S2","Estadistico","Conclusión")
pruebasComo se puede observar, no es posible agrupar las supervivencias, cada grupo está bien diferenciado de los demás.
Medias
m1<-Media(Surv_m1_fit)
m2<-Media(Surv_m2_fit)
m3<-Media(Surv_m3_fit)
medias<-data.frame(c("Casado","Soltero","Divorciado"),
c(m1[1],m2[1],m3[1]),
c(m1[2],m2[2],m3[2]),
c(m1[3],m2[3],m3[3]),
c(m1[4],m2[4],m3[4]))
kable(medias , caption = "Media de la duración de la llamada por estado civil"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Estado","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Estado | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| Casado | 313.6 | 23.74 | 304.1 | 323.2 |
| Soltero | 366.8 | 83.88 | 348.9 | 384.8 |
| Divorciado | 342.9 | 117.39 | 321.7 | 364.1 |
El grupo que posee un mayor valor medio de la duracion es el de los clientes solteros
Cuantiles
La siguiente tabla muetras los cuantiles \(0.25,0.5\) y \(0.75\) para las duraciones de las llamadas de acuerdo a cada tipo de empleo.
c1<-Cuantiles(Surv_m1_fit)
c2<-Cuantiles(Surv_m2_fit)
c3<-Cuantiles(Surv_m3_fit)
cuantiles<-data.frame(c1$p,c1$c,c2$c,c3$c)
kable(cuantiles , caption = "Cuantiles de las duraciones"
, align = c('c', 'c', 'c', 'c'), col.names=c("Porcentaje","Casado","Soltero","Divorciado")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | Casado | Soltero | Divorciado |
|---|---|---|---|
| 25% | 102 | 107 | 105 |
| 50% | 181 | 196 | 183 |
| 75% | 333 | 372 | 345 |
De igual manera, el grupo de los clientes solteros logra acumular mas segundos de duracion que los casados y divorciados.
Como las medias son mayores a las medianas en los tres casos, las supervivencias tienen sesgo hacia la derecha.
Funciones de Riesgo
Funciones de riesgo \(\hat h(t)\) y de riesgo acumulado \(\hat H(t)\) por estado civil.
ht1<-Surv_m1_fit$n.event/Surv_m1_fit$n.risk
ht2<-Surv_m2_fit$n.event/Surv_m2_fit$n.risk
ht3<-Surv_m3_fit$n.event/Surv_m3_fit$n.risk
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=Surv_m1_fit$time,y=ht1,colour="Casado"))+
geom_line(aes(x=Surv_m2_fit$time,y=ht2,colour="Soltero"))+
geom_line(aes(x=Surv_m3_fit$time,y=ht3,colour="Divorciado"))+
theme_bw() +
labs(title="h(t) por estado civil", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p2<-ggplot()+
geom_step(aes(x=Surv_m1_fit$time,y=Surv_m1_fit$cumhaz, colour="Casado"))+
geom_step(aes(x=Surv_m2_fit$time,y=Surv_m2_fit$cumhaz, colour="Soltero"))+
geom_step(aes(x=Surv_m3_fit$time,y=Surv_m3_fit$cumhaz, colour="Divorciado"))+
theme_bw()+
labs(title="H(t) por estado civil", x="Duración", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p1;p2
#ggarrange(p1,, ncol = 2, nrow = 1)Para los casados el riesgo acumulado es significativamente mayor que para los otros dos grupos. Esto indica que a medida que aumenta la duracion del contacto, son los clientes casados los mas propensos a no realizar sus depósitos y a seguir siendo molestados en futuras campañas.
Máxima Segregación
Comparamos las supervivencias de Casado y Soltero, sobre una probabilidad mayor a \(0.05\)
#estado civil
maxima<-Maxima_segregacion(Surv_m1_fit,Surv_m2_fit,"Casado","Soltero",0.05)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.0501 | 1019 | 1375 | 356.3 |
Nivel Educativo
Los niveles registrados para los clientes son:
- Primaria
- Secundaria
- Terciaria(Curso profesional/título universitario)
- Desconocido
#separamos los datos de acuerdo al nivel educativo
e1=subset(d,education=="primary", select = c(duration,Status))
e2=subset(d,education=="secondary", select = c(duration,Status))
e3=subset(d,education=="tertiary", select = c(duration,Status))
e4=subset(d,education=="unknown", select = c(duration,Status))
#Supervivencias
Surv_e1_fit <- survfit(Surv(time = e1$duration, event = e1$Status)~1)
a1 <-data.frame(Tiempo = Surv_e1_fit$time, Surv = Surv_e1_fit$surv)
Surv_e2_fit <- survfit(Surv(time = e2$duration, event = e2$Status)~1)
a2 <-data.frame(Tiempo = Surv_e2_fit$time, Surv = Surv_e2_fit$surv)
Surv_e3_fit <- survfit(Surv(time = e3$duration, event = e3$Status)~1)
a3 <-data.frame(Tiempo = Surv_e3_fit$time, Surv = Surv_e3_fit$surv)
Surv_e4_fit <- survfit(Surv(time = e4$duration, event = e4$Status)~1)
a4 <-data.frame(Tiempo = Surv_e4_fit$time, Surv = Surv_e4_fit$surv)
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="Primaria"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Secundaria"))+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="Terciaria"))+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="Desconocida"))+
labs(title="Supervivencia por nivel educativo", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Pruebas de Hipótesis
prueba1<-LogRank(Surv_e1_fit,Surv_e2_fit)
prueba2<-LogRank(Surv_e1_fit,Surv_e3_fit)
prueba3<-LogRank(Surv_e1_fit,Surv_e4_fit)
prueba4<-LogRank(Surv_e2_fit,Surv_e3_fit)
prueba5<-LogRank(Surv_e2_fit,Surv_e4_fit)
prueba6<-LogRank(Surv_e3_fit,Surv_e4_fit)
s1<-c("Primaria","Primaria","Primaria", "Secundaria","Secundaria","Terciaria")
s2<-c("Secundaria","Terciaria","Desconocida","Terciaria","Desconocida","Desconocida")
pruebas<-data.frame(s1,s2,c(prueba1[1],prueba2[1],prueba3[1],prueba4[1],prueba5[1],prueba6[1]),
c(prueba1[2],prueba2[2],prueba3[2],prueba4[2],prueba5[2],prueba6[2]))
colnames(pruebas)<-c("S1","S2","Estadistico","Conclusión")
pruebasComo se puede observar, podemos agrupar a las personas que tienen hasta el nivel de educación primario con los que se desconoce su nivel.
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="Primaria"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Secundaria"))+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="Terciaria"))+
labs(title="Supervivencia por nivel educativo", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Se puede ver que para las personas con educacion terciaria la supervivencia es mayor que para los demás niveles.
Medias
m1<-Media(Surv_e1_fit)
m2<-Media(Surv_e2_fit)
m3<-Media(Surv_e3_fit)
medias<-data.frame(c("Primaria","Secundaria","Terciaria"),
c(m1[1],m2[1],m3[1]),
c(m1[2],m2[2],m3[2]),
c(m1[3],m2[3],m3[3]),
c(m1[4],m2[4],m3[4]))
kable(medias , caption = "Media de la duración de la llamada por nivel educativo"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Educacion","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Educacion | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| Primaria | 311.9 | 52.55 | 297.7 | 326.1 |
| Secundaria | 325.6 | 28.62 | 315.2 | 336.1 |
| Terciaria | 352.0 | 84.51 | 334.0 | 370.0 |
Las personas con educacion terciaria poseen una mayor media de duracion del contacto que los demas niveles.
Cuantiles
c1<-Cuantiles(Surv_e1_fit)
c2<-Cuantiles(Surv_e2_fit)
c3<-Cuantiles(Surv_e3_fit)
cuantiles<-data.frame(c1$p,c1$c,c2$c,c3$c)
kable(cuantiles , caption = "Cuantiles de las duraciones"
, align = c('c', 'c', 'c', 'c'), col.names=c("Porcentaje","Primaria","Secundaria","Terciaria")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | Primaria | Secundaria | Terciaria |
|---|---|---|---|
| 25% | 103 | 106 | 101 |
| 50% | 180 | 188 | 184 |
| 75% | 322 | 343 | 362 |
Funciones de Riesgo
ht1<-Surv_e1_fit$n.event/Surv_e1_fit$n.risk
ht2<-Surv_e2_fit$n.event/Surv_e2_fit$n.risk
ht3<-Surv_e3_fit$n.event/Surv_e3_fit$n.risk
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=Surv_e1_fit$time,y=ht1,colour="Primaria"))+
geom_line(aes(x=Surv_e2_fit$time,y=ht2,colour="Secundaria"))+
geom_line(aes(x=Surv_e3_fit$time,y=ht3,colour="Terciaria"))+
theme_bw() +
labs(title="h(t) por nivel educativo", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p2<-ggplot()+
geom_step(aes(x=Surv_e1_fit$time,y=Surv_e1_fit$cumhaz, colour="Primaria"))+
geom_step(aes(x=Surv_e2_fit$time,y=Surv_e2_fit$cumhaz, colour="Secundaria"))+
geom_step(aes(x=Surv_e3_fit$time,y=Surv_e3_fit$cumhaz, colour="Terciaria"))+
theme_bw()+
labs(title="H(t) por nivel educativo", x="Duración", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p1;p2
#ggarrange(p1,, ncol = 2, nrow = 1)La funcion de riesgo acumulado de educacion terciaria es la que menos riesgo acumula, haciendo sentido que sea el grupo con la supervivencia mas alta, por ende, las personas con educacion terciaria o mayor son los que más realizan los depositos al momento de ser contactados.
Máxima Segregación
Contrastamos las supervivencias de nivel Pirmaria y nivel Terciaria
#nivel educativo
maxima<-Maxima_segregacion(Surv_e1_fit,Surv_e3_fit,"Primario","Terciario",0.05)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.0503 | 1009 | 1330 | 321 |
Crédito en Mora
Se separan los clientes de acuerdo a si tienen o no crédito en mora.
- No
- Si
de1=subset(d,default=="no", select = c(duration,Status))
de2=subset(d,default=="yes", select = c(duration,Status))
#Supervivencias
Surv_de1_fit <- survfit(Surv(time = de1$duration, event = de1$Status)~1)
a1 <-data.frame(Tiempo = Surv_de1_fit$time, Surv = Surv_de1_fit$surv)
Surv_de2_fit <- survfit(Surv(time = de2$duration, event = de2$Status)~1)
a2 <-data.frame(Tiempo = Surv_de2_fit$time, Surv = Surv_de2_fit$surv)
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="No"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Si"))+
labs(title="Supervivencia por credito en mora", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
La supervivencia para los que tienen credito en mora es bastante mas baja y corta que para aquellos que no cuentan con el credito.
Pruebas de Hipótesis
prueba1<-LogRank(Surv_de1_fit,Surv_de2_fit)
s1<-c("No")
s2<-c("Si")
pruebas<-data.frame(s1,s2,c(prueba1[1]),c(prueba1[2]))
colnames(pruebas)<-c("S1","S2","Estadistico","Conclusión")
pruebasComo se puede observar, no hay forma de agrupar las supervivencias.
Medias
m1<-Media(Surv_de1_fit)
m2<-Media(Surv_de2_fit)
medias<-data.frame(c("No","Si"),
c(m1[1],m2[1]),
c(m1[2],m2[2]),
c(m1[3],m2[3]),
c(m1[4],m2[4]))
kable(medias , caption = "Media de la duración de la llamada por estado del credito"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Educacion","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Educacion | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| No | 334.0 | 22.08 | 324.8 | 343.2 |
| Si | 259.4 | 111.55 | 238.7 | 280.1 |
La media de duracion es mayor para los que no tiene credito en mora.
Cuantiles
c1<-Cuantiles(Surv_de1_fit)
c2<-Cuantiles(Surv_de2_fit)
cuantiles<-data.frame(c1$p,c1$c,c2$c)
kable(cuantiles , caption = "Cuantiles de las duraciones"
, align = c('c', 'c', 'c'), col.names=c("Porcentaje","No","Si")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | No | Si |
|---|---|---|
| 25% | 104 | 97 |
| 50% | 185 | 172 |
| 75% | 346 | 305 |
Para los cuantiles pasa lo mismo, aquellos que tienen credito en mora duran menos en la llamada que los que no lo tienen.
Funciones de Riesgo
ht1<-Surv_de1_fit$n.event/Surv_de1_fit$n.risk
ht2<-Surv_de2_fit$n.event/Surv_de2_fit$n.risk
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=Surv_de1_fit$time,y=ht1,colour="No"))+
geom_line(aes(x=Surv_de2_fit$time,y=ht2,colour="Si"))+
theme_bw() +
labs(title="h(t) por estado de credito en mora", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p2<-ggplot()+
geom_step(aes(x=Surv_de1_fit$time,y=Surv_de1_fit$cumhaz, colour="No"))+
geom_step(aes(x=Surv_de2_fit$time,y=Surv_de2_fit$cumhaz, colour="Si"))+
theme_bw()+
labs(title="H(t) por estado de credito en mora", x="Duración", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
#p1;p2
ggarrange(p1,p2, ncol = 2, nrow = 1)
Al largo plazo las personas sin credito en mora acumulan mayor riesgo pero como se observa, aquellos que poseen esta caracteristica, acumulan el riesgo de forma más apresurada hasta que la supervivencia llega a 0 antes que para el otro grupo.
Lo anterior indica que aquellos con creditos en mora, será más dificil que realicen los depositos.
Máxima Segregación
Veamos la máxima segregacion sobre una probabilidad de \(0\)
#credito en mora
maxima<-Maxima_segregacion(Surv_de1_fit,Surv_de2_fit,"No","Si",0)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.0077 | 3284 | 1501 | 1783 |
Préstamo para vivienda
Categorias:
- No
- Si
h1=subset(d,housing=="no", select = c(duration,Status))
h2=subset(d,housing=="yes", select = c(duration,Status))
#Supervivencias
Surv_h1_fit <- survfit(Surv(time = h1$duration, event = h1$Status)~1)
a1 <-data.frame(Tiempo = Surv_h1_fit$time, Surv = Surv_h1_fit$surv)
Surv_h2_fit <- survfit(Surv(time = h2$duration, event = h2$Status)~1)
a2 <-data.frame(Tiempo = Surv_h2_fit$time, Surv = Surv_h2_fit$surv)
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="No"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Si"))+
labs(title="Supervivencia por prestamo para vivienda", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Para aquellos que no tienen prestamo para la vivienda, su supervivencia es mayor que para los que si lo tienen.
Pruebas de Hipótesis
Veamos si es posible agrupar las supervivencias
prueba1<-LogRank(Surv_h1_fit,Surv_h2_fit)
s1<-c("No")
s2<-c("Si")
pruebas<-data.frame(s1,s2,c(prueba1[1]),c(prueba1[2]))
colnames(pruebas)<-c("S1","S2","Estadistico","Conclusión")
pruebasComo se puede observar, no hay forma de agrupar las supervivencias.
Medias
m1<-Media(Surv_h1_fit)
m2<-Media(Surv_h2_fit)
medias<-data.frame(c("No","Si"),
c(m1[1],m2[1]),
c(m1[2],m2[2]),
c(m1[3],m2[3]),
c(m1[4],m2[4]))
kable(medias , caption = "Media de la duración de la llamada por prestamo para la vivienda"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Educacion","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Educacion | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| No | 361 | 44.95 | 347.9 | 374.2 |
| Si | 314 | 41.76 | 301.3 | 326.7 |
la media es mayor para los que no pidieron prestamos para la vivienda.
Cuantiles
c1<-Cuantiles(Surv_h1_fit)
c2<-Cuantiles(Surv_h2_fit)
cuantiles<-data.frame(c1$p,c1$c,c2$c)
kable(cuantiles , caption = "Cuantiles de las duraciones"
, align = c('c', 'c', 'c'), col.names=c("Porcentaje","No","Si")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | No | Si |
|---|---|---|
| 25% | 104 | 104 |
| 50% | 187 | 184 |
| 75% | 373 | 328 |
Para los cuantiles se sigue teniendo el mismo comportamiento, aquellas personas que no poseen creditos para la vivienda, alcanzan una mayor duracion en el contacto.
Funciones de Riesgo
ht1<-Surv_h1_fit$n.event/Surv_h1_fit$n.risk
ht2<-Surv_h2_fit$n.event/Surv_h2_fit$n.risk
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=Surv_h1_fit$time,y=ht1,colour="No"))+
geom_line(aes(x=Surv_h2_fit$time,y=ht2,colour="Si"))+
theme_bw() +
labs(title="h(t) por prestamo para vivienda", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p2<-ggplot()+
geom_step(aes(x=Surv_h1_fit$time,y=Surv_h1_fit$cumhaz, colour="No"))+
geom_step(aes(x=Surv_h2_fit$time,y=Surv_h2_fit$cumhaz, colour="Si"))+
theme_bw()+
labs(title="H(t) por prestamo para vivienda", x="Duración", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
#p1;p2
ggarrange(p1,p2, ncol = 2, nrow = 1)
Como era de esperarse aquellas personas que cuentan con un credito para la vivienda acumulan mayor riesgo. Esto quiere decir que para personas que hayan pedido este crédito, les será más dificil poder realizar el depósito que se solicita en la campaña.
Máxima Segregación
Veamos la maxima segregacion para una probabilidad mayor a \(0\)
#vivienda
maxima<-Maxima_segregacion(Surv_h1_fit,Surv_h2_fit,"No","Si",0)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.0114 | 3104 | 2241 | 863.1 |
| 0.0114 | 3104 | 2241 | 863.1 |
Prestamo personal
Categorias:
- No
- Si
l1=subset(d,housing=="no", select = c(duration,Status))
l2=subset(d,housing=="yes", select = c(duration,Status))
#Supervivencias
Surv_l1_fit <- survfit(Surv(time = l1$duration, event = l1$Status)~1)
a1 <-data.frame(Tiempo = Surv_l1_fit$time, Surv = Surv_l1_fit$surv)
Surv_l2_fit <- survfit(Surv(time = l2$duration, event = l2$Status)~1)
a2 <-data.frame(Tiempo = Surv_l2_fit$time, Surv = Surv_l2_fit$surv)
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="No"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Si"))+
labs(title="Supervivencia por prestamo ersonal", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
La supervivencia para personas sin prestamos personales es mayor que para aquellas que si cuentan con un prestamo
Pruebas de Hipótesis
prueba1<-LogRank(Surv_l1_fit,Surv_l2_fit)
s1<-c("No")
s2<-c("Si")
pruebas<-data.frame(s1,s2,c(prueba1[1]),c(prueba1[2]))
colnames(pruebas)<-c("S1","S2","Estadistico","Conclusión")
pruebasComo se puede observar, no hay forma de agrupar las supervivencias.
Medias
m1<-Media(Surv_l1_fit)
m2<-Media(Surv_l2_fit)
medias<-data.frame(c("No","Si"),
c(m1[1],m2[1]),
c(m1[2],m2[2]),
c(m1[3],m2[3]),
c(m1[4],m2[4]))
kable(medias , caption = "Media de la duración de la llamada por prestamo personal"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Educacion","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Educacion | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| No | 361 | 44.95 | 347.9 | 374.2 |
| Si | 314 | 41.76 | 301.3 | 326.7 |
La media es mayor para las personas que no cuentan con prestamos personales.
Cuantiles
c1<-Cuantiles(Surv_l1_fit)
c2<-Cuantiles(Surv_l2_fit)
cuantiles<-data.frame(c1$p,c1$c,c2$c)
kable(cuantiles , caption = "Cuantiles de las duraciones"
, align = c('c', 'c', 'c'), col.names=c("Porcentaje","No","Si")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | No | Si |
|---|---|---|
| 25% | 104 | 104 |
| 50% | 187 | 184 |
| 75% | 373 | 328 |
Como se viene observando, las personas que no tienen prestamos personales acumulan mas segundos de duracion en los contactos.
Funciones de Riesgo
ht1<-Surv_l1_fit$n.event/Surv_l1_fit$n.risk
ht2<-Surv_l2_fit$n.event/Surv_l2_fit$n.risk
#grafica de h(t)
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=Surv_l1_fit$time,y=ht1,colour="No"))+
geom_line(aes(x=Surv_l2_fit$time,y=ht2,colour="Si"))+
theme_bw() +
labs(title="h(t) por prestamo personal", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p2<-ggplot()+
geom_step(aes(x=Surv_l1_fit$time,y=Surv_l1_fit$cumhaz, colour="No"))+
geom_step(aes(x=Surv_l2_fit$time,y=Surv_l2_fit$cumhaz, colour="Si"))+
theme_bw()+
labs(title="H(t) por prestamo personal", x="Duración", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
#p1;p2
ggarrange(p1,p2, ncol = 2, nrow = 1)
Las personas que cuentan con prestamos logran acumular más riesgo. Esto se interpreta como que para aquellos clientes que cuenten con un prestamo personal, les será más dificil atender el llamado y realizar el deposito dado que deben realizar otros pagos tambien.
Máxima Segregación
Veamos la máxima distancia entre las dos supervivencia para una probabilidad mayor a \(0\)
#personal
maxima<-Maxima_segregacion(Surv_l1_fit,Surv_l2_fit,"No","Si",0)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.0114 | 3104 | 2241 | 863.1 |
| 0.0114 | 3104 | 2241 | 863.1 |
Contacto
Categorias:
Celular
Telefono
Desconocido
Las funciones de supervivencia son las siguientes:
c1=subset(d,contact=="cellular", select = c(duration,Status))
c2=subset(d,contact=="telephone", select = c(duration,Status))
c3=subset(d,contact=="unknown", select = c(duration,Status))
#Supervivencias
Surv_c1_fit <- survfit(Surv(time = c1$duration, event = c1$Status)~1)
a1 <-data.frame(Tiempo = Surv_c1_fit$time, Surv = Surv_c1_fit$surv)
Surv_c2_fit <- survfit(Surv(time = c2$duration, event = c2$Status)~1)
a2 <-data.frame(Tiempo = Surv_c2_fit$time, Surv = Surv_c2_fit$surv)
Surv_c3_fit <- survfit(Surv(time = c3$duration, event = c3$Status)~1)
a3 <-data.frame(Tiempo = Surv_c3_fit$time, Surv = Surv_c3_fit$surv)
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="Celular"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Telefono"))+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="Desconocido"))+
labs(title="Supervivencia por metodo de contacto", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Pruebas de Hipótesis
Veamos si es posible agrupar algunas supervivencias.
prueba1<-LogRank(Surv_c1_fit,Surv_c2_fit)
prueba2<-LogRank(Surv_c1_fit,Surv_c3_fit)
prueba2<-LogRank(Surv_c2_fit,Surv_c3_fit)
s1<-c("Celular","Celular","Telefono")
s2<-c("Telefono","Desconocido","Desconocido")
pruebas<-data.frame(s1,s2,c(prueba1[1],prueba2[1],prueba3[1]),
c(prueba1[2],prueba2[2],prueba3[2]))
colnames(pruebas)<-c("S1","S2","Estadistico","Conclusión")
pruebasComo se puede observar, podemos agrupar las supervivencias de Telefono y Desconocido.
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="Celular"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Telefono"))+
labs(title="Supervivencia por metodo de contacto", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Se puede observar que en un inicio la supervivencia para celular es mayor pero pasados los 1000 segundos, los papeles se invierten y ahora la supervivencia de telefono es la mayor, esto se mantiene hasta que ambas llegan a \(0\).
Medias
m1<-Media(Surv_c1_fit)
m2<-Media(Surv_c2_fit)
medias<-data.frame(c("Celular","Telefono"),
c(m1[1],m2[1]),
c(m1[2],m2[2]),
c(m1[3],m2[3]),
c(m1[4],m2[4]))
kable(medias , caption = "Media de la duración de la llamada por metodo de contacto"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Estado","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Estado | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| Celular | 353.2 | 26.42 | 343.2 | 363.3 |
| Telefono | 342.0 | 459.07 | 300.0 | 383.9 |
Para celular la media es mayor que para telefono, se puede observar que la varianza de la media para telefono es mucho mayor, lo que representa mayor volatilidad.
Cuantiles
La siguiente tabla muetras los cuantiles \(0.25,0.5\) y \(0.75\) para las duraciones de las llamadas de acuerdo al metodo de contacto.
c1<-Cuantiles(Surv_c1_fit)
c2<-Cuantiles(Surv_c2_fit)
cuantiles<-data.frame(c1$p,c1$c,c2$c)
kable(cuantiles , caption = "Cuantiles de las duraciones"
, align = c('c', 'c', 'c'), col.names=c("Porcentaje","Celular","Telefono")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | Celular | Telefono |
|---|---|---|
| 25% | 107 | 75 |
| 50% | 192 | 158 |
| 75% | 367 | 313 |
Los cuantiles son bastante mayores para celular con respecto a telefono.
Funciones de Riesgo
Funciones de riesgo \(\hat h(t)\) y de riesgo acumulado \(\hat H(t)\) por estado civil.
ht1<-Surv_c1_fit$n.event/Surv_c1_fit$n.risk
ht2<-Surv_c2_fit$n.event/Surv_c2_fit$n.risk
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=Surv_c1_fit$time,y=ht1,colour="Celular"))+
geom_line(aes(x=Surv_c2_fit$time,y=ht2,colour="Telefono"))+
theme_bw() +
labs(title="h(t) por metodo de contacto", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p2<-ggplot()+
geom_step(aes(x=Surv_c1_fit$time,y=Surv_c1_fit$cumhaz, colour="Celular"))+
geom_step(aes(x=Surv_c2_fit$time,y=Surv_c2_fit$cumhaz, colour="Telefono"))+
theme_bw()+
labs(title="H(t) por metodo de contacto", x="Duración", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
#p1;p2
ggarrange(p1,p2, ncol = 2, nrow = 1)
Se observa el comportamiento previsto, en un inicio Telefono acumula mas riesgo, pero pasados los \(1000\) segundos celular empieza a acumular más riesgo y más rapido hasta que la supervivencia llega a \(0\) mas rapido que para Telefono.
Máxima Segregación
Veamos cual es la maxima distancia entre las supervivencias de Celular y Telefono para una probabilidad mayor a \(0.02\)
#contacto
maxima<-Maxima_segregacion(Surv_c1_fit,Surv_c2_fit,"Celular","Telefono",0.02)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.0214 | 1990 | 2264 | 274.1 |
| 0.0214 | 1990 | 2264 | 274.1 |
| 0.0214 | 1990 | 2264 | 274.1 |
| 0.0214 | 1990 | 2264 | 274.1 |
| 0.0214 | 1990 | 2264 | 274.1 |
| 0.0214 | 1990 | 2264 | 274.1 |
Campaña Anterior
Resultado de la campaña de marketing anterior, cuyas categorias son:
- fracaso
- inexistente
- éxito
- otro
Las funciones de supervivencia son las siguientes:
o1=subset(d,poutcome=="failure", select = c(duration,Status))
o2=subset(d,poutcome=="other", select = c(duration,Status))
o3=subset(d,poutcome=="success", select = c(duration,Status))
o4=subset(d,poutcome=="unknown", select = c(duration,Status))
#Supervivencias
Surv_o1_fit <- survfit(Surv(time = o1$duration, event = o1$Status)~1)
a1 <-data.frame(Tiempo = Surv_o1_fit$time, Surv = Surv_o1_fit$surv)
Surv_o2_fit <- survfit(Surv(time = o2$duration, event = o2$Status)~1)
a2 <-data.frame(Tiempo = Surv_o2_fit$time, Surv = Surv_o2_fit$surv)
Surv_o3_fit <- survfit(Surv(time = o3$duration, event = o3$Status)~1)
a3 <-data.frame(Tiempo = Surv_o3_fit$time, Surv = Surv_o3_fit$surv)
Surv_o4_fit <- survfit(Surv(time = o4$duration, event = o4$Status)~1)
a4 <-data.frame(Tiempo = Surv_o4_fit$time, Surv = Surv_o4_fit$surv)
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="Fallo"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Otro"))+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="Exito"))+
geom_step(mapping = aes(x=a4$Tiempo, y=a4$Surv, color="Desconocido"))+
labs(title="Supervivencia por resultado de campaña anterior", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Pruebas de Hipótesis
Veamos si es posible agrupar algunas supervivencias.
prueba1<-LogRank(Surv_o1_fit,Surv_o2_fit)
prueba2<-LogRank(Surv_o1_fit,Surv_o4_fit)
prueba2<-LogRank(Surv_o2_fit,Surv_o4_fit)
s1<-c("Fallo","Fallo","Otro")
s2<-c("Otro","Desconocido","Desconocido")
pruebas<-data.frame(s1,s2,c(prueba1[1],prueba2[1],prueba3[1]),
c(prueba1[2],prueba2[2],prueba3[2]))
colnames(pruebas)<-c("S1","S2","Estadistico","Conclusión")
pruebasComo se puede observar, podemos agrupar fallo con otro y otro con desconocido, dado que el estadistico es menor para Otro-Desconocido agrupamos estos 2.
ggplot()+
geom_step(mapping = aes(x=a1$Tiempo, y=a1$Surv, color="Fallo"))+
geom_step(mapping = aes(x=a2$Tiempo, y=a2$Surv, color="Otro"))+
geom_step(mapping = aes(x=a3$Tiempo, y=a3$Surv, color="Exito"))+
labs(title="Supervivencia por resultado de campaña anterior", x="Duracion", y="S(t)")+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal",
plot.title = element_text(hjust = 0.5))
Veamos que exito posee una mayor supervivencia que los demas resultados.
Medias
m1<-Media(Surv_o1_fit)
m2<-Media(Surv_o2_fit)
m3<-Media(Surv_o3_fit)
medias<-data.frame(c("Fallo","Otro","Exito"),
c(m1[1],m2[1],m3[1]),
c(m1[2],m2[2],m3[2]),
c(m1[3],m2[3],m3[3]),
c(m1[4],m2[4],m3[4]))
kable(medias , caption = "Media de la duración de la llamada por campaña anterior"
, align = c('c', 'c', 'c', 'c','c'), col.names=c("Resultado","Media","Varianza de la media","Limite Inferior","Limite Superior")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Resultado | Media | Varianza de la media | Limite Inferior | Limite Superior |
|---|---|---|---|---|
| Fallo | 293.8 | 41.2 | 281.2 | 306.4 |
| Otro | 327.1 | 135.0 | 304.3 | 349.8 |
| Exito | 902.9 | 2574.9 | 803.5 | 1002.4 |
Sumado a lo anterior, la media de exito tambien es superior y por mucho con respecto a las medias de los otros resultados.
Cuantiles
La siguiente tabla muetras los cuantiles \(0.25,0.5\) y \(0.75\) para las duraciones de las llamadas de acuerdo a cada resultado de la campaña anterior.
c1<-Cuantiles(Surv_o1_fit)
c2<-Cuantiles(Surv_o2_fit)
c3<-Cuantiles(Surv_o3_fit)
cuantiles<-data.frame(c1$p,c1$c,c2$c,c3$c)
kable(cuantiles , caption = "Cuantiles de las duraciones"
, align = c('c', 'c', 'c', 'c'), col.names=c("Porcentaje","Fracaso","Otro","Exito")) %>%
kable_styling(full_width = F, bootstrap_options = c("condensed"))| Porcentaje | Fracaso | Otro | Exito |
|---|---|---|---|
| 25% | 106 | 100 | 239 |
| 50% | 183 | 192 | 638 |
| 75% | 329 | 383 | 1563 |
Como se esperaba, exito acumula más segundos de llamada en los cuantiles \(0.25,0.5\) y \(0.75\) con respecto a los otros resultados.
Funciones de Riesgo
Funciones de riesgo \(\hat h(t)\) y de riesgo acumulado \(\hat H(t)\) por estado civil.
ht1<-Surv_o1_fit$n.event/Surv_o1_fit$n.risk
ht2<-Surv_o2_fit$n.event/Surv_o2_fit$n.risk
ht3<-Surv_o3_fit$n.event/Surv_o3_fit$n.risk
#grafica de h(t)
p1<-ggplot() +
geom_line(aes(x=Surv_o1_fit$time,y=ht1,colour="Fracaso"))+
geom_line(aes(x=Surv_o2_fit$time,y=ht2,colour="Otro"))+
geom_line(aes(x=Surv_o3_fit$time,y=ht3,colour="Exito"))+
theme_bw() +
labs(title="h(t) por resultado anterior", x="Tiempo", y=TeX("$\\hat{h}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
p2<-ggplot()+
geom_step(aes(x=Surv_o1_fit$time,y=Surv_o1_fit$cumhaz, colour="Fracaso"))+
geom_step(aes(x=Surv_o2_fit$time,y=Surv_o2_fit$cumhaz, colour="Otro"))+
geom_step(aes(x=Surv_o3_fit$time,y=Surv_o3_fit$cumhaz, colour="Exito"))+
theme_bw()+
labs(title="H(t) por resultado anterior", x="Duración", y=TeX("$\\hat{H}(t)$"))+
theme(plot.title = element_text(hjust = 0.5))+
scale_colour_discrete("")+
theme(legend.position = "bottom", legend.direction = "horizontal")
#p1;p2
ggarrange(p1,p2, ncol = 2, nrow = 1)
Como era de esperarse exito acumula la menor cantidad de riesgo, en cambio fracaso es el grupo que mas riesgo acumula.
La interpretacion de esto es que si en una campaña anterior se tuvo exito con un cliente, se espera que esto vuelva a ocurrir, generando asi un buen historial. Dado que esta es la primer campaña para una gran cantidad de clientes, aun se necesita recopilar mas informacion para conocer de antemano si un cliente podrá realizar el depósito en nuevas campañas.
Máxima Segregación
Comparemos las supervivencias de Fallo y Exito para una probabilidad mayor a \(0.25\)
#anterior
maxima<-Maxima_segregacion(Surv_o1_fit,Surv_o3_fit,"Fallo","Exito",0.25)
maxima$grafica
maxima$tabla| Probabilidad | Tiempo 1 | Tiempo 2 | Distancia maxima |
|---|---|---|---|
| 0.2506 | 328 | 1596 | 1268 |
Riesgos Proporcionales
Modelos
A continuación buscaremos covariables a partir de un modelo de Cox, esto con el fin de ver si alguna de estas es capaz de aportar mas informacion al modelo de supervivencia.
d4<-d
#Agrupamos los datos de acuerdo a los resultados de las pruebas de hipotesis
for(i in 1:length(d4$duration)){
if(d4$job[i]=="admin." | d4$job[i]== "technician"){
d4$job[i]= "blue-collar"
}
if(d4$job[i]=="self-employed" | d4$job[i]== "student" | d4$job[i]== "unemployed"){
d4$job[i]= "services"
}
if(d4$education[i]=="unknown"){
d4$education[i]= "primary"
}
if(d4$contact[i]=="unknown"){
d4$contact[i]= "telephone"
}
if(d4$poutcome[i]=="unknown"){
d4$poutcome[i]= "other"
}
}
#generamos el modelo de Cox
proporcional <- coxph(Surv(duration, Status) ~ ., d4)
summary(proporcional)## Call:
## coxph(formula = Surv(duration, Status) ~ ., data = d4)
##
## n= 45211, number of events= 39922
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 1.85e-03 1.00e+00 6.13e-04 3.02 0.00255 **
## jobentrepreneur -1.35e-02 9.87e-01 2.86e-02 -0.47 0.63597
## jobhousemaid -7.20e-03 9.93e-01 3.13e-02 -0.23 0.81795
## jobmanagement 9.96e-03 1.01e+00 1.66e-02 0.60 0.54940
## jobretired -1.12e-01 8.94e-01 2.78e-02 -4.01 6.1e-05 ***
## jobservices -1.14e-02 9.89e-01 1.41e-02 -0.81 0.41651
## jobunknown -5.47e-03 9.95e-01 6.39e-02 -0.09 0.93168
## maritalmarried -1.72e-02 9.83e-01 1.63e-02 -1.06 0.28952
## maritalsingle -1.58e-03 9.98e-01 1.88e-02 -0.08 0.93333
## educationsecondary -8.66e-03 9.91e-01 1.37e-02 -0.63 0.52866
## educationtertiary 2.85e-03 1.00e+00 1.77e-02 0.16 0.87220
## defaultyes 2.71e-02 1.03e+00 3.68e-02 0.73 0.46281
## balance -6.01e-06 1.00e+00 1.86e-06 -3.24 0.00121 **
## housingyes -3.09e-03 9.97e-01 1.23e-02 -0.25 0.80207
## loanyes 2.72e-02 1.03e+00 1.37e-02 1.98 0.04750 *
## contacttelephone 1.99e-03 1.00e+00 1.51e-02 0.13 0.89506
## day 6.80e-03 1.01e+00 7.16e-04 9.49 < 2e-16 ***
## monthaug 3.05e-01 1.36e+00 2.69e-02 11.35 < 2e-16 ***
## monthdec 4.34e-02 1.04e+00 9.66e-02 0.45 0.65367
## monthfeb 3.03e-01 1.35e+00 3.13e-02 9.68 < 2e-16 ***
## monthjan 3.14e-02 1.03e+00 3.65e-02 0.86 0.38926
## monthjul 1.21e-01 1.13e+00 2.55e-02 4.74 2.2e-06 ***
## monthjun 2.90e-01 1.34e+00 2.90e-02 9.99 < 2e-16 ***
## monthmar 4.00e-01 1.49e+00 6.97e-02 5.73 1.0e-08 ***
## monthmay 1.32e-01 1.14e+00 2.42e-02 5.46 4.7e-08 ***
## monthnov 2.02e-01 1.22e+00 2.72e-02 7.42 1.2e-13 ***
## monthoct 1.80e-01 1.20e+00 5.43e-02 3.31 0.00094 ***
## monthsep 2.13e-01 1.24e+00 6.14e-02 3.46 0.00054 ***
## campaign 3.64e-02 1.04e+00 1.56e-03 23.29 < 2e-16 ***
## pdays 3.27e-04 1.00e+00 8.69e-05 3.77 0.00016 ***
## previous 1.12e-03 1.00e+00 1.92e-03 0.58 0.56031
## poutcomeother -4.40e-02 9.57e-01 2.55e-02 -1.73 0.08415 .
## poutcomesuccess -1.20e-01 8.87e-01 4.67e-02 -2.57 0.01008 *
## yyes -1.86e+01 8.53e-09 9.42e+01 -0.20 0.84356
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.00e+00 9.98e-01 1.00e+00 1.00e+00
## jobentrepreneur 9.87e-01 1.01e+00 9.33e-01 1.04e+00
## jobhousemaid 9.93e-01 1.01e+00 9.34e-01 1.06e+00
## jobmanagement 1.01e+00 9.90e-01 9.78e-01 1.04e+00
## jobretired 8.94e-01 1.12e+00 8.47e-01 9.45e-01
## jobservices 9.89e-01 1.01e+00 9.62e-01 1.02e+00
## jobunknown 9.95e-01 1.01e+00 8.78e-01 1.13e+00
## maritalmarried 9.83e-01 1.02e+00 9.52e-01 1.01e+00
## maritalsingle 9.98e-01 1.00e+00 9.62e-01 1.04e+00
## educationsecondary 9.91e-01 1.01e+00 9.65e-01 1.02e+00
## educationtertiary 1.00e+00 9.97e-01 9.69e-01 1.04e+00
## defaultyes 1.03e+00 9.73e-01 9.56e-01 1.10e+00
## balance 1.00e+00 1.00e+00 1.00e+00 1.00e+00
## housingyes 9.97e-01 1.00e+00 9.73e-01 1.02e+00
## loanyes 1.03e+00 9.73e-01 1.00e+00 1.06e+00
## contacttelephone 1.00e+00 9.98e-01 9.73e-01 1.03e+00
## day 1.01e+00 9.93e-01 1.01e+00 1.01e+00
## monthaug 1.36e+00 7.37e-01 1.29e+00 1.43e+00
## monthdec 1.04e+00 9.58e-01 8.64e-01 1.26e+00
## monthfeb 1.35e+00 7.39e-01 1.27e+00 1.44e+00
## monthjan 1.03e+00 9.69e-01 9.61e-01 1.11e+00
## monthjul 1.13e+00 8.86e-01 1.07e+00 1.19e+00
## monthjun 1.34e+00 7.48e-01 1.26e+00 1.41e+00
## monthmar 1.49e+00 6.71e-01 1.30e+00 1.71e+00
## monthmay 1.14e+00 8.76e-01 1.09e+00 1.20e+00
## monthnov 1.22e+00 8.17e-01 1.16e+00 1.29e+00
## monthoct 1.20e+00 8.36e-01 1.08e+00 1.33e+00
## monthsep 1.24e+00 8.08e-01 1.10e+00 1.40e+00
## campaign 1.04e+00 9.64e-01 1.03e+00 1.04e+00
## pdays 1.00e+00 1.00e+00 1.00e+00 1.00e+00
## previous 1.00e+00 9.99e-01 9.97e-01 1.00e+00
## poutcomeother 9.57e-01 1.04e+00 9.10e-01 1.01e+00
## poutcomesuccess 8.87e-01 1.13e+00 8.09e-01 9.72e-01
## yyes 8.53e-09 1.17e+08 6.15e-89 1.18e+72
##
## Concordance= 0.638 (se = 0.002 )
## Likelihood ratio test= 22649 on 34 df, p=<2e-16
## Wald test = 1133 on 34 df, p=<2e-16
## Score (logrank) test = 13818 on 34 df, p=<2e-16Como podemos observar a partir del nivel de significancia de las variables, las que son estadisticamente significativas para el modelo son:
age
job: solo para el tipo retired
balance
loan: para yes
day
month
campaign
pdays
poutcome
Los resultados de Likelihood ratio test, Wald test y Score test permiten rechazar la hipótesis nula de que todos los coeficientes son cero. La concordancia es del \(63.8\%\) lo que indica que los datos son descritos de buena manera por el modelo.
Procedemos a realizar un siguiente modelo tomando en cuenta solo las variables significativas:
proporcional2 <- coxph(Surv(duration, Status) ~
age + job + balance + loan + day + month + campaign + pdays + poutcome, d4)
summary(proporcional2)## Call:
## coxph(formula = Surv(duration, Status) ~ age + job + balance +
## loan + day + month + campaign + pdays + poutcome, data = d4)
##
## n= 45211, number of events= 39922
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 3.16e-03 1.00e+00 5.39e-04 5.87 4.5e-09 ***
## jobentrepreneur -1.90e-02 9.81e-01 2.81e-02 -0.68 0.49883
## jobhousemaid 3.76e-02 1.04e+00 3.10e-02 1.21 0.22455
## jobmanagement -1.63e-02 9.84e-01 1.33e-02 -1.22 0.22255
## jobretired -2.28e-01 7.96e-01 2.78e-02 -8.21 2.2e-16 ***
## jobservices -4.27e-02 9.58e-01 1.40e-02 -3.05 0.00226 **
## jobunknown 2.96e-02 1.03e+00 6.35e-02 0.47 0.64075
## balance -1.13e-05 1.00e+00 1.91e-06 -5.92 3.2e-09 ***
## loanyes 7.47e-02 1.08e+00 1.37e-02 5.46 4.6e-08 ***
## day 5.72e-03 1.01e+00 6.98e-04 8.19 2.6e-16 ***
## monthaug 3.21e-01 1.38e+00 2.54e-02 12.62 < 2e-16 ***
## monthdec -3.19e-01 7.27e-01 9.61e-02 -3.32 0.00091 ***
## monthfeb 3.07e-01 1.36e+00 3.08e-02 9.97 < 2e-16 ***
## monthjan 1.83e-01 1.20e+00 3.58e-02 5.11 3.3e-07 ***
## monthjul 1.71e-01 1.19e+00 2.50e-02 6.84 8.2e-12 ***
## monthjun 3.36e-01 1.40e+00 2.62e-02 12.85 < 2e-16 ***
## monthmar -2.26e-01 7.98e-01 6.94e-02 -3.26 0.00112 **
## monthmay 2.61e-01 1.30e+00 2.27e-02 11.51 < 2e-16 ***
## monthnov 3.02e-01 1.35e+00 2.69e-02 11.21 < 2e-16 ***
## monthoct -1.55e-01 8.56e-01 5.36e-02 -2.89 0.00380 **
## monthsep -2.01e-01 8.18e-01 6.08e-02 -3.30 0.00097 ***
## campaign 4.09e-02 1.04e+00 1.55e-03 26.36 < 2e-16 ***
## pdays 1.04e-04 1.00e+00 8.11e-05 1.28 0.19964
## poutcomeother -6.64e-02 9.36e-01 2.50e-02 -2.66 0.00793 **
## poutcomesuccess -1.12e+00 3.25e-01 4.66e-02 -24.12 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.003 0.997 1.002 1.004
## jobentrepreneur 0.981 1.019 0.929 1.037
## jobhousemaid 1.038 0.963 0.977 1.103
## jobmanagement 0.984 1.016 0.959 1.010
## jobretired 0.796 1.256 0.754 0.841
## jobservices 0.958 1.044 0.932 0.985
## jobunknown 1.030 0.971 0.910 1.166
## balance 1.000 1.000 1.000 1.000
## loanyes 1.078 0.928 1.049 1.107
## day 1.006 0.994 1.004 1.007
## monthaug 1.378 0.726 1.311 1.449
## monthdec 0.727 1.376 0.602 0.878
## monthfeb 1.359 0.736 1.279 1.443
## monthjan 1.200 0.833 1.119 1.287
## monthjul 1.186 0.843 1.129 1.246
## monthjun 1.400 0.714 1.330 1.474
## monthmar 0.798 1.254 0.696 0.914
## monthmay 1.298 0.770 1.242 1.357
## monthnov 1.352 0.740 1.282 1.425
## monthoct 0.856 1.168 0.771 0.951
## monthsep 0.818 1.222 0.726 0.922
## campaign 1.042 0.960 1.039 1.045
## pdays 1.000 1.000 1.000 1.000
## poutcomeother 0.936 1.069 0.891 0.983
## poutcomesuccess 0.325 3.075 0.297 0.356
##
## Concordance= 0.575 (se = 0.002 )
## Likelihood ratio test= 2602 on 25 df, p=<2e-16
## Wald test = 2311 on 25 df, p=<2e-16
## Score (logrank) test = 2409 on 25 df, p=<2e-16Para este modelo todas las variables son significativas, se reduce la concordancia a un \(58\%\), los resultados siguen señalando que se rechaza la hipotesis nula en la que los coeficientes son iguales a cero.
Comparemos los \(AIC\) de ambos modelos:
mod<-data.frame(c("Modelo 1","Modelo 2"),
c(AIC(proporcional),AIC(proporcional2)))
colnames(mod)<-c("Modelo","AIC")
modAunque el \(AIC\) sea mayor en el modelo 2, nos quedamos con este pues para el siguiente paso necesitamos que el resultado del modelo sea no singular (cosa que pasó con el modelo 1)
Lo siguiente es verificar los supuestos del modelo de Cox lo cual se puede hacer mediante pruebas estadísticas basadas en los residuos de Schoenfeld escalados.
Usamos la función cox.zph(). Para cada covariable, la función cox.zph() correlaciona el conjunto correspondiente de residuos de Schoenfeld escalados con el tiempo, para probar la independencia entre los residuos y el tiempo. Además, realiza una prueba global para el modelo en su conjunto.
test <- cox.zph(proporcional2)
test## chisq df p
## age 15.418 1 8.6e-05
## job 38.886 6 7.5e-07
## balance 5.050 1 0.02463
## loan 11.005 1 0.00091
## day 120.134 1 < 2e-16
## month 264.983 11 < 2e-16
## campaign 789.465 1 < 2e-16
## pdays 0.728 1 0.39361
## poutcome 14.225 2 0.00081
## GLOBAL 1200.802 25 < 2e-16Se observa que el modelo no es adecuado pues no cumple el supuesto principal para el modelo de Cox, se observa tambien que la unica variable con las que se acepta el supuesto es pdays, veamos que sucede al usar un modelo con unicamente esta variable.
proporcional3 <- coxph(Surv(duration, Status) ~ pdays, d4)
summary(proporcional3)## Call:
## coxph(formula = Surv(duration, Status) ~ pdays, data = d4)
##
## n= 45211, number of events= 39922
##
## coef exp(coef) se(coef) z Pr(>|z|)
## pdays -4.75e-04 1.00e+00 5.31e-05 -8.94 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## pdays 1 1 0.999 1
##
## Concordance= 0.514 (se = 0.001 )
## Likelihood ratio test= 83.3 on 1 df, p=<2e-16
## Wald test = 79.9 on 1 df, p=<2e-16
## Score (logrank) test = 80 on 1 df, p=<2e-16test <- cox.zph(proporcional3)
test## chisq df p
## pdays 0.0324 1 0.86
## GLOBAL 0.0324 1 0.86Podemos ver que la variable es significativa, se sigue rechazando la hipotesis nula de que el coeficiente es cero. la concordancia nos dice que se rechaza el ajuste del modelo por parte de los datos.
Este modelo cumple con el riesgo proporcional.
mod<-data.frame(c("Modelo 1","Modelo 2","Modelo 3"),
c(AIC(proporcional),AIC(proporcional2),AIC(proporcional3)))
colnames(mod)<-c("Modelo","AIC")
modPese a que el \(AIC\) del modelo actual sea el mayor, esto se contrasta con que es el unico modelo que cumple el supuesto de proporcionalidad por lo que nos quedamos con este.