Tutorial Análise de sobrevivencia

1 Introdução

Em diversas áreas do conhecimento o objeto de estudo é o tempo até a ocorrência de um determinado evento de interesse – por exemplo, a morte de uma pessoa ou a falha de uma máquina. Situações desse tipo são comuns também em análise de crédito. Frequentemente queremos estimar o tempo até que um contrato de financiamento com parcelas atrasadas se torne adimplente ou o tempo até que ele se torne inadimplente pela primeira vez. Note que, em ambas as situações descritas, o evento esperado pode nunca ocorrer. Quando é esse o caso, a observação é considerada censurada. As técnicas estatísticas utilizadas para se analisar dados de duração na presença de censura fazem parte de uma subárea conhecida como Análise de Sobrevivência.

Dado um contrato de financiamento qualquer, nesta nota vamos utilizar a Análise de Sobrevivência para examinar três situações:

1. O tempo até recuperação de uma parcela que entrou em inadimplência em uma determinada data;

2. O tempo até recuperação de um contrato que entrou em indimplência em uma determinada data;

3. O tempo tempo até a primeira inadimplência de um contrato.

Para cada um dos casos elencados acima iremos calcular curvas de sobrevivência e estimar probabilidades condicionais de sobrevida. Todas as análises descritas neste texto serão demonstradas com dados da base de contratos pro-soluto da Direcional no período que vai de 2017 a 2022. Na próxima seção faremo uma breve introdução aos principais tópicos de análise de sobrevivência. As seções seguintes mostram as aplicações dos modelos desenvolvidos à análise de crédito.

2 Análise de Sobrevivência

Um modelo de análise de sobrevivência, tipicamente, contém três elementos básicos descritos a seguir:

• Tempo: tempo desde o início do período de observação (que pode ser o dia de assinatura de um contrato, por exemplo) até o evento de interesse ou até a censura. É representado como uma variável contínua que toma valores em \(R^{+}\).
• Evento: desfecho de interesse. Pode ser, por exemplo, o momento em que um cliente entra em inadimplência pela primeira vez. É representado por uma variável dicotômica.
• Censura: A censura ocorre quando a informação a respeito do desfecho de determinado indivíduo não é conhecida. Ela acontece, por exemplo, em um experimento que termina antes que um dos participantes apresente o evento de interesse. Existem basicamente três tipos de censura:
  • Censura à direita: é o tipo mais comum de censura. Ocorre quando o sujeito deixa o experimento no tempo \(t=T_{1}\), porém o evento de interesse ocorre em \(t>T_{1}\).
  • Censura à esquerda: ocorre quando o evento de interesse já aconteceu antes do início do período de observação. Exemplo: um estudo para se determinar o tempo até alfabetização no qual algumas crianças já começam o experimento sabendo ler e escrever.
  • Censura Intervalar: ocorre quando o instante exato do evento de interesse não é conhecido. Sabe-se apenas que ele aconteceu entre dois momentos de observação.

Nota: Há uma diferença entre observações censuradas e truncadas. Em certos experimentos os participantes são incluídos apenas se seus tempos de vida estiverem acima (truncagem à esquerda) ou abaixo (truncagem à direita) de um limiar. Por exemplo, um ensaio clínico para testar determinada droga pode incluir somente pessoas com idades acima de um certo patamar. Indivíduos truncados não fazem parte do estudo e, portanto, não contribuem com nenhuma informação para os pesquisadores. A truncagem é um efeito do design do experimento.

2.1 Função de Sobrevivência

A Função de Sobrevivência (ou confiabilidade) \(S(t)\) dá a probabilidade de que o objeto de estudo (um indivíduo ou uma máquina, por exemplo) sobreviva além de um tempo \(t\). Seja \(T\) uma variável aleatória contínua positiva que denota o tempo até o evento de interesse, a Função de Sobrevivência pode ser representada matematicamente como o complementar da função de probabilidade acumulada \(F(t)=P\{t<T\}\), ou seja: \[S(t)=P(\{T>t\})=1-F(t)=\int_{t}^{\infty}f(x)dx\]

2.2 Função Risco

A Função Risco (ou taxa de falha) representa a taxa instantânea de ocorrência do evento de interesse. Ela é dada por:

\[\lambda (t)= \lim_{dt \rightarrow 0}\frac{\text{Pr}\{t \leq T+dt|T\geq t\}}{dt}\] O numerador da expressão acima representa a probabilidade condicional de falha entre \(t\) e \(t+dt\) dado que a falha não ocorreu até o instante \(t\). Essa probabilidade condicional pode ser reescrita como a razão entre a probabilidade de que o evento ocorra no intervalo \([t, t+dt)\) e a probabilidade de que ele ocorra em \(T>t\). Note que a primeira probabilidade é a função densidade de probabilidade dos tempos de falha \(f(t)\) multiplicada pelo intervalo infinitesimal \(dt\). Já a segunda probabilidade é simplesmente a Função de sobrevivência \(S(t)\). Sendo assim, a função risco pode ser expressa como:

\[\lambda(t)=\frac{f(t)}{S(t)}. \] Como \(f(t)=-\frac{dS(t)}{dt}\), a equação acima é equivalente a \[\lambda(t)=-\frac{d}{dt}log(S(t)).\]

Importante: A função risco não é uma função densidade de probabilidade.

Exemplo:

Seja uma Função de Sobrevivência Exponencial, ou seja:

\[S(t)=\exp(-\lambda t).\]

A função densidade de probabilidade correspondente será dada por \(f(t)=\lambda \exp(-\lambda t)\). Dessa forma, temos:

\[\lambda(t)=\frac{f(t)}{S(t)}=\frac{\lambda \exp(-\lambda t)}{\exp(-\lambda t)}=\lambda. \] Portanto, o risco associado à curva de sobrevivência exponencial é constante. Para calcular o tempo esperado até o evento podemos fazer:

\[E(t)=\int_{0}^{\infty}tf(t)dt=\frac{1}{\lambda}.\] Alternativamente, usando integração por partes é possível provar que

\[E(t)=\int_{0}^{\infty}S(t)dt.\] Ou seja: o tempo médio até o evento é área sob a curva de sobrevivência.

2.3 O estimador de Kaplan Meier

Na maior parte das vezes, a distribuição dos tempos até o evento de interesse não é conhecida. Nessas situações, a curva de sobrevivência pode ser estimada empiricamente através do estimador de Kaplan-Meier, também conhecido como estimador de produto-limite, conforme a equação abaixo:

\[\hat{S(t)}=\prod_{i:t_{i}<t}\left(1-\frac{d_{i}}{n_{i}}\right),\] em que \(t_{i}\) representa um instante no qual pelo menos um evento ocorreu; \(d_{i}\) o número de eventos ocorridos em \(t_{i}\) e \(n_{i}\) o número de indíviduos em risco – ou seja, que sobreviveram até \(t_{i}\) e não foram censurados. Quando a função de sobrevivência estimada \(\hat{S(t)}\) é traçada contra o tempo obtemos a curva de Kaplan-Meier.

A utilização do estimador de Kaplan-Meier requer que as premissas a seguir sejam observadas:

• A probabilidade de ocorrência do evento de interesse é a mesma para observações censuradas e não censuradas;
• O momento de entrada do indivíduo no experimento não altera a probabilidade de ocorrência do evento de interesse;
• Os indivíduos são independentes;
• Os tempos de censura são independentes dos tempos de sobrevivência;
• A probabilidade de sobrevivência é constante entre dois tempos de sobrevivência;
• Não há excesso de observações censuradas.

A violação de uma ou mais dessas condições pode comprometer a interpretação da curva de Kaplan-Meier. Especialmente se a amostra for pequena.

2.4 Modelo de Aceleração de Falha

Até o momento consideramos que os tempos até o evento de interesse de todos os indivíduos do experimento são determinados pela mesma função de sobrevivência. Em dados reais essa suposição dificilmente será verdadeira uma vez que, na maior parte dos casos, os tempos de sobrevivência estarão condicionados a um vetor de covariáveis que pode acelerar ou atrasar a ocorrência do evento. Por esse motivo, modelos paramétricos de regressão capazes de capturar tais efeitos são chamados de modelos de aceleração de falha ou AFT (Accelerated Failures Times).

Seja \(T_{i}\) uma variável aleatória que representa o tempo até falha associado ao indivíduo \(i\). Como \(T_{i}>0\), podemos modelar \(T_{i}\) usando a especificação abaixo:

\[\log(T_{i})=\mathbf{x}_{i}^{\prime}\beta+\epsilon_{i},\] onde \(\mathbf{x}_{i}^{\prime}\) é uma matriz de covariáveis e \(\epsilon_{i}\) uma erro de média zero e variância \(\sigma^{2}\). Exponenciando ambos os lados da equação temos:

\[T_{i}=\exp(\mathbf{x}_{i}^{\prime}\beta)\lambda_{0}\] onde \(\lambda_{0}=\exp(\epsilon_{i})\). Para interpretar o modelo suponha uma matriz de covariáveis composta apenas por uma variável explicativa dicotômica (i.e. que assume valores zero ou um conforme o grupo ao qual determinada observação pertence). Se o coeficiente correspondente for, por exemplo, \(\log(2)\), então o tempo até falha no grupo em que \(X_{i}=1\) será do dobro do observado no grupo \(X_{i}=0\).

2.5 Modelo de Cox

Uma abordagem alternativa ao AFT é a classe de modelos de riscos proporcionais proposta por Cox (1972). Neste tipo de modelo, o efeito das covariáveis é multiplicativo com relação à função risco. Seja \(x_{i}\) um vetor de covariáveis associadas ao indivíduo \(i\), de acordo com o modelo de Cox, a taxa de risco do sujeito no tempo \(t\) pode ser escrita como:

\[\lambda_{i}(t \mid X_{i})=\lambda_{0}(t)\exp(x_{i}^{\prime}\beta),\] onde \(\lambda_{0}(t)\) é uma taxa de risco de referência que não depende do indivíduo \(i\). O termo \(\exp(x_{i}^{\prime}\beta)\) funciona como um fator de escala que indica se as covariáveis aumentam o diminuem o risco do sujeito \(i\).

Uma das vantagens do modelo de Cox sobre o AFT é que ele não requer se conheça de antemão a distribuição dos tempos de sobrevivência \(T\). A estimação da taxa de risco de referência é feita de maneira não paramétrica através do estimador de Breslow. Já a estimação dos parâmetros \(\beta\) é feita via maximização da função de verossimilhança parcial. Dessa forma, o modelo de Cox é semiparamétrico.

A interpretação dos coeficientes em uma regressão de Cox é bastante direta. Seja um modelo com apenas uma covariável, o efeito de um incremento de uma unidade da variável explicativa \(X\) será dado por:

\[\begin{aligned} \lambda(t \mid X+1) &=\lambda_{0}(t)\exp(\beta(x+1))\\ \lambda(t \mid X+1)&=\lambda_{0}(t)\exp(\beta x+\beta)\\ \lambda(t \mid X+1)&=\lambda_{0}(t)\exp(\beta x)\exp(\beta)\\ \lambda(t \mid X+1)&=\lambda(t\mid X)\exp(\beta)\\ \frac{\lambda(t \mid X+1)}{\lambda(t\mid X)} &=\exp(\beta) \end{aligned}\]

Ou seja: um incremento de uma unidade em uma variável explicativa \(j\) produz um aumento (ou diminuição) no risco de um indivíduo da ordem de \(\exp(\beta_{j})\).

2.6 Probabilidade de sobrevida

Em algumas situações é interessante saber a probabilidade condicional de ocorrência do evento de interesse dado que o indivíduo já sobreviveu por um determinado período. Essa probabilidade de sobrevida pode ser calculada através da expressão abaixo

\[ CS(x \mid T) = \frac{S(x+T)}{S(T)},\] em que \(x\) é o tempo adicional até a ocorrência do evento de interesse e \(T\) o tempo já decorrido sem que o evento tenha acontecido. As funções de sobrevivência do lado direito da equação podem ser estimadas por métodos bem estabelecidos em análise de sobrevivência como o estimador de Kaplan Meier ou o Modelo de Cox.

3 Tempo até inadimplência

Os conceitos introduzidos na seção anterior podem ser aplicados para se estimar o tempo até inadimplência em uma carteira de créditos imobiliários. Para essa análise usamos a premissa de que todos os contratos (os indivíduos do modelo) vão entrar em default ao menos uma vez. Se esse evento acontecer antes que todas as parcelas tenham sido quitadas, então dizemos que houve um default real; já se o evento não ocorrer ao longo do tempo de vida do contrato, o default é dito virtual. Do ponto de vista da análise de dados, os defaults virtuais são observações censuradas. Nas próximas seções vamos apresentar a modelagem em detalhes utilizando como exemplo a base de dados Pro-Soluto da Direcional Engenharia e o pacote Survival do R.

3.1 Preparação dos dados

Vamos utilizar os seguintes pacotes nesta seção.

library(data.table)
library(magrittr)
library(lubridate)
library(janitor)
library(moments)
library(arrow)
library(dplyr)
library(glue)
library(purrr)
library(zoo)
library(stringr)
library(openxlsx)
library(survival)
library(survminer)
library(locfit)
library(ggfortify)
library(condSURV)
library(plotly)
library(tidyr)
library(readxl)
library(genderBR)
library(ciTools)
library(arules)
library(knitr)
library(casebase)
library(visreg)
library(splines)
library(kableExtra)
library(hence.aws.s3)
library(hence.plots)

Ao longo de toda a análise usaremos a base de dados Pro Soluto da Direcional. O código abaixo carrega os dados da AWS.

source("../R/global-s3.R")
dados <-
  hence_s3_read_file(
    p_bucket = g_bucket,
    p_group = g_group,
    p_project = g_project,
    p_dat_folder = g_dir_raw,
    p_file_name = "direcional_20220321-carteira_completa.parquet",
    p_method = "paws",
    p_file_type = "parquet"
  ) %>%
  setDT() %>%
  clean_names()

Inicialmente, iremos selecionar as variáveis de interesse para a modelagem e criar dummies referentes às safras dos contratos.

cart <- dados %>%
  select("id_venda", "uf_empreendimento", "data_venda",
         "data_vencimento", "data_habite_se", "data_recebimento",
         "categoria_produto",
 "valor_venda", "valor_nominal_das_parcelas", "tipo_venda", "nome_cliente") %>%
  mutate(safra = as.yearqtr(data_venda)) %>%
  mutate(safra = as.character(safra))

A tabela a seguir resume algumas características da carteira.

Empresa:	Direcional
Data de Referencia:	2022-03-17
Data minima de venda	2017-01-01
Data maxima de venda	2022-03-09
Número de Contratos	50.867

Após o pré-processamento inicial, a base deve ter a estrutura abaixo. Para cada id_venda podemos visualizar todas as parcelas – quitadas ou em aberto – correspondentes.

head(cart)

        id_venda uf_empreendimento data_venda data_vencimento data_habite_se
          <char>            <char>     <Date>          <Date>         <Date>
1: 103-273VA-111                AM 2017-07-28      2017-09-15     2019-02-22
2: 103-273VA-111                AM 2017-07-28      2018-06-15     2019-02-22
3: 103-273VA-111                AM 2017-07-28      2018-07-20     2019-02-22
4: 103-273VA-111                AM 2017-07-28      2018-02-15     2019-02-22
5: 103-273VA-111                AM 2017-07-28      2017-10-15     2019-02-22
6: 103-273VA-111                AM 2017-07-28      2017-11-15     2019-02-22
   data_recebimento categoria_produto valor_venda valor_nominal_das_parcelas
             <Date>            <char>       <num>                      <num>
1:       2017-09-13 MCMV Faixas 2 e 3      179000                    1950.08
2:       2018-06-15 MCMV Faixas 2 e 3      179000                    1000.00
3:       2018-09-06 MCMV Faixas 2 e 3      179000                     842.93
4:       2018-03-07 MCMV Faixas 2 e 3      179000                     829.41
5:       2017-10-06 MCMV Faixas 2 e 3      179000                     829.41
6:       2017-11-06 MCMV Faixas 2 e 3      179000                     829.41
            tipo_venda              nome_cliente   safra
                <char>                    <char>  <char>
1: Financiamento Banco EDINEY FRANÇA DE OLIVEIRA 2017 Q3
2: Financiamento Banco EDINEY FRANÇA DE OLIVEIRA 2017 Q3
3: Financiamento Banco EDINEY FRANÇA DE OLIVEIRA 2017 Q3
4: Financiamento Banco EDINEY FRANÇA DE OLIVEIRA 2017 Q3
5: Financiamento Banco EDINEY FRANÇA DE OLIVEIRA 2017 Q3
6: Financiamento Banco EDINEY FRANÇA DE OLIVEIRA 2017 Q3

No entanto, o interesse está em modelar o tempo até a primeira inadimplência de cada contrato. Sendo assim, queremos um conjunto de dados em que cada id_venda apareça associado ao tempo até a primeira inadimplência ou até a censura. O código abaixo faz exatamente isso: para cada contrato calculamos a diferença, em dias, entre a data de venda e a data em que houve o primeiro atraso no pagamento. Caso determinado contrato não tenha atrasado nenhuma vez até a data de referência, ele será considerado uma observação censurada e, então, calcularemos a diferença em dias entre a data de referência e a data de venda. Note ainda que a definição de inadimplência pode variar – em algumas situações um atraso de um dia será considerado inadimplência e, em outras, um atraso de até 90 dias pode ser tolerado. Adicionalmente, para a modelagem de sobrevivência precisamos de covariáveis e de uma variável binária que atribui 0 para observações censuradas e 1 caso contrário. Neste exemplo, as covariáveis utilizadas são UF, Categoria do Produto, ltv, genero do cliente e safra do contrato.

data_ref <- max(cart$data_recebimento, na.rm = TRUE)

N <- 1 # Número de dias em atraso para que a parcela seja considerada inadimplente

cart_inad <- cart %>%
  .[order(data_vencimento), .SD, by = "id_venda"] %>%
  .[, dias_atraso := ifelse(!is.na(data_recebimento),
    data_recebimento - data_vencimento,
    data_ref - data_vencimento
  )] %>%
  .[, gender := get_gender(nome_cliente)] %>%
  group_by(id_venda) %>%
  filter(!all(!is.na(data_recebimento) & dias_atraso <= 0)) %>%
  mutate(data_inad = ifelse(dias_atraso <= N, NA, data_vencimento + N)) %>%
  summarise(
    data_inad = ifelse(all(dias_atraso <= N), data_ref,
      min(data_inad, na.rm = T)
    ),
    data_venda = first(data_venda),
    uf_empreendimento = first(uf_empreendimento),
    categoria_produto = first(categoria_produto),
    ltv = sum(valor_nominal_das_parcelas) / first(valor_venda),
    tipo_venda = first(tipo_venda),
    gender = first(gender),
    safra = first(safra)
  ) %>%
  mutate(data_inad = (as_date(data_inad, origin = lubridate::origin))) %>%
  mutate(tempo_ate_inad = as.numeric(data_inad - data_venda)) %>%
  mutate(censura = ifelse(data_inad == data_ref, 0, 1)) %>%
  mutate(ltv_discr = discretize(ltv, breaks = 4)) %>%
  filter(!categoria_produto %in% c("Kit", "Comercial", "Outros") & !uf_empreendimento == "RO") %>%
  drop_na()

cart_inad <- cart_inad %>%
  filter(ltv < 1) %>%
  setDT()

head(cart_inad)

        id_venda  data_inad data_venda uf_empreendimento categoria_produto
          <char>     <Date>     <Date>            <char>            <char>
1: 103-273VA-100 2017-09-02 2017-07-01                AM MCMV Faixas 2 e 3
2: 103-273VA-107 2017-10-29 2017-07-23                AM MCMV Faixas 2 e 3
3: 103-273VA-108 2017-10-01 2017-07-18                AM MCMV Faixas 2 e 3
4: 103-273VA-111 2018-02-16 2017-07-28                AM MCMV Faixas 2 e 3
5: 103-273VA-112 2017-12-01 2017-07-31                AM MCMV Faixas 2 e 3
6: 103-273VA-113 2017-11-01 2017-07-26                AM MCMV Faixas 2 e 3
          ltv          tipo_venda gender   safra tempo_ate_inad censura
        <num>              <char> <char>  <char>          <num>   <num>
1: 0.08070344 Financiamento Banco   Male 2017 Q3             63       1
2: 0.08603352 Financiamento Banco Female 2017 Q3             98       1
3: 0.06694486 Financiamento Banco Female 2017 Q3             75       1
4: 0.11929598 Financiamento Banco   Male 2017 Q3            203       1
5: 0.16972369 Financiamento Banco   Male 2017 Q3            123       1
6: 0.14987089 Financiamento Banco Female 2017 Q3             98       1
           ltv_discr
              <fctr>
1:    [0.0732,0.121)
2:    [0.0732,0.121)
3: [0.000197,0.0732)
4:    [0.0732,0.121)
5:      [0.167,63.9]
6:     [0.121,0.167)

Por fim, para utilizar as funções do pacote ‘Survival’, precisamos criar um objeto da classe survival.

surv_object <-
  survival::Surv(time = cart_inad$tempo_ate_inad, event = cart_inad$censura)

3.2 Curva de Sobrevivência

Para calcular a Curva de Sobrevivência vamos utilizar o estimador de Kaplan-Meier. Inicialmente, vamos considerar a base completa.

km_geral <-
  survival::survfit(surv_object ~ 1, data = cart_inad, type = "kaplan-meier")

A curva de sobrevivência pode ser traçada utilizando funções próprias do pacote ‘Survival’ conforme código abaixo.

km_geral_plot <- ggsurvplot(
  fit = km_geral,
  title = "Tempo até inadimplência",
  xlab = "Tempo em dias",
  ylab = "Probabilidade de sobrevivência",
  conf.int = T,
  conf.int.alpha = c(0.1),
  legend = "none",
  censor = FALSE
)
km_geral_plot

Ignoring unknown labels:
• fill : "Strata"
Ignoring unknown labels:
• fill : "Strata"

A figura acima mostra a probabilidade de um indivíduo sobreviver ao evento de interesse – neste caso entrar em inadimplência – até o tempo \(t\). Exemplo: a probabilidade de um cliente chegar ao dia 500 do contrato sem ter entrado em inadimplência é de, aproximadamente, 25%. O tempo mediano até o evento fica armazenado no objeto survfit criado. No caso desse exemplo ele é de 235 dias.

km_geral

Call: survfit(formula = surv_object ~ 1, data = cart_inad, type = "kaplan-meier")

         n events median 0.95LCL 0.95UCL
[1,] 40618  30282    235     230     237

Também é possível visualizar a curva de Kaplan-Meier considerando subconjuntos dos dados. Neste caso, precisamos que as covariáveis sejam categóricas ou contínuas discretizadas (como é o caso do LTV neste exemplo). Antes de traçar as curvas, é necessário que sejam criados os objetos ‘survival’ correspondentes. Em seguida podemos obter as curvas utilzando a função ‘ggsurvplot’.

km_uf <- survival::survfit(surv_object ~ uf_empreendimento,
  data = cart_inad,
  type = "kaplan-meier"
)

km_categoria <- survival::survfit(
  surv_object ~
    categoria_produto,
  data = cart_inad,
  type = "kaplan-meier"
)

km_tipo_venda <- survival::survfit(
  surv_object ~
    tipo_venda,
  data = cart_inad,
  type = "kaplan-meier"
)

km_gender <- survival::survfit(
  surv_object ~
    gender,
  data = cart_inad,
  type = "kaplan-meier"
)


km_safra <- survival::survfit(
  surv_object ~
    safra,
  data = cart_inad,
  type = "kaplan-meier"
)

km_LTV <- survival::survfit(
  surv_object ~
    ltv_discr,
  data = cart_inad,
  type = "kaplan-meier"
)

Kaplan-Meier por UF

km_uf_plot <- ggsurvplot(
  fit = km_uf,
  title = "Tempo até inadimplência",
  xlab = "Tempo em dias",
  ylab = "Probabilidade de sobrevivência",
  conf.int = F,
  conf.int.alpha = c(0.1),
  censor = FALSE
)
km_uf_plot

Kaplan Meier por categoria do produto

km_categoria_plot <- ggsurvplot(
  fit = km_categoria,
  title = "Tempo até inadimplência",
  xlab = "Tempo em dias",
  ylab = "Probabilidade de sobrevivência",
  conf.int = F,
  conf.int.alpha = c(0.1),
  censor = FALSE
)
km_categoria_plot

Kaplan-Meier por tipo de venda

km_tipovenda_plot <- ggsurvplot(
  fit = km_tipo_venda,
  title = "Tempo até inadimplência",
  xlab = "Tempo em dias",
  ylab = "Probabilidade de sobrevivência",
  conf.int = F,
  conf.int.alpha = c(0.1),
  censor = FALSE
)

km_tipovenda_plot

Kaplan-Meier por sexo

km_gender <-
  survival::survfit(
    surv_object ~ gender, data = cart_inad, type = "kaplan-meier")

km_gender_plot <- ggsurvplot(
  fit = km_gender,
  title = "Tempo até inadimplência",
  xlab = "Tempo em dias",
  ylab = "Probabilidade de sobrevivência",
  conf.int = F,
  conf.int.alpha = c(0.1),
  censor = FALSE
)

km_gender_plot

Kaplan-Meier por safra

km_safra_plot <- ggsurvplot(
  fit = km_safra,
  title = "Tempo até inadimplência",
  xlab = "Tempo em dias",
  ylab = "Probabilidade de sobrevivência",
  conf.int = F,
  conf.int.alpha = c(0.1),
  censor = FALSE
)

km_safra_plot

Kaplan-Meier por LTV

km_LTV_plot <- ggsurvplot(
  fit = km_LTV,
  title = "Tempo até inadimplência",
  xlab = "Tempo em dias",
  ylab = "Probabilidade de sobrevivência",
  conf.int = F,
  conf.int.alpha = c(0.1),
  censor = FALSE
)

km_LTV_plot

O exame das curvas de Kaplan-Meier em subamostras permite a identificação de categorias nas quais as taxas de sobrevivência são maiores ou menores. No caso da análise por UF, por exemplo, é possível identificar que contratos originados no Pará tendem a entrar em inadimplência mais rapidamente enquanto que contratos de Minas Gerais e São Paulo demoram mais, na média.

3.3 Curva de Risco

A curva de risco pode ser compreendida como uma taxa instantânea de ocorrência do evento de interesse. Ela permite ao usuário identificar períodos nos quais os indivíduos estão mais, ou menos, susceptíveis a falhas. Nesta nota optamos por estimá-la através da implementação disponível no pacote ‘casebase’. No exemplo abaixo, queremos traçar as curvas de risco para cada uma das safras disponíveis em nosso banco de dados. O código pode demorar alguns minutos para rodar.

hazard_safra <-
  fitSmoothHazard(censura ~ ns(log(tempo_ate_inad), df = 3) + safra,
  data = cart_inad,
  time = "tempo_ate_inad"
)

plot(hazard_safra,
  hazard.params = list(
    xvar = "tempo_ate_inad",
    by = "safra",
    ylab = "Hazard"
  )
)

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

Error in `plot()`:
! `...` must be empty.
✖ Problematic argument:
• ylab = "Hazard"

Note que, para traçar as curvas por UF ou categoria, bastaria mudar a covariável do primeiro argumento da função ‘fitSmoothHazard’. A forma das curvas depende fortemente da parametrização utilizada na função. Por esse motivo, a interpretação dos gráficos deve ser feita com cautela e, preferencialmente, de maneira qualitativa. No caso da figura acima, podemos identificar claramente que as safras mais antigas apresentam maior risco de inadimplência enquanto que as safras de 2020 e 2021 apresentam risco menor.

3.4 Modelo AFT

Podemos utilizar as covariáveis empregadas no cálculo das curvas de sobrevivência via Kaplan-Meier para estimar um modelo AFT. A função ‘survreg’ do pacote ‘survival’ requer três inputs básicos: a forma funcional do modelo, a base de dados e a distribuição assumida para os tempos de sobrevivência. Note que, na maior parte dos casos, é impossível saber, a priori, o processo gerador da variável de interesse \(T\). A função ‘survreg’ aceita as seguintes distribuições: Exponencia, Weibull, Normal, LogNormal, Logística e Log-Logística. Usualmente, recomenda-se estimar o modelo utilizando mais de uma distribuição e selecionar a melhor especificação através de critérios de informação como o AIC (Akaike Information Criterion) e o BIC (Bayesian Information Criterion). No exemplo abaixo, utilizamos a distribuição LogNormal já que é aquela que resulta no modelo com o menor AIC.

AFT <-
  survreg(
    surv_object ~ uf_empreendimento + categoria_produto + ltv + tipo_venda +
      gender + safra,
    data = cart_inad,
    dist = "lognormal"
  )

summary(AFT)

Call:
survreg(formula = surv_object ~ uf_empreendimento + categoria_produto + 
    ltv + tipo_venda + gender + safra, data = cart_inad, dist = "lognormal")
                                      Value Std. Error      z       p
(Intercept)                         5.06281    0.05013 100.99 < 2e-16
uf_empreendimentoCE                -0.01003    0.02208  -0.45 0.64954
uf_empreendimentoDF                 0.06672    0.02014   3.31 0.00092
uf_empreendimentoGO                -0.18210    0.02741  -6.64 3.1e-11
uf_empreendimentoMG                 0.07791    0.01866   4.17 3.0e-05
uf_empreendimentoPA                -0.13590    0.22966  -0.59 0.55401
uf_empreendimentoPE                -0.01877    0.04749  -0.40 0.69271
uf_empreendimentoRJ                -0.10736    0.01720  -6.24 4.3e-10
uf_empreendimentoSP                -0.06525    0.01745  -3.74 0.00018
categoria_produtoMCMV Faixas 2 e 3  0.14557    0.01200  12.14 < 2e-16
categoria_produtoMédio              0.16641    0.03013   5.52 3.3e-08
categoria_produtoMédio-Alto        -0.03802    0.07343  -0.52 0.60465
ltv                                -0.58775    0.08114  -7.24 4.4e-13
tipo_vendaRecurso Próprio           1.16618    0.01156 100.92 < 2e-16
genderMale                         -0.05128    0.01022  -5.02 5.2e-07
safra2017 Q2                        0.02012    0.06001   0.34 0.73746
safra2017 Q3                       -0.21747    0.05901  -3.69 0.00023
safra2017 Q4                       -0.20510    0.05741  -3.57 0.00035
safra2018 Q1                        0.07370    0.05443   1.35 0.17568
safra2018 Q2                       -0.11434    0.05259  -2.17 0.02969
safra2018 Q3                       -0.04301    0.05302  -0.81 0.41723
safra2018 Q4                        0.04587    0.05157   0.89 0.37376
safra2019 Q1                        0.04661    0.05261   0.89 0.37561
safra2019 Q2                        0.11479    0.05121   2.24 0.02498
safra2019 Q3                        0.14539    0.05122   2.84 0.00453
safra2019 Q4                        0.16717    0.05070   3.30 0.00098
safra2020 Q1                       -0.03028    0.05168  -0.59 0.55790
safra2020 Q2                        0.03115    0.05086   0.61 0.54024
safra2020 Q3                       -0.15488    0.05088  -3.04 0.00234
safra2020 Q4                        0.09868    0.05081   1.94 0.05209
safra2021 Q1                        0.08001    0.05059   1.58 0.11377
safra2021 Q2                       -0.02329    0.05048  -0.46 0.64453
safra2021 Q3                       -0.24804    0.05121  -4.84 1.3e-06
safra2021 Q4                        0.06756    0.05399   1.25 0.21082
safra2022 Q1                        0.29657    0.08711   3.40 0.00066
Log(scale)                         -0.03370    0.00421  -8.00 1.2e-15

Scale= 0.967 

Log Normal distribution
Loglik(model)= -203304   Loglik(intercept only)= -209109.1
    Chisq= 11610.22 on 34 degrees of freedom, p= 0 
Number of Newton-Raphson Iterations: 4 
n= 40618 

extractAIC(AFT)[2]

[1] 406679.9

Como todas as variáveis explicativas são categóricas, devemos interpretar os resultados em relação à um elemento de cada grupo que é tomado como base de comparação. No caso da variável UF, o elemento de comparação é o estado do Amazonas. Sendo assim, contratos originados em estados com coeficientes negativos na regressão apresentam, em média, tempo até inadimplência menor que no Amazonas. Esse efeito é mais evidente conforme o coeficiente de regressão fica mais negativo. Dessa forma, tudo o mais constante, o estado de Minas Gerais é aquele que apresentam maior tempo até inadimplência. Análise semelhante pode ser feita em relação às outras variáveis do modelo.

3.5 Modelo de Cox

A implementação da regressão de Cox pode ser feita através da função ‘coxph’ do pacote survival. Os argumentos da função são a especificação do modelo e o data frame com os dados, conforme código abaixo:

cox <-
  coxph(
    surv_object ~ uf_empreendimento + categoria_produto + ltv_discr + 
      tipo_venda + gender + safra,
    data = cart_inad
  )

summary(cox)

Call:
coxph(formula = surv_object ~ uf_empreendimento + categoria_produto + 
    ltv_discr + tipo_venda + gender + safra, data = cart_inad)

  n= 40618, number of events= 30282 

                                        coef exp(coef)  se(coef)       z
uf_empreendimentoCE                -0.065361  0.936729  0.025348  -2.579
uf_empreendimentoDF                -0.130350  0.877788  0.022671  -5.750
uf_empreendimentoGO                 0.146272  1.157511  0.030484   4.798
uf_empreendimentoMG                -0.144835  0.865165  0.021259  -6.813
uf_empreendimentoPA                 0.186609  1.205156  0.244203   0.764
uf_empreendimentoPE                 0.004239  1.004248  0.056440   0.075
uf_empreendimentoRJ                 0.054658  1.056179  0.019130   2.857
uf_empreendimentoSP                -0.077142  0.925758  0.019787  -3.899
categoria_produtoMCMV Faixas 2 e 3 -0.145661  0.864451  0.013640 -10.679
categoria_produtoMédio             -0.192968  0.824509  0.035414  -5.449
categoria_produtoMédio-Alto         0.033072  1.033625  0.083929   0.394
ltv_discr[0.0732,0.121)             0.007857  1.007888  0.016063   0.489
ltv_discr[0.121,0.167)              0.033576  1.034146  0.016754   2.004
ltv_discr[0.167,63.9]               0.136328  1.146058  0.018570   7.341
tipo_vendaRecurso Próprio          -1.003948  0.366430  0.013651 -73.546
genderMale                          0.050450  1.051745  0.011641   4.334
safra2017 Q2                        0.008999  1.009039  0.063769   0.141
safra2017 Q3                        0.244253  1.276667  0.062579   3.903
safra2017 Q4                        0.141467  1.151962  0.060991   2.319
safra2018 Q1                       -0.135200  0.873542  0.058059  -2.329
safra2018 Q2                        0.144284  1.155212  0.055907   2.581
safra2018 Q3                        0.107450  1.113435  0.056453   1.903
safra2018 Q4                        0.031921  1.032436  0.055113   0.579
safra2019 Q1                        0.015698  1.015822  0.056318   0.279
safra2019 Q2                       -0.017132  0.983014  0.054866  -0.312
safra2019 Q3                       -0.033537  0.967019  0.055096  -0.609
safra2019 Q4                       -0.046442  0.954620  0.054503  -0.852
safra2020 Q1                        0.125779  1.134032  0.055782   2.255
safra2020 Q2                        0.019545  1.019737  0.054864   0.356
safra2020 Q3                        0.053085  1.054520  0.055110   0.963
safra2020 Q4                       -0.031744  0.968754  0.055270  -0.574
safra2021 Q1                        0.006200  1.006219  0.055279   0.112
safra2021 Q2                        0.132843  1.142071  0.055260   2.404
safra2021 Q3                        0.339579  1.404356  0.056711   5.988
safra2021 Q4                       -0.043723  0.957219  0.064746  -0.675
safra2022 Q1                       -0.480131  0.618703  0.163055  -2.945
                                   Pr(>|z|)    
uf_empreendimentoCE                 0.00992 ** 
uf_empreendimentoDF                8.94e-09 ***
uf_empreendimentoGO                1.60e-06 ***
uf_empreendimentoMG                9.56e-12 ***
uf_empreendimentoPA                 0.44477    
uf_empreendimentoPE                 0.94013    
uf_empreendimentoRJ                 0.00427 ** 
uf_empreendimentoSP                9.67e-05 ***
categoria_produtoMCMV Faixas 2 e 3  < 2e-16 ***
categoria_produtoMédio             5.07e-08 ***
categoria_produtoMédio-Alto         0.69354    
ltv_discr[0.0732,0.121)             0.62474    
ltv_discr[0.121,0.167)              0.04506 *  
ltv_discr[0.167,63.9]              2.12e-13 ***
tipo_vendaRecurso Próprio           < 2e-16 ***
genderMale                         1.46e-05 ***
safra2017 Q2                        0.88778    
safra2017 Q3                       9.50e-05 ***
safra2017 Q4                        0.02037 *  
safra2018 Q1                        0.01988 *  
safra2018 Q2                        0.00986 ** 
safra2018 Q3                        0.05699 .  
safra2018 Q4                        0.56246    
safra2019 Q1                        0.78044    
safra2019 Q2                        0.75485    
safra2019 Q3                        0.54272    
safra2019 Q4                        0.39416    
safra2020 Q1                        0.02414 *  
safra2020 Q2                        0.72166    
safra2020 Q3                        0.33541    
safra2020 Q4                        0.56573    
safra2021 Q1                        0.91070    
safra2021 Q2                        0.01622 *  
safra2021 Q3                       2.13e-09 ***
safra2021 Q4                        0.49948    
safra2022 Q1                        0.00323 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

                                   exp(coef) exp(-coef) lower .95 upper .95
uf_empreendimentoCE                   0.9367     1.0675    0.8913    0.9844
uf_empreendimentoDF                   0.8778     1.1392    0.8396    0.9177
uf_empreendimentoGO                   1.1575     0.8639    1.0904    1.2288
uf_empreendimentoMG                   0.8652     1.1558    0.8299    0.9020
uf_empreendimentoPA                   1.2052     0.8298    0.7468    1.9450
uf_empreendimentoPE                   1.0042     0.9958    0.8991    1.1217
uf_empreendimentoRJ                   1.0562     0.9468    1.0173    1.0965
uf_empreendimentoSP                   0.9258     1.0802    0.8905    0.9624
categoria_produtoMCMV Faixas 2 e 3    0.8645     1.1568    0.8416    0.8879
categoria_produtoMédio                0.8245     1.2128    0.7692    0.8838
categoria_produtoMédio-Alto           1.0336     0.9675    0.8768    1.2184
ltv_discr[0.0732,0.121)               1.0079     0.9922    0.9767    1.0401
ltv_discr[0.121,0.167)                1.0341     0.9670    1.0007    1.0687
ltv_discr[0.167,63.9]                 1.1461     0.8726    1.1051    1.1885
tipo_vendaRecurso Próprio             0.3664     2.7290    0.3568    0.3764
genderMale                            1.0517     0.9508    1.0280    1.0760
safra2017 Q2                          1.0090     0.9910    0.8905    1.1434
safra2017 Q3                          1.2767     0.7833    1.1293    1.4433
safra2017 Q4                          1.1520     0.8681    1.0222    1.2982
safra2018 Q1                          0.8735     1.1448    0.7796    0.9788
safra2018 Q2                          1.1552     0.8656    1.0353    1.2890
safra2018 Q3                          1.1134     0.8981    0.9968    1.2437
safra2018 Q4                          1.0324     0.9686    0.9267    1.1502
safra2019 Q1                          1.0158     0.9844    0.9097    1.1344
safra2019 Q2                          0.9830     1.0173    0.8828    1.0946
safra2019 Q3                          0.9670     1.0341    0.8680    1.0773
safra2019 Q4                          0.9546     1.0475    0.8579    1.0622
safra2020 Q1                          1.1340     0.8818    1.0166    1.2650
safra2020 Q2                          1.0197     0.9806    0.9158    1.1355
safra2020 Q3                          1.0545     0.9483    0.9466    1.1748
safra2020 Q4                          0.9688     1.0323    0.8693    1.0796
safra2021 Q1                          1.0062     0.9938    0.9029    1.1214
safra2021 Q2                          1.1421     0.8756    1.0248    1.2727
safra2021 Q3                          1.4044     0.7121    1.2566    1.5695
safra2021 Q4                          0.9572     1.0447    0.8431    1.0867
safra2022 Q1                          0.6187     1.6163    0.4495    0.8517

Concordance= 0.686  (se = 0.002 )
Likelihood ratio test= 7819  on 36 df,   p=<2e-16
Wald test            = 7233  on 36 df,   p=<2e-16
Score (logrank) test = 7769  on 36 df,   p=<2e-16

Os coeficientes estimados exponenciados fornecem a razão de risco (hazard ratio). Razões acima de um indicam que a covariável está associada à um aumento de risco, já razões abaixo estão relacionadas a atributos protetivos, isto é, que diminuem o risco de de mortalidade. A probabilidade de ocorrência de um evento num dado período de tempo é igual ao risco multiplicado pela duração do intervalo.

No nosso exemplo, a razão de risco de Minas Gerais é de 0.87. Isso significa que, tudo mais constante, um contrato originado em MG possui um risco de inadimplência aproximadamente 13% menor em relação à um contrato proveniente do Amazonas (estado tomado como referência na análise).

4 Recuperação de Crédito

4.1 Objetivo

O objetivo deste exercício é estimar a probabilidade de que um contrato permaneça em inadimplência por \(x\) dias adicionais dado que está em default há \(T\) dias

4.2 Preparação dos dados

A preparação da base de dados envolve os seguintes passos:

1. Filtrar apenas contratos que possuam parcelas cuja data de recebimento é maior do que a data de vencimento ou que possuam parcelas em aberto na data em que a base de dados foi gerada.
2. Calcular o tempo até a recuperação do contrato – ou seja: o tempo, em dias, decorrido do momento em que o contrato entrou em inadimplência até o momento em que todas as parcelas em atraso estão quitadas. Note que, por esse critério, o tempo de inadimplência do contrato pode ser maior que o tempo máximo em default de suas parcelas em aberto.
   • Exemplo: considere um cliente que atrasou duas parcelas consecutivas – a primeira com vencimento em 01/01 e a segunda em 01/02. Suponha que ele conseguiu quitar a parcela referente a janeiro em 15/02, porém não foi capaz de pagar a parcela atrasada de fevereiro. Nesse caso, o tempo em inadimplência do contrato é de 45 dias e não de 14 dias.
3. Caso o contrato continue em default na data de referência da análise, ele entrará no modelo como uma observação censurada.

A função abaixo calcula os dias em atraso de cada contrato de acordo com os passo (2) e (3). Ela será aplicada a todos os contratos filtrados de acordo com o passo (1).

retorna_dias_atraso <- function(df) {
  if (!"data.frame" %in% class(df)) {
    stop("parametro df nao é data frame")
  }
  dias_atraso <- df %>%
    mutate(start_inad = ifelse(
      index == 1, data_vencimento,
      ifelse(data_vencimento < lag(data_recebimento), NA, data_vencimento)
    )) %>%
    fill(start_inad, .direction = "down") %>%
    mutate(start_inad = as_date(start_inad, origin = lubridate::origin)) %>%
    mutate(end_inad = ifelse(is.na(data_recebimento) & index == nrow(.),
      data_ref, ifelse(data_recebimento <
        lead(data_vencimento) &
        !is.na(data_recebimento), data_recebimento, NA)
    )) %>%
    mutate(end_inad = (as_date(end_inad, origin = lubridate::origin))) %>%
    select("start_inad", "end_inad") %>%
    drop_na(.) %>%
    mutate(tempo_ate_evento = as.numeric(end_inad - start_inad)) %>%
    mutate(evento = ifelse(end_inad == data_ref, 0, 1))
  return(dias_atraso)
}

Finalmente, o código abaixo gera a base de dados final que será utilizada no cálculo das probabilidades condicionais.

data_ref <- max(dados$data_recebimento, na.rm = TRUE)

cart_inad_rec <- dados %>%
  filter(data_recebimento > data_vencimento | (data_vencimento < data_ref &
    is.na(data_recebimento))) %>%
  select(
    "id_venda", "cpf_cnpj", "uf_empreendimento", "categoria_produto",
    "valor_venda", "data_vencimento", "data_recebimento",
    "valor_nominal_das_parcelas"
  ) %>%
  mutate(dias_atraso = ifelse(!is.na(data_recebimento),
    data_recebimento - data_vencimento,
    data_ref - data_vencimento
  )) %>%
  select("id_venda", "data_vencimento", "data_recebimento", "dias_atraso") %>%
  group_by(id_venda) %>%
  arrange(data_vencimento, .by_group = T) %>%
  mutate(index = row_number(id_venda)) %>%
  group_modify(~ retorna_dias_atraso(.x))

head(cart_inad_rec)

    id_venda start_inad   end_inad tempo_ate_evento evento
      <char>     <Date>     <Date>            <num>  <num>
1: 1-CVBI-11 2019-12-20 2019-12-26                6      1
2: 1-CVBI-11 2020-01-20 2020-01-23                3      1
3: 1-CVBI-11 2020-02-20 2020-02-26                6      1
4: 1-CVBI-11 2020-05-20 2020-05-22                2      1
5: 1-CVBI-11 2020-09-20 2020-09-21                1      1
6: 1-CVBI-11 2020-11-20 2020-12-03               13      1

Para estimar as funções de sobrevivência necessárias para computar as probabilidades condicionais vamos utilizar o estimador de Kaplan Meier. Sendo assim, como no exemplo anterior envolvendo o cálculo do tempo até a primeira inadimplência, devemos transformar a base de dados em um objeto do pacote ‘survival’.

surv_object <- survival::Surv(
  time = cart_inad_rec$tempo_ate_evento,
  event = cart_inad_rec$evento
)

4.3 Resultados

Agora podemos traçar a curva de sobrevivência associada aos dados de tempo até recuperação. Ela mostra a probabilidade de um contrato com pelo menos um dia de atraso permanecer inadimplente por \(n\) dias. Note que a inflexão da curva ocorre perto do dia 180, sinalizando que, depois de 6 meses de atraso, o contrato se torna praticamente irrecuperável.

f1 <- survival::survfit(surv_object ~ 1,
  data = cart_inad_rec,
  type = "kaplan-meier"
)

kmplot_tempo_recuperacao <- ggsurvplot(
  fit = f1,
  title = "Probabilidade do contrato permanecer inadimplente \n dado que está com pelo menos um dia de atraso",
  xlab = "Tempo em dias",
  ylab = "Probabilidade de sobrevivência",
  conf.int = T,
  conf.int.alpha = c(0.1),
  legend = "none",
  censor = FALSE
)

kmplot_tempo_recuperacao

Ignoring unknown labels:
• fill : "Strata"
Ignoring unknown labels:
• fill : "Strata"

Observe que a inflexão da curva de Kaplan Meier ocorre perto do dia 180, sinalizando que, depois de 6 meses de atraso, o contrato se torna praticamente irrecuperável. Sendo assim, vamos usar a equação de probabilidade de sobrevida para calcular a probabilidade de que um contrato inadimplente a \(T\) dias alcance \(x=180\) dias de atraso. O código abaixo faz isso.

prob_cond <- NULL
for (i in 1:179) {
  prob_cond[i] <- f1$surv[180] / f1$surv[i]
}


prob_cond <- as.data.frame(prob_cond) %>%
  mutate(dias = row_number()) %>%
  setDT()

hence_line_chart(
  p_x = dias,
  p_y = prob_cond,
  p_dados = prob_cond,
  p_xlab = "Dias em atraso",
  p_ylab = "Probabilidade",
  p_title = "Probabilidade de que o contrato chegue \n a pelo menos 180 dias de atraso"
) %>%
  ggplotly()