Modelo de Regressão de Cox

Pacotes Necessários

library(survival)
library(survminer)
library(ggplot2)
library(ggExtra)
library(ggfortify)

Os dados a seguir , disponibilizados no site da FioCruz, temos propostas de modelos para pacientes em diálise, onde é removido os resíduos e os líquidos do seu corpo que os seus rins não são capazes de remover. Assim, também se pretende manter o equilíbrio do seu corpo, através da correção dos níveis das várias substâncias tóxicas presentes no sangue.

Modelo 1 : idade
Modelo 2 : idade + cdiab + congenita + crim
Modelo 3 : idade + cdiab + congenita + crim + grande

setwd("~/Área de Trabalho/Rmardown/Sobrevivência")
dialise=read.table("dialise.csv",header=T,sep=",")
attach(dialise)

Primeiramente é feito o gráfico de Kaplan-Meier separados por cada covariável categórica. O objetivo do gráfico é verificar visualmente se há não proporcionalidade, que como o nome do modelo já deixa implícito, inviabiliza o ajuste do modelo de Cox.

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")

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")

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")

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. A seguir, temos os ajustes dos modelos propostos e suas interpretações.

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 )
## Likelihood ratio test= 435.1  on 1 df,   p=<2e-16
## Wald test            = 408.7  on 1 df,   p=<2e-16
## Score (logrank) test = 415.2  on 1 df,   p=<2e-16

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.007 )
## Likelihood ratio test= 473.3  on 4 df,   p=<2e-16
## Wald test            = 438.5  on 4 df,   p=<2e-16
## Score (logrank) test = 451.3  on 4 df,   p=<2e-16

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.007 )
## Likelihood ratio test= 481.4  on 5 df,   p=<2e-16
## Wald test            = 447.4  on 5 df,   p=<2e-16
## Score (logrank) test = 460.8  on 5 df,   p=<2e-16

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.

Comparando os modelos pelo teste da razão de verossimilhanças e escolhemos o modelo com melhor ajuste. Os modelos 2 e 3 foram bem ajustados, com a inclusão da covariável grande, o modelo 3 foi o melhor.

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

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\). A suposição de riscos proporcionais (PH) pode ser verificada usando testes estatísticos e diagnósticos gráficos com base nos resíduos de Schoenfeld escalados. Para testar a suposição de riscos proporcionais (PH),temos:

test.ph <- cox.zph(modelo3);test.ph

##                    chisq df      p
## idade              4.837  1 0.0278
## factor(cdiab)      4.161  1 0.0414
## factor(crim)       0.357  1 0.5503
## factor(congenita)  3.679  1 0.0551
## factor(grande)     7.187  1 0.0073
## GLOBAL            20.187  5 0.0012

O teste não é estatisticamente significativo para algumas das variáveis uma das covariáveis e o teste global também é estatisticamente significativo. Portanto, não podemos assumir os riscos proporcionais.

Em princípio, os resíduos de Schoenfeld são independentes do tempo. Um gráfico que mostra um padrão não aleatório em relação ao tempo é evidência de violação da suposição de PH. Na figura abaixo, a linha central é um ajuste de spline de suavização ao gráfico, com as linhas tracejadas representando um erro padrão +/- 2 em torno do ajuste.

testph<-cox.zph(modelo3)
ggcoxzph(testph)

Para testar observações influentes ou outliers, podemos visualizar os resíduos de desvio ou os valores dfbeta que fornecem uma verificação de observações influentes. Especificando o tipo de argumento = “dfbeta”, temos as mudanças estimadas nos coeficientes de regressão ao excluir cada observação por vez; da mesma forma, type = “dfbetas” produz as mudanças estimadas nos coeficientes divididas por seus erros padrão.

ggcoxdiagnostics(modelo3, type = "dfbeta",
                 linear.predictions = FALSE, ggtheme = theme_bw())

## `geom_smooth()` using formula 'y ~ x'

Os gráficos de índice acima mostram que a comparação das magnitudes dos maiores valores dfbeta com os coeficientes de regressão sugere que algumas observações é influente individualmente.

Também é possível verificar outliers visualizando os resíduos de desvio. O resíduo de desvio é uma transformação normalizada do resíduo de martingale. Esses resíduos devem ser distribuídos simetricamente aproximadamente em torno de zero com um desvio padrão de 1.

ggcoxdiagnostics(modelo3, type = "deviance",
                 linear.predictions = FALSE, ggtheme = theme_bw())

## `geom_smooth()` using formula 'y ~ x'

Como os resíduos não pertercem ao intervalo (-1;1) e tendo muitos valores outliers, tanto positivos quanto negativos, assim obtendo uma predição fraca pelo modelo.

Modelo de Regressão de Cox

Leandro Valter

Tiago Almeida de Oliveira

Modelo de Regressão de Cox

Pacotes Necessários