ANÁLISE DE SOBREVIVÊNCIA E CONFIABILIDADE
AULA 1: ESTIMAÇÃO NÃO-PARAMÉTRICA
\[\\[0.05in]\]1 PACOTES NECESSÁRIOS
Para esta primeira aula prática de Análise de Sobrevivência e Confiabilidade, precisaremos dos seguintes pacotes:
require(ggplot2)
require(survival)
require(survminer)Pacote
ggplot2: Pacote fundamental para construção de gráficos noRPacote
survival: É o pacote básico para o estudo de Análise de Sobrevivência e Confiabilidade. Contém as principais rotinas, incluindo definição de objetosSurv, curvas Kaplan-Meier, modelos Cox e modelos paramétricos.Pacote
survminer: O pacotesurvivaldoRajusta e traça curvas de sobrevivência (confiabilidade) usando gráficos básicos doR. Os gráficos do pacotesurvminersão baseados na estrutura gráfica do pacoteggplot2:
2 BASES DE DADOS
Esta seção é dedicada à descrição das bases de dados utilizadas nesta aula.
2.1 HEPATITE
Um estudo clínico aleatorizado com o objetivo de investigar o efeito da terapia com esteroide no tratamento de hepatite viral aguda. Vinte e nove pacientes com esta doença foram aleatorizados para receber um placebo ou o tratamento com esteroide. Cada paciente foi acompanhado por 16 semanas ou até a morte (evento de interesse) ou até a perda de acompanhamento. Os dados são retratados a seguir:
2.2 REINCIDÊNCIA DE TUMOR SÓLIDO
Considere este outro exemplo em que se deseja avaliar os tempos de reincidência de 10 pacientes com tumor sólido. Dos 10 pacientes, seis deles apresentaram reincidência aos 3; 6,5; 6,5; 10; 12 e 15 meses de seus respectivos ingressos no estudo; um deles não re tornou após 8,4 meses de acompanhamento e três deles permaneceram sem reincidência após 4; 5,7 e 10 meses de acompanhamento. Os esquemas que ilustram hipoteticamente o acompanhamento dos pacientes deste estudo são apresentados na figura abaixo:
Do esquema (a), apresentado nesta figura, observa-se que o experimento foi planejado para durar 18 meses e teve início com três pacientes. Após ter decorrido um mês do início do experimento, ocorreu o ingresso do quarto paciente e assim sucessivamente, até o décimo paciente, que ingressou após decorridos 14 meses de andamento do experimento.
O esquema apresentado em (b) mostra, por sua vez, quanto tempo cada paciente permaneceu no estudo. Note que o uso do referencial “zero” neste último esquema possibilita que o tempo até a ocorrência do evento de interesse ou da censura de cada paciente sob estudo seja observado de maneira mais facil e direta do que no esquema (a).
2.3 ISOLANTE ELÉTRICO
Os dados mostrados a seguir representam o tempo até a ruptura de um tipo de isolante elétrico sujeito a uma tensão de estresse de 35 Kvolts. O teste consistiu em deixar 25 destes isolantes funcionando até que 15 deles falhassem (censura do tipo II), obtendo-se os seguintes resultados (em minutos):
0,19 0,78 0,96 1,31 2,78 3,16 4,67 4,85 6,50 7,35 8,27 12,07 32,52 33,91 36,71
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. Para comparação da sobrevivência, vamos ver como
o sexo impacta a sobrevivência do paciente. Essa base de dados
está armazenada no objeto DadosAIDS.csv
3 ESTIMADOR DE KAPLAN-MEIER PARA \(S(t)\)
Proposto por Kaplan e Meier em 1958 para estimar a função de sobrevivência (confiabilidade) em dados censurados. Também é chamado de estimador limite-produto. Esse estimador é, sem dúvidas, o mais utilizado em estudos clínicos e vem ganhando cada vez mais espaço em estudos de confiabilidade. É uma adaptação do estimador da função de sobrevivência (confiabilidade) para dados não censurados.
Como foi visto em sala, o estimador de Kaplan-Meier para a função de sobrevivência (confiabilidade) depende das quantidades \(n_j\) e \(d_j\) e é calculado recursivamente. A fórmula geral para o seu cálculo é:
\[\hat{S}_{KM}(t)=\prod_{j:t_j\leq t}\left(1-\frac{d_j}{n_j}\right),\]
- \(t_1 < t_2 < · · · < t_k\) são os \(k\) tempos distintos de ocorrência do evento de interesse.
- \(d_j\) é o número de ocorrências do evento em \(t_j\) .
- \(n_j\) é o número de elementos sob risco no tempo \(t_j\) . Ou seja, o número de elementos que não apresentaram o evento e nem foram censuradas até o instante imediatamente anterior a \(t_j\).
A variância do estimador de Kaplan-Meier é dada pela fórmula de Greenwood:
\[Var\left(\hat{S}_{KM}(t)\right)=\left[\hat{S}_{KM}(t)\right]^2\sum_{j:t_j<t}\left(\frac{d_j}{n_j(n_j-d_j)}\right)\]
Como \(\hat{S}_{KM}(t)\), para \(t\) fixo, tem distribuição normal assintoticamente, o intervalo aproximado de \(100(1-\alpha)\%\) de confiança para \(S(t)\) é dado por
\[\hat{S}_{KM}(t) \pm z_{\alpha/2}\sqrt{Var\left(\hat{S}_{KM}(t)\right)},\]
em que \(\alpha/2\) denota o \(\alpha/2-\)percentil da distribuição Normal padrão.
Entretanto, para valores extremos de \(t\), o intervalo de confiança anterior pode apresentar um limite inferior negativo ou um limite superior maior que 1. Isso não faz sentido, uma vez que \(S(t)\) é uma probabilidade. Para corrigir isso, utiliza-se a seguinte transformação para \(\hat{S}_{KM}(t)\):
\[\hat{U}_{KM}(t)=ln\left[-ln\left(\hat{S}_{KM}(t)\right)\right]\]
A variância dessa variável transformada é:
\[Var\left(\hat{U}_{KM}(t)\right)=\frac{\sum_{j:t_j<t}\frac{d_j}{n_j(n_j-d_j)}}{\left[ln\left(\hat{S}_{KM}(t)\right)\right]^2}\]
Assim um outro intervalo aproximado de \(100(1-\alpha)\%\) de confiança para \(S(t)\) é dado por
\[\left[\hat{S}_{KM}(t)\right]^{exp\left\{\pm z_{\alpha/2}\sqrt{Var(\hat{U}_{KM}(t))}\right\}}\]
Esse intervalo assume valores no suporte [0,1]. É importante ressaltar que há também a transformação log para construção de intervalos de confiança, que veremos nessa aula prática
Para exemplificar, vamos utilizar os dados de Hepatite. Lembre-se que sempre observamos uma amostra de \(n\) elementos. Os dados observados de sobrevivência (confiabilidade) para o elemento \(i\), com \(i = 1, 2, 3, ..., n\), são representados, em geral, pelo par
\[(t_i , \delta_i ),\]
em que:
- \(t_i\) é o tempo observado até a ocorrência do evento de interesse ou o tempo observado de censura para o elemento \(i\).
- \(\delta_i\) é uma função indicadora de ocorrência do evento:
\(\delta_i=\begin{cases} 1, \quad \textrm{se } t_i \textrm{ é um tempo de ocorrência do evento de interesse}\\ 0, \quad \textrm{se } t_i \textrm{ é um tempo censurado} \end{cases}\)
Como a base é pequena, vamos criar os dados manualmente:
tempos <- c(1,2,3,3,3,5,5,16,16,16,16,16,16,16,16,1,1,1,1,4,5,7,8,10,10,12,16,16,16)
cens <- c(0,0,1,1,0,0,0,0,0,0,0,0,0,0,0,1,1,1,0,0,1,1,1,1,0,0,0,0,0)
grupos <- c(rep("Controle",15),rep("Esteroide",14))
Dados1 <- data.frame(tempos,cens,grupos)- O objeto
temposguarda os valores dos tempos \(t_i\), para \(i=1,2,...,29\) - O objeto
censguarda a indicadora \(\delta_i\), para \(i=1,2,...,29\) - o objeto
gruposindica o grupo o qual os pacientes estão submetidos: Controle ou Esteroide
As estimativas de Kaplan-Meier para os grupos
Controle e Esteroide, associados aos dados de Hepatite, podem
ser obtidas a partir da função survfit(). Essa
função possui os seguintes argumentos
básicos:
data: A base de dados.formula: Essa fórmula quase sempre tem essa estrutura:Surv(tempos,cens)~gruposSe você não tem grupos em sua base de dados, utilizeSurv(tempos,cens)~1.conf.int: Nível de confiança do intervalo.conf.type: Para o intervalo de confiança básico, utilizeconf.type = "plain". Para o intervalo de confiança com a variável transformada \(\hat{U}_{KM}(t)=ln\left(\hat{S}_{KM}(t)\right)\) useconf.type = "log"(default doR). Para o intervalo de confiança com a variável transformada \(\hat{U}_{KM}(t)=ln\left[-ln\left(\hat{S}_{KM}(t)\right)\right]\), utilizeconf.type = "log-log".
Assim, vamos estimar a função de sobrevivência pelo estimador de Kaplan-Meier a partir dos dados de Hepatite utilizando um intervalo com 95% de confiança sem transformação:
ekm<- survfit(data=Dados1,formula = Surv(tempos,cens)~grupos,
conf.type = "plain", conf.int=0.95)
summary(ekm)## Call: survfit(formula = Surv(tempos, cens) ~ grupos, data = Dados1,
## conf.type = "plain", conf.int = 0.95)
##
## grupos=Controle
## time n.risk n.event survival std.err lower 95% CI
## 3.000 13.000 2.000 0.846 0.100 0.650
## upper 95% CI
## 1.000
##
## grupos=Esteroide
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 14 3 0.786 0.110 0.571 1.000
## 5 9 1 0.698 0.128 0.448 0.948
## 7 8 1 0.611 0.138 0.340 0.882
## 8 7 1 0.524 0.143 0.243 0.805
## 10 6 1 0.437 0.144 0.155 0.718Agora vamos estimar a função de sobrevivência pelo estimador de Kaplan-Meier a partir dos dados de Hepatite utilizando um intervalo com 95% de confiança com a transformação \(\hat{U}_{KM}(t)=ln\left(\hat{S}_{KM}(t)\right)\):
ekm2<- survfit(data=Dados1,formula = Surv(tempos,cens)~grupos,
conf.type = "log", conf.int=0.95)
summary(ekm2)## Call: survfit(formula = Surv(tempos, cens) ~ grupos, data = Dados1,
## conf.type = "log", conf.int = 0.95)
##
## grupos=Controle
## time n.risk n.event survival std.err lower 95% CI
## 3.000 13.000 2.000 0.846 0.100 0.671
## upper 95% CI
## 1.000
##
## grupos=Esteroide
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 14 3 0.786 0.110 0.598 1.000
## 5 9 1 0.698 0.128 0.488 0.999
## 7 8 1 0.611 0.138 0.392 0.952
## 8 7 1 0.524 0.143 0.306 0.896
## 10 6 1 0.437 0.144 0.229 0.832Agora vamos estimar a função de sobrevivência pelo estimador de Kaplan-Meier a partir dos dados de Hepatite utilizando um intervalo com 95% de confiança com a transformação \(\hat{U}_{KM}(t)=ln\left[-ln\left(\hat{S}_{KM}(t)\right)\right]\):
ekm3<- survfit(data=Dados1,formula = Surv(tempos,cens)~grupos,
conf.type = "log-log", conf.int=0.95)
summary(ekm3)## Call: survfit(formula = Surv(tempos, cens) ~ grupos, data = Dados1,
## conf.type = "log-log", conf.int = 0.95)
##
## grupos=Controle
## time n.risk n.event survival std.err lower 95% CI
## 3.000 13.000 2.000 0.846 0.100 0.512
## upper 95% CI
## 0.959
##
## grupos=Esteroide
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 14 3 0.786 0.110 0.472 0.925
## 5 9 1 0.698 0.128 0.378 0.876
## 7 8 1 0.611 0.138 0.298 0.819
## 8 7 1 0.524 0.143 0.227 0.754
## 10 6 1 0.437 0.144 0.164 0.683Para todas as saídas, observe que:
time: Os tempos distintos de ocorrência do evento (\(t_j\))n.risk: \(n_j\)n.event: \(d_j\)survival: Estimativa de Kaplan-Meier para \(S(t_j)\)std.error: Desvio padrão associado à estimativa de Kaplan-Meier (raiz quadrada da variância de Greenwood)lower CIeupper CI: limites do intervalo de confiança solicitado e com a confiança especificada.
Para visualizar o gráfico com a função de sobrevivência
estimada, vamos usar o objeto ekm3. O gráfico é
feito com a função ggsurvplot() do pacote
survminer com os comandos:
ggsurvplot(ekm3,legend.title = " ",
legend.labs = c("Controle", "Esteroide"),
legend = c(0.15, 0.2))+
xlab("Tempo (semanas)")+
ylab("Sobrevivência estimada") 
Para incluir o intervalo de confiança, basta
acrescentar o comando conf.int=T
ggsurvplot(ekm3,legend.title = " ",
legend.labs = c("Controle", "Esteroide"),
legend = c(0.15, 0.4),
conf.int=T)+
xlab("Tempo (semanas)")+
ylab("Sobrevivência estimada") 
Para mais detalhes sobre o a função
ggsurvplot, basta consultar o site ggsurvplot
Assim, um estimador de Kaplan-Meier para \(\Lambda(t)\) é dado por:
\[\hat{\Lambda}_{KM}(t)=-ln(\hat{S}_{KM}(t))\]
No R, os valores dessa função podem ser obtidos a partir
dos seguintes comandos:
racum_controle<- -log(unique(ekm$surv))[1:2]
racum_controle## [1] 0.0000000 0.1670541racum_esteroide<- -log(unique(ekm$surv))[c(1,3:7)]
racum_esteroide## [1] 0.0000000 0.2411621 0.3589451 0.4924765 0.6466272 0.8289487O gráfico com as estimativas da função de risco acumulada
pode ser feito com o acrescimo do comando fun="cumhaz", aos
gráficos anteriores:
ggsurvplot(ekm3,legend.title = " ",
legend.labs = c("Controle", "Esteroide"),
legend = c(0.15, 0.8),
fun="cumhaz")+
xlab("Tempo (semanas)")+
ylab("Risco acumulado") 
Também é possível acrescentar os intervalos de confiança:
ggsurvplot(ekm3,legend.title = " ",
legend.labs = c("Controle", "Esteroide"),
legend = c(0.15, 0.8),
fun="cumhaz",
conf.int=T)+
xlab("Tempo (semanas)")+
ylab("Risco acumulado") 
4 ESTIMADOR DE NELSON-AALEN PARA \(S(t)\)
Este estimador é mais recente que o estimador de Kaplan-Meier e ele se diferencia por começar o processo de estimação a partir da função de risco acumulada e é dada por:
\[\begin{equation} \hat{\Lambda}_{NA}(t)=\sum_{j:t_j\leq t}\left(\frac{d_j}{n_j}\right) \end{equation}\]
A notação é a mesma utilizada no estimador de Kaplan-Meier. A ideia por trás dessa fórmula reside no fato de que a função de risco acumulado no tempo \(t\) traz exatamente o risco acumulado de apresentar o evento de interesse no tempo \(t\). Assim, é bastante intuitivo somar as probabilidades \(\left(\frac{d_j}{n_j}\right)\) que retratam, cada uma, o risco de ocorrência do evento de interesse no tempo \(t_j\).
A partir do estimador de Nelson-Aalen para a função de risco acumulado, pode-se estimar a função de sobrevivência (confiabilidade) \(S(t)\), pois sabemos que \(S(t)=exp(-\Lambda(t))\)Assim, um estimador de Nelson-Aalen para \(S(t)\) é dado por: \[\hat{S}_{NA}(t)=exp(-\hat{\Lambda}_{NA}(t))\]
A variância do estimador de Nelson Aalen é dada pela fórmula: \[Var\left(\hat{S}_{NA}(t)\right)=\left[\hat{S}_{NA}(t)\right]^2\sum_{j:t_j<t}\left(\frac{d_j}{n_j^2}\right)\]. Como \(\hat{S}_{NA}(t)\), para \(t\) fixo, tem distribuição normal assintoticamente, o intervalo aproximado de \(100(1-\alpha)\%\) de confiança para \(S(t)\) é dado por
\[\hat{S}_{NA}(t) \pm z_{\alpha/2}\sqrt{Var\left(\hat{S}_{NA}(t)\right)},\]
em que \(\alpha/2\) denota o \(\alpha/2-\)percentil da distribuição Normal padrão.
No R, as estimativas de Nelson-Aalen para a
função de sobrevivência associada aos dados de Hepatite podem ser
obtidas também a partir da função survfit() com a inclusão
do argumento stype=2. Para um intervalo de
confiança construído sem transformação, seguem os comandos:
na<-survfit(data=Dados1,formula = Surv(tempos,cens)~grupos,
stype=2, conf.type = "plain", conf.int=0.95)
summary(na)## Call: survfit(formula = Surv(tempos, cens) ~ grupos, data = Dados1,
## stype = 2, conf.type = "plain", conf.int = 0.95)
##
## grupos=Controle
## time n.risk n.event survival std.err lower 95% CI
## 3.0000 13.0000 2.0000 0.8574 0.0933 0.6746
## upper 95% CI
## 1.0000
##
## grupos=Esteroide
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 14 3 0.807 0.0999 0.611 1.000
## 5 9 1 0.722 0.1201 0.487 0.958
## 7 8 1 0.637 0.1326 0.377 0.897
## 8 7 1 0.553 0.1394 0.279 0.826
## 10 6 1 0.468 0.1414 0.190 0.745Para um intervalo de confiança com a transformação log:
na2<-survfit(data=Dados1,formula = Surv(tempos,cens)~grupos,
stype=2, conf.type = "log", conf.int=0.95)
summary(na2)## Call: survfit(formula = Surv(tempos, cens) ~ grupos, data = Dados1,
## stype = 2, conf.type = "log", conf.int = 0.95)
##
## grupos=Controle
## time n.risk n.event survival std.err lower 95% CI
## 3.0000 13.0000 2.0000 0.8574 0.0933 0.6928
## upper 95% CI
## 1.0000
##
## grupos=Esteroide
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 14 3 0.807 0.0999 0.633 1.000
## 5 9 1 0.722 0.1201 0.521 1.000
## 7 8 1 0.637 0.1326 0.424 0.958
## 8 7 1 0.553 0.1394 0.337 0.906
## 10 6 1 0.468 0.1414 0.259 0.846Para um intervalo de confiança com a transformação log-log:
na3<-survfit(data=Dados1,formula = Surv(tempos,cens)~grupos,
stype=2, conf.type = "log-log", conf.int=0.95)
summary(na3)## Call: survfit(formula = Surv(tempos, cens) ~ grupos, data = Dados1,
## stype = 2, conf.type = "log-log", conf.int = 0.95)
##
## grupos=Controle
## time n.risk n.event survival std.err lower 95% CI
## 3.0000 13.0000 2.0000 0.8574 0.0933 0.5406
## upper 95% CI
## 0.9623
##
## grupos=Esteroide
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 14 3 0.807 0.0999 0.515 0.933
## 5 9 1 0.722 0.1201 0.412 0.887
## 7 8 1 0.637 0.1326 0.328 0.833
## 8 7 1 0.553 0.1394 0.255 0.773
## 10 6 1 0.468 0.1414 0.191 0.706Para fazer o gráfico com a função de sobrevivência
estimada, basta usar os comandos já utilizados do pacote
survminer (vamos considerar o objeto na3):
ggsurvplot(na3,legend.title = " ",
legend.labs = c("Controle", "Esteroide"),
legend = c(0.15, 0.2))+
xlab("Tempo (semanas)")+
ylab("Sobrevivência estimada") 
Mais uma vez, para acrescentar os intervalos, acrescente o
comando conf.int=T.
ggsurvplot(na3,legend.title = " ",
legend.labs = c("Controle", "Esteroide"),
legend = c(0.15, 0.2),
conf.int=T)+
xlab("Tempo (semanas)")+
ylab("Sobrevivência estimada") 
Assim como antes, no R, os valores das estimativas da
função de risco acumulado podem ser obtidos a partir dos seguintes
comandos:
racum_controle<- -log(unique(na$surv))[1:2]
racum_controle## [1] 0.0000000 0.1538462racum_esteroide<- -log(unique(na$surv))[c(1,3:7)]
racum_esteroide## [1] 0.0000000 0.2142857 0.3253968 0.4503968 0.5932540 0.7599206O gráfico com as estimativas da função de risco acumulada
pode ser feito com o acrescimo do comando fun="cumhaz", aos
gráficos anteriores:
ggsurvplot(na3,legend.title = " ",
legend.labs = c("Controle", "Esteroide"),
legend = c(0.15, 0.8),
fun="cumhaz")+
xlab("Tempo (semanas)")+
ylab("Risco acumulado") 
Também é possível acrescentar os intervalos de confiança:
ggsurvplot(na3,legend.title = " ",
legend.labs = c("Controle", "Esteroide"),
legend = c(0.15, 0.8),
fun="cumhaz",
conf.int=T)+
xlab("Tempo (semanas)")+
ylab("Risco acumulado") 
5 EXERCÍCIOS
Utilize a transformação log-log, quando necessário
1- Base de dados Tumor Sólido:
A partir da descrição dos dados, construa a base de dados e guarde em um objeto chamado
Dados2.Obtenha as estimativas para as funções de sobrevivência a partir do estimador de Kaplan-Meier. Faça um gráfico com o intervalo construído.
Obtenha as estimativas para as funções de sobrevivência a partir do estimador de Nelson-Aalen. Faça um gráfico com o intervalo construído.
Faça o gráfico da função de risco acumulada para o ajuste da letra b).
Estime o tempo mediano de ocorrência do evento de interesse (usando interpolação linear)
2- Base de dados Isolantes Elétricos.
A partir da descrição dos dados, construa a base de dados e guarde em um objeto chamado
Dados3.EObtenha as estimativas para as funções de sobrevivência a partir do estimador de Kaplan-Meier. Faça um gráfico com o intervalo construído.
Faça o gráfico da função de risco acumulada para o ajuste escolhido para a letra b).
Estime o tempo mediano de ocorrência do evento de interesse (usando interpolação linear)
3- Base de dados AIDS.
Leia a base de dados e guarde em um objeto
Dados4.Considerando os dados completos (sem estratificar por sexo), obtenha as estimativas para a função de sobrevivência a partir do estimador de Kaplan-Meier e faça um gráfico com o intervalo construído
Faça o gráfico da função de risco acumulada para o ajuste da letra b).
Considerando o sexo, obtenha as estimativas para as funções de sobrevivência a partir do estimador de Kaplan-Meier e faça um gráfico com o intervalo construído. Qual sexo tende a sobreviver mais, descritivamente?
Faça o gráfico da função de risco acumulada para o ajuste da letra d).
Com a interpolação linear em cada sexo, estime:
- O tempo mediano de ocorrência do evento de interesse.
- O tempo para o qual 25% das pessoas morrem por AIDS.