library(rmarkdown); library(knitr); library(readxl); library(survival); library(survminer); library(coin)

Introdução

A Análise de Sobrevivência é uma das áreas da estatística que mais cresceu nas últimas duas décadas do século passado. A razão deste crescimento é o desenvolvimento e aprimoramento de técnicas estatísticas combinado com computadores cada vez mais velozes.

A análise de sobrevivência centra-se na descrição, para um determinado indivíduo ou grupo de indivíduos, de um ponto definido de evento chamado a falha (ocorrência de uma doença, cura de uma doença, óbito, recaída após resposta ao tratamento…) que ocorre após um período de tempo chamado tempo de falha (ou tempo de seguimento em estudos baseados em coorte/população) durante o qual os indivíduos são observados. Para determinar o tempo de falha, é então necessário definir um tempo de origem (que pode ser a data de inclusão, a data de diagnóstico…).

A análise de sobrevivência é um conjunto de procedimentos estatísticos para análise de dados em que a variável de resultado de interesse é o tempo até a ocorrência de um evento, ou seja, é a área que se debruça sobre estudos de tempos até à ocorrência de determinado evento de interesse.

Evento define-se como um qualquer acontecimento cuja ocorrência em função do tempo, que passou entre a observação inicial e o tempo em que esse evento ocorre, interessa estudar. Este está frequentemente definido como “morte” e o tempo até à ocorrência do evento (terminal) é frequentemente definido como tempo de vida ou tempo de sobrevivência.

O nome “Análise de Sobrevivência” advém desta técnica ter sido inicialmente desenvolvida para aplicações nas ciências biomédicas relacionadas com os tempos de sobrevivência de pacientes com doenças terminais (eg.. cancro) ou com condições de risco de vida (eg, ocorrência de ataque cardíaco ou tempo de falência de outros órgãos).

É normalmente expressa através da probabilidade de sobrevivência que é a probabilidade de o evento de interesse não ter ocorrido por uma duração $t$.

Este tempo é denominado tempo de falha, podendo ser o tempo até a morte do paciente, bem como até a cura ou recidiva de uma doença. Em estudos de câncer, é usual o registro das datas correspondentes ao diagnóstico da doença, à remissão (após o tratamento, o paciente fica livre dos sintomas da doença), à recorrência da doença (recidiva) e à morte do paciente. O tempo de falha pode ser, por exemplo, do diagnóstico até a morte ou da remissão até a recidiva.

A principal característica de dados de sobrevivência é a presença de censura, que é a observação parcial da resposta. Isto se refere a situações em que, por alguma razão, o acompanhamento do paciente foi interrompido, seja porque o paciente mudou de cidade, o estudo terminou para a análise dos dados ou, o paciente morreu de causa diferente da estudada. Isto significa que toda informação referente à resposta se resume ao conhecimento de que o tempo de falha é superior àquele observado. Sem a presença de censura, as técnicas estatísticas clássicas, como análise de regressão e planejamento de experimentos, poderiam ser utilizadas na análise deste tipo de dados, provavelmente usando uma transformação para a resposta.

Suponha, por exemplo, que o interesse seja o de comparar o tempo médio de vida de três grupos de pacientes. Se não houver censuras, pode-se usar as técnicas usuais de análise de variância para se fazer tal comparação.

No entanto, se houver censuras, que é o mais provável, tais técnicas não podem ser utilizadas pois elas necessitam de todos os tempos de fálha.

Contextos de Utilização

O termo análise de sobrevivência refere-se basicamente a situações médicas envolvendo dados censurados. Entretanto, condições similares ocorrem em outras áreas em que se usam as mesmas técnicas de análise de dados.

Em contextos médicos, o evento utilizado é frequentemente a morte do paciente. Nesse contexto, o teste pode ser usado, por exemplo, para determinar o tempo médio de vida de uma pessoa com uma determinada doença. O conceito de evento não precisa, no entanto, de ser a morte ou a doença. Ele pode ser algo positivo como a remissão ou a recuperação.

A análise de sobrevivência é usada para vários fins e diversos eventos dentro da área médica, como, por exemplo, para estimar:

O tempo de vida que têm os portadores de uma determinada doença;
O tempo que demora até a reincidência ou reaparecimento de uma determinada doença;
O tempo que demora uma determinada doença a desenvolver-se. ou
O tempo que demora um paciente a passar do estado de vigilância de uma determinada doença até iniciar o tratamento contra essa doença.

Mas esta análise de sobrevivência pode ser usada em muitas outras áreas e em eventos do dia-a-dia.

Por exemplo, na área da economia e da gestão, o teste pode ser usado para:

Determinar o tempo médio que os licenciados e os não licenciados demoram a encontrar um emprego depois de ficarem desempregados ou determinar o tempo médio até os recém-licenciados encontrarem o primeiro emprego; Nestes casos, encontrar emprego é o evento.
Determinar o tempo médio que um fica em stock sem ser vendido. Aqui, o evento é o produto ser vendido.

Na área da engenharia, o teste é usado por exemplo para:

Determinar o tempo médio de duração de um material; ou
Estimar o tempo de produção de uma máquina até à próxima avaria;

Existem ainda muitas outras aplicações desta análise noutras áreas. Por exemplo:

Na banca, esta análise pode ser utilizada para analisar os prazos de reembolso total de um empréstimo.
Nos seguros, a análise de sobrevivência é usada para avaliar a rentabilidade das apólices e portanto determinar os prémios dos seguros.
Na sociologia, a análise de sobrevivência é usada para avaliar o tempo médio que um adulto demora a casar-se ou a determinar quantos anos mantêm o casamento. Neste caso, o divórcio é o evento que está a ser avaliado.

A Análise de sobrevivência permite responder às seguintes perguntas/questões:

Qual é e Probabilidade de um paciente sobreviver por mais de dois anos após o diagnóstico da AIDS?
Qual é o risco de um paciente diagnosticado com AIDS vir a óbito em três anos após o diagnóstico?
Qual seria o número esperado de óbitos em pacientes acompanhado por cinco anos?
Qual é o tempo mediano de sobrevivência?

Tipologia dos Estudos

No presente texto, o interesse envolve situações em que a variável resposta é o tempo até a ocorrência de um evento de interesse.

Estudos Clínicos

Os estudos clínicos são investigações científicas realizadas com o objetivo de verificar uma determinada hipótese de interesse. Estas investigações são conduzidas coletando dados e analisando-os por meio de métodos estatísticos. Em geral, estes estudos podem ser divididos nas seguintes três etapas:

Formulação da hipótese de interesse.
Planejamento e coleta dos dados.
Análise estatística dos dados para testar a hipótese formulada.

Estas etapas são comuns em qualquer estudo envolvendo análise estatística de dados.

A primeira etapa de um estudo clínico é gerada pela curiosidade científica do pesquisador. Identificar fatores de risco ou prognóstico para uma doença é um objetivo que aparece com freqüência em estudos clínicos. A comparação de drogas ou diferentes opções terapêuticas é outro objetivo usualmente encontrado neste tipo de estudo.

Já na segunda etapa, o pesquisador define o desenho do estudo, seleciona a população-alvo, estabelece critérios de inclusão e exclusão, determina o tamanho da amostra e escolhe os métodos adequados para a coleta dos dados. Nesta fase, é fundamental garantir que os dados sejam coletados de forma padronizada, ética e confiável, de modo a minimizar vieses e erros que possam comprometer os resultados do estudo.

Enquanto que na terceira etapa, que corresponde à análise estatística dos dados, são aplicadas técnicas estatísticas apropriadas para resumir, explorar e interpretar os dados coletados. Nesta fase, a hipótese inicialmente formulada é testada, permitindo ao pesquisador avaliar se há evidências estatísticas suficientes para aceitá-la ou rejeitá-la. Os resultados obtidos são então interpretados à luz do conhecimento científico existente, possibilitando a formulação de conclusões e recomendações clínicas.

Em análise de sobrevivência, a resposta é por natureza longitudinal. O delineamento de estudos com respostas dessa natureza pode ser observacional ou experimental, assim como ele pode ser etrospectivo ou prospectivo.

As quatro formas básicas de estudos clínicos são: descritivo, casocontrole, coorte e clínico aleatorizado.

Os três primeiros são observacionais e o quarto é experimental, pois existe a intervenção do pesquisador ao alocar, de forma aleatória, tratamento ao paciente.

O uso das técnicas de análise de sobrevivência é mais frequente nos estudos de coorte e ensaios clínicos. Entretanto, o seu uso é também possível nos demais estudos, desde que os tempos até a ocorrência do evento de interesse possam ser claramente definidos e obtidos.

O estudo envolvendo somente uma amostra, usualmente de doentes, é descritivo, pois não existe um grupo de comparação. Nestes estudos o objetivo é freqüentemente a identificação de fatores de prognóstico para a doença em estudo.

O estudo caso-controle é usualmente retrospectivo. Dois grupos, um de doentes (casos) e outro de não-doentes (controles), são comparados em relação à exposição a um ou mais fatores de interesse. Este estudo é simples, de baixo custo e rápido, pois a informação já se encontra disponível.

No entanto, ele sofre algumas limitações por estar sujeito a alguns tipos de vícios. Esses vícios estão relacionados à informação disponível sobre a história da exposição, assim como a incerteza sobre a escolha do grupo controle.

As limitações dos estudos caso-controle podem ser vencidas pelos estudos conhecidos por coorte. Em um estudo de coorte, dois grupos, um exposto e outro não-exposto ao fator de interesse, são acompanhados por um período de tempo registrando-se a ocorrência da doença ou evento de interesse.

A vantagem deste estudo sobre o caso-controle é poder avaliar a comparabilidade dos grupos no início do estudo e identificar as variáveis de interesse a serem medidas. Por outro lado, é um estudo longo e mais caro, pois os indivíduos são acompanhados por um período de tempo muitas vezes superior a um ano. Também não é um estudo indicado para doenças consideradas raras.

A forma mais consagrada de pesquisa clínica é o estudo clínico aleatorizado, que é importante por ser experimental. Isto significa que existe a intervenção direta do pesquisador ao alocar, de forma aleatória, tratamento ao paciente. Este procedimento garante a comparabilidade dos grupos. Este estudo é bastante analisado na literatura e pode-se citar os seguintes livros, entre outros, Pocock (1983) e Friedman et al. (1998).

Caracterizando os Dados

Os conjuntos de dados de sobrevivência são caracterizados pelos tempos de falha e, muito freqüentemente, pelas censuras. Estes dois componentes constituem a resposta. Em estudos clínicos, um conjunto de covariáveis é também, geralmente, medido em cada paciente.

Os seguintes elementos constituem o tempo de falha: o tempo inicial, a escala de medida e o evento de interesse (falha).

Тempo de Falha

O tempo de início do estudo deve ser precisamente definido. Os indivíduos devem ser comparáveis na origem do estudo, com exceção de diferenças medidas pelas covariáveis. Em um estudo clínico aleatorizado, a data da aleatorização é a escolha natural para a origem do estudo. A data do início do tratamento de doenças ou do diagnóstico também são outras escolhas possíveis.

A escala de medida é quase sempre o tempo real ou “de relógio”, apesar de existirem outras alternativas. Em testes de engenharia podem surgir outras escalas de medida, como o número de ciclos, a quilometragem de um carro ou qualquer outra medida de carga.

O terceiro elemento é o evento de interesse. Estes eventos são, na maioria dos casos, indesejáveis e, como já mencionado, chamados de falha.

É importante, em estudos de sobrevivência, definir de forma clara e precisa o que vem a ser a falha. Em algumas situações, a definição de falha já é clara, tais como morte ou recidiva, mas em outras pode assumir termos ambíguos. Por exemplo, fabricantes de produtos alimentícios desejam saber informações sobre o tempo de vida de seus produtos expostos em balcões frigoríficos de supermercados.

O tempo de falha vai do tempo inicial de exposição (chegada ao supermercado) até o produto ficar "inapropriado ao consumo". Este evento deve ser claramente definido antes do estudo ter seu início. Por exemplo, o produto fica inapropriado para o consumo quando atingir mais do que uma determinada concentração de microorganismos por $mm^2$ de área do produto.

O evento de interesse (falha) pode ainda ocorrer devido a uma única causa ou devido a duas ou mais.

Censura e Dados Truncados

As dificuldades específicas relacionadas à análise de sobrevivência decorrem, em grande parte, do fato de que apenas alguns indivíduos vivenciaram o evento e, consequentemente, os tempos de sobrevivência serão desconhecidos para um subconjunto do grupo de estudo.

A censura ocorre quando no final do seguimento, alguns dos indivíduos não tiveram o evento de interesse, e assim o seu verdadeiro tempo para o evento é desconhecido.

Censuraé a observação parcial da resposta que foiinterrompida por alguma razão, não permitindo a observação completa do tempo de falha.

Esses tempos de sobrevivência censurados subestimam o tempo real (mas desconhecido) até o evento.

Um outro fenómeno que é geralmente confundido com a censura é o truncamento , embora seja um conceito totalmente diferente. Truncamento é caracterizado por uma condição que exclui certos indivíduos do estudo. Nestes estudos, os pacientes não são acompanhados a partir do tempo inicial, mas somente após experimentarem um certo evento.

Por exemplo, isto acontece se, para estimação da distribuição do tempo de vida dos moradores de uma certa localidade, for usada uma amostra retirada do banco de dados da previdência local.

Desta forma, somente moradores que atingiram a aposentadoria fazem parte da amostra. Estas observações são conhecidas por truncadas à esquerda.

No estudo da taxa de sobrevivência a determinada doença em determinada localidade, se eu recorrer aos dados do recenseamento eleitoral na junta de freguesia, apenas considero os sujeitos maiores de idade. Dessa forma, existe uma truncatura à esquerda, uma vez que os menores de idade são excluídos do estudo.

Em estudos de AIDS, a data da infecção é uma origem de tempo bastante utilizada e o evento de interesse pode ser o desenvolvimento da AIDS. Neste caso, o número de pacientes infectados é desconhecido. Então, indivíduos já infectados e que ainda não desenvolveram a doença são desconhecidos para o pesquisador e não são incluídos na amostra. Neste caso, somente pacientes que têm comprovada a doença fazem parte da amostra. Estas observações são chamadas de truncadas à direita.

As Observações Censuradas

São consideradas observações censuradas os casos em que:

um paciente ainda não vivenciou o desfecho relevante, como recidiva ou óbito, até o encerramento do estudo, ou seja, o evento ainda não ocorreu;
O evento ocorreu mas de causas não relacionadas (por exemplo, a pessoa morreu mas não da doença que estava a ser estudada).
um paciente é perdido para o acompanhamento durante o período do estudo, ou seja, foram feitas outras intervenções (ou seja, foram introduzidos outros factores); ou
um paciente vivencia um evento diferente que impossibilita o acompanhamento posterior.

Os dados censurados geralmente são representados por uma estrela, asterisco ou símbolo ‘mais’.

Todas as observações da amostra devem ser classificadas em censuradas ou não censuradas. Os casos censurados devem ser codificados com 0 e os não censurados com 1. Por exemplo, num caso de um estudo médico, o 1 seria usado para codificar os casos de pessoas que tinham morrido efectivamente da doença e o 0 os restantes casos.

A percentagem de casos censurados geralmente é um indicador de qualidade destes testes. Idealmente, a análise de sobrevivência deve ser feita quando se garante que cerca de metade dos dados são censurados e a outra metade atingiram o evento.

Ressalta-se o fato de que, mesmo censurados, todos os resultados provenientes de um estudo de obrevivência devem ser usados na análise estatística. Duas razões justificam tal procedimento:

mesmo sendo incompletas, as observações censuradas fornecem informações sobre o tempo de vida de pacientes;
omissão das censuras no cálculo das estatísticas de interesse pode acarretar conclusões viciadas.

Tipos das Censuras

Visualizando o processo de sobrevivência de um indivíduo como uma linha do tempo, seu evento (caso ocorra) está além do final do período de acompanhamento. Essa situação é frequentemente chamada de censura à direita (tempo de falha é superior ao tempo registrado)

A censura também pode ocorrer se observarmos a presença de um estado ou condição, mas não soubermos onde ele começou. Por exemplo, considere um estudo que investiga o tempo até a recorrência de um câncer após a remoção cirúrgica do tumor primário. Se os pacientes fossem examinados 3 meses após a cirurgia para determinar a recorrência, aqueles que apresentassem recorrência teriam um tempo de sobrevida censurado à esquerda, pois o momento real da recorrência ocorreu menos de 3 meses após a cirurgia, ou seja, ocorre quando o tempo registrado é maior do que o tempo de falha.

Os dados de tempo de evento também podem ser censurados por intervalo, o que significa que os indivíduos entram e saem do período de observação. Se considerarmos o exemplo anterior e os pacientes também forem examinados aos 6 meses, aqueles que estiverem livres da doença aos 3 meses e forem perdidos para o acompanhamento entre 3 e 6 meses serão considerados censurados por intervalo.

Mecanismo de Censuras

Alguns mecanismos de censura são diferenciados em estudos clínicos. Censura do tipo I é aquela em que o estudo será terminado após um período pré-estabelecido de tempo. As observações cujo evento de interesse não foi observado até este tempo são ditas censuradas.

Censura do tipo II é aquela em que o estudo será terminado após ter ocorrido o evento de interesse em um número pré-estabelecido de indivíduos/observações.

Um terceiro mecanismo de censura, o do tipo aleatório, é o que mais ocorre na prática médica. Isto acontece quando um paciente é retirado no decorrer do estudo sem ter ocorrido a falha, ou também, por exemplo, se o paciente morrer por uma razão diferente da estudada.

A intervalar é um tipo mais geral de censura que acontece, por exemplo, em estudos em que os pacientes são acompanhados em visitas periódicas e é conhecido somente que o evento de interesse ocorreu em um certo intervalo de tempo. Pelo fato de o tempo de falha $T$ não ser conhecido exatamente, mas sim pertencer a um intervalo, isto é, $T \in (L_{inf}, L_{sup}]$, tais dados são denominados por sobrevivência intervalar ou, mais usualmente, por dados de censura intervalar.

Figura 1.a: Ilustração de Mecanismo de Censura

(a) todos os pacientes experimentaram o evento antes do final do estudo;
(b) alguns pacientes não experimentaram o evento até o final do estudo;

Figura 1.b: Ilustração de Mecanismo de Censura

(c) o estudo foi finalizado após a ocorrência de um número pré-estabelecido de falhas; e
(d) o acompanhamento de alguns pacientes foi interrompido por alguma razão e alguns pacientes não experimentaram o evento até o final do estudo.

PROBLEMA AO TRABALHAR COM DADOS DE SOBREVIVÊNCIA:

Nem todos os indivíduos da amostra são observados até os seus respectivos tempos de falha.
Em situações onde o tempo entre o início da observação e a falha é muito longo, os dados são frequentemente analisados antes que o evento de interesse tenha ocorrido em todos os pacientes.
Pode também haver desistências ou perdas de seguimento.

Representação dos Dados

Os dados de sobrevivência para o indivíduo i $(i = 1,...,n)$ sob estudo são representados, em geral, pelo par $(t_i, \delta_i)$ sendo $t_i$ o tempo de falha ou de censura e $\delta_i$ a variável indicadora de falha ou censura, isto é,

\[ \delta_i = \begin{cases} 1, & \text{se } t_i \text{ é um tempo de falha}, \\ 0, & \text{se } t_i \text{ é um tempo censurado}. \end{cases} \]

Desta forma, a variável aleatória resposta em análise de sobrevivência é representada por duas colunas no banco de dados.

Na presença de covariáveis medidas no $i-ésimo$ indivíduo, tais como, dentre outras, x; = (sexo;, idade;, tratamento;), os dados ficam representados por (ti, i,x;).

No caso particular de dados de sobrevivência intervalar, tem-se, ainda, a representação $(\ell_i, u_i, \delta_i, \text{x}_i)$ em que $\ell_i$ e $u_i$ são, respectivamente, os limites inferior e superior do intervalo observado para o $i-ésimo$ indivíduo.

Covariáveis

Na análise de sobrevivência, a semelhanºa de outras análises, é importante recolher, para cada sujeito, os valores de outras covariáveis de interesse, que se espera influenciar o tempo de sobrevivência.

Ao incluir covariáveis, obtemos diferenças na função de sobrevivência consoante o valor da covariável e assim podemos testar se existem diferenças significativas entre valores diferente e grupos de interesse.

No contexto de estudos de sobrevivência a doenças, salientamos as seguintes covariáveis, que têm geralmente um efeito significativo e que portanto devem ser consideradas (Serranho, P., 2015):

Coorte: Em estudos de sobrevivência, a função de sobrevivência é fortemente afetada pela época temporal em que é feito o estudo. É evidente que a taxa de sobrevivência ao cancro tem vindo a aumentar ao longo dos anos, devido aos avanços da medicina.

Assim, caso tenhamos dados transversais a várias épocas, é importante introduzir no modelo uma variável coorte. Por exemplo, se tivermos dados de sobrevivência a cancro recolhidos em períodos de observação entre 1981 e 1984, depois entre 1990 e 1992 e entre 1996 e 1999, será indicado considerar no modelo uma variável coorte com valores 0, 1 e 2, consoante o período em que o sujeito foi observado, uma vez que é esperado que a função de sobrevivência seja diferente consoante o coorte considerado. Os Coortes podem ser décadas, séculos ou períodos definidos no tempo e permitem ao modelo dar estimativas diferentes consoante o coorte considerado.

Grupo de Controlo vs. Grupo de Estudo: A análise de sobrevivência tem aplicações para mostrar que determinado(s) tratamento(s) inovador(es) tem melhores resultados que o tratamento convencional. Assim, é importante recolher dados de sujeitos em vários grupos diferentes, como o grupo de controlo (sujeitos sem tratamento tomando um placebo, ou com o tratamento convencional) e o grupo de estudo (sujeito ao tratamento A, B e/ou C). Assim, a variável Grupo deve ser considerada no modelo com valores 0 para o grupo de controlo e 1,2,3, … para os grupos sujeitos ao tratamento A, B, C, etc …

Outras covariáveis: Consoante a doença e o objetivo do estudo importa recolher dados de outros fatores que possam influenciar o tempo de sobrevivência.

Por exemplo, em geral, a idade do paciente na data de diagnóstico é importante para o modelo. A idade do sujeito influencia a progressão da maioria das doenças. O género do sujeito é uma variável de interesse na maioria dos casos, uma vez que influencia não só a taxa de incidência, mas também a progressão da doença. Em estudos de cancro, os estadio do cancro no momento do diagnóstico influencia fortemente a função de sobrevivência, pelo que o estadio deve ser considerado no modelo.

Tempo de Sobrevivência

A variável aleatória não-negativa $T$, usualmente contínua, que representa o tempo até a ocorrência de um evento (falha, morte, etc.), é geralmente especificada em análise de sobrevivência pela sua função de sobrevivência ou pela função de taxa de falha (ou risco).

O evento de interesse é geralmente denominado neste contexto por falha, sendo que a variável de interesse é portanto geralmente denominada por tempo de sobrevivência que, pode ser dada em minutos, horas, dias, meses ou anos.

Tabela de sobrevivência

A tabela de sobrevivência (tabela de mortalidade ou tabela actuarial) é um dos métodos mais antigos (remontando os estudos de demografia e mortalidade apresentados à cidade de Londres por J. Graunt na segunda metade do séc. XVII) apresentando dados de sobrevivência de uma amostra de sujeitos seguidos ao longo do tempo (coorte).

Este tipo de tabela apresenta, para cada idade ou tempo, a probabilidade de uma pessoa daquela idade, ou naquele tempo, experienciar o evento de interesse (e.g. morte; conversão; …) antes da próxima avaliação temporal.

A partir destes dados é possível estimar a probabilidade condicional de um indivíduo sobreviver a um determinado tempo ou idade, dada a probabilidade de sobrevivência no início do estudo, bem como estimar a probabilidade de sobrevivência ou esperança de vida para diferentes tempos ou idades.

Função de Sobrevivência

Esta é uma das principais funções probabilísticas usadas para descrever estudos de sobrevivência. A função de sobrevivência é definida como a probabilidade de uma observação não falhar até um certo tempo t, ou seja, a probabilidade de uma observação sobreviver ao tempo t. Em termos probabilísticos é definida como:

\[ S(t)=P(T \geq t) \]

Em consequência, a função de distribuição acumulada é definida como a probabilidade de uma observação não sobreviver ao tempo t, isto é,

\[ F(t)=P(T \leq t)=1-S(t) \]

em que $T$ denota a variável aleatória contínua e não-negativa que representa o tempo decorrido até o evento de interesse; $𝑆(𝑡)$ é uma função monótona não crescente. A sua representação gráfica é chamada de curva de sobrevivência.

Na Figura 2 pode ser observada a forma típica de duas funções de sobrevivência. Estas curvas que, nesse caso, representam as funções de sobrevivência de dois grupos de pacientes, o grupo 1, tratado com a droga A e o grupo 2, com a droga B, fornecem informações importantes.

Note, por exemplo, que o tempo de vida dos pacientes do grupo 1 é superior ao dos pacientes do grupo 2 na maior parte do tempo de acompanhamento. Para os pacientes do grupo 1, o tempo para o qual cerca de 50% (tempo mediano) deles morrem é de 20 anos, enquanto que, para os pacientes do grupo 2, este tempo é menor (10 anos).

Outra informação importante e possível de ser retirada desta figura é o percentual de pacientes que ainda estão vivos até um determinado tempo de interesse. Por exemplo, para os pacientes do grupo 1, cerca de 90% deles ainda estão vivos após 10 anos do início do estudo, enquanto que, para os do grupo 2, apenas 50%.

Figura 2: Funções de sobrevivência para dois grupos de pacientes.

Da mesma forma, no caso de funções suaves e bem comportadas, como a função densidade de probabilidade f satisfaz por definição

\[ F(t) = \int_{-\infty}^{t} f(u) \, du \]

que no caso da variável tempo $T$ (por assumir apenas valores positivos) se reduz a

\[ F(t) = \int_{0}^{t} f(u) \, du \]

temos, diferenciando em t de ambos os lados, que

\[ F'(t) = f(t)=-S'(t) \]

Taxa de Falha - Hazard rate

A taxa de variação da falha num intervalo de tempo é a razão entre a probabilidade de uma falha ocorrer nesse intervalo de tempo e a amplitude do intervalo, condicionada ao facto da falha não ter acontecido antes.

Considerando um instante $t$, a taxa de falha $\lambda(t)$ nesse instante é obtida fazendo o limite do intervalo anterior em torno de $t$ para uma amplitude nula, ou seja, é definida por

\[ \lambda(t) = \lim_{h \to 0} \frac{P(t \le T \le t + h \mid T \ge t)}{h} \]

Esta quantidade mede a probabilidade instantânea de ocorrência do evento em $t$, dado que o indivíduo sobreviveu até esse instante.

Pela definição de probabilidade condicional, temos:

\[ P(t \le T \le t + h \mid T \ge t) = \frac{P(t \le T \le t + h)}{P(T \ge t)} \]

Utilizando a função de distribuição acumulada, temos:

\[ \small F(t + h) - F(t) = (1 - S(t + h)) - (1 - S(t)) \]

Substituindo no limite anterior:

\[ \lambda(t) = \lim_{h \to 0} \frac{S(t) - S(t + h)}{h\, S(t)} \]

Assumindo que $S(t)$ é diferenciável e, como a densidade de $T$ é dada por $f(t) = -S'(t)$, obtemos a forma clássica:

\[ \boxed{ \lambda(t) = \frac{f(t)}{S(t)} } \]

A Figura seguinte mostra três funções de taxa de falha. A função crescente indica que a taxa de falha do paciente aumenta com o transcorrer do tempo. Este comportamento mostra um efeito gradual de envelhecimento. A função constante indica que a taxa de falha não se altera com o passar do tempo. A função decrescente mostra que a taxa de falha diminui à medida que o tempo passa.

Figura 3: Funções de Taxa de Falha

A taxa de falha para o tempo de vida de seres humanos é uma combinação das curvas apresentadas em diferentes períodos de tempo. Ela é conhecida como curva da banheira e tem uma taxa de falha decrescente no período inicial, representando a mortalidade infantil, constante na faixa intermediária e crescente na porção final.

Uma representação desta curva é mostrada na Figura seguinte. A função de taxa de falha é mais informativa do que a função de sobrevivência.

Figura 4: Representação da função de taxa de falha conhecida por Curva da Banheira.

Diferentes funções de sobrevivência podem ter formas semelhantes, enquanto as respectivas funções de taxa de falha podem diferir drasticamente. Desta forma, a modelagem da função de taxa de falha é um importante método para dados de sobrevivência.

Taxa de Falha Acumulada

A função de taxa de falha acumulada, $\Lambda(t)$, não tem uma interpretação direta, mas pode ser útil na avaliação da função de maior interesse que é a de taxa de falha, $\lambda(t)$.

Isto acontece essencialmente na estimação não-paramétrica em que $\Lambda(t)$ apresenta um estimador com propriedades ótimas e $\lambda(t)$ é difícil de ser estimada.

\[ \Lambda(t) =\int_{0}^{t} \lambda(u) \, du =\int_{0}^{t} \frac{-S'(t)} {S(t)}\, du =-[\text{ln} \, S(u)]_{0}^{t} \]

e como $S(0) = 1$ (ou seja, no instante inicial todos os sujeitos estão vivos), temos $\text{ln} \, S(0) = 0$, logo

\[ \Lambda(t) =\int_{0}^{t} \lambda(u) \, du =- \text{ln} \, S(t) \]

Esta função reflecte assim o risco acumulado que um sujeito tem de experimentar o evento no tempo $t$, considerando que sobreviveu até esse tempo (ou simplesmente a soma dos riscos de um sujeito entre o tempo $0$ e o tempo $t$ no caso da integral acima se referindo a tempos discretos).

Técnicas Não-Paramétricas

Os objetivos de uma análise estatística envolvendo dados de sobrevivência estão geralmente relacionados, em medicina, à identificação de fatores de prognóstico para uma certa doença ou à comparação de tratamentos em um estudo clínico enquanto controlado por outros fatores. Vários exemplos podem ser encontrados na literatura médica.

No estudo de leucemia pediátrica (Colosimo, E. A & Giolo, S. R (2024)), leucometria registrada no diagnóstico (contagem de células brancas) e idade são conhecidos fatores de prognóstico para o tempo de vida de crianças com leucemia.

O estudo e análise de dados de sobrevivência frequentemente envolvem situações em que algumas observações são censuradas, ou seja, o evento de interesse (como uma falha ou morte) não ocorre durante o período de observação, mas o tempo até a falha é conhecido apenas parcialmente. A presença de censura em dados é um desafio para as técnicas tradicionais de análise descritiva, como médias, desvios padrão e gráficos (como histogramas e box-plots), que assumem a disponibilidade de informações completas.

Um exemplo típico ocorre ao tentar construir um histograma para ilustrar o tempo até a falha de um evento. Quando não há censura nos dados, o histograma é construído dividindo o eixo do tempo em intervalos e contando o número de falhas que ocorrem em cada intervalo. No entanto, quando há observações censuradas, a frequência exata de falhas em cada intervalo se torna desconhecida, tornando impossível construir um histograma tradicional.

Mesmo assim, algumas técnicas podem ser usadas com precaução. Um exemplo é o gráfico de dispersão, onde se analisa a relação entre uma covariável contínua (como a leucometria) e o tempo de sobrevivência. Embora a censura dificulte a interpretação do gráfico, ele ainda pode oferecer informações valiosas sobre a relação entre as variáveis.

O estudo de um conjunto de dados de leucemia pediátrica ilustra esse ponto. Nesse estudo, o tempo de sobrevivência (medido pelo tempo entre remissão e recidiva) é correlacionado com a raiz quadrada da leucometria ao diagnóstico (número de células brancas no sangue). A transformação de raiz quadrada é comum em variáveis com ampla escala de medida, como a leucometria. O gráfico resultante mostra uma associação negativa entre a leucometria e o tempo de sobrevivência, isto é, pacientes com maior contagem de células brancas tendem a ter tempos de sobrevivência mais curtos. Além disso, o gráfico revela que, para os pacientes com contagem de leucometria mais baixa, há uma maior proporção de observações censuradas, o que indica que esses pacientes não falharam dentro do período de estudo.

Portanto, apesar dos desafios impostos pelas observações censuradas, técnicas como gráficos de dispersão podem ainda ser úteis para visualizar e entender relações entre variáveis, desde que sejam interpretadas com cuidado e levando em consideração a natureza da censura.

Figura 2.1: Gráfico de dispersão do tempo de sobrevivência versus araiz quadrada da contagem inicial de leucócitos para os dados de leucemia pediátrica.

Nos textos básicos de estatística, uma análise descritiva consiste essencialmente em encontrar medidas de tendência central e variabilidade. Como a presença de censuras invalida este tipo de tratamento aos dados de sobrevivência, o principal componente da análise descritiva envolvendo dados de tempo de vila é a função de sobrevivência.

Nesta situação, o procedimento inicial para dados com censura é encontrar a estimativa para a função de sobrevivência e a partir dela estimar as estatísticas de interesse. Para realizar essa estimativa podem ser considerados três estimadores: Kaplan‐Meier, Nelson Aalen e Tabela de vida.

Os estimadores não paramétricos não assumem qualquer modelo de distribuição para a associação entre as covariáveis de interesse e a taxa de falha acumulada. Assim, estes são aplicáveis na maioria dos casos em que não existe informação para o modelo próprio da função de sobrevivência.

Estimação na Ausência de Censura

Nesta seção é tratada a estimação das funções de sobrevivência e de taxa de falha em uma situação sem censura. A função de taxa de falha é difícil de ser estimada em termos não-paramétricos. A dificuldade é a mesma de se estimar a função de densidade.

Alguns textos apresentam uma estimativa para esta função como sendo a variação da função de taxa de falha acumulada. No entanto, esta estimativa não é boa, principalmente para amostras de tamanho pequeno.

A Figura seguinte apresenta um histograma que mostra a distribuição dos tempos de falha associados a um certo grupo de indivíduos. Este histograma foi obtido a partir de uma amostra de 54 observações não-censuradas.

Figura 5: Histograma dos tempos de falha de um certo grupo de indivíduos.

Uma estimativa para a taxa de falha no período compreendido entre 400 e 500 horas é dada por:

\[ \hat{\lambda}([400,500)) = \frac{\text{nº de falhas no período } [400,500)} {\text{nº de indivíduos que não falharam até } t = 400} \]

\[ \hat{\lambda}([400,500)) = \frac{9}{21} = 0{,}429 \]

Em palavras, a taxa de falha é de 42,9% durante o período de 100 horas, compreendido entre 400 e 500 horas a partir do início do estudo.

Isto significa que, entre 100 indivíduos que sobreviveram até 400 horas, espera-se que 57 sobrevivam mais 100 horas.

A probabilidade de sobrevivência no tempo t = 400 horas é, por sua vez, estimada por:

\[ \small \hat{S}(400) = \frac{\text{nº de indivíduos que não falharam até o t = 400}} {\text{nº de indivíduos em estudo}} \]

\[ \hat{S}(400) = \frac{21}{54} = 0{,}389 \]

Em palavras, este número significa que 39% destes indivíduos sobrevivem mais do que 400 horas.

Usando o mesmo tipo de cálculo para os outros intervalos de tempo, obtêm-se os resultados mostrados na Tabela seguinte:

Tabela 1: Estimativas das taxas de falha e das probabilidades de sobrevivência para os dados do histograma mostrado na Figura 5

\[ \begin{array}{c c c} \text{Intervalo} & \text{Taxa de Falha (%/hora)} & \text{Sobrevivência (%)} \\ 0-100 & 0,037 & 100,0 \\ 100-200 & 0,096 & 96,3 \\ 200-300 & 0,213 & 87,0 \\ 300-400 & 0,432 & 68,5 \\ 400-500 & 0,429 & 38,9 \\ 500-600 & 0,583 & 22,2 \\ 600-700 & 0,800 & 9,3 \\ 700-800 & 1,000 & 1,9 \\ \end{array} \]

A forma utilizada para calcular as estimativas para as taxas de falha foi bastante intuitiva. Uma forma alternativa consiste em apresentar a taxa de falha em termos da função de sobrevivência.

\[ \hat{\lambda}([400,500)) = \frac{\hat{S}(400)-\hat{S}(500)}{(500-400)\hat{S}(400)}=0{,}0043/\text{hora} \]

Estimador de Kaplan-Meier

O estimador não-paramétrico de Kaplan-Meier é a estimativa de máxima verosimilhança da função de sobrevivência, proposto por Kaplan e Meier (1958) para estimar a função de sobrevivência, é também chamado de estimador limite-produto. É uma função em escada em que o início dos degraus da escada coincídem com os instantes onde ocorrem falhas.

Assim, o primeiro passo para a estimativa de Kaplan-Meier é dividir o intervalo em n sub-intervalos, com limites nos pontos de falha. A partir daí, começamos por estimar a função de sobrevivência no primeiro instante de falha $t_1$ como a razão entre o número de sujeitos em risco no instante $t_2$ sobre o número se sujeitos inicialmente no estudo.

No segundo instante de falha $t_2$ consideramos a relação

\[ \small P(T \geq t_2) = P(T \geq t_1, T \geq t_2) = P(T \geq t_1) \cdot P(T \geq t_2 \mid t \geq 1) \]

ou seja, a probabilidade de um sujeito sobreviver até $t_2$ é igual à probabilidade de sobreviver até $t_1$ vezes a probabilidade de sobreviver até $t_2$ sabendo que sobreviveu até $t_1$.

Este processo é depois repetido sucessivamente para se ter as estimativas sucessivas nos vários instantes de falha.

`Exemplo 1:`

Consideramos os dados recolhidos e ordenados de tempos de falha seguintes, em que o sinal + indica uma censura em vez de uma falha:

\[ 1, 3+, 4,5,5+,6+,7+,8,8,9+,10,11+,13+ \]

Começamos por estabelecer os extremos dos intervalos $t_i$ a considerar para a estimativa de Kaplan-Meier, ou seja,

\[ \small t_0 \, \, \text{(instante inicial)}, \, \, t_1=1, \, \, t_2=4, \, \, t_3=5, \, \, t_4=8, \, \, t_5=10 \]

Em seguida, para cada intervalo $[t_i, t_i+1[$ é preciso definir o número $n_i$ de sujeitos em risco (ou seja, vivos) e o número $d_i$ de falhas no intervalo.

A probabilidade de sobrevivência inicial é 1, logo consideramos que \[ S(0^+)=P(T\geq0)=1 \] ou seja, a probabilidade de sobreviver ao instante 0 é 1.

De seguida, temos que no instante $t = 1$ ocorre 1 falha em 13 sujeitos em risco. Assim, a estimativa para o sujeito não sobreviver ao instante $t = 1$ é $1/13$ (número de falhas sobre número de sujeitos em risco), pelo que a estimativa para sobreviver ao instante 1 é dada por

\[ \hat{S}(1^+)=1-\frac{1}{13}=0{,}92308 \]

No instante seguinte $t = 4$, ocorre 1 falha em 11 sujeitos em risco. A estimativa para o sujeito não sobreviver até ao instante $t = 4$, condicionada ao sujeito estar vivo no instante anterior $t = 1$ é então dada por:

\[ P(T \leq 4 | T \geq 1)=\frac{1}{11}=0{,}090909 \]

pelo que

\[ P(T \geq 4 \mid T \geq 1)=1-\frac{1}{11}=0{,}90909 \]

Assim, a estimativa da função de sobrevivência é dada por

\[ \small \hat{S}(4^+) = \underbrace{P(T \ge 1)}_{\hat{S}(1^+)} \cdot P(T \ge 4 \mid T \ge 1) = 0.92308 \times 0.90909 \]

De forma semelhante, a estimativa para a probabilidade de sobrevivência ao instante $t = 5$, condicionada ao sujeito estar vivo no instante $t = 4$ é

\[ P(T \geq 5 \mid T \geq 4)=1-\frac{1}{10}=0{,}9 \]

pelo que a estimativa para o sujeito estar vivo no instante t = 5 é dada por

\[ \small \hat{S}(5^+) = \hat{S}(4^+) \cdot P(T \ge 5 \mid T \ge 4) = 0.83916 \times 0.9 = 0.75524 \]

Repetindo o processo, obtem-se a estimativa de Kaplan-Meier para a função de sobrevivência para $\hat{S}(8)$ e $\hat{S}(10)$.

t_1 = c(1,3,4,5,5,6,7,8,8,9,10,11,13);    cens_1 = c(1,0,1,1,0,0,0,1,1,0,1,0,0)

ekm_1 = survfit(Surv(t_1, cens_1)~1);     summary(ekm_1)  # `Surv(t_1, cens_1)´ = Temp. Sobrev.

## Call: survfit(formula = Surv(t_1, cens_1) ~ 1)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1     13       1    0.923  0.0739        0.789        1.000
##     4     11       1    0.839  0.1045        0.657        1.000
##     5     10       1    0.755  0.1232        0.549        1.000
##     8      6       2    0.503  0.1669        0.263        0.964
##    10      3       1    0.336  0.1765        0.120        0.941

Finalmente, com o comando seguinte, obtemos o gráfico da função de sobrevivência estimada pelo estimador de Kapla-Meier:

plot(ekm_1, mark.time = TRUE, cex = 1.2,  xlab = "Tempo", ylab = "Sobrevivência" )

De notar que no gráfico da estimativa da função de sobrevivência, os instantes em que ocorrem censuras são assinalados automaticamente com o símbolo +.

`Exemplo 2:`

Observamos que o procedimento para se obter a estimativa de Kaplan-Meier envolve uma seqüência de passos, em que o próximo depende do anterior. Isto significa, por exemplo, que:

\[ \small \hat{S}(5^+) = P(T \ge 5) = P(T \ge 5, \, \, \,T \ge 1) = P(T \ge 1) \cdot P(T \ge 5\mid T \ge 1) \]

Desta forma, para o indivíduo sobreviver por 5 semanas, ele vai precisar sobreviver, em um primeiro passo, à primeira semana e depois sobreviver à quinta semana, sabendo-se que ele sobreviveu à primeira. Os tempos de 1 e 5 semanas foram tomados por serem os dois primeiros tempos distintos de falha nos dados dod grupo esteróide.

Os passos são gerados a partir de intervalos definidos pela ordenação dos tempos de falha de forma que cada um deles começa em um tempo observado e termina no próximo tempo.

A Tabela 2.2 apresenta os tempos ordenados e mostra a existência de 6 intervalos, iniciando com $[0,1)$, até o sexto intervalo que é $[10, 16)$. O limite superior deste último intervalo é definido como sendo 16 por ser este o maior tempo de acompanhamento do estudo.

Tabela 2: Estimativas de Kaplan-Meier para o grupo esteróide.

\[ \begin{array}{ccccc} \hline t_j & Intervalos & d_j & n_j & \hat{S}(t_j^+) \\ \hline 0 & [0, 1) & 0 & 14 & 1.000 \\ 1 & [1, 5) & 3 & 14 & 0.786 \\ 5 & [5, 7) & 1 & 9 & 0.698 \\ 7 & [7, 8) & 1 & 8 & 0.611 \\ 8 & [8, 10) & 1 & 7 & 0.524 \\ 10 & [10, 16) & 1 & 6 & 0.437 \\ \hline \end{array} \]

Todos os indivíduos estavam vivos em $t = 0$ e se mantêm até a primeira morte que ocorre em $t = 1$ semana. Então, a estimativa de $S(t)$ deve ser 1 neste intervalo compreendido entre 0 e 1 semana. No valor correspondente a 1 semana, a estimativa deve cair devido a três mortes que ocorrem neste tempo.

No segundo intervalo, $[1,5)$, existem, então, 14 indivíduos que estavam vivos (sob risco) antes de $t = 1$ e 3 morrem.

Desta forma, a estimativa da probabilidade condicional de morte neste intervalo é $3/14$ e a probabilidade de sobreviver é $1- 3/14$. Isto pode ser escrito como:

\[ \small \hat{S}(1^+) = P(T \ge 0 ) \cdot P(T \ge 1 \mid T \ge 0) = 1 \times \frac{11}{14} = 0.786 \]

Assim, sucessivamente, para qualquer $t$, $S(t)$ pode ser escrito em termos de probabilidades condicionais.

Suponha que existam $n$ pacientes no estudo e $k≤n$ falhas distintas nos tempos $t_1 < t_2 <... < t_k$. Considerando $S(t)$ uma função discreta com probabilidade maior que zero somente nos tempos de falha $t_j$, $j = 1,..., k$, tem-se que:

\[ S(t_j) = (1-q_1) \cdot (1-q_2) \cdot \,... \, \cdot (1-q_j) \]

em que $q_j$ é a probabilidade de um indivíduo morrer no intervalo $[t_{j-1},tj)$ sabendo que ele não morreu até $t_{j-1}$ e considerando $t_0 = 0$. Ou seja, pode-se escrever q; como:

\[ q_j = P(T \in [t_{j-1}, \, t_j) \mid T \ge t_{j-1}) \]

Desta forma, a expressão geral de $S(t)$ é escrita em termos de probabilidades condicionais. O estimador de Kaplan-Meier se reduz, então, a estimar $q_j$ que é dado por:

\[ \hat{q}(t) = \frac{\text{Nº de falhas em } t_{j-1}} {\text{Nº de observação sob ricsco em } t_{j-1}} \]

A expressão geral do estimador de Kaplan-Meier pode ser apresentada após estas considerações preliminares. Formalmente, considere:

$t_1 <t_2<t_k$, os $k$ tempos distintos e ordenados de falha, ou seja, instantes em que foram observadas falhas na amostra;
$d_j$ o número de falhas em $t_j$, $j = 1,..., k$; e
$n_j$ o número de indivíduos sob risco em $t_j$, ou seja, os indivíduos que não falharam e não foram censurados até o instante imediatamente anterior a $t_j$.
Além disso, podemos definir $c_1$ como o número de censuras que ocorrem para estas observações com tempo de monitoramento também exatamente igual a $t$.
Por fim, $t_0$ corresponde ao menor tempo de monitoramento entre todos os monitoramentos realizados na amostra.

Então, o estimador de Kaplan-Meier para a função de sobrevivência é dado pela função em escada

\[ \hat S(t) = \prod_{j:t_j \le t} \left(\frac{n_j-d_j}{n_j}\right)= \prod_{j:t_j \le t} \left(1 - \frac{d_j}{n_j}\right) \]

Naturalmente, o estimador de Kaplan-Meier se reduz à função de sobrevivência empírica se não existirem censuras.

Condições e Pressupostos

O teste de Kaplan-Meier tem vários pressupostos que têm de ser assegurados (Laerd, 2015).

O estado do evento consiste em apenas dois estados mutuamente exclusivos: o evento e os censurados. O evento pode ser a morte, a operação, a reincidência da doença, a avaria da máquina, o encontrar emprego, etc.

Os eventos são também colectivamente exaustivos. Isso significa que todos os casos têm de ser classificados num destes estados. Por exemplo, se se estiver a fazer um estudo de 10 anos na área do cancro, é necessário classificar todos os casos em censurados ou morte.
O tempo para o evento ou a censura (também conhecida como tempo de sobrevivência) deve estar definido e deve ser medido com precisão. Isto significa por exemplo nos estudos clínicos, que a data da última consulta deve ser registada em termos de dia e ano assim como a data do óbito, no caso do sujeito ter falecido.
A censura deve ser independente do evento, ou seja, quando se regista um sujeito como censurado deve ser porque

a) o sujeito desistiu do estudo ou

b) na data de fim do estudo o evento não tinha ocorrido. Por exemplo, ao fim de 5 anos do estudo, o sujeito continuava vivo.
Não devem existir alterações estruturais no estudo, conhecidas como secular trends ou secular changes. Uma característica comum dos estudos em que se usa a análise de sobrevivência é que não só decorre um longo período entre o início e o fim da experiência mas os participantes também não iniciaram a experiência - ou o tratamento - ao mesmo tempo e muitas vezes até se demora bastante tempo a recrutar ou seleccionar os sujeitos que integrarão o estudo. Os participantes no estudo não terão sido diagnosticados com a doença ao mesmo tempo, pelo que a data de início do estudo será a data de diagnóstico mais antiga entre todos os participantes. Pode acontecer existirem factores que afectem a probabilidade de ocorrência do evento e isso causa desvios no estudo. Por exemplo, a taxa de morte pode ter-se reduzido em virtude da introdução de novos medicamentos ou da alteração do protocolo de tratamento, o que aumentaria a taxa de sobrevivência dos sujeitos que fossem adicionados ao estudo mais tarde (i.e. reduzindo o right-censoring, ou censura a juzante). Também pode ter sido introduzidas alterações por exemplo nos planos de diagnóstico que tenham provocado diagnósticos mais rápidos e mais precoses, o que pioria as taxas de sobrevivência dos participantes que se teriam adicionado ao estudo mais cedo, antes dessas alterações (i.e. reduzindo o left-censoring, ou censura a montante). Estes factores (novos tratamentos, melhor diagnóstico) são exemplos de secular trends que podem criar desvios nos resultados.
Deve existir uma quantidade e padrão semelhante de casos censurados por grupo. Isto porque o teste de Kaplan-Meier baseia-se na ideia de que o a censura é igual em todos os grupos testado, em termos de quantidade e padrão.

Quando este critério não é assegurado, podem ser retiradas conclusões falsas sobre as distribuições de sobrevivência. Para garantir o cumprimento deste pressuposto, deve-se calcular a percentagem de casos censorados por grupo, para determinar a quantidade de censura e produzir um gráfico de pontos (scatterplot) para verificar o padrão da censura.

Propriedades do Estimador

As principais propriedades do estimador de Kaplan-Meier são basicamente as seguintes:

é não-viciado pára amostras grandes,
é fracamente consistente,
converge assintoticamente para um processo gaussiano e
é estimador de máxima verossimilhança de $S(t)$.

A partir das estimativas apresentadas na Tabela 2, o mais prático é construir um gráfico, por meio do qual é possível responder a possíveis perguntas de interesse. Este gráfico é construído mantendo o valor de $S(t)$ constante entre os tempos de falha. A forma gráfica do estimador de Kaplan-Meier para os grupos esteróide e controle, é apresentada na Figura 2.3.

As estimativas para o grupo controle são de simples cálculo, pois neste grupo existe somente um tempo distinto de falha. Em ambos os grupos, $S(t)$ não atinge o valor zero. Isto sempre acontece quando o maior tempo observado na amostra for uma censura.

As estimativas para o grupo esteróide apresentadas na Tabela 2, bem como as estimativas para o grupo controle e suas respectivas representações gráficas apresentadas, podem ser obtidas no pacote estatístico R por meio dos comandos apresentados a seguir:

 t_2 = 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_2 = 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)
 
 gr_2 = c(rep(1,15),rep(2,14));  ekm_2 = survfit(Surv(t_2, cens_2)~gr_2); summary(ekm_2)

## Call: survfit(formula = Surv(t_2, cens_2) ~ gr_2)
## 
##                 gr_2=1 
##         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 
## 
##                 gr_2=2 
##  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.832

 plot(ekm_2, mark.time = TRUE, lty=c(2,1), xlab = "Tempo (semanas)", ylab = "S(t) estimada")
 
 legend(1, 0.3, lty = c(2,1), c("Controle","Esteróide"), lwd = 1, bty = "n")

Figura 6: Estimativas de Kaplan-Meier para os grupos controle e esteróide dos dados de hepatite

`Exemplo 3:`

Consideremos os dados apresentados na tabela seguinte, de um exemplo hipotético ou ilustrativo.

\[ \begin{array}{ccc} \text{Observação} & \it{Status} & \text{Tempo de Monitoramento} \\ \hline 1 & \text{Censura} & 2 \\ 2 & \text{Evento} & 3 \\ 3 & \text{Evento} & 1 \\ 4 & \text{Evento} & 4 \\ 5 & \text{Censura} & 2 \\ 6 & \text{Evento} & 3 \\ \hline \end{array} \]

Realizando os cálculos das probabilidades de sobrevivência para os tempos de monitoramento, obtem-se:

Tabela 4.a: Estimativas de Kaplan-Meier para o grupo esteróide.

\[ \begin{array}{p{5cm}|p{5cm}|p{5cm}|p{5cm}|p{5cm}} \hline t & n_t & d_t & c_t & \hat{S}(t)\\ \hline 1 & 6 & 1 & 0 & \hat{S}(1^+) = \left( \frac{6-1}{6} \right) = 0,833\\ 2 & 5 & 0 & 2 & \hat{S}(2^+) = \left( \frac{6-1}{6} \right) \cdot \left( \frac{5-0}{5} \right) = 0,833 \\ 3 & 3 & 2 & 0 & \hat{S}(3^+) = \left( \frac{6-1}{6} \right) \cdot \left( \frac{5-0}{5} \right) \cdot \left( \frac{3-2}{3} \right) = 0,277 \\ 4 & 1 & 1 & 0 & \hat{S}(4^+) = \left( \frac{6-1}{6} \right) \cdot \left( \frac{5-0}{5} \right) \cdot \left( \frac{3-2}{3} \right) \cdot \left( \frac{1-1}{1} \right) = 0,000\\ \hline \end{array} \]

É importante que não haja confusão entre o tempo de monitoramento e o instante eni que se dá o início do monitoramento de cada observação.

É o primeiro que nos interessa, já que o nosso intuito é calcular as probabilidades de sobrevivência ao evento para cada período de monitoramento, independentemente de quando se inicia.

Podemos inicialmente observar que o’s tempos de monitoramento foram dispostos de forma crescente, mesmo que isto não tenha sido verificado. Assim, podemos verificar que, para um tempo de monitoramento menor do que 1, nenhuma observação apresentou evento ou censura $(n_1 = 6)$, porém uma delas apresentou evento exatamente no tempo $t = 1$ $(e_1 = 1)$.

Já para um tempo de monitoramento menor do que 2, verificamos que cinco observações ainda não apresentaram evento ou censura $(n_2 = 5)$, porém duas delas apresentaram censura exatamente no tempo $t = 2$ $(c_2 = 2)$. Como não ocorreu nenhum evento no tempo de monitoramento $t = 2$ $(e_2 =0)$, o cálculo da probabilidade não sofre nenhuma alteração $( S (1) = S (2) =0,833 )$.

Por outro lado, as duas censuras que ocorreram em $t = 2$ fazem com que apenas três observações não tenham apresentado evento ou censura para um tempo de monitoramento menor do que 3 $(n_3 = 3)$ e, como mais duas apresentaram evento em $t = 3$ $(e_3 = 2)$, isso precisa ser levado em consideração para o cálculo da probabilidade de sobrevivência ao evento para um tempo de monitoramento maior do que $t = 3$ $( S ( 3) )$ .

Por fim, como apenas uma observação ainda não apresentou evento ou censura para um tempo de monitoramento menor do que 4 $(n_4 = 1)$, porém esta mesma observação sofre evento em $t = 4$ $(e_4 = 1)$, a probabilidade de sobrevivência ao evento para um tempo de monitoramento maior do que $t = 4$ é igual a zero $(S(4) = 0)$. Obviamente, a probabilidade de sobrevivência ao tempo máximo de monitoramento é sempre igual a zero $(S(t = \mathrm{máx}) =0)$ e a probabilidade de sobrevivência a um tempo nulo de monito ramento é sempre igual a 1 $(S(0)=1)$.

A quantidade de censuras que ocorrem para um determinado tempo de monitoramento t não interfere no cálculo da probabilidade de sobrevivência para o tempo de monitoramento maior do que t.

Entretanto, caso ocorram censuras em $t$, este fato influenciará no cálculo das probabilidades de sobrevivência ao evento para tempos de monitoramento maiores do que $t+1$.

Com base nos cálculos das probabilidades de sobrevivência ao evento para os diferentes tempos de monito ramento, podemos elaborar a curva da função de sobrevivência ao evento, também conhecida por curva de probabilidades de sobrevivência de Kaplan-Meier.

$Figura 7: Curva de probabilidades de sobrevivência de Kaplan-Meier \widehat S(t).$

Figura 7: Curva de probabilidades de sobrevivência de Kaplan-Meier $\widehat S(t)$.

As curvas de probabilidades de falha de Kaplan-Meier tipicamente também apresentam a forma de degraus, porém agora ascendentes, já que as probabilidades de ocorrência do evento para tempos de monitoramento maiores tendem a ser mais elevadas.

$Figura 7: Curva de probabilidades de falha de Kaplan-Meier \lambda(t)$

Figura 7: Curva de probabilidades de falha de Kaplan-Meier $\lambda(t)$

Realizando os cálculos das taxas de falha para os tempos de monitoramento, temos:

Tabela 4.b: Estimativas de Kaplan-Meier para o grupo esteróide.

\[ \begin{array}{c|c|c} \hline t & \hat{S}(t) & \lambda(t) & \Lambda(t) \\ \hline 1 & \hat{S}(1^+) = 0,833 & \lambda(1) = \frac{\hat{S}(0) - \hat{S}(1)}{(1) \cdot \hat{S}(0)} = \frac{1,000 - 0,833}{1} = 0,167 & 0,167 \\ 2 & \hat{S}(2^+) = 0,833 & \lambda(2) = \frac{\hat{S}(1) - \hat{S}(2)}{(1) \cdot \hat{S}(1)} = \frac{0,833 - 0,833}{0,833} = 0,000 & 0,167 \\ 3 & \hat{S}(3^+) = 0,277 & \lambda(3) = \frac{\hat{S}(2) - \hat{S}(3)}{(1) \cdot \hat{S}(2)} = \frac{0,833 - 0,277}{0,833} = 0,666 & 0,833 \\ 4 & \hat{S}(4^+) = 0,000 & \lambda(4) = \frac{\hat{S}(3) - \hat{S}(4)}{(1) \cdot \hat{S}(3)} = \frac{0,277 - 0,000}{0,277} = 1,000 & 1,833 \\ \hline \end{array} \]

Assim, a taxa de risco de ocorrência do evento para o tempo de monitoramento t = 1 é igual a 0,167, visto que apenas uma observação apresentou evento em t = 1 entre as seis que começaram a ser monitoradas $(t = 0)$.

Já para $t =2$, a taxa de falha é igual a 0,000, uma vez que, das cinco observações que foram monitoradas por um período de tempo maior do que 1, nenhuma apresentou evento em $t = 2$ (apenas censuras).

Para o tempo de monitoramento $t = 3$, a taxa de risco de ocorrência do evento é igual a $0,666$, já que duas observações apresentaram evento em $t = 3$ entre as três que foram monitoradas por um período maior do que 2.

Por fim, para o tempo de monitoramento $t = 4$, a taxa de falha é igual a 1,000, uma vez que apenas uma observação foi monitorada por um período de tempo maior do que três, tendo esta apresentado evento em $t = 4$.

Em outras palavras, o risco de ser evento para um período máximo de monitoramen!o é igual a 1 (100%). Além disso, a última coluna da Tabela, que apresenta os valores acumulados de $\Lambda(t)$ ao longo dos tempos de monitoramento.

Tempo Médio e Mediano

A utilização direta da curva de Kaplan-Meier nos informa a probabilidade estimada de sobrevivência para um determinado tempo. Um exemplo é a probabilidade do paciente sobreviver a 12 semanas de tratamento. A estimativa de Kaplan-Meier para este valor é diretamente obtida da Figura 7 e é igual a $44\%$.

Se o valor do tempo de interesse estiver ao longo de um degrau da curva de Kaplan-Meier, pode-se também utilizar uma interpolação linear. Por exemplo, como havia sido observado na Seção 7.2, a probabilidade estimada de um paciente do grupo esteróide sobreviver a 6 semanas obtida diretamente da curva de Kaplan-Meier é de $0,698$.

No entanto, se a interpolação linear for utilizada, obtém-se:

\[ \frac{7-5}{0,611-0,698} = \frac{6-5}{\widehat S(6)-0,698} \]

cuja solução é a estimativa de $0,655$. Esta última estimativa deve ser preferida (Colosimo et al., 2002). A partir da curva de Kaplan-Meier também é possível obter estimativas de percentis. Uma informação muito útil é o tempo mediano de vida, ou seja, instante de tempo em que a função de sobrevivência é igual a 50%. Como a curva de Kaplan-Meier é uma função escada, a estimativa mais adequada para a mediana é novamente obtida por meio de uma interpolação linear. Isto é,

\[ \frac{10-8}{0,437-0,524} = \frac{\text{MED}-8}{0,5-0,524} \]

cuja solução é a estimativa de $8,55$ semanas.

Esta forma de estimar estes valores é equivalente a conectar por retas as estimativas de Kaplan-Meier, em vez de se utilizar $S(t)$ na forma de escada. Esta forma usualmente gera uma melhor representação da distribuição contínua dos tempos de vida (Colosimo et al., 2002). De forma análoga, pode-se obter estimativas de outros percentis da distribuição dos tempos de vida dos pacientes.

Marôco, J (2021) sugere/propõe outra forma de outra forma de determinar os percentis $p \times 100$ dos tempos de sobrevivência como:

\[ t_p = \min \, \{t_j:S(t_j) \le p \} \text{ com } i = 1, \, ..., \,k \]

Relembrando que, por exemplo, a mediana é o percentil 50 $(i.e., \, \, p = 0,5)$

A variância assintótica do estimador de percentis $(\widehat t_p)$ é expressa por:

\[ \mathrm{Var}(\widehat t_p) = \frac{\mathrm{Var} \left( \widehat S(\widehat t_p) \right)}{\left[ f(\widehat t_p) \right]^2} \]

A dificuldade em se obter uma estimativa para $f(\widehat t_p)$ inviabiliza a utilização desta expressão. Brookmeyer e Crowley (1982) propõem um estimador alternativo para a mediana invertendo a região de rejeição de um teste não-paramétrico que não necessita estimar $f(\widehat t_p)$.

Outra quantidade que pode ser de interesse é o tempo médio de vida do paciente. Esta quantidade, no entanto, nem sempre é estimada adequadamente utilizando estimadores não-paramétricos em estudos incluindo censuras. Pode ser mostrado, por argumentos probabilísticos, que o tempo médio de vida é dado pela área (integral) sob a função de sobrevivência.

Uma estimativa para o tempo médio é então obtida calculando-se a área sob a curva de Kaplan-Meier estimada. Como esta curva é uma função escada, esta integral é simplesmente a soma de áreas de retângulos, isto é,

\[ \hat t_m = t_1 + \sum_{j=1}^{k-1} \widehat S(t_j)(t_{j+1}-t_j) \]

em que $t_1 <... < t_k$ são os k tempos distintos e ordenados de falha.

Entretanto, surge um problema se o maior tempo observado for uma censura. Isto acontece com freqüência em estudos clínicos, como é o caso dos dados de hepatite. Neste caso, a curva de Kaplan-Meier não atinge o valor zero e o valor do tempo médio de vida fica subestimado.

Nesses casos, tal estimativa deve ser interpretada com bastante cuidado ou talvez até mesmo evitada. Uma alternativa é utilizar a mediana em vez do tempo médio de vida. Ambas são medidas de tendência central, representando um valor típico da distribuição do tempo de vida da população sob estudo. A mediana, no entanto, pode ser extraída facilmente da função de sobrevivência que foi estimada anteriormente para os pacientes do grupo esteróide.

Nesta forma, utilizam-se os modelos paramétricos para dados de sobrevivência. Kaplan e Meier (1958) mostraram que a variância assintótica de $\widehat t_m$ pode ser estimada por:

\[ \widehat{\mathrm{Var}}(\widehat t_m) = \frac{r}{r-1} \left[ \sum_{j=1}^{r-1} \frac{\left( A_j \right)^2}{n_j \cdot (n_j-d_j)}\right] \]

\[ A_j = \widehat S(t_j)(t_{j+1}-t_j) \, + \, ... \, + \, \widehat S(t_{r-1})(t_r-t_{r-1}) \]

com $r$ é o número de observações não censuradas, isto é, o número de falhas. Observe que $r$ é igual ao número de falhas e não ao número de falhas distintas.

Outra quantidade possivelmente de interesse é o tempo médio restante de vida daqueles pacientes que se encontram livres do evento em um determinado tempo $t$, ou seja, o tempo de vida médio restante para os sujeitos com idade $t$, denotada por $\mathrm{vmr(t)}$. Este tempo é estimado pela área sob a curva de sobrevivência à direita de t dividido por $\widehat S(t)$, isto é,

\[ \widehat{\mathrm{vmr}}(t) = \frac{\text{área sob a curva }\widehat S(t) \text{ à direita de } t}{\widehat S(t)} \]

Este estimador apresenta as mesmas limitações de $\widehat t_m$.

`Exemplo 5`: Tumor Sólido

Para os dados desse exemplo (Reincidência de Tumor Sólido, as estimativas da função de sobrevivência $S(t)$ e seus respectivos intervalos de 95% de confiança, obtidos utilizando-se o estimador de Kaplan-Meier, encontram-se na Tabela 2.5. Para obtenção das estimativas pontuais e intervalares utilizou-se, no R, os comandos:

t_4 = c(3,4,5.7,6.5,6.5,8.4,10,10,12,15)   ;       cens_4 = c(1,0,0,1,1,0,1,0,1,1)

ekm_4 = survfit(Surv(t_4, cens_4)~1, conf.type = "plain")   ;      summary(ekm_4)

## Call: survfit(formula = Surv(t_4, cens_4) ~ 1, conf.type = "plain")
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   3.0     10       1    0.900  0.0949        0.714        1.000
##   6.5      7       2    0.643  0.1679        0.314        0.972
##  10.0      4       1    0.482  0.1877        0.114        0.850
##  12.0      2       1    0.241  0.1946        0.000        0.622
##  15.0      1       1    0.000     NaN          NaN          NaN

A representação gráfica de $\widehat S(t)$, com os respectivos intervalos de $95%$ de confiança para todo $t$ tal que $0 ≤t ≤ 15$ é mostrada na Figura 2.5 e foi obtida no R por meio dos comandos:

plot(ekm_4, mark.time = TRUE, conf.int = T,  xlab = "Tempo (em meses)", ylab = "S(t) estimada", bty = "n")

Figura 2.5: Sobrevivência e respectivos intervalos de 95% de confiança estimados a partir do estimador de Kaplan-Meier para os dados de tumor sólido.

Note, a partir da Tabela 2.5 e também da Figura 2.5, que os intervalos de confiança obtidos para $\widehat S(t)$ são relativamente amplos. Isto se deve, em particular, ao tamanho amostral relativamente pequeno $(n = 10)$, incluindo 4 censuras.

Para o tempo mediano, obtido por meio de uma interpolação linear, tem-se que:

\[ \frac{10-6,5}{0,482-0,643} = \frac{\text{MED}-6,5}{0,5-0,643} = 9,6\text{ meses} \]

Assim, $9,6$ meses é uma estimativa do tempo em que 50% dos pacientes permanecem vivos. Tem-se, ainda, para os pacientes deste exemplo, um tempo médio de vida estimado de $\widehat t_m = 10,1\text{ meses}$. Tal estimativa pode ser obtida no R por meio dos comandos:

t_41 = t_4[cens_4 == 1];       tj = c(0,as.numeric(levels(as.factor(t_41))))    

surv = c(1, as.numeric(levels(as.factor(ekm_4$surv))))    

surv_1 = sort(surv, decreasing = T);    k = length(tj)-1;    

prod = matrix(0,k,1)       
  for(j in 1:k){      
    prod[j] = (tj[j+1] - tj[j])*surv_1[j]   
    }    

tm_1 = sum(prod);  tm_1

## [1] 10.0875

Observe, neste exemplo, que o tempo médio apresenta-se bem estimado, uma vez que o maior tempo observado trata-se de uma falha. O mesmo não seria verdade, como discutido anteriormente, se o referido tempo correspondesse a uma censura.

A variância estimada de tm foi também obtida e esta resultou em:

\[ \widehat{\mathrm{Var}}(\widehat t_m) = \frac{6}{5} \left[ \frac{\left( A_1 \right)^2}{9 \cdot 10}+\frac{\left( A_2 \right)^2}{6 \cdot 7}+\frac{\left( A_3 \right)^2}{5 \cdot 6}+\frac{\left( A_4 \right)^2}{3 \cdot 4}+\frac{\left( A_5 \right)^2}{1 \cdot 2}\right] \]

sendo:

\[ A_5 = 0,24 \cdot (15-12)=0,723 \ \, \mid \ \, A_4 = 0,482 \cdot (12-10)+A_5=1,687 \]

\[ A_3 = A_2 = 0,643 \cdot (10-6,5)+ A_4 = 3,938 \ \, \mid \ \, A_1 = 0,9\cdot (6,5-3)+A_2 = 7,088 \]

Para pacientes que sobreviverem até, por exemplo, o tempo $t = 10 \text{ meses}$, estima-se, também, que os mesmos tenham um tempo médio de vida restante de:

\[ \widehat{\mathrm{vmr}}(t) = \frac{\text{área sob a curva }\widehat S(t) \text{ à direita de } t=10}{\widehat S(10)} = 3,5\text{ meses} \]

Intervalo de Confiança

Para que se possa construir intervalos de confiança e testar hipóteses para $S(t)$, é necessário, no entanto, avaliar a precisão do estimador de Kaplan-Meier. Este estimador, assim como outros, está sujeito a variações que devem ser descritas em termos de estimações intervalares.

A expressão para a variância assintótica do estimador de Kaplan-Meier é dada por: \[ \widehat {Var} (\hat S(t)) = [\hat S(t)]^2 \cdot \sum_{j:t_j < t} \left(\frac{d_j}{n_j\cdot(n_j-d_j}\right) \]

Esta expressão é conhecida como fórmula de Greenwood e pode ser obtida a partir de propriedades do estimador de máxima verossimilhança.

A estimativa da variância de S(5), para o exemplo considerado, é, então, dada por:

\[ \small \widehat {Var} (\hat S(5)) = (0,698)^2 \cdot \left[ \frac{3}{\left(14 \right)\left(11 \right)} + \frac{1}{\left(9 \right)\left(8 \right)} \right]^2 =0,0163 \]

Como $S(t)$, para $t$ fixo, tem distribuição assintótica Normal, segue que um intervalo aproximado de $100(1- \alpha)\%$ de confiança para $S(t)$ é dado por:\[ \hat S(t)) \pm z_{\alpha/2} \cdot \sqrt{\widehat {Var} (\hat S(t))} \]

O intervalo de 95% de confiança para $S(6)$ é $0,698 ±1,96/0,0163$, ou seja, $(0,45; 0,95)$.

Para obtenção deste intervalo no R deve-se usar:

ekm_21 = survfit(Surv(t_2, cens_2)~gr_2, conf.type = "plain", conf.int = 0.95) ; summary(ekm_21)

## Call: survfit(formula = Surv(t_2, cens_2) ~ gr_2, conf.type = "plain", 
##     conf.int = 0.95)
## 
##                 gr_2=1 
##         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 
## 
##                 gr_2=2 
##  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.718

Entretanto, para valores extremos de $t$, este intervalo de confiança pode apresentar limite inferior negativo ou limite superior maior do que 1.

Nesses casos, o problema é resolvido utilizando-se uma transformação para S(t) como, por exemplo, $\hat{U}(t) = \text{log}[- \text{log}(\hat{S}(t))]$, sugerida por Kalbfleish e Prentice (1980), que tem variância assintótica estimada por:

\[ \widehat{\mathrm{Var}}(\hat{U}(t)) = \frac{ \sum_{j : t_j < t} \frac{d_j}{n_j (n_j - d_j)} }{ \left[ \sum_{j : t_j < t} \log \left( \frac{n_j - d_j}{n_j} \right) \right]^2 } = \frac{ \sum_{j : t_j < t} \frac{d_j}{n_j (n_j - d_j)} }{ \left[ \log \,(\hat{S}(t)) \right]^2 } \]

Assim, um intervalo aproximado de $100(1-\alpha)\%$ de confiança para $S(t)$ é dado por:

\[ \widehat{S}(t) ^ {\exp \left\{ \pm z_{\alpha/2} \sqrt{\widehat{\mathrm{Var}}(\hat{U}(t))} \right\}} \]que assume valores no intervalo $[0,1]$ e resulta no intervalo $(0,38; 0,88)$ de 95% de confiança para $S(6)$ no exemplo dos dados de esteróide. Este intervalo é obtido no R por meio dos comandos:

ekm_22 = survfit(Surv(t_2, cens_2)~gr_2, conf.type="log-log");      summary(ekm_22)

## Call: survfit(formula = Surv(t_2, cens_2) ~ gr_2, conf.type = "log-log")
## 
##                 gr_2=1 
##         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 
## 
##                 gr_2=2 
##  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.683

Pode-se confirmar consultando o help do R para mais informações sobre os tipos de intervalos disponíveis. Os intervalos produzidos por default neste pacote usam a transformação Û(t) = log{S(t)] sendo, os mesmos, obtidos equivalentemente por uma das duas linhas de comandos apresentadas a seguir:

ekm_23 = survfit(Surv(t_2, cens_2)~gr_2, conf.type="log") ;      summary(ekm_23)

## Call: survfit(formula = Surv(t_2, cens_2) ~ gr_2, conf.type = "log")
## 
##                 gr_2=1 
##         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 
## 
##                 gr_2=2 
##  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.832

ekm_24 = survfit(Surv(t_2, cens_2)~gr_2);                        summary(ekm_24)

## Call: survfit(formula = Surv(t_2, cens_2) ~ gr_2)
## 
##                 gr_2=1 
##         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 
## 
##                 gr_2=2 
##  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.832

Outros Estimadores N-P

Como constataram, o estimador de Kaplan-Meier é, sem dúvida, o mais utilizado para se estimar $S(t)$ em análise de sobrevivência. Existem muitos pacotes estatísticos que disponibilizam este estimador e ele também é apresentado em vários textos de estatística básica.

Entretanto, outros dois estimadores de $S(t)$ têm importância na literatura mais especializada desta área. Eles são: o estimador de Nelson-Aalen e o estimador da tabela de vida.

O primeiro, de Nelson-Aalen, é mais recente que o de Kaplan-Meier e apresenta, aparentemente, propriedades similares ao deste último.

O segundo tem uma importância histórica, pois foi utilizado em informações provenientes de censos demográficos para, essencialmente, estimar características associadas ao tempo de vida dos seres humanos. Este estimador foi proposto por demógrafos e atuários no final do século XIX e utilizado basicamente em grandes amostras.

Estimador de Nelson-Aalen

Este estimador, como mencionado anteriormente, é mais recente do que de Kaplan-Meier e baseia-se na função de sobrevivência expressa por:

\[ S(t) = \exp \left\{ -\Lambda(t)\right\} \]

em que $\Lambda (t)$ é a função de risco acumulado.

Um estimador para $\Lambda (t)$ foi inicialmente proposto por Nelson (1972) e retomado por Aalen (1978), que provou suas propriedades assintóticas usando processos de contagem. Este estimador é denominado na literatura por Nelson-Aalen e tem a seguinte forma:

\[ \widetilde{\Lambda}(t) = \sum_{j : \, t_j < t} \left( \frac{d_j}{n_j} \right) \]

em que $d_j$ e $n_j$ são definidos como no estimador de Kaplan-Meier. A variância deste estimador, proposta por Aalen (1978), é dada por:

\[ \widehat{\mathrm{Var}}(\widetilde{\Lambda}(t) ) = \sum_{j : \, t_j < t} \left( \frac{d_j}{n_j^2} \right) \]

Um estimador alternativo para a variância de ${\widetilde \Lambda(t)}$ proposto por Klein (1991) é:

\[ \widehat{\mathrm{Var}}(\widetilde{\Lambda}(t) ) = \sum_{j : \, t_j < t} \left( \frac{(n_j-d_j)\cdot d_j}{n_j^3} \right) \]

mas, por apresentar menor vício, o estimador proposto por Aalen é preferível a este último.

Desse modo, e com base no estimador de Nelson-Aalen, um estimador para a função de sobrevivência é expresso por (conhecido em alguma bibliografia por Estimador de Fleming-Harrington):

\[ \widetilde{S}(t) = \exp \left\{ -\widetilde{\Lambda}(t) \right\} \]

A variância deste estimador, devido a Aalen e Johansen (1978), pode ser obtida por:

\[ \widehat{\mathrm{Var}}(\widetilde S(t)) = \left[ \widetilde S(t) \right]^2 \cdot \sum_{j: \, t_j < t} \left(\frac{d_j}{n_j^2}\right) \]

O estimador $\widetilde S(t)$ e o de Kaplan-Meier apresentam, na maioria das vezes, estimativas muito próximas para $S(t)$. Bohoris (1994) mostrou que $\widetilde S(t) \geq \widehat S(t)$ para todo $t$, ou seja, as estimativas obtidas por meio do estimador de Nelson-Aalen são maiores ou iguais às obtidas por meio do estimador de Kaplan-Meier.

A Tabela seguinte apresenta as estimativas de Nelson-Aalen para o grupo esteróide dos dados de hepatite. Observe que estas estimativas são bem próximas das de Kaplan-Meier mostradas na Tabela anterior, mesmo neste caso em que a amostra é relativamente pequena.

Tabela 5: Estimativas de Kaplan-Meier para o grupo esteróide.

\[ \begin{array}{ccccccc} t_j & n_j & d_j & \widetilde \Lambda (t_j) & \widetilde S (t_j^+) & \text{e. p } \left( \widetilde S (t_j^+) \right) & \text{I. C } \left( \widetilde S (t_j^+) \right)_{95\text{%}} \\ \hline 0 & 14 & 0 & 0 & 1 & --- & --- \\ 1 & 14 & 3 & 0,214 & 0,807 & 0,0999 & (0,633; \,\, 1,000) \\ 5 & 9 & 1 & 0,325 & 0,722 & 0,1201 & (0,521; \,\, 1,000) \\ 7 & 8 & 1 & 0,450 & 0,637 & 0,1326 & (0,424; \,\, 0,958) \\ 8 & 7 & 1 & 0,593 & 0,553 & 0,1394 & (0,337; \,\, 0,906) \\ 10 & 6 & 1 & 0,760 & 0,468 & 0,1414 & (0,259; \,\, 0,846) \\ \hline \end{array} \]

O estimador de Kaplan-Meier tem a vantagem de estar disponível em vários pacotes estatísticos, o que não acontece em geral com o de NelsonAalen.

No pacote estatístico R, por exemplo, as estimativas de Nelson Aalen para o grupo esteróide dos dados de hepatite apresentadas na Tabela acima podem ser obtidas por meio dos comandos:

 ena_1 = survfit(coxph(Surv(t_2[gr_2 == 2],cens_2[gr_2 == 2])~1, method = "breslow")) ; summary(ena_1)

## Call: survfit(formula = coxph(Surv(t_2[gr_2 == 2], cens_2[gr_2 == 2]) ~ 
##     1, method = "breslow"))
## 
##  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.846

Cabe aqui uma observação final sobre a função $\Lambda(t)$. Ela não tem interpretação probabilística mas tem utilidade na seleção de modelos. O gráfico da estimativa desta função, em papéis especiais, é utilizado para verificar a adequação de modelos paramétricos.

Se a partir da função de sobrevivência quisermos obter o estimador de Nelson Aalen para a taxa de falha acumulada, podemos fazê-lo aplicando $− \ln$ à função de sobrevivência dada pelo comando ena$surv.

 racum = -log(ena_1$surv);     cbind(ena_1$time, racum)

##             racum
## [1,]  1 0.2142857
## [2,]  4 0.2142857
## [3,]  5 0.3253968
## [4,]  7 0.4503968
## [5,]  8 0.5932540
## [6,] 10 0.7599206
## [7,] 12 0.7599206
## [8,] 16 0.7599206

`Exemplo 4.a): Hipotético`

Pergunta: Determine a estimativa para a função de sobrevivência com base no estimador de Nelson-Aalen.

Tal como nos exemplos anteriores, necessitamos de introduzir os tempos e as censuras, o que pode ser feito com os comandos:

t_3 = c(5,8,8,9,10,10,14,4,5,5,6,8,9,9,10,11,12,5,5,6,6,7,8,8,10)

cens_3 = c(1,1,0,1,1,0,1,1,1,0,1,1,1,0,0,1,0,0,1,1,0,1,1,0,0)

De seguida, introduzimos (pela ordem corresponde aos tempos e censuras) a variável grupos, por exemplo, como

gr_3 = c(rep(0,7),rep(1,10),rep(2,8))

O estimador de Nelson-Aalen pode ser obtido através dos comandos seguintes (ver figura seguinte):

ena_2 = survfit(Surv(t_3,cens_3)~gr_3,type="fh2",conf.type="plain", error="t"); summary(ena_2)

## Call: survfit(formula = Surv(t_3, cens_3) ~ gr_3, error = "t", type = "fh2", 
##     conf.type = "plain")
## 
##                 gr_3=0 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5      7       1    0.867   0.124       0.6242        1.000
##     8      6       1    0.734   0.161       0.4181        1.000
##     9      4       1    0.571   0.190       0.1988        0.944
##    10      3       1    0.409   0.193       0.0315        0.787
##    14      1       1    0.151   0.167       0.0000        0.477
## 
##                 gr_3=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     10       1    0.905  0.0905        0.727        1.000
##     5      9       1    0.810  0.1210        0.572        1.000
##     6      7       1    0.702  0.1451        0.417        0.986
##     8      6       1    0.594  0.1578        0.285        0.903
##     9      5       1    0.486  0.1617        0.169        0.803
##    11      2       1    0.295  0.1772        0.000        0.642
## 
##                 gr_3=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5      8       1    0.882   0.110       0.6663        1.000
##     6      6       1    0.747   0.156       0.4420        1.000
##     7      4       1    0.582   0.189       0.2107        0.953
##     8      3       1    0.417   0.194       0.0362        0.797

plot(ena_2, mark.time = TRUE, lty = c(1,2,3), xlab = "Tempo", ylab = "S(t) por N-A")

legend(1,0.3,lty = c(1,2,3), c("Controlo","Tratamento 1", "Tratamento 2"), bty="n")

Os estimadores de Kaplan-Meier e de Nelson-Aalen para a função de sobrevivência dão geralmente estimativas muito próximas, como se pode verificar por comparação das figuras 2.3 e 2.4. Ambos geram estimativas que são funções em escada nos mesmos intervalos, que são definidos pelos instantes de falha.

Como é visível, nestes modelos não paramétricos só é permitido o estudo da influência de covariáveis categórias. Sempre que existir uma covariável contínua, esta deve ser considerada por classes, por forma a permitir considerar e comparar as estimativas obtidas em cada classe. Para a incorporação de covariáveis contínuas precisamos de outro tipo de métodos, nomeadamente paramétricos ou semiparamétricos.

Estimador da Tabela Atuarial

A construção de uma tabela de vida consiste em dividir o eixo do tempo em um certo número de intervalos.

Suponha que o eixo do tempo seja dividido em s intervalos definidos pelos pontos de corte, $t_1, \, t_2, ... , \, t_s$. Isto é, $I_j = [t_{j-1}, \, t_j)$, para $j = 1, ..., s$, em que $t_ = 0$ e $t_s = +\infty$.

O estimador da tabela de vida apresenta a forma do estimador de Kaplan-Meier, mas utiliza um estimador ligeiramente diferente, uma vez que, neste caso, tem-se para $d_j$ e $n_j$ que:

$d_j = \, \text{n.º de falhas no intervalo} [t_{j-1}, \, t_j)$ e
$\small n_j = \left[ \text{n.º sob risco em } \, t_{j-1} \right] - \left[ \frac{1}{2} \times \text{n.º de censuras em } \, [t_{j-1}, \, t_j) \right]$

Assim, a estimativa para $q_j$ na tabela de vida é dada por:

\[ \widehat{q}_j=\frac{d_j}{n_j} \]

A explicação para o $\left[ \frac{1}{2} \times \text{n.º de censuras em } \, [t_{j-1}, \, t_j) \right]$ é que observações para as quais a censura ocorreu no intervalo $[t_{j-1}, \, t_j)$ são tratadas como se estivessem sob risco durante a metade do intervalo considerado.

O estimador da tabela de vida fica expresso por:

\[ \widehat S(t) = \prod_{i=1}^j \left(1-\widehat q_{i-1}\right) \]

para $j = 1,..., s$ e $\widehat q_0 = 0$. A representação gráfica da função de sobrevivência é uma escada, com valor constante em cada intervalo de tempo.

Suponha o exemplo da hepatite com os dados do grupo esteróide divididos em 4 intervalos: $[0; 5)$, $[5, 10)$, $[10, 15)$ e $[15, 16)$. A estimativa de $\widehat q_2$ correspondente ao intervalo $[5, 10)$ é:

\[ \widehat{q}_2=\frac{3}{9}=0,33 \]

Isto significa que a probabilidade de morte até a $\text{10ª}$ semana da terapia com esteróide para aqueles que sobreviveram à $\text{5ª }$semana é de 33,3%. О cálculo pode ser estendido da mesma forma para os outros intervalos e estes valores são mostrados na Tabela seguinte.

A estimativa para a função de sobrevivência no tempo $t = 10$ semanas é:

\[ \widehat{S}(10)=(1-0,231) \cdot (1-0,333)=0,513 \]

Isto significa que um paciente no grupo esteróide tem uma probabilidade de 51,3% de sobreviver a 10 semanas de tratamento.

Na Tabela seguinte estão também apresentados os valores estimados da função de sobrevivência para os outros intervalos de tempo.

Tabela 6: Estimativas de Kaplan-Meier para o grupo esteróide.

\[ \begin{array}{cccccc} \text{Intervalo} & \text{N.º sob} & \text{N.º de} & \text{N.º de} & \widehat q_j & \widehat S(t_{I_j}) \\ I_j & \text{Risco} & \text{Falhas} & \text{Censuras} & \text{(%)} & \text{(%)} \\ \hline [0, \, \, 5) & 14 & 3 & 2 & 23,1 & 100,0 \\ [5, \, \, 10) & 9 & 3 & 0 & 33,3 & 76,9 \\ [10, \, \, 15) & 6 & 1 & 2 & 20,0 & 51,3 \\ [15, \, \, \infty) & 3 & 0 & 3 & 0,0 & 41,0 \\ \hline \end{array} \]

A variância assintótica estimada para $\widehat S(t_{I_j})$ é, neste caso, obtida por:

\[ \widehat{\mathrm{Var}}(\widehat S (I_j)) \cong \left[ \widehat S(I_j) \right]^2 \cdot \sum_{\ell=1}^j \left(\frac{\widehat q_\ell}{n_\ell \cdot (1-\widehat q_\ell)}\right) \,, \, \, \, \, \, j=2,\,...,\,s \]

Vantagens e Desvantagens

Fácil de entender;
Suposições fracas (não impõe distribuição para $T$).
Pouco eficientes;
Díficil de incluir covariáveis na análise estatística.

Comparação dos Estimadores

A grande diferença entre os estimadores de $S(t)$ está no número de intervalos utilizados para a construção de cada um deles. O estimador de Kaplan-Meier e o de Nelson-Aalen são sempre baseados em um número de intervalos igual ao número de tempos de falha distintos, enquanto que, na tabela de vida, os tempos de falha são agrupados em intervalos de forma arbitrária. Isto faz com que a estimativa obtida pelo estimador de KaplanMeier seja baseada freqüentemente em um número de intervalos maior que a estimativa obtida por meio da tabela de vida.

No exemplo discutido nesta seção, o eixo do tempo foi dividido em cinco intervalos de tempo, correspondendo a cada falha distinta, para o estimador de Kaplan-Meier, enquanto que, no estimador da tabela de vida foram utilizados quatro intervalos de tempo. É natural esperar que quanto maior o número de intervalos, melhor será a aproximação para a verdadeira distribuição do tempo de falha.

Pode-se então perguntar: por que não usar cinco ou mais intervalos para o cálculo do estimador da tabela de vida? Isto poderia ser feito. No entanto, observa-se, na prática, que isto não acontece devido às suas origens. A justificativa reside no fato deste estimador ter sido proposto por demógrafos e atuários no século passado e usado sempre em grandes amostras (por exemplo, proveniente de censos demográficos).

A divisão em um número arbitrário e grande de intervalos é justificada por ser a amostra muito grande, o que não acontece em resultados provenientes de estudos clínicos ou ensaios de confiabilidade.

O uso da tabela de vida, considerando um número igual ou maior de intervalos que o do estimador de Kaplan-Meier, gera estimativas exatamente iguais às estimativas de Kaplan-Meier, se o mecanismo de censura for do tipo I ou do tipo II. Entretanto, se o mecanismo de censura for do tipo aleatório, as estimativas são próximas, mas não necessariamente coincidentes.

Nesta última situação, alguns autores estudaram as propriedades assintóticas dos dois estimadores. Estes estudos mostraram a superioridade do estimador de Kaplan-Meier. Ele é um estimador não-viciado para a função de sobrevivência eme grandes amostras, enquanto o estimador da tabela de vida não o é, com um vício que fica pequeno à medida que o comprimento dos intervalos diminui.

Com amostras de pequeno ou médio porte, existe também alguma evidência empírica da superioridade do estimador de Kaplan-Meier. Desta forma, o mais indicado é, então, usar o estimador de Kaplan-Meier ou eventualmente o de Nelson-Aalen, em vez daquele da tabela de vida, quando o interesse do pesquisador se concentrar em informações provenientes da função de sobrevivência.

Comparação de Curvas de $\small S(t)$

Uma das aplicaçõe principais da análise de sobrevivência é a comparação entre grupos, isto é, determinar se existem diferenças estatisticamente significativas entre as funções de sobrevivência de dois ou mais grupos distintos da população. Os grupos podem ser várias regiões do país, várias faixas etárias, grupos de controlo e de tratamento, entre outros.

O teste de Kaplan-Meier também pode ser usado, por exemplo, para comparar a eficácia de dois ou mais tipos de tratamento.

Por exemplo, podem ser definidos como objectivos do estudo:

Determinar se um determinado medicamento ou protocolo médico tem efeito na taxa de sobrevivência dos doentes e se é mais eficaz que outro.
Determinar se o tempo (ou a distribuição do tempo) que decorre até ocorrer uma determinada doença difere do tipo de estilo de vida da pessoa, da sua idade e de outros factores.

Nestes casos, vai haver um factor de comparação. No primeiro, o factor vai ser os tipos de medicamentos que foi utilizada. O teste vai comparar a distribuição de sobrevivência entre, por exemplo, 4 tipos de medicamentos e determinar se a sua eficácia é igual. No caso da sua eficácia não ser semelhante, o teste dá indicações que um determinado medicamento é mais eficaz que outro.

No caso de se utilizar apenas um medicamento mas com doses diferentes, a dosagem seria o factor de comparação.

No segundo caso, o factor estilo de vida poderia ter três grupos de pessoas: os desportistas, os semi-activos e os sedentários. O teste indicaria se as distribuições de sobrevivência são diferentes para cada estilo de vida e que diferenças existem.

`Exemplo 4.b): Hipotético`

Pergunta: Determine a estimativa de Kaplan-Meier da função de sobrevivência para cada um dos três grupos.

Temos então a estimativa de Kaplan-Meier para a função de sobrevivência em cada grupo e os respetivo gráfico dados pelos comandos:

ekm_3 = survfit(Surv(t_3, cens_3)~gr_3);                    summary(ekm_3)

## Call: survfit(formula = Surv(t_3, cens_3) ~ gr_3)
## 
##                 gr_3=0 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5      7       1    0.857   0.132        0.633            1
##     8      6       1    0.714   0.171        0.447            1
##     9      4       1    0.536   0.201        0.257            1
##    10      3       1    0.357   0.198        0.121            1
##    14      1       1    0.000     NaN           NA           NA
## 
##                 gr_3=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     4     10       1    0.900  0.0949       0.7320        1.000
##     5      9       1    0.800  0.1265       0.5868        1.000
##     6      7       1    0.686  0.1515       0.4447        1.000
##     8      6       1    0.571  0.1638       0.3258        1.000
##     9      5       1    0.457  0.1662       0.2242        0.932
##    11      2       1    0.229  0.1817       0.0481        1.000
## 
##                 gr_3=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     5      8       1    0.875   0.117        0.673            1
##     6      6       1    0.729   0.165        0.468            1
##     7      4       1    0.547   0.201        0.266            1
##     8      3       1    0.365   0.200        0.124            1

plot(ekm_3, mark.time = TRUE, lty = c(1,2,3), xlab = "Tempo", ylab = "S(t) por Kaplan-Meier")

legend(10,0.95, lty = c(1,2,3), c("Controlo", "Tratamento 1", "Tratamento 2"))

A comparação na figura ilustra-nos que aparentemente não existe grande diferença entre os vários grupos. No entanto, é necessário fazer um teste estatístico adequado para a verificação desta hipótese.

Teste Log-Rank

Um procedimento natural usaria os resultados assintóticos de $\widehat S(t)$, apresentados na seção anterior, para testar a igualdade de funções de sobrevivência em um determinado tempo $t$. Esta forma, no entanto, não faria uso eficiente dos dados disponíveis, pois não se estaria usando todo o período do estudo. Estatísticas mais comumente usadas podem ser vistas como generalizações para dados censurados, de conhecidos testes não-paramétricos. O teste Log-Rank (Mantel, 1966) é o mais usado em análise de sobrevivência.

Gehan (1965) propôs uma generalização para a estatística de Wilcoxon. Outras generalizações foram propostas por Peto e Peto (1972) e Prentice (1978), entre outros. Latta (1981) fez uso de simulações de Monte Carlo para comparar vários testes não-paramétricos.

Neste texto, ênfase será dada ao teste Log-Rank. Este teste é muito utilizado em análise de sobrevivência e é particularmente apropriado quando a razão das funções de risco dos grupos a serem comparados é aproximadamente constante. Isto é, as populações têm a propriedade de riscos proporcionais. A estatística deste teste é a diferença entre o número observado de falhas em cada grupo e uma quantidade que, para muitos propósitos, pode ser pensada como o correspondente número esperado de falhas sob a hipótese nula, ou seja, compara a distribuição de ocorência dos eventos observados em cada estrato com a distribuição que seria esperada se a incidência fosse igual em todos os estratos.

A expressão do teste Log-Rank é obtida de forma similar a do conhecido teste de Mantel-Haenszel (1959), para combinar tabelas de contingência. O teste Log-Rank tem, também, a mesma expressão do teste escore para o modelo de regressão de Cox.

Considere, inicialmente, o teste de igualdade de duas funções de sobrevivência $S_1(t)$ e $S_2(t)$. Sejam $t_1 < t_2 < ... <t_k$ os tempos de falha distintos da amostra formada pela combinação das duas amostras individuais.

Suponha que no tempo $t_j$ aconteçam $d_j$ falhas e que $n_j$ indivíduos estejam sob risco em um tempo imediatamente inferior a $t_j$ na amostra combinada e, respectivamente, $d_{ij}$ e $n_{ij}$ na amostra $i$; $i = 1, \,2$ e $ј = 1,..., k$.

Em cada tempo de falha $t$, os dados podem ser dispostos em forma de uma tabela de contingência $2 \times 2$ com $d_{ij}$ e $n_{ij}$ falhas e e $n_{ij}-d_{ij}$ sobreviventes na coluna $i$. Isto é mostrado na Tabela 7.

Tabela 7: Tabela de contingência gerada no tempo $t_j$.

\[ \begin{array}{c|cc|c|c} & \text{Grupo 1} & \text{Grupo 2} & --- & \text{Grupo }p & \text{Total} \\ \hline \text{Falha} & d_{1j} & d_{2j} & --- & d_{pj} & d_{j} \\ \text{Não Falha} & n_{1j}-d_{1j} & n_{2j}-d_{2j} & --- & n_{pj}-d_{pj} & n_{j}-d_{j} \\ \hline \text{Total} & n_{1j} & n_{2j} & --- & n_{pj} & n_{j} \\ \end{array} \]

Condicional à experiência de falha e censura até o tempo $t_j$ (fixando as marginais de coluna) e ao número de falhas no tempo $t_j$ (fixando as marginais de linha), a distribuição de $d_{2j}$ é, então, uma hipergeométrica:

\[ \frac{ \binom{n_{1j}}{d_{1j}} \binom{n_{2j}}{d_{2j}}}{\binom{n_j}{d_j}} \]

A média de $d_{2j}$ é $w_{2j} = n_{2j}d_jn_j^{-1}$, o que equivale a dizer que, se não houver diferença entre as duas populações no tempo $t_j$, o número total de falhas $(d_j)$ pode ser dividido entre as duas amostras de acordo com a razão entre o número de indivíduos sob risco em cada amostra e o número total sob risco.

A variância de $d_{2j}$ obtida a partir da distribuição hipergeométrica é:

\[ \left( V_j\right)_2 = n_{2j}\cdot(n_j-n_{2j})\cdot d_j \cdot(n_j-d_j)\cdot n_j^{-2} \cdot(n_j-1)^{-1} \]

Então, a estatística $d_{2j}-w_{2j}$ tem média zero e variância $(V_j)_2$. Se as $k$ tabelas de contingência forem independentes, um teste aproximado para a igualdade das duas funções de sobrevivência pode ser baseado na estatística:

\[ T = \frac{ \left[ \sum_{j=1}^{k} \left( d_{2j} - w_{2j}\right) \right]^2}{\sum_{j=1}^{k} \left( V_j\right)_2} \, \, \, \, \sim \, \, \, \, \chi_1^2 \]

que, sob a hipótese nula $H_0: S_1(t) = S_2(t)$ para todo $t$ no período de acompanhamento, tem uma distribuição qui-quadrado com 1 grau de liberdade para grandes amostras.

O objetivo principal do estudo dos dados de hepatite é comparar a terapia com esteróide e o grupo controle. As curvas de Kaplan-Meier para os dois grupos apresentadas na Figura 7 indicam que, possivelmente, a terapia com esteróide não é um tratamento adequado para pacientes com hepatite viral aguda.

No entanto, é necessária uma evidência quantitativa deste fato e, sendo assim, foi utilizado um teste de significância. O valor do teste Log-Rank para a comparação entre os dois grupos resultou em $T = 3,67$, o que implica em um valor $p = 0,055$, indicando uma diferença entre as duas curvas de sobrevivência. O valor deste teste e seu correspondente valor p podem ser obtidos no pacote estatístico R por meio dos comandos:

 t_5 = 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_5 = 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)
 
 gr_5 = c(rep(1,15), rep(2,14));          survdiff(Surv(t_5, cens_5)~gr_5, rho = 0)

## Call:
## survdiff(formula = Surv(t_5, cens_5) ~ gr_5, rho = 0)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## gr_5=1 15        2     4.81      1.64      3.67
## gr_5=2 14        7     4.19      1.89      3.67
## 
##  Chisq= 3.7  on 1 degrees of freedom, p= 0.06

A generalização do testeLog-Rankpara a igualdade de r > 2 funções de sobrevivência $S_i(t), ..., S_r(t)$ não é complicada. Considere a mesma notação anterior, com o índice $i$ variando, agora, entre $1$ e $r$. Desta forma, os dados podem ser arranjados em forma de uma tabela de contingência $2 \times r$ com $d_{ij}$ falhas e $n_{ij} - d_{ij}$ sobreviventes na coluna $i$. Ou seja, a Tabela 7 passaria a ter $r$ colunas em vez de simplesmente duas.

Condicional à experiência de falha e censura até o tempo $t$; e ao número de falhas no tempo $t_j$, a distribuição conjunta de $d_{2j}, ..., d_{rj}$ é, então, uma hipergeométrica multivariada, isto é,

\[ \frac{ \prod_{i=1}^r \binom{n_{1j}}{d_{1j}} }{\binom{n_j}{d_j}} \]

A média de $d_{ij}$ é $w_{ij} = n_{ij}d_jn_j^{-1}$, bem como a variância de $d_{ij}$ e a covariância de $d_{ij}$ e $d_{lj}$ são, respectivamente, e

\[ \left( V_j\right)_{ii} = n_{ij}\cdot(n_j-n_{ij})\cdot d_j \cdot(n_j-d_j)\cdot n_j^{-2} \cdot(n_j-1)^{-1} \] e \[ \left( V_j\right)_{il} = - n_{ij}\cdot n_{lj}\cdot d_j \cdot(n_j-d_j)\cdot n_j^{-2} \cdot(n_j-1)^{-1} \]

Então, a estatística $v'_j = (d_{2j} - w_{2j}, ..., d_{rj} - w_{rj})$ tem média zero e matriz de variância-covariância $V_j$ de dimensão $r - 1$, com $(V_j)_ii$ na diagonal principal e os elementos $(V_j)_ii$, $i = 2, ..., r$ fora da diagonal principal.

Pode-se, então, formar a estatística $v$, somando sobre todos os tempos distintos de falha, isto é,

\[ v = \sum_j^k v_j \]

com $v$ um vetor de dimensão $(r-1) \times 1$, cujos elementos são as diferenças entre os totais observados e esperados de falha.

Considerando, novamente, a suposição de que as k tabelas de contingência são independentes, a variância da estatística $v$ será

\[V = V_1 + V_2 \, \, + \, \,,..., \, \,V_k\]

Um teste aproximado para a igualdade das r funções de sobrevivência pode ser baseado na estatística:

\[ T=v' \cdot V^{-1} \cdot v \, \, \, \, \sim \, \, \, \, \chi_{r-1}^2 \]

que, sob $H_0$ (igualdade das curvas), tem uma distribuição qui-quadrado com $r-1$ graus de liberdade para amostras grandes. Os graus de liberdade são $r-1$ e não $r$, pois os elementos de $v$ somam zero.

`Exemplo 5: Dados`

O valor da estatística Log-Rank que, sob a hipótese de igualdade das curvas de sobrevivência, tem uma distribuição qui-quadrado com dois graus de liberdade, resultou em $T = 12,6$. Isto gera um $\text{valor } p = 0,0019$, o que indica a existência de diferenças entre os grupos.

Constatada a presença de diferenças entre os grupos, existe, então, a necessidade de identificar quais curvas diferem entre si. Isto é usualmente chamado de comparações múltiplas.

Em planeamento de experimentos em que é assumido um modelo linear com resposta normal, existem vários métodos disponíveis para a realização de tais comparações. O mesmo não acontece com dados de sobrevivência.

De forma a encontrar as diferenças entre os grupos, uma possibilidade é fazer comparações dos grupos, dois a dois, controlando o erro do tipo I pelo método de Bonferroni. Como existem três grupos, três testes dois a dois entre os grupos são possíveis.

O método de Bonferroni utiliza um nível de significância de $0,05/3 = 0,017$ para cada um dos testes, de forma a garantir uma conclusão geral a um nível de no máximo 0,05.

A Tabela 2.7 mostra os resultados dos testes Log-Rank realizados para as comparações dos grupos dois a dois.

t_6 = c(7,8,8,8,8,12,12,17,18,22,30,30,30,30,30,30, 8,8,9,10,10,14, 15,15,18,19,21,22,22,23,25,8,8,
        8,8,8,8,9,10,10,10,11,17,19)

 cens_6 = c(rep(1,10), rep(0,6), rep(1,15), rep(1,13)); gr_6 = c(rep(1,16), rep(2,15), rep(3,13))
 
 ekm_6 = survfit(Surv(t_6, cens_6)~gr_6);               summary(ekm_6)

## Call: survfit(formula = Surv(t_6, cens_6) ~ gr_6)
## 
##                 gr_6=1 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     7     16       1    0.938  0.0605        0.826        1.000
##     8     15       4    0.688  0.1159        0.494        0.957
##    12     11       2    0.562  0.1240        0.365        0.867
##    17      9       1    0.500  0.1250        0.306        0.816
##    18      8       1    0.438  0.1240        0.251        0.763
##    22      7       1    0.375  0.1210        0.199        0.706
## 
##                 gr_6=2 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     8     15       2   0.8667  0.0878       0.7106        1.000
##     9     13       1   0.8000  0.1033       0.6212        1.000
##    10     12       2   0.6667  0.1217       0.4661        0.953
##    14     10       1   0.6000  0.1265       0.3969        0.907
##    15      9       2   0.4667  0.1288       0.2717        0.802
##    18      7       1   0.4000  0.1265       0.2152        0.743
##    19      6       1   0.3333  0.1217       0.1630        0.682
##    21      5       1   0.2667  0.1142       0.1152        0.617
##    22      4       2   0.1333  0.0878       0.0367        0.484
##    23      2       1   0.0667  0.0644       0.0100        0.443
##    25      1       1   0.0000     NaN           NA           NA
## 
##                 gr_6=3 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     8     13       6   0.5385  0.1383       0.3255        0.891
##     9      7       1   0.4615  0.1383       0.2566        0.830
##    10      6       3   0.2308  0.1169       0.0855        0.623
##    11      3       1   0.1538  0.1001       0.0430        0.550
##    17      2       1   0.0769  0.0739       0.0117        0.506
##    19      1       1   0.0000     NaN           NA           NA

 plot(ekm_6, lty = c(1,4,2), mark.time = TRUE, xlab = "Tempo", ylab = "S(t) estimada")
 
 legend(1, 0.3, lty = c(1,4,2), c("Grupo 1","Grupo 2","Grupo 3"), lwd = 1, bty = "n", cex = 0.8)

Dos resultados apresentados na Tabela 2.7 pode-se concluir pela existência de diferenças significativas entre os grupos 1 e 3 e entre os grupos 2 e 3.

Entre os grupos 1 e 2, não foram encontradas evidências de diferenças.

A. diferença entre os grupos 1 e 3 atesta a eficácia da imunização pela malária na presença de infecções pela malária e pela esquistossomose.

Por outro lado, a diferença entre os grupos 2 e 3 mostra o impacto na mortalidade dos camundongos devido à infecção pela esquistossomose.

Tabela 2.7: Resultados dos testes Log-Rank utilizados para as comparações dos grupos, dois a dois, considerados no estudo da malária.

\[ \begin{array}{c|c|c} \text{Grupos comparados} & \text{Estatística de Teste} & \text{Valor }p \\ \hline 1 \times 2 & 2,5 & 0,112 \\ 2 \times 3 & 8,0 & 0,005 \\ 1 \times 3 & 7,9 & 0,005 \\ \hline \end{array} \]

Esta análise foi realizada no R por meio dos comandos a seguir:

 survdiff(Surv(t_6, cens_6)~gr_6, rho = 0)

## Call:
## survdiff(formula = Surv(t_6, cens_6) ~ gr_6, rho = 0)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## gr_6=1 16       10    17.00    2.8816    6.4111
## gr_6=2 15       15    14.51    0.0167    0.0317
## gr_6=3 13       13     6.49    6.5190   10.4447
## 
##  Chisq= 12.6  on 2 degrees of freedom, p= 0.002

 survdiff(Surv(t_6[1:31], cens_6[1:31])~gr_6[1:31], rho = 0)

## Call:
## survdiff(formula = Surv(t_6[1:31], cens_6[1:31]) ~ gr_6[1:31], 
##     rho = 0)
## 
##               N Observed Expected (O-E)^2/E (O-E)^2/V
## gr_6[1:31]=1 16       10     13.7      1.01      2.53
## gr_6[1:31]=2 15       15     11.3      1.23      2.53
## 
##  Chisq= 2.5  on 1 degrees of freedom, p= 0.1

 survdiff(Surv(t_6[17:44], cens_6[17:44])~gr_6[17:44], rho = 0)

## Call:
## survdiff(formula = Surv(t_6[17:44], cens_6[17:44]) ~ gr_6[17:44], 
##     rho = 0)
## 
##                N Observed Expected (O-E)^2/E (O-E)^2/V
## gr_6[17:44]=2 15       15    20.53      1.49      7.98
## gr_6[17:44]=3 13       13     7.47      4.08      7.98
## 
##  Chisq= 8  on 1 degrees of freedom, p= 0.005

 survdiff(Surv(c(t_6[1:16],t_6[32:44]),c(cens_6[1:16],cens_6[32:44]))~c(gr_6[1:16], gr_6[32:44]),rho=0)

## Call:
## survdiff(formula = Surv(c(t_6[1:16], t_6[32:44]), c(cens_6[1:16], 
##     cens_6[32:44])) ~ c(gr_6[1:16], gr_6[32:44]), rho = 0)
## 
##                               N Observed Expected (O-E)^2/E (O-E)^2/V
## c(gr_6[1:16], gr_6[32:44])=1 16       10    15.34      1.86      7.86
## c(gr_6[1:16], gr_6[32:44])=3 13       13     7.66      3.72      7.86
## 
##  Chisq= 7.9  on 1 degrees of freedom, p= 0.005

 # pode-se, alternativamente, usar a opção "subset"
 
 data = as.data.frame(cbind(t_6, cens_6, gr_6)) 
 
 # grupos 1 e 2
 data_1 = subset(data, subset = gr_6 %in% c(1,2)); survdiff(Surv(t_6, cens_6) ~ gr_6, rho = 0, data_1)

## Call:
## survdiff(formula = Surv(t_6, cens_6) ~ gr_6, data = data_1, rho = 0)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## gr_6=1 16       10     13.7      1.01      2.53
## gr_6=2 15       15     11.3      1.23      2.53
## 
##  Chisq= 2.5  on 1 degrees of freedom, p= 0.1

 # grupos 2 e 3
 data_2 = subset(data, subset = gr_6 %in% c(2,3)); survdiff(Surv(t_6, cens_6) ~ gr_6, rho = 0, data_2)

## Call:
## survdiff(formula = Surv(t_6, cens_6) ~ gr_6, data = data_2, rho = 0)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## gr_6=2 15       15    20.53      1.49      7.98
## gr_6=3 13       13     7.47      4.08      7.98
## 
##  Chisq= 8  on 1 degrees of freedom, p= 0.005

 # grupos 1 e 3 
 data_3 = subset(data, subset = gr_6 %in% c(1,3)); survdiff(Surv(t_6, cens_6) ~ gr_6, rho = 0, data_3)

## Call:
## survdiff(formula = Surv(t_6, cens_6) ~ gr_6, data = data_3, rho = 0)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## gr_6=1 16       10    15.34      1.86      7.86
## gr_6=3 13       13     7.66      3.72      7.86
## 
##  Chisq= 7.9  on 1 degrees of freedom, p= 0.005

Aplicar o teste de Log-Rank, para o Exemplo 4

survdiff(Surv(t_3, cens_3)~gr_3, rho = 0)

## Call:
## survdiff(formula = Surv(t_3, cens_3) ~ gr_3, rho = 0)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## gr_3=0  7        5     5.92  0.143973    0.2920
## gr_3=1 10        6     5.93  0.000778    0.0015
## gr_3=2  8        4     3.14  0.232789    0.3460
## 
##  Chisq= 0.5  on 2 degrees of freedom, p= 0.8

Como se verifica, temos um $\text{valor-} p \text{ de }0.795$, pelo que para uma significância de $5%$ não existe evidência estatística de diferenças entre os três grupos. Assim, não só não existe diferenças entre a administração do medicamento em 2 ou 3 doses, nem existe diferença entre a administração ou não do medicamento. Desta forma, concluímos que não existe evidência estatística que o efeito do medicamento seja eficaz.

Teste de Wilcoxon

No caso particular da comparação de duas funções de sobrevivência, a seguinte forma geral inclui os testes mais importantes na literatura e generaliza a estatística $T$, dada por:

\[ S = \frac{ \left[ \sum_{j=1}^{k} u_j \cdot \left( d_{2j} - w_{2j}\right) \right]^2}{\sum_{j=1}^{k} u_j^2 \cdot \left( V_j\right)_2} \]

com $u_j$ os pesos que especificam os testes.

Sob a hipótese nula de que as funções de sobrevivência não diferem, a estatística $S$ tem distribuição quiquadrado com 1 grau de liberdade para amostras grandes. O teste Log-Rank é obtido tomando-se $u_j= 1$, para $j = 1, \, ..., \, \, k$.

Outro teste bastante utilizado na prática é o de Wilcoxon obtido quando se toma $u_j = n_j$. Este teste foi adaptado para dados censurados a partir do conhecido teste não-paramétrico de Wilcoxon (Gehan, 1965, Breslow, 1970).

Teste de Tarone e Ware

O teste de Tarone e Ware (1977) propõe peso $u_j = \sqrt{n_j}$, que fica entre os pesos dos testes Log-Rank e de Wilcoxon.

Tabela 2.8: Testes não-paramétricos para comparação das curvas de sobrevivência obtidas para os grupos esteróide e controle dos dados de hepatite.

\[ \begin{array}{lcc} \text{Testes} & \text{Estatística de Teste} & \text{Valor }p \\ \hline Log-Rank \text{ (Mantel-Haenszel)} & 3,67 & 0,055 \\ \text{Wilcoxon-Gehan (Breslow)} & 3,19 & 0,074 \\ \text{Tarono-Ware} & 3,43 & 0,064 \\ \hline \end{array} \]

A escolha do peso na expressão direciona o tipo de diferença a ser detectado nas funções de sobrevivência. O teste de Wilcoxon, que utiliza peso igual ao número de indivíduos sob risco, coloca mais peso na porção inicial do eixo do tempo.

No início do estudo, todos os indivíduos estão sob risco e saindo do estado “sob risco” à medida que falham ou são censurados. O teste Log-Rank, por outro lado, coloca mesmo peso para todo o eixo do tempo, o que reforça o enfoque nos tempos maiores quando comparado ao teste de Wilcoxon. O teste de Tarone-Ware se localiza em uma situação intermediária.

Portante, segue o exemplo de aplicação dos testes acima (a partir de outro conjunto de dados) com maior exactidão no r:

 t_7 = 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_7 = 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)
 
 gr_7 = c(rep("Controle",15),rep("Esteroide",14))
 
 logrank_test(Surv(t_7,cens_7)~factor(gr_7), type = "logrank")

## 
##  Asymptotic Two-Sample Logrank Test
## 
## data:  Surv(t_7, cens_7) by factor(gr_7) (Controle, Esteroide)
## Z = 1.9245, p-value = 0.05429
## alternative hypothesis: true theta is not equal to 1

 logrank_test(Surv(t_7,cens_7)~factor(gr_7), type = "Gehan-Breslow")

## 
##  Asymptotic Two-Sample Gehan-Breslow Test
## 
## data:  Surv(t_7, cens_7) by factor(gr_7) (Controle, Esteroide)
## Z = 1.7897, p-value = 0.07349
## alternative hypothesis: true theta is not equal to 1

 logrank_test(Surv(t_7,cens_7)~factor(gr_7), type = "Tarone-Ware")

## 
##  Asymptotic Two-Sample Tarone-Ware Test
## 
## data:  Surv(t_7, cens_7) by factor(gr_7) (Controle, Esteroide)
## Z = 1.859, p-value = 0.06303
## alternative hypothesis: true theta is not equal to 1

Outros testes

Peto e Peto (1972) e Prentice (1978) sugerem utilizar uuma função do peso que depende diretamente da experiência passada de sobrevivência observada das duas amostras combinadas. A função do peso é uma modificação do estimador de Kaplan-Meier e é definido de tal forma que seu valor é conhecido antes da falha ocorrer. O estimador modificado da função de sobrevivência é:

\[ \widetilde{S}(t) = \prod_{j : \, t_j < t} \left( \frac{n_j+1-d_j}{n_j+1} \right) \]

e os pesos utilizados são:

\[ u_j = \widetilde{S}(t_{j-1}) \cdot \frac{n_j}{n_j+1} \]

Este estimador é conhecido por Peto-Prentice.

Outra classe de pesos para a expressão foi proposta por Harrington-Fleming (1982) e é dada por:

\[ u_j = \left[ \widetilde{S}(t_{j-1}) \right]^p \]

Se $\rho = 0$, obtém-se $u_j = 1$ e tem-se, então, o teste Log-Rank.

Entretanto, se $\rho = 1$, o peso é o Kaplan-Meier no tempo de falha anterior e, nesse caso, tem-se um teste similar ao de Wilcoxon.

A principal vantagem dos testes de Peto-Prentice e Harrington-Fleming é que a ponderação é feita relativa à experiência de sobrevivência anterior. Isto não acontece com o teste Log-Rank.

O teste de Wilcoxon, em particular, pondera pelo número de indivíduos sob risco que depende da experiência de sobrevivência, assim como da. de censura.

Se o padrão de censura é nitidamente diferente nos dois grupos, então o teste pode rejeitar ou não rejeitar, não somente com base nas diferenças das sobrevivências entre os grupos mas, também, devido ao padrão de censura.

 logrank_test(Surv(t_7,cens_7)~factor(gr_7), type = "Peto-Peto")

## 
##  Asymptotic Two-Sample Peto-Peto Test
## 
## data:  Surv(t_7, cens_7) by factor(gr_7) (Controle, Esteroide)
## Z = 1.8571, p-value = 0.06329
## alternative hypothesis: true theta is not equal to 1

 logrank_test(Surv(t_7,cens_7)~factor(gr_7), type = "Fleming-Harrington")

## 
##  Asymptotic Two-Sample Fleming-Harrington Test
## 
## data:  Surv(t_7, cens_7) by factor(gr_7) (Controle, Esteroide)
## Z = 1.9245, p-value = 0.05429
## alternative hypothesis: true theta is not equal to 1

Tópicos de Análise de Sobrevivência