Introdução à Análise de sobrevivência
Em diversas situações práticas em análise de dados lidamos com o problema de analisar o tempo até a ocorrência de um evento de interesse. Perguntas como “Em que momento o cliente deixará de ser cliente?” ou “Em que momento o colaborador deixará a empresa?” são relevantes para uma companhia e podem ser respondidas com técnicas de análise de sobrevivência. Nesse artigo iremos aplicar a metodologia na prática em um exemplo com dados abertos.
Modelos de sobrevivência são bastante utilizados em outras áreas, especialmente em medicina, como por exemplo para avaliar o efeito de tratamentos no tempo de vida de pacientes; e em engenharia, com aplicações na análise do tempo até a falha de algum dispositivo. Entretanto, o uso de tais modelos fora da academia ainda não é tão difundido como outras técnicas de regressão e classificação.
Censura
Ao iniciarmos qualquer estudo, precisamos definir uma data de início e fim para a coleta dos dados. Como estamos analisando o tempo até a ocorrência de um evento, em muitos casos pode acontecer do evento ainda não ter ocorrido. Por exemplo, ao avaliar o tempo até um colaborador deixar o seu cargo, teremos muitos casos de colaboradores que ainda não fizeram isso, logo não temos esse tempo disponível. Assim, dizemos que o dado é “censurado” ou “incompleto”. A única coisa que podemos garantir é que o evento irá ocorrer após esse tempo coletado. Um outro caso que pode acarretar em censura é a perda de acompanhamento ou desistência do estudo.
Há diferentes tipos de censura, sendo a censura à direita a mais comum, e que iremos utilizar nos exemplos, que é o caso em que o evento irá ocorrer após o fim do estudo. Também temos a censura à esquerda, aonde o evento ocorreu antes do estudo começar (mais incomum) e a censura intervalar, caso em que não sabemos exatamente quanto o evento ocorreu, mas sabemos que com certeza ele ocorreu dentro de um intervalo de tempo conhecido.
Como as técnicas de análise de sobrevivência lidam com censura, não há perda de informação ou viés ao analisarmos os tempos, coisa que aconteceria se descartássemos os dados censurados.
Portanto, para modelar o problema através de análise de sobrevivência devemos ter a informação do tempo e um indicador de que o tempo coletado é completo ou censurado, aonde geralmente usa-se uma variável binária (1 para o caso em que o tempo não é censurado, 0 para a censura).
Curvas de Sobrevivência
Uma análise comum a ser feita para esse tipo de problema é a construção de curvas de sobrevivência. A curva de sobrevivência descreve a probabilidade do evento ocorrer após um determinado período, e nos dá um bom resumo do comportamento do tempo até a ocorrência do evento. Podemos, por exemplo, comparar diferentes curvas para diferentes grupos, podendo testar o efeito de uma variável no tempo de sobrevivência.
A forma mais simples de estimar uma curva de sobrevivência é através do estimador de Kaplan-Meier (KM), que gera uma curva de sobrevivência não paramétrica (não é necessário fazer nenhuma suposição sobre a distribuição de probabilidade utilizada para modelar os dados).
É possível gerar curvas de sobrevivência estimadas via KM com a biblioteca do python chamada lifelines. No R temos a opção do pacote survival, que já vem instalado de forma padrão.
Abaixo temos um exemplo de código com um banco de dados com o tempo até o desligamento de um colaborador e com outras informações como idade, gênero, setor da empresa, meio de transporte utilizado para ir ao trabalho, etc. A base de dados está disponível no Kaggle (https://www.kaggle.com/datasets/davinwijaya/employee-turnover?resource=download) e contém 1129 observações.
Devemos passar o vetor com os tempos T e o indicador de censura E:
from lifelines.datasets import load_rossi
from lifelines import KaplanMeierFitterfrom lifelines import KaplanMeierFitterT = df['stag']
E = df['event']kmf = KaplanMeierFitter()
kmf.fit(T, event_observed=E)
kmf.plot_survival_function()
Abaixo a curva de sobrevivência com intervalo de confiança de 95%:
Uma quantidade de interesse é o tempo mediano (e o seu intervalo de confiança), que é o momento em 50% da amostra passou pelo evento, sendo o ponto da curva em que a probabilidade chega a 0.5. O tempo mediano pode ser recuperado com o seguinte comando:
from lifelines.utils import median_survival_times
median_ = kmf.median_survival_time_
ci = kmf.confidence_interval_
median_confidence_interval_ = median_survival_times(ci)No exemplo, o valor da mediana é 50.72 e o intervalo de confiança de 95% é (45.53, 54.01).
Comparando Curvas de Sobrevivência
Avaliar o impacto de alguma variável nos tempos de sobrevivência é outra análise recorrente. Podemos fazer isso comparando as curvas geradas pelo método KM para os grupos. Também é possível realizar um teste estatístico em que a hipótese nula é de que as curvas são iguais (o fator não tem influência no tempo de sobrevivência), assim, para um baixo p-valor, rejeitamos a hipótese a nula de igualdade entre as curvas. Abaixo o código para gerar duas curvas comparando os tempos entre os colaboradores que vão a pé para o trabalho e os que vão de carro ou ônibus:
ax = plt.subplot(111)kmf.fit(T[hg], event_observed=E[hg], label="Feminino")
kmf.plot_survival_function(ax=ax)kmf.fit(T[~hg], event_observed=E[~hg], label="Masculino")
kmf.plot_survival_function(ax=ax)
Visualmente parece existir uma diferença entre as curvas. Nota-se que os colaboradores que não vão a pé ao trabalho tendem a se desligar mais rápido, dado que a probabilidade cai de forma mais acentuada.
Para aplicar o teste log-rank, que tem como hipótese nula a igualdade entre as curvas de sobrevivência, precisamos passar as duas colunas com tempos e as duas colunas com indicador de censura:
from lifelines.statistics import logrank_testresults = logrank_test(T[wa], T[~wa], E[wa], E[~wa])
results.print_summary()
O p-valor calculado foi de 0.02, logo temos indício de que essa variável pode ser relevante.
Modelos de Sobrevivência
Se quisermos avaliar o efeito de diversas variáveis ao mesmo tempo, devemos utilizar modelos mais complexos. Entre as opções disponíveis estão os modelos paramétricos, aonde é necessário escolher uma distribuição de probabilidade para os dados, com um número de parâmetros bem definidos. Modelos de tempo de falha acelerado (com a sigla em inglês AFT) são exemplos de modelos de sobrevivência paramétricos em que o tempo de sobrevivência passa a depender de um conjunto de features (ou covariáveis).
Com modelos paramétricos podemos fazer predições e também avaliar o efeito de cada feature. Nesses modelos temos coeficientes com a influência de cada variável e e um p-valor. O p-valor está associado ao teste em que a hipótese nula é a de que o coeficiente é zero, caso em que a feature não influencia no tempo de sobrevivência.
Nesse cenário o tempo de sobrevivência deve seguir uma distribuição de probabilidade pré-definida. Abaixo o exemplo do ajuste do modelo que assume os tempos segue uma distribuição Weibull. Precisamos definir a coluna com os tempos (duration_col) e a coluna com o indicador de censura (event_col):
from lifelines import WeibullAFTFitter#transformando variáveis do tipo 'object' em categórica
for i in range(len(df.dtypes)):
if df.dtypes[i] == 'object':
df[df.dtypes.index[i]] = pd.Categorical(df[df.dtypes.index[i]])#criando dummies para as variáveis categóricas
mdf = pd.get_dummies(df.drop(['traffic'], axis=1), drop_first = True)aft = WeibullAFTFitter()
aft.fit(mdf, duration_col='stag', event_col='event')print(aft.print_summary(3))
O modelo retorna os coeficientes, intervalo de confiança e p-valor associado a cada variável. Na imagem abaixo temos tais valores para algumas variáveis. Nota-se, por exemplo, que a idade parece ser fator relevante (baixo p-valor), enquanto o gênero não parece interferir no tempo até o desligamento. Para variável relacionada ao meio de transporte o p-valor associado foi 0.005, corroborando com o resultado obtido pelo KM.
Diversas alternativas de previsão estão disponíveis, como a mediana do tempo, a curva de sobrevivência como um todo ou algum percentil específico.
X = mdf.loc[:10]
aft.predict_survival_function(X)
aft.predict_median(X)
aft.predict_percentile(X, p=0.9)C-index
Uma forma de avaliar a qualidade da predição de modelos de sobrevivência é através do C-index, que é a proporção de vezes em que a ordenação de um par de observações é concordante com o que aconteceu na realidade. Por exemplo, se o tempo de sobrevivência de um par de observações for (10, 12) e o valor predito para esses dois indivíduos for (8, 15), temos um par concordante. Um modelo completamente aleatório apresenta um C-index médio de 0.5, logo modelos úteis devem apresentar valores maiores do que 0.5. Vale ressaltar que o C-index é bom em avaliar a ordenação das predições, mas não o erro entre o valor predito e observado. Assim, modelos com maiores C-index ordenam melhor as observações em relação ao tempo de sobrevivência.
O C-index para o exemplo em questão (isso considerando toda a base de dados, não houve a separação do treino e teste), é de 0.64. Calculando o c-index através de validação cruzada com 10 folds, temos que a média da métrica é de 0.60, valor relativamente baixo. Para melhorar tal métrica, poderíamos alterar a distribuição escolhida para o modelagem, pois talvez a Weibull não está bem ajustada aos dados.
from lifelines.utils import k_fold_cross_validationscores = k_fold_cross_validation(aft, mdf, 'stag', event_col='event', k=10, scoring_method='concordance_index')
print(scores)
Análise de Sobrevivência com xgboost
Também há a alternativa de ajustar um modelo de sobrevivência com xgboost, biblioteca de boosting muito popular. Com o xgboost podemos escalar para problemas com alto volume de dados, mas perdemos a possibilidade de testar a significância das features, bem como construir curvas de sobrevivência de forma mais fácil.
A indicação de censura no xgboost é um pouco diferente. O modelo aceita os 3 tipos de censura, sendo necessário definir um limite inferior e superior para os dados. Em caso de não censura esses limites são iguais, para o caso de censura à direita, devemos colocar o valor ‘Inf’ para o limite superior, ou seja, o evento ocorreu entre o tempo coletado e o “infinito”.
import numpy as np
import xgboost as xgb#features
X = mdf.drop(['stag', 'event'], axis=1)
dtrain = xgb.DMatrix(X)#para o caso de censura à direita o limite inferior é sempre o valor do target
y_lower_bound = mdf['stag']
#em caso de censura, limite superior é o infinito
y_upper_bound = [np.Inf if event == 0 else time for time, event in zip(mdf['stag'], mdf['event'])]dtrain.set_float_info('label_lower_bound', y_lower_bound)
dtrain.set_float_info('label_upper_bound', y_upper_bound)
É necessário escolher a distribuição de probabilidade que o modelo irá se basear, alterar a função objetivo para “survival:aft” e definir um parâmetro de escala da distribuição escolhida (aft_loss_distribution_scale), que idealmente deve ser escolhido por métodos de validação cruzada. O modelo gera apenas estimativas para o tempo mediano com método predict.
params = {'objective': 'survival:aft',
'eval_metric': 'aft-nloglik',
'aft_loss_distribution': 'normal',
'aft_loss_distribution_scale': 1,
'tree_method': 'hist', 'learning_rate': 0.05,
'max_depth': 6}
bst = xgb.train(params, dtrain, num_boost_round=100)
pred = bst.predict(dtrain)Também é possível calcular o C-index com 10 folds, mas o código é um pouco mais elaborado, já que não há um método ou atributo específico no xgboost para o cálculo dessa métrica em um contexto de análise de sobrevivência. A média do C-index foi de 0.57, valor baixo, fica a sugestão de alterar a distribuição ou o valor da escala. O banco de dados também é pequeno, ficando difícil construir um preditor bom com uma quantidade tão baixa de observações.
from lifelines.utils import concordance_index
from sklearn.model_selection import KFold#criando objeto com indices dos folds
kf = KFold(n_splits=10)
kf.get_n_splits(X)#features
X = mdf.drop(['stag', 'event'], axis=1)
#tempos
y = mdf['stag']
#indicador de evento
ind = mdf['event']
#lista para guardar o c-index de cada fold
c_index_list = []
c_index_train_list = []#iterando em cada fold
for train_index, test_index in kf.split(X):
#dmatrix com treino e test
dtrain = xgb.DMatrix(X.iloc[train_index, :])
dtest = xgb.DMatrix(X.iloc[test_index, :])
#tempos e indicadores de evento para treino e teste
y_train, ind_train = y[train_index], ind[train_index]
y_test, ind_test = y[test_index], ind[test_index]
#para o caso do treino, adicionando a informação de censura
y_lower_bound = y_train
#em caso de censura, limite superior é o infinito
y_upper_bound = [np.Inf if event == 0 else time for time, event in zip(y_train, ind_train)]dtrain.set_float_info('label_lower_bound', y_lower_bound)
dtrain.set_float_info('label_upper_bound', y_upper_bound)
#treinando o modelo
params = {'objective': 'survival:aft',
'eval_metric': 'aft-nloglik',
'aft_loss_distribution': 'normal',
'aft_loss_distribution_scale': 1,
'tree_method': 'hist', 'learning_rate': 0.05,
'max_depth': 4}
bst = xgb.train(params, dtrain, num_boost_round=100)
#predição
pred_test = bst.predict(dtest)
pred_train = bst.predict(dtrain)
#calculando c-index
c_index = concordance_index(y_test, pred_test, ind_test)
#adicionando na lista
c_index_list.append(c_index)
Conclusão e algumas referências
Por fim, há uma vasta literatura e modelos de sobrevivência disponíveis, como por exemplo o modelo de Cox, bastante utilizado na área médica, que modela o risco estatístico do evento acontecer em um instante qualquer. Muita pesquisa ainda é feita na área de sobrevivência e é uma classe de modelo muito útil, com diversas vantagens, especialmente por lidar muito bem com dados censurados.
Algumas referências:
Colosimo Enrico Antônio, & Giolo, S. R. (2006). Análise de sobrevivência aplicada. Edgard Blücher.
