ANÁLISE DE SOBREVIVÊNCIA E CONFIABILIDADE

AULA 7: MODELO DE REGRESSÃO DE COX

\[\\[0.05in]\]

1 PACOTES NECESSÁRIOS

Para esta sétima aula prática de Análise de Sobrevivência e Confiabilidade, precisaremos do seguinte pacote:

require(survival)
  • Pacote survival: É o pacote básico para o estudo de Análise de Sobrevivência e Confiabilidade. Contém as principais rotinas, incluindo definição de objetos Surv, curvas Kaplan-Meier, modelos Cox e modelos paramétricos.

2 BASES DE DADOS

2.1 CÂNCER DE LARINGE

Neste exemplo, os dados considerados referem-se a um estudo, descrito em Klein e Moeschberger (2003), realizado com 90 pacientes do sexo masculino diagnosticados no perÍodo com cancer de laringe. Para cada paciente, foram registrados, no diagnóstico, a idade (em anos) e o estágio da doença (I = tumor primário, II = envolvimento de nódulos, III = metástases e IV = combinações dos 3 estágios anteriores), bem como seus respectivos tempos de morte ou censura (em meses). Os estágios encontram-se ordenados pelo grau de seriedade da doença (menos sério para mais sério). A base de dados está armazenada no objeto laringe.txt.

2.2 CANCER DE OVÁRIO

Um estudo realizado para comparar dois tratamentos pós-cirúrgicos de câncer de ovário, envolveu o acompanhamento de 26 mulheres após a cirurgia de remoção do tumor. A resposta foi o tempo (em dias), contado a partir do início do tratamento (aleatorização) até a morte do paciente. As seguintes covariáveis foram registradas: X1 = tratamento, X2 = idade, X3 = resíduo (1 se o resíduo da doença foi parcialmente removido e 2 se foi completamente removido) e status (1 se a condição do doente no início do estudo era boa e 2 se ruim). Os dados encontram-se no arquivo ovario.txt

3 INTRODUÇÃO

O modelo de regressão de Cox permite a análise de dados em que a resposta é o tempo até a ocorrência de um evento de interesse, ajustando por covariáveis. De forma genérica, suponha a existência de \(p\) covariáveis, de modo que \(\pmb{x}\) seja um vetor com os componentes \(\pmb{x}=(x_1,x_2,...,x_p)^t\). A expressão geral do modelo de regressão de Cox é:

\[\lambda(t|\pmb{x})=\lambda_0(t)exp\{\pmb{\beta}\pmb{x}\}\]

Este modelo é composto pelo produto de dois componentes, um não paramétrico e outro paramétrico.

Conforme mencionado em sala, o modelo de Regressão de Cox também é chamado de modelo de riscos proporcionais. Os riscos proporcionais constituem uma suposição para este modelo. Isto é, a razão entre as funções de risco de dois indivíduos com covariáveis em valores diferentes é constante no tempo.

Para começar a nossa análise dos dados de câncer de laringe, vamos carregar a base de dados:

dados<-read.table("laringe.txt",header=T)

4 ESTIMAÇÃO DOS PARÂMETROS \(\pmb{\beta}\)

O modelo de Cox é caracterizado pelos coeficientes \(\pmb{\beta}\), que medem os efeitos das covariáveis sobre a função de risco. Para que o modelo fique determinado, estas quantidades devem ser estimadas a partir dos dados amostrais. Assim, consideramos, novamente, a situação em que \(n\) elementos foram observados e, em cada um deles, foi coletado o vetor:

\[(t_i,\delta_i,x_{i1},x_{i2},...,x_{ip})^t\]

em que, para o elemento \(i\):

  • \(t_i\) é o tempo observado.
  • \(\delta_i\) é a função indicadora de ocorrência do evento.
  • \((x_{i1},x_{i2},...,x_{ip})\) são as covariáveis.

A solução proposta por Cox consiste em condicionar a construção da função de verossimilhança ao conhecimento da história passada de falhas e censuras para eliminar a função de risco básica. Este método é chamado de método de máxima verossimilhança parcial. Partindo do pressuposto de proporcionalidade, é possível estimar os efeitos das covariáveis sem ter que fazer qualquer suposição a respeito da distribuição do tempo de vida.

Considere que em uma amostra de \(n\) indivíduos existam \(k\leq n\) falhas distintas nos tempos \(t_1, t_2, ..., t_k\). A ideia básica deste método é considerar a probabilidade condicional da \(i\)-ésima observação vir a falhar no tempo \(t_i\) conhecendo quais observações estão sob risco em \(t_i\)}. Esta probabilidade condicional, que é a razão entre o risco do indivíduo falhar em \(t_i\) e a soma dos riscos de falha de todos os indivíduos em risco, é a contribuição de cada indivíduo no tempo de falha \(t_i\).

Então a verossimilhança individual \(L_i\) será:

\[L_i=\frac{\lambda_i(t|\pmb{x_i})}{\sum_{j\in R(t_i)}\lambda_j(t|\pmb{x_j})}=\frac{\lambda_0(t)exp\{\pmb{\beta}\pmb{x_i}\}}{\sum_{j\in R(t_i)}\lambda_0(t)exp\{\pmb{\beta}\pmb{x_j}\}}=\frac{exp\{\pmb{\beta}\pmb{x_i}\}}{\sum_{j\in R(t_i)}exp\{\pmb{\beta}\pmb{x_j}\}}\]

em que \(R(t_i)\) é o conjunto dos índices das observações sob risco no tempo \(t_i\).

Assim, condicional a história de falhas e censuras até o tempo \(t_i\), o componente não paramétrico desaparece da expressão de verossimilhança. A função de verossimilhança parcial é dada por:

\[L(\pmb{\beta})=\prod_{i=1}^k\frac{exp\{\pmb{\beta}\pmb{x_i}\}}{\sum_{j\in R(t_i)}exp\{\pmb{\beta}\pmb{x_j}\}}=\prod_{i=1}^n\left[\frac{exp\{\pmb{\beta}\pmb{x_i}\}}{\sum_{j\in R(t_i)}exp\{\pmb{\beta}\pmb{x_j}\}}\right]^{\delta_i}\]

Os estimadores de máxima verossimilhança de \(\pmb{\beta}\) são obtidos a partir da verossimilhança parcial \(L(\pmb{\beta})\).

Alguns autores provaram que os estimadores de máxima verossimilhança para o modelo de Cox são consistentes e assintoticamente normais sob certas condições de regularidade. No R o ajuste do modelo de regressão de COX será feito a partir da função coxph do pacote survival. A sintaxe é muito parecida com outros modelos que ja ajustamos. Como temos duas covariáveis, precisamos fazer uma seleção das covariáveis para o modelo final. Inicialmente, vamos ajustar todos os modelos possíveis:

mod1 <- coxph(Surv(tempos,cens)~idade, data=dados,method="breslow")
mod2 <- coxph(Surv(tempos,cens)~estagio, data=dados, method="breslow")
mod3 <- coxph(Surv(tempos,cens)~idade+estagio, data=dados, method="breslow")
mod4 <- coxph(Surv(tempos,cens)~idade+estagio+idade*estagio, data=dados, method="breslow")

Para escolher qual destes é o melhor modelo, vamos usar como critério o AIC e BIC:

AIC<-c(AIC(mod1),AIC(mod2),AIC(mod3),AIC(mod4)) 
BIC<-c(BIC(mod1),BIC(mod2),BIC(mod3),BIC(mod4)) 
AIC_BIC<-cbind(AIC,BIC)
rownames(AIC_BIC)<-c("Idade","Estágio","Idade + Estágio", "Idade + Estágio + Interação")
colnames(AIC_BIC)<-c("AIC","BIC")
AIC_BIC
##                                  AIC      BIC
## Idade                       393.8118 395.7238
## Estágio                     384.1625 389.8986
## Idade + Estágio             384.3589 392.0070
## Idade + Estágio + Interação 384.1550 397.5391

Após uma análise dos critérios, eles apontam para os modelos mod2 e mod4. Recomenda-se darmos sequência a ambos os modelos para fazermos, futuramente, a escolha entre o melhor deles. Os comandos a seguir, trazem um resumo de ambos os modelos:

mod2
## Call:
## coxph(formula = Surv(tempos, cens) ~ estagio, data = dados, method = "breslow")
## 
##               coef exp(coef) se(coef)     z        p
## estagioII  0.06576   1.06797  0.45844 0.143   0.8859
## estagioIII 0.61206   1.84423  0.35520 1.723   0.0849
## estagioIV  1.72284   5.60040  0.41966 4.105 4.04e-05
## 
## Likelihood ratio test=16.26  on 3 df, p=0.001001
## n= 90, number of events= 50
mod4
## Call:
## coxph(formula = Surv(tempos, cens) ~ idade + estagio + idade * 
##     estagio, data = dados, method = "breslow")
## 
##                       coef exp(coef)  se(coef)      z      p
## idade            -0.002559  0.997444  0.026051 -0.098 0.9218
## estagioII        -7.946142  0.000354  3.678209 -2.160 0.0307
## estagioIII       -0.122500  0.884706  2.468331 -0.050 0.9604
## estagioIV         0.846986  2.332605  2.425717  0.349 0.7270
## idade:estagioII   0.120254  1.127783  0.052307  2.299 0.0215
## idade:estagioIII  0.011351  1.011416  0.037449  0.303 0.7618
## idade:estagioIV   0.013673  1.013767  0.035967  0.380 0.7038
## 
## Likelihood ratio test=24.27  on 7 df, p=0.001021
## n= 90, number of events= 50

As saídas representam:

  • coef: Estimativa dos parâmetros \(\pmb{\beta}\).
  • exp(coef): Exponencial das estimativas dos parâmetros \(\pmb{\beta}\) (útil para interpretação).
  • se(coef): Erro padrão associado às estimativas dos parâmetros \(\pmb{\beta}\).
  • z: Estatística de teste de Wald.
  • p: P-Valor associado à estatística de teste de Wald.

Observe como não há estimativa para o intercepto. A função de risco de base absorve esse componente.

5 ADEQUAÇÃO DO AJUSTE

Considere que se o \(i\)-ésimo indivíduo com vetor de covariáveis \(\pmb{x_i}=(x_{1i},...,x_{pi})^t\) é observado falhar. Tem-se para este indivíduo um vetor de resíduos de Schoenfeld \(\pmb{r_i} = (r_{i1},...r_{ip})^t\) dado por:

\[r_{i1}=x_{iq}-\frac{\sum_{j \in R(t_i)}x_{jq}exp\{\hat{\pmb{\beta}}\pmb{x_j}\}}{\sum_{j \in R(t_i)}exp\{\hat{\pmb{\beta}}\pmb{x_j}\}}\]

Estes resíduos são interpretados como a diferença entre os valores observados de covariáveis de um indivíduo com tempo de ocorrência do evento \(t_i\) e os valores esperados em \(t_i\) dado o grupo de risco \(R(t_i)\)}.

Estes resíduos são \textbf{definidos apenas nos tempos de ocorrência do evento de interesse. O número de vetores de resíduos é igual ao número de covariáveis ajustadas no modelo. Dessa forma, através do gráfico dos resíduos padronizados de Schoenfeld contra o tempo, é possível verificar a existência ou não de proporcionalidade. Isto é, se a suposição de riscos proporcionais for satisfeita não deverá existir nenhuma tendência sistemática no gráfico. Para auxiliar na detecção de uma possível falha na suposição de proporcionalidade, uma curva suavizada com bandas de confiança é adicionada neste gráfico.

Temos dois modelos candidatos: mod2 somente com a variável estagio e mod4 com as variáveis idade, estagio e a interação.

Vamos, individualmente para cada modelo, obter os gráficos dos resíduos de Schoenfeld e os testes de correlação para estes resíduos. Os resíduos de Schoenfeld são obtidos a partir da função resid e a avaliação da proporcionalidade a partir dos gráficos é feito com o comando cox.zph:

### MODELO COM ESTAGIO
Schoenfeld1<-resid(mod2,type="scaledsch")
par(mfrow=c(1,3))
plot(cox.zph(mod2, transform="identity",terms=F))

### MODELO COM ESTAGIO + IDADE + INTERAÇÃO
Schoenfeld2<-resid(mod4,type="scaledsch")
par(mfrow=c(2,4))
plot(cox.zph(mod4, transform="identity",terms=F))

A obtenção de medidas estatísticas, bem como a realização de testes de hipóteses são, desse modo, de grande utilidade nessas situações. O coeficientes de correlação de Pearson \((\rho)\) entre os resíduos de Schoenfeld e \(t\), para cada covariável, é uma dessas medidas. Valores de \(\rho\) próximos de 0 mostram não haver evidências para a rejeição da suposição de riscos proporcionais.

Os comandos a seguir fazem os testes de correlação nos resíduos contra o tempo para ambos os modelos:

cox.zph(mod2, transform="identity")
##         chisq df    p
## estagio   3.1  3 0.38
## GLOBAL    3.1  3 0.38
cox.zph(mod4, transform="identity")
##               chisq df    p
## idade          1.34  1 0.25
## estagio        4.75  3 0.19
## idade:estagio  4.48  3 0.21
## GLOBAL         6.73  7 0.46

Note, das análises anteriores, que, embora não significativo ao nível de 5%, o estágio III é marginalmente significativo no mod2 (p-valor = 0,074). Comparando-se então, os resultados de ambos os diagnósticos para os dois modelos, decidiu-se pelo uso do modelo de Cox com a presença da interação.

6 INTERPRETAÇÃO DOS COEFICIENTES

A interpretação dos coeficientes no Modelo de Cox é feita considerando a razão de riscos. Suponha que a razão entre as funções de risco para os indivíduos distintos \(i\) e \(j\) seja dada por:

\[\frac{\lambda(t|\pmb{x_i})}{\lambda(t|\pmb{x_j})}=\frac{\lambda_0(t)exp\{\pmb{\beta}\pmb{x_i}\}}{\lambda_0(t)exp\{\pmb{\beta}\pmb{x_j}\}}=exp\{\pmb{\beta}\pmb{x_i}-\pmb{\beta}\pmb{x_j}\},\]

que não depende do tempo.

De acordo com o modelo escolhido, temos a seguinte constante de proporcionalidade:

\[\exp\{\hat{\beta_1}idade+\hat{\beta_2}estagioII+\hat{\beta_3}estagioIII+\hat{\beta_4}estagioIV+\\\hat{\beta_5}Idade*estagioII+\hat{\beta_6}Idade*estagioIII+\hat{\beta_7}Idade*estagioIV\}\}\]

Por exemplo, para os pacientes \(i\) e \(j\), em que ambos encontram-se no estágio II da doença, mas um deles apresenta idade de 65 anos e o outro apresenta uma idade de 50 anos, tem-se que a razão dos riscos entre eles é:

\[\frac{\hat{\lambda}(t|\pmb{x_i})}{\hat{\lambda}(t|\pmb{x_j})}=\frac{exp\{\hat{\beta_1}*65+\hat{\beta_2}+\hat{\beta_5}*65\}}{exp\{\hat{\beta_1}*50+\hat{\beta_2}+\hat{\beta_5}*50\}}\\=\frac{exp\{-0.002*65-7.946+0.120*65\}}{exp\{-0.002*50-7.946+0.120*50\}}=5.87\]

Assim, o risco de morte de um paciente com 65 anos de idade e no estágio II da doença é aproximadamente 6 vezes o risco de morte de um paciente com 50 anos e no mesmo estágio.

Por outro lado, considere os indivíduos \(k\) e \(l\), em que ambos tem 50 anos de idade, sendo que um deles se concontra no estágio IV da doença e o outro no estágio III. A razão dos riscos de morte é:

\[\frac{\hat{\lambda}(t|\pmb{x_k})}{\hat{\lambda}(t|\pmb{x_l})}=\frac{exp\{\hat{\beta_1}*50+\hat{\beta_4}+\hat{\beta_7}*50\}}{exp\{\hat{\beta_1}*50+\hat{\beta_3}+\hat{\beta_6}*50\}}\\=\frac{exp\{-0.002*50+0.847+0.014*50\}}{exp\{-0.002*50-0.123+0.011*50\}}=2.99\]

Assim, o risco de morte do paciente com 50 anos e no estágio IV da doença é de aproximadamente 3 vezes o risco de morte do paciente com 50 anos de idade e no estágio III da doença.

7 ESTIMAÇÃO DE FUNÇÕES RELACIONADAS A \(\lambda_0(t)\)

Os coeficientes de regressão \(\pmb{\beta}\) são as quantidades de maior interesse na modelagem estatística dos dados. Entretanto, funções relacionadas a \(\lambda_0(t)\) são também importantes no modelo de Cox. Estas funções se referem basicamente à função de risco acumulada de base \[\Lambda_0(t)=\int_0^t\lambda_0(u)du\] e à correspondentes função de sobrevivência de base: \[S_0(t)=exp\{-\Lambda_0(t)\}\]

Por exemplo, a função de sobrevivência \[S(t|x)=S_0(t)^{exp\{\pmb{\beta}\pmb{x}\}}\] é extremamente útil quando se deseja concluir a análise em termos de percentis associados a grupos de indivíduos.

Uma estimativa simples para \(\Lambda_0(t)\) é uma função escada com saltos nos tempos distintos de falha e expressa por:

\[\hat{\Lambda}_0(t)\sum_{j:t_j<t}\frac{d_j}{\sum_{l\in R_j}exp\{\hat{\pmb{\beta x_l}}\}}\]

em que \(d_j\) é o número de falhas em \(t_j\).

O estimador \(\hat{\Lambda}_0(t)\) é uma versão do estimador de Nelson-Aalen em um contexto com covariáveis.

Consequentemente, as funções de sobrevivência \(S_0(t)\) e \(S(t|\pmb{x})\) podem ser estimadas por, respectivamente,

\[\hat{S}_0(t)=exp\{-\hat{\Lambda}_0(t)\}\]

e

\[\hat{S}(t|\pmb{x})=\left[\hat{S}_0(t)\right]^{exp\{\hat{\pmb{\beta x}}\}}\]

Nas figuras abaixo estão representadas as curvas de sobrevivência estimadas para pacientes com idades de 50 e 65 anos, em cada um dos 4 estágios da doença. Desta figura, pode-se observar que as curvas de sobrevivência estimadas para os estágios I, III e IV nio apresentam diferenças muito acentuadas, quando comparadas entre as idades de 50 e 65 anos. No estágio II, contudo, observa-se um decréscimo expressivo desta curva para os pacientes de 65 anos de idnde, quando comparados aos de 50 anos. Este fato justifica, assim, a presença da interação entre o estágio e idade no modelo, em especial entre idade e o estágio II.

plot(survfit(mod4,newdata=list(idade=50,estagio="I")),
     conf.int = F,col="blue",xlab="Tempos (meses)", ylab="Sobrevivência Estimada",main="50 anos")
lines(survfit(mod4,newdata=list(idade=50,estagio="II")),
      conf.int = F,col="red")
lines(survfit(mod4,newdata=list(idade=50,estagio="III")),
      conf.int = F,col="green")
lines(survfit(mod4,newdata=list(idade=50,estagio="IV")),
      conf.int = F,col="orange")
legend(0,0.35,legend=c("Estagio I", "Estagio II", "Estagio III", "Estagio IV"),
       col=c("blue", "red", "green", "orange"),lty=1)

plot(survfit(mod4,newdata=list(idade=50,estagio="I")),
     conf.int = F,col="blue",xlab="Tempos (meses)", ylab="Sobrevivência Estimada",main="65 anos")
lines(survfit(mod4,newdata=list(idade=65,estagio="II")),
      conf.int = F,col="red")
lines(survfit(mod4,newdata=list(idade=65,estagio="III")),
      conf.int = F,col="green")
lines(survfit(mod4,newdata=list(idade=65,estagio="IV")),
      conf.int = F,col="orange")
legend(0,0.35,legend=c("Estagio I", "Estagio II", "Estagio III", "Estagio IV"),
       col=c("blue", "red", "green", "orange"),lty=1)

8 EXERCÍCIO

1- Para os dados de câncer de ovário:

  1. Ajuste um modelo de Cox a estes dados fazendo seleção de covariáveis.

  2. No modelo ajustado em a):

  • Faça a análise de resíduos
  • Faça a interpretação dos coeficientes
  • Estime a sobrevivência para indivíduos distintos de acordo com a(s) covariável(is) do modelo.