ANÁLISE DE SOBREVIVÊNCIA E CONFIABILIDADE

AULA 4: ESCOLHA DO MODELO PROBABILÍSTICO

\[\\[0.05in]\]

1 PACOTES NECESSÁRIOS

Para esta quarta aula prática de Análise de Sobrevivência e Confiabilidade, precisaremos dos seguintes pacotes:

require(flexsurv)
require(reslife)
require(dplyr)

Pacote flexsurv: Responsável pelo ajuste de modelos paramétricos flexíveis aos dados de tempo até a ocorrência de um evento de interesse.
Pacote reslife: Responsável por fornecer uma estimativa para o tempo médio residual de ocorrência de um evento de interesse.
Pacote dplyr: Responsável por facilitar a manipulação de base de dados.

2 BASES DE DADOS

Esta seção é dedicada à descrição das bases de dados utilizadas nesta aula.

2.1 CÂNCER DE BEXIGA

Essa base de dados é proveniente de assessorias estatísticas realizadas no Departamento de Estatística da UFPR. Os dados contém os tempos de reincidência, em meses, de um grupo de 20 pacientes com câncer de bexiga que foram submetidos a um procedimento cirúrgico realizado por laser. Os tempos obtidos foram:

3, 5, 6, 7, 8, 9, 10, 10+, 12, 15, 15+, 18, 19, 20, 22, 25, 28, 30, 40, 45+,

em que o simbolo + indica censura.

2.2 TRATAMENTO QUIMIOTERÁPICO

No estudo analisado neste exemplo, são apresentados os tempos, em dias até a ocorrência dos primeiros sinais de alterações indesejadas no estado geral de saúde de 45 pacientes de ambos os sexos que receberam tratamento quimioterápico após terem sido submetidos à cirurgia de intestino. Houve um acompanhamento total de 250 dias desde a entrada do primeiro paciente até o término do estudo.

7, 8, 10, 12, 13, 14+, 19, 23, 25+, 26, 27, 31, 31+, 49, 59+, 64+, 87, 89, 107, 117, 119, 230+, 233+, 130, 148, 153, 156, 159, 191, 222, 200+, 203+, 210+, 220+, 220+, 228+, 235+, 240+ , 240+, 240+, 241+, 245+, 247+, 248+, 250+,

em que o símbolo + indica a censura.

2.3 ISOLANTE ELÉTRICO

0 fabricante de um tipo de isolante elétrico quer conhecer o comportamento de seu produto funcionando a uma temperatura de 200Cº. Um teste de vida foi realizado nestas condições uasando-se 60 isoladores elétricos. O teste terminou quando 45 deles havia falhado (censura do tipo II). As 15 unidades que não haviam falhado ao final do teste foram, desta forma, censuradas no tempo t= 2729 horas:

151, 164, 336, 365, 403, 454, 455, 473, 538, 577, 592, 628, 632, 647, 675, 727, 785, 801, 811, 816, 867, 893, 930, 937, 976, 1008, 1040, 1051, 1060, 1183, 1329, 1334, 1379, 1380, 1633, 1769, 1827, 1831, 1849, 2016, 2282, 2415, 2430, 2686, 2729, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+, 2729+,

em que o símbolo + indica a censura.

2.4 AIDS

O estudo é constituído de pacientes portadores de HIV atendidos entre 1986 e 2000 no Instituto de Pesquisa Clínica Evandro Chagas (Ipec/Fiocruz). Dessa coorte, obteve-se uma amostra de 193 indivíduos que foram diagnosticados como portadores de AIDS (critério CDC 1993) durante o período de acompanhamento. O objetivo é mensurar o tempo, em dias, do diagnóstico até a morte dos pacientes. Essa base de dados está armazenada no objeto DadosAIDS.csv

3 INTRODUÇÃO

A escolha do modelo a ser utilizado é um tópico extremamente importante na análise paramétrica de dados de tempo de ocorrência de um evento de interesse. Toda a inferência a partir do método de máxima verossimilhança para as quantidades de interesse só poderá ser devidamente interpretada após ter sido definido um modelo probabilístico adequado para os dados.

Por exemplo, somente após ter definido que o modelo log-normal se ajusta bem aos dados que podemos fazer as devidas interpretações de quantidades de interesse estimadas pelo método de máxima verossimilhança. Entretanto, se o modelo log-normal for usado inadequadamente para um certo conjunto de dados, toda a análise estatística fica comprometida e consequentemente as respostas às perguntas de interesse ficam distorcidas.

Mas por que usar o modelo log-normal e não o de Weibull nesse exemplo hipotético? Algumas vezes existiam evidências provenientes de testes realizados no passado de que um certo modelo se ajusta bem a um certo tipo de dado. No entanto, em muitas situações, este tipo de informação não se encontra disponível. A solução para estas situações é basicamente empírica.

Sabe-se que as distribuições apresentadas aqui são típicas para dados de tempos de ocorrência de um evento de interesse. A proposta empírica consiste em ajustar os modelos probabilísticos apresentados (exponencial, de Weibull etc.) e, com base na comparação entre medidas de qualidade de ajuste, decidir qual deles “melhor” explica os dados amostrais. Nesta aula, veremos as seguintes técnicas para a escolha do modelo probabilístico:

Métodos gráficos
AIC e BIC

Antes de seguir com a escolha do modelo probabilístico, vamos ler nossa base de dados do câncer de bexiga:

 tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
 cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)
 dados1<-data.frame(tempos,cens)

Na aula passada, aprendemos a ajustar diversos modelos probabilísticos para estes dados. Vamos ajustar todos eles:

# Exponencial
ajuste_exp<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "exp", data=dados1)
# Weibull
ajuste_wei<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "weibull", data=dados1)
# Log-Normal
ajuste_lognorm<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "lnorm", data=dados1)
# Log-Logístico
ajuste_loglogi<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "llogi", data=dados1)
# Gama
ajuste_gama<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "gamma", data=dados1)
# Gama Generalizada
ajuste_gengama<-flexsurvreg(Surv(tempos, cens) ~ 1, dist = "gengamma.orig", data=dados1)

3.1 MÉTODOS GRÁFICOS

O método gráfico mais utilizado consiste na* comparação da função de sobrevivência (confiabilidade) do modelo proposto com o estimador de Kaplan-Meier. Neste procedimento, ajustam-se os modelos paramétricos propostos ao conjunto de dados e, a partir das estimativas dos parâmetros de cada modelo, estimam-se suas respectivas funções de sobrevivência (confiabilidade). Para o conjunto de dados, obtém-se, também, a estimativa de Kaplan-Meier para a função de sobrevivência (confiabilidade).

Finalmente, comparam-se graficamente as funções de sobrevivência (confiabilidade) estimadas para cada modelo proposto com as funções de sobrevivência (confiabilidade) estimadas pelo método de Kaplan-Meier. O(s) modelo(s) adequado(s) é(são) aquele(s) em que sua curva de sobrevivência (confiabilidade) mais se aproximar daquela do estimador de Kaplan-Meier. Em outras palavras, o(s) melhor(es) modelo(s) será(ão) aquele(s) cujos pontos no gráfico estiverem mais próximos da reta y=x.

A figura abaixo traz um exemplo da utilização desse método gráfico considerando o ajuste dos modelos Exponencial, Weibull e Log-normal, respectivamente, para um conjunto de dados:

Percebe-se que o modelo Exponencial possui os seus pontos mais distantes da reta y=x, o que pode indicar um ajuste ruim para estes dados.

Para o conjunto de dados de câncer de bexiga, podemos replicar esse gráfico para os 6 modelos ajustados e comparar com a estimativa obtida pelo método Kaplan-Meier. Primeiro, vamos obter as estimativas de Kaplan-Meier para a função de sobrevivência e guardar em um objeto:

ekm <- survfit(Surv(tempos,cens)~1, data=dados1)
Surv_ekm <- data.frame(ekm$time,ekm$surv)
colnames(Surv_ekm) <- c("Tempo", "Surv_EKM")

Similarmente, para os demais ajustes:

# Exponencial
Surv_pam <- data.frame(summary(ajuste_exp)[[1]][,c(1,2)],
                       summary(ajuste_wei)[[1]][,2],
                       summary(ajuste_lognorm)[[1]][,2],
                       summary(ajuste_loglogi)[[1]][,2],
                       summary(ajuste_gama)[[1]][,2],
                       summary(ajuste_gengama)[[1]][,2])
colnames(Surv_pam) <- c("Tempo", "Surv_exp", "Surv_wei", "Surv_lognorm", "Surv_loglogi", "Surv_gama", "Surv_gengama")

O comando a seguir junta todas as estimativas:

# Juntando as estimativas paramétricas com as estimativas não paramétricas
Surv_Final <- merge(Surv_pam,Surv_ekm,by="Tempo",all.x=T)

Enfim os gráficos que comparam as funções de sobrevivência estimadas de forma paramétrica com a função de sobrevivência estimada de forma não paramétrica são feitos com os seguintes comandos:

par(mfrow=c(2,3))
plot(Surv_Final$Surv_EKM, Surv_Final$Surv_exp, 
     xlab= "S(t): Kaplan-Meier", ylab = "S(t): Exponencial")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(Surv_Final$Surv_EKM, Surv_Final$Surv_wei, 
     xlab= "S(t): Kaplan-Meier", ylab = "S(t): Weibull" )
lines(c(0,1), c(0,1), type="l", lty=1)
plot(Surv_Final$Surv_EKM, Surv_Final$Surv_lognorm, 
     xlab= "S(t): Kaplan-Meier", ylab = "S(t): LogNormal")
lines(c(0,1), c(0,1), type="l", lty=1)
plot(Surv_Final$Surv_EKM, Surv_Final$Surv_loglogi, 
     xlab= "S(t): Kaplan-Meier", ylab = "S(t): LogLogistica" )
lines(c(0,1), c(0,1), type="l", lty=1)
plot(Surv_Final$Surv_EKM, Surv_Final$Surv_gama, 
     xlab= "S(t): Kaplan-Meier", ylab = "S(t): Gama" )
lines(c(0,1), c(0,1), type="l", lty=1)
plot(Surv_Final$Surv_EKM, Surv_Final$Surv_gengama, 
     xlab= "S(t): Kaplan-Meier", ylab = "S(t): Gama Generalizada" )
lines(c(0,1), c(0,1), type="l", lty=1)

Nota-se que, diferentemente dos demais, o modelo exponencial possi um maior desvio com respeito a linha reta. Isso pode ser um indicativo de que esse modelo pode não representar bem os dados.

Uma outra forma de comparar as funções é colocar, em um mesmo gráfico, todas as curvas de sobrevivência (confiabilidade) estimadas versus o o tempo \(t\). Alguns autores, por exemplo Nelson (1990), sugerem a comparação da função de risco acumulada estimada de forma não paramétrica com a função de risco acumulada que foi estimada pelos modelos probabilísticos.

3.2 AIC e BIC

AIC, o Critério de Informação de Akaike, é geralmente considerado como o primeiro critério de seleção do modelo. Hoje, o AIC continua a ser a mais amplamente conhecida e usada ferramenta de seleção de modelos entre os profissionais. O AIC foi introduzido por Hirotugu Akaike em 1973.

É uma métrica que mensura a qualidade de um modelo estatístico visando também a sua simplicidade (parcimônia). O AIC leva em consideração e penaliza a complexidade dos modelos e tende a favorecer a escolha de modelos mais simples.

Considere um determinado modelo estatístico, ajustado de acordo com dados observados, \(k\) o número de parâmetros de tal modelo e \(\hat{L}\) o valor máximo da função de verossimilhança. Então, o valor de AIC do modelo considerado é dado por}:

\[AIC = 2K - 2ln(\hat{L})\]

Dado uma coleção de modelos candidatos para os dados, o modelo com menor AIC é o escolhido de acordo com este critério. Assim, o AIC bonifica a qualidade de ajuste (altos valores para a função de verossimilhança) e, por outro lado, penaliza a quantidade de parâmetros do modelo. O AIC fornece, portanto, uma métrica para comparação e seleção de modelos, em que menores valores de AIC representam uma maior qualidade e simplicidade, segundo este critério. Note que o AIC não fornece uma medida de qualidade do modelo global, apenas relativa no que diz respeito à comparação entre modelos candidatos. Dessa forma, se todos os modelos propostos se ajustam mal aos dados, o AIC não explicita tal fato.

O Bayesian Information Criterion (BIC), também conhecido como Critério de Informação Bayesiano, é uma métrica estatística utilizada para selecionar o melhor modelo estatístico em uma análise de dados. Ele é baseado na teoria da probabilidade bayesiana e é amplamente utilizado em áreas como machine learning, deep learning e inteligência artificial.

O BIC é uma medida que busca equilibrar a complexidade do modelo e a qualidade do ajuste aos dados. Ele leva em consideração dois componentes principais: o número de parâmetros do modelo e a verossimilhança dos dados. O objetivo é encontrar o modelo que maximize a verossimilhança dos dados, mas penalize modelos mais complexos, evitando o overfitting.

Considere um determinado modelo estatístico, ajustado de acordo com \(n\) dados observados, \(k\) o número de parâmetros de tal modelo e \(\hat{L}\) o valor máximo da função de verossimilhança. Então, o valor de AIC do modelo considerado é dado por}:

\[BIC = k*log(n) - 2ln(\hat{L})\]

Dado uma coleção de modelos candidatos para os dados, o modelo com menor BIC é o escolhido de acordo com este critério.

No R, a extração do AIC e do BIC de cada modelo é feita de forma simples, pelos seguintes comandos:

AIC<-c(ajuste_exp$AIC,
       ajuste_wei$AIC,
       ajuste_lognorm$AIC,
       ajuste_loglogi$AIC,
       ajuste_gama$AIC,
       ajuste_gengama$AIC)
BIC<-c(ajuste_exp$BIC,
       ajuste_wei$BIC,
       ajuste_lognorm$BIC,
       ajuste_loglogi$BIC,
       ajuste_gama$BIC,
       ajuste_gengama$BIC)
Modelos<-c("Modelo Exponencial", "Modelo Weibull",
           "Modelo Log-Normal", "Modelo Log-Logístico",
           "Modelo Gama", "Modelo Gama Generalizado")
Comparacao<-data.frame(Modelos,AIC,BIC)
Comparacao

##                    Modelos      AIC      BIC
## 1       Modelo Exponencial 138.5478 139.5435
## 2           Modelo Weibull 136.2667 138.2582
## 3        Modelo Log-Normal 135.4798 137.4713
## 4     Modelo Log-Logístico 136.0611 138.0525
## 5              Modelo Gama 135.6951 137.6865
## 6 Modelo Gama Generalizado 137.3870 140.3742

Assim, segundo ambos os critérios, o modelo Log-Normal é o que mais se ajusta bem aos dados. Com a escolha deste modelo, podemos retormar as suas estimativas do tempo médio de ocorrência do evento de interesse, o tempo mediano de ocorrência do evento de interesse, o percentil 80 do tempo de ocorrência do evento de interesse e tempo médio residual de ocorrência do evento de interesse no instante t=20 meses.

Mean_lognorm<-summary(ajuste_lognorm,type="mean",cl=0.95,se=T)
Mean_lognorm

##  
##        est      lcl      ucl       se
## 1 20.28026 13.66221 30.83926 4.649566

Median_lognorm<-summary(ajuste_lognorm,type="median",cl=0.95,se=T)
Median_lognorm

##  
##        est     lcl      ucl       se
## 1 15.13751 10.5531 21.81799 2.824792

quantile80_lognorm<-summary(ajuste_lognorm,type="quantile",
                        quantiles=0.8,cl=0.95,se=T)
quantile80_lognorm

##  
##   quantile    est      lcl      ucl       se
## 1      0.8 28.814 19.14536 44.62491 6.734768

Resid20_lognorm<-reslifefsr(ajuste_lognorm,20)
Resid20_lognorm

## [1] 17.15706

Além de considerar medidas estatísticas de qualidade de ajuste, pode ser útil uma boa conversa com o médico/engenheiro para obter informações extra que auxiliem na escolha do melhor modelo que vai ser ajustado aos dados.

3.3 EXERCÍCIOS

1- Tratamento Quimioterápico:

tempos<-c(7, 8, 10, 12, 13, 14, 19, 23, 25, 26, 27, 
          31, 31, 49, 59, 64, 87, 89, 107, 117, 119, 
          230, 233, 130, 148, 153, 156, 159, 191, 222, 
          200, 203, 210, 220, 220, 228, 235, 240, 240,
          240, 241, 245, 247, 248, 250)
cens<-c(1,1,1,1,1,0,1,1,0,1,1,1,0,1,0,0,1,1,1,1,1,0,
        0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0)
dados1<-data.frame(tempos,cens)

Ajuste todos os 6 modelos paramétricos e obtenha também as estimativas de Kaplan-Meier.
Faça o método gráfico comparando Kaplan-Meier x Modelos Paramétricos. Qual a conclusão?
Obtenha o AIC e BIC dos modelos ajustados. Qual a conclusão?
Com o modelo ajustado, estime o tempo mediano até a ocorrência dos primeiros sinais de alterações indesejadas.

2- Isolante Elétrico:

tempos<-c(151, 164, 336, 365, 403, 454, 455, 473, 538, 
          577, 592, 628, 632, 647, 675, 727, 785, 801, 
          811, 816, 867, 893, 930, 937, 976, 1008, 1040, 
          1051, 1060, 1183, 1329, 1334, 1379, 1380, 1633, 
          1769, 1827, 1831, 1849, 2016, 2282, 2415, 2430, 
          2686, 2729, rep(2729,15))
cens<-c(rep(1,45),rep(0,15))
dados1<-data.frame(tempos,cens)

Ajuste todos os 6 modelos paramétricos e obtenha também as estimativas de Kaplan-Meier.
Faça o método gráfico comparando Kaplan-Meier x Modelos Paramétricos. Qual a conclusão?
Obtenha o AIC e BIC dos modelos ajustados. Qual a conclusão?
Com o modelo ajustado, estime o tempo G de garantia o qual 10% dos isolantes falharam até o tempo G.

3- AIDS

Ajuste todos os 6 modelos paramétricos e obtenha também as estimativas de Kaplan-Meier.
Faça o método gráfico comparando Kaplan-Meier x Modelos Paramétricos. Qual a conclusão?
Obtenha o AIC e BIC dos modelos ajustados. Qual a conclusão?
Com o modelo ajustado, estime o tempo mediano até a morte desses pacientes.