Modelo de Cox

Análise de Sobrevivência

19 de maio de 2018

Leandro Valter
leandro.vvalter@gmail.com
Tiago Almeida de Oliveira
tiagoestatistico@gmail.com
Universidade Estadual da Paraíba
Centro de Ciência e Tecnologia
Departamento de Estatística

Algumas questões dos livros: Análise de Sobrevivência: teoria e aplicações em saúde da e Análise de sobrevivência aplicada.

Carregando pacotes necessários para análise.

library(survival)
library(dplyr)
library(ggplot2)
library(lme4)
library(tibble)
library(ggExtra)
library(ggfortify)
library(survminer)

Colosimo

1.(5.9)

a)

setwd("D:/UEPB/Sobrevivência")
dados=read.table("ovario.txt",header=T)
attach(dados)
head(dados)
##   tempos censura tumor
## 1     28       1     0
## 2     89       1     0
## 3    175       1     0
## 4    195       1     0
## 5    309       1     0
## 6    377       0     0
modelo1<-coxph(Surv(tempos,censura)~factor(tumor), x = T, method="breslow", dados)
summary(modelo1)
## Call:
## coxph(formula = Surv(tempos, censura) ~ factor(tumor), data = dados, 
##     x = T, method = "breslow")
## 
##   n= 35, number of events= 22 
## 
##                 coef exp(coef) se(coef)     z Pr(>|z|)  
## factor(tumor)1 1.119     3.061    0.497 2.251   0.0244 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## factor(tumor)1     3.061     0.3267     1.155     8.107
## 
## Concordance= 0.588  (se = 0.059 )
## Rsquare= 0.151   (max possible= 0.978 )
## Likelihood ratio test= 5.73  on 1 df,   p=0.01672
## Wald test            = 5.07  on 1 df,   p=0.02441
## Score (logrank) test = 5.51  on 1 df,   p=0.01889
anova(modelo1)
## Analysis of Deviance Table
##  Cox model: response is Surv(tempos, censura)
## Terms added sequentially (first to last)
## 
##                loglik  Chisq Df Pr(>|Chi|)  
## NULL          -67.168                       
## factor(tumor) -64.305 5.7258  1    0.01672 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Com um p-valor significativo (0,0167), o modelo foi melhor que o modelo nulo.

logrank = 5,51

FioCruz

6.2)

  • Modelo 1 : idade
  • Modelo 2 : idade + cdiab + congenita + crim
  • Modelo 3 : idade + cdiab + congenita + crim + grande
setwd("D:/UEPB/Sobrevivência")
dialise=read.table("dialise.csv",header=T,sep=",")
attach(dialise)

Kaplan-Meier para covariável diabetes.

diab=survfit(Surv(tempo,status) ~ cdiab, data = dialise)
autoplot(diab,conf.int = F)+
  labs(x="Tempo de Sobrevivência (Dias)",y="Probabilidade de Sobrevivência",
    title="Kaplan-Meier DIABETES",colour="Diabetes")+theme_bw()

Kaplan-Meier para covariável congenita.

congenita=survfit(Surv(tempo,status) ~ congenita, data = dialise)
autoplot(congenita,conf.int = F)+
  labs(x="Tempo de Sobrevivência (Dias)",y="Probabilidade de Sobrevivência",
    title="Kaplan-Meier Congenita",colour="Congenita")+theme_bw()

Kaplan-Meier para covariável crim.

crim=survfit(Surv(tempo,status) ~ crim, data = dialise)
autoplot(crim,conf.int = F)+
  labs(x="Tempo de Sobrevivência (Dias)",y="Probabilidade de Sobrevivência",
    title="Kaplan-Meier Causa renal",colour="crim")+theme_bw()

Kaplan-Meier para covariável grande.

grande=survfit(Surv(tempo,status) ~ grande, data = dialise)
autoplot(grande,conf.int = F)+
  labs(x="Tempo de Sobrevivência (Dias)",y="Probabilidade de Sobrevivência",
    title="Kaplan-Meier Unidade de tratamento",colour="grande")+theme_bw()

A partir da análise gráfica, pode-se dizer que existe evidência de proporcionalidade. A covariável congênita apresenta um padrão diferente, mas ainda mantém a tendência.

b)Ajuste cada modelo explicativo acima utilizando o modelo de riscos proporcionais de Cox, tomando o cuidado de interpretar os parâmetros a cada modelo.

modelo1=coxph(Surv(tempo,status)~idade, data=dialise, x=T)
summary(modelo1)
## Call:
## coxph(formula = Surv(tempo, status) ~ idade, data = dialise, 
##     x = T)
## 
##   n= 6805, number of events= 1603 
## 
##           coef exp(coef) se(coef)     z Pr(>|z|)    
## idade 0.034903  1.035520 0.001727 20.22   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## idade     1.036     0.9657     1.032     1.039
## 
## Concordance= 0.644  (se = 0.008 )
## Rsquare= 0.062   (max possible= 0.98 )
## Likelihood ratio test= 435.1  on 1 df,   p=0
## Wald test            = 408.7  on 1 df,   p=0
## Score (logrank) test = 415.2  on 1 df,   p=0

No modelo 1, a idade é um fator de risco de 1,035, ou seja, cada ano de idade a mais no início do tratamento aumenta o risco de óbito em 3,5%.

modelo2=coxph(Surv(tempo,status)~idade + factor(cdiab) + factor(crim) + factor(congenita), data = dialise, x = T)
summary(modelo2)
## Call:
## coxph(formula = Surv(tempo, status) ~ idade + factor(cdiab) + 
##     factor(crim) + factor(congenita), data = dialise, x = T)
## 
##   n= 6805, number of events= 1603 
## 
##                         coef exp(coef)  se(coef)      z Pr(>|z|)    
## idade               0.034448  1.035048  0.001751 19.677  < 2e-16 ***
## factor(cdiab)1      0.285606  1.330568  0.060308  4.736 2.18e-06 ***
## factor(crim)1       0.030331  1.030796  0.066975  0.453  0.65064    
## factor(congenita)1 -0.715243  0.489073  0.226181 -3.162  0.00157 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## idade                 1.0350     0.9661    1.0315    1.0386
## factor(cdiab)1        1.3306     0.7516    1.1822    1.4975
## factor(crim)1         1.0308     0.9701    0.9040    1.1754
## factor(congenita)1    0.4891     2.0447    0.3139    0.7619
## 
## Concordance= 0.649  (se = 0.008 )
## Rsquare= 0.067   (max possible= 0.98 )
## Likelihood ratio test= 473.3  on 4 df,   p=0
## Wald test            = 438.5  on 4 df,   p=0
## Score (logrank) test = 451.3  on 4 df,   p=0

No modelo 2, a idade sendo um fator de risco. A diabetes se mostrou um importante fator de risco, com um sobrerisco de 33%. A causa renal não foi significativa, não sendo considerada um fator de risco. A doença não sendo congênita se mostrou um fator risco, ou seja, têm risco duas vezes maior de ir a óbito do que as com causa congênita.

modelo3=coxph(Surv(tempo,status)~idade + factor(cdiab) + factor(crim) + factor(congenita)+factor(grande), data = dialise, x = T)
summary(modelo3)
## Call:
## coxph(formula = Surv(tempo, status) ~ idade + factor(cdiab) + 
##     factor(crim) + factor(congenita) + factor(grande), data = dialise, 
##     x = T)
## 
##   n= 6805, number of events= 1603 
## 
##                         coef exp(coef)  se(coef)      z Pr(>|z|)    
## idade               0.033995  1.034580  0.001754 19.386  < 2e-16 ***
## factor(cdiab)1      0.295484  1.343777  0.060398  4.892 9.97e-07 ***
## factor(crim)1       0.019132  1.019316  0.067083  0.285  0.77550    
## factor(congenita)1 -0.727439  0.483145  0.226222 -3.216  0.00130 ** 
## factor(grande)1     0.176994  1.193624  0.063390  2.792  0.00524 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                    exp(coef) exp(-coef) lower .95 upper .95
## idade                 1.0346     0.9666    1.0310    1.0381
## factor(cdiab)1        1.3438     0.7442    1.1938    1.5126
## factor(crim)1         1.0193     0.9811    0.8937    1.1625
## factor(congenita)1    0.4831     2.0698    0.3101    0.7527
## factor(grande)1       1.1936     0.8378    1.0542    1.3515
## 
## Concordance= 0.652  (se = 0.008 )
## Rsquare= 0.068   (max possible= 0.98 )
## Likelihood ratio test= 481.4  on 5 df,   p=0
## Wald test            = 447.4  on 5 df,   p=0
## Score (logrank) test = 460.8  on 5 df,   p=0

No modelo 3, a inclusão da covariável grande manteve a significância das outras covariáveis. O tamanho da unidade se mostrou significativo, os pacientes atendidos em unidades grandes tem o risco aumentado em 19% de ir a óbito dos que são atendidos em unidades menores.

c)Compare os modelos usando o teste da razão de verossimilhanças e escolha o modelo com melhor ajuste.

anova(modelo1, modelo2, modelo3, test = "Chisq")
## Analysis of Deviance Table
##  Cox model: response is  Surv(tempo, status)
##  Model 1: ~ idade
##  Model 2: ~ idade + factor(cdiab) + factor(crim) + factor(congenita)
##  Model 3: ~ idade + factor(cdiab) + factor(crim) + factor(congenita) + factor(grande)
##   loglik   Chisq Df P(>|Chi|)    
## 1 -13037                         
## 2 -13017 38.2457  3 2.507e-08 ***
## 3 -13013  8.0653  1  0.004512 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Os modelos 2 e 3, foram bem ajustados, com a inclusão da covariável grande, o modelo 3 foi o melhor

d)Qual o poder explicativo do modelo escolhido? Calcule a razão entre o R2 do modelo escolhido e o R2 máximo.

O modelo 3 tem o poder explicativo de 6,8%, apresentando um concordância de 65,2%, ou seja, \(R^2/R^2_{máx}=0,068/0.98=0.652\).

e)Faça o gráfico dos índices de prognóstico.

6.3) Faça uma análise do banco de dados ipec.csv utilizando o modelo de Cox, considerando três modelos explicativos aninhados.

  • Modelo 1 : idade + sexo
  • Modelo 2 : idade + sexo + acompanhamento
  • Modelo 3 : idade + sexo + acompanhamento + tratamento

Carregando o banco de dados.

setwd("D:/UEPB/Sobrevivência")
ipec=read.table("ipec.txt",h=T)
attach(ipec)

a)Faça gráficos de Kaplan-Meier estratificados por cada variável categórica para verificar o pressuposto do modelo de Cox.

Kaplan-Meier para covariável sexo.

sexo =survfit(Surv(tempo, status) ~ sexo, data = ipec)
autoplot(sexo,conf.int = F)+
  labs(x="Tempo de Sobrevivência (Dias)",y="Probabilidade de Sobrevivência",
    title="Kaplan-Meier SEXO",colour="Sexo")+theme_bw()

Kaplan-Meier para covariável escolaridade.

esc = survfit(Surv(tempo, status) ~ escola, data = ipec)
autoplot(esc,conf.int = F)+
  labs(x="Tempo de Sobrevivência (Dias)",y="Probabilidade de Sobrevivência",
       title="Kaplan-Meier ESCOLARIDADE",colour="Escolaridade")+theme_bw()

Kaplan-Meier para covariável acompanhamento.

acompan = survfit(Surv(tempo, status) ~ acompan, data = ipec)
autoplot(acompan,conf.int = F)+
  labs(x="Tempo de Sobrevivência (Dias)",y="Probabilidade de Sobrevivência",
       title="Kaplan-Meier ACOMPANHAMENTO",colour="Acompanhamento")+theme_bw()

Kaplan-Meier para covariável tratamento.

trat = survfit(Surv(tempo, status) ~ tratam, data = ipec)
autoplot(trat,conf.int = F)+
  labs(x="Tempo de Sobrevivência (Dias)",y="Probabilidade de Sobrevivência",
       title="Kaplan-Meier TRATAMENTO",colour="Tratamento")+theme_bw()

b)Ajuste cada modelo acima utilizando o modelo de riscos proporcionais de Cox, tomando o cuidado de interpretar os parâmetros de cada modelo.

Ajuste do modelo de Cox tendo idade e sexo como covariáveis.

modelo1 = coxph(Surv(tempo, status) ~ idade + sexo, data = ipec, x = T)
summary(modelo1)
## Call:
## coxph(formula = Surv(tempo, status) ~ idade + sexo, data = ipec, 
##     x = T)
## 
##   n= 193, number of events= 90 
## 
##           coef exp(coef) se(coef)      z Pr(>|z|)  
## idade -0.01274   0.98734  0.01162 -1.096    0.273  
## sexoM  0.55622   1.74407  0.27611  2.015    0.044 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##       exp(coef) exp(-coef) lower .95 upper .95
## idade    0.9873     1.0128    0.9651     1.010
## sexoM    1.7441     0.5734    1.0152     2.996
## 
## Concordance= 0.556  (se = 0.033 )
## Rsquare= 0.029   (max possible= 0.988 )
## Likelihood ratio test= 5.64  on 2 df,   p=0.05965
## Wald test            = 5.09  on 2 df,   p=0.07834
## Score (logrank) test = 5.18  on 2 df,   p=0.07485

No modelo 1, o sexo masculino foi significativo, sendo assim um fator de risco, com uma razão de risco de 1,74%, ou seja, pelo fato do paciente ser do sexo masculino tem quase duas vezes mais chance de ocorrência do evento. Porém, o modelo 1 aprensenta um \(R^{2}\) de 0,2%, indicando um ajuste ruim, já que tem um poder explicativo muito baixo. Além disso, o teste de Wald teve um p-valor de 0,078, denunciando que o modelo 1 não é melhor que o modelo nulo.

Ajuste do modelo de Cox tendo idade, sexo e acompanhamento como covariáveis.

modelo2 = coxph(Surv(tempo, status) ~ idade + sexo + factor(acompan), data = ipec, x= T)
summary(modelo2)
## Call:
## coxph(formula = Surv(tempo, status) ~ idade + sexo + factor(acompan), 
##     data = ipec, x = T)
## 
##   n= 193, number of events= 90 
## 
##                       coef exp(coef)  se(coef)      z Pr(>|z|)    
## idade            -0.001659  0.998342  0.012013 -0.138    0.890    
## sexoM             0.272185  1.312830  0.281849  0.966    0.334    
## factor(acompan)1  1.707324  5.514184  0.401021  4.257 2.07e-05 ***
## factor(acompan)2  2.517632 12.399195  0.445248  5.654 1.56e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## idade               0.9983    1.00166    0.9751     1.022
## sexoM               1.3128    0.76171    0.7556     2.281
## factor(acompan)1    5.5142    0.18135    2.5126    12.101
## factor(acompan)2   12.3992    0.08065    5.1808    29.675
## 
## Concordance= 0.711  (se = 0.033 )
## Rsquare= 0.231   (max possible= 0.988 )
## Likelihood ratio test= 50.74  on 4 df,   p=2.535e-10
## Wald test            = 37.01  on 4 df,   p=1.795e-07
## Score (logrank) test = 49  on 4 df,   p=5.844e-10

O modelo 2, tem um melhor ajuste que o modelo 1 devido a inclusão da covariável acompanhamento, com um \(R^{2}\) de 23%, juntamente com o teste de Wald, com um p-valor significativo, indicando que o modelo 2 também é melhor que o modelo nulo. A covariável acompanhamento foi significativa, sendo um importante fator de risco, sinalizando que os pacientes que tiveram um acompanhamento tiveram 5,5 vezes menos chance de ocorrência do evento caso exista um internação posterior e 12,39 vezes menos chance de ocorrência do evento caso seja internado imediatamente.

Ajuste do modelo de Cox tendo idade, sexo, acompanhamento e tratamento como covariáveis.

modelo3 = coxph(Surv(tempo, status) ~ idade + sexo + factor(acompan) + factor(tratam), data = ipec, x = T)
summary(modelo3)
## Call:
## coxph(formula = Surv(tempo, status) ~ idade + sexo + factor(acompan) + 
##     factor(tratam), data = ipec, x = T)
## 
##   n= 193, number of events= 90 
## 
##                       coef exp(coef)  se(coef)      z Pr(>|z|)    
## idade             0.001428  1.001429  0.012120  0.118  0.90621    
## sexoM             0.074236  1.077061  0.285759  0.260  0.79503    
## factor(acompan)1  1.676179  5.345096  0.408363  4.105 4.05e-05 ***
## factor(acompan)2  2.153003  8.610675  0.467205  4.608 4.06e-06 ***
## factor(tratam)1  -1.241920  0.288829  0.301134 -4.124 3.72e-05 ***
## factor(tratam)2  -2.096740  0.122856  0.470499 -4.456 8.33e-06 ***
## factor(tratam)3  -2.945019  0.052601  1.018814 -2.891  0.00384 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                  exp(coef) exp(-coef) lower .95 upper .95
## idade               1.0014     0.9986  0.977921    1.0255
## sexoM               1.0771     0.9285  0.615179    1.8857
## factor(acompan)1    5.3451     0.1871  2.400800   11.9002
## factor(acompan)2    8.6107     0.1161  3.446289   21.5141
## factor(tratam)1     0.2888     3.4623  0.160072    0.5212
## factor(tratam)2     0.1229     8.1396  0.048855    0.3089
## factor(tratam)3     0.0526    19.0110  0.007141    0.3874
## 
## Concordance= 0.789  (se = 0.033 )
## Rsquare= 0.372   (max possible= 0.988 )
## Likelihood ratio test= 89.84  on 7 df,   p=1.11e-16
## Wald test            = 76.87  on 7 df,   p=5.973e-14
## Score (logrank) test = 100.1  on 7 df,   p=0

O modelo 3, com a inclusão da covariável tratamento, teve o melhor ajuste, com um \(R^{2}\) de 37,2%, tendo um bom poder explicativo, com o teste de Wald indicando que o modelo 3 é melhor que o modelo nulo (p-valor=0). As covariáveis acompanhamento e tratamento foram significativas, sendo importantes fatores de risco, sinalizando que os pacientes que tiveram um acompanhamento com internação posterior tem 5,3 vezes menos chance de ocorrência do evento e pacientes com internação imediata tem 8,6 vezes menos chando de ocorrência do evento. Já os pacientes que tiveram uma terapia antirretroviral do tipo mono tem 4,46 vezes menos chance de ocorrência de evento, já do tipo combinada tem 8,1 vezes menos chance de ocorrência e os pacientes que passaram pelo tratamento antirretroviral do tipo potente tem 19 vezes menos chance de ocrrência do evento.

c)Compare os modelos usando o teste da razão de verossimilhanças e o gráfico dos índices de prognóstico.

anova(modelo1, modelo2, modelo3, test = "Chisq")
## Analysis of Deviance Table
##  Cox model: response is  Surv(tempo, status)
##  Model 1: ~ idade + sexo
##  Model 2: ~ idade + sexo + factor(acompan)
##  Model 3: ~ idade + sexo + factor(acompan) + factor(tratam)
##    loglik  Chisq Df P(>|Chi|)    
## 1 -422.30                        
## 2 -399.75 45.097  2 1.612e-10 ***
## 3 -380.20 39.106  3 1.648e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

d) Qual o poder explicativo do modelo escolhido? Calcule a razão entre o \(R^{2}\) do modelo escolhido e o \(R^{2}\) máximo.

O modelo 3 tem o poder explicativo de 37,2%, apresentando um concordância de 78,8%, ou seja, \(R^2/R^2_{máx}=0,357/0.988=0.788\).