Com pacotes extras para fazer gráfico e o pacote survival para sobrevivência

require(rjags)
require(coda)
require(DT)
require(tidyverse)
require(car)
#install.packages("MCMCpack")
require(MCMCpack)
#install.packages("bayesplot")
require(bayesplot)
#install.packages("ggrepel")
require(grid)
require(survival)

#para fazer o ggsurvplot: gráficos de sobrevivência
#install.packages("survminer") 
require(survminer)

Aplicação de Bayesiana em Análise de sobrevivência

Análise de sobrevivência no R:
- material .
Artigos:
- artigo I .
- Artigo estudo sobre papas .
Exemplos introdutórios e materiais didáticos
- Material I
- Material II
Referência principal (modelo & conjunto de dados):
- Referência principal (modelo & conjunto de dados)
- aplicações em R

Exemplo 1: Dados de transplante renal

Tempo de sobrevida em meses (t) de 148 pacientes submetidos ao transplante renal;
- Os dados contém uma uma covariável, x: o número total de incompatibilidades de antígeno HLA-B ou DR entre doador e receptor;
- Os dados contém censura (variável “status”) - característica principal inerente aos dados de análise de sobrevivência: 1: quando o tempo t é observado, 0: quando o paciente é censurado no tempo t (quando ele abandona o estudo por exemplo);
- A censura é a direita;

dados=read.table("http://www.mas.ncl.ac.uk/~nmf16/teaching/mas8391/renal.txt",header=TRUE)
datatable(dados, cap="Tempos de vida de pacientes após transplante renal",options = list(pageLength = 5))

Estimador de Kaplan-Meier sem ou com extratos

O estimador de Kaplan-Meier (conhecido como estimador do limite-produto) é um estimador não paramétrico para a função de sobrevivência (ou confiabilidade): \[ S(t)=\prod_{t_i \leq t} \left(\frac{n_i-d_i}{n_i}\right) \] onde \(n_i\) é o número de indivíduos em risco no tempo \(t_i\), \(d_i\) é o número total de óbitos no tempo \(t_i\).
Este estimador é utilizado largamente;
Possui propriedades assinótica - converge para a distribuição normal quando \(n\) (tamanho amostral) cresce.
Observamos na análise sem extratos (sem considerar a variável X, que o gráfico de KM para a sobrevida destes pacientes se estabiliza em 0.75);
Na análise por extratos, verificamos que a covariável \(X:\) número de incompatibilidades parece diminuir a estimativa de sobrevida: quanto mais incompatibilidades, menor é a sobrevida do paciente ao longo do tempo.

fit = survfit(Surv(t, status) ~ 1, data = dados,type="kaplan-meier")
ggsurvplot(fit, censor.shape="|", censor.size = 4,main="KM")

fit2 = survfit(Surv(t, status) ~ x, data = dados,type="kaplan-meier")
ggsurvplot(fit2, censor.shape="|", censor.size = 4,main="KM por grupos")

Modelo paramétrico - distribuição exponencial para os tempos de vida

Presença de covariável: X utilizada como fator - filtra-se os conjuntos de dados e aplica-se o modelo de sobrevivência exponencial.
A função survreg implementa , ou modelo de tempo de falha acelerada, em que a covariável X é um fator que “acelera” o tempo até o óbito do indivíduo. Coeficientes negativos implicam que a sobrevida dos indivíduos com incompatibilidades é menor - as incompatibilidades são prejudiciais para a sobrevida dos indivíduos.
PENDENTE verificar pressupostos do modelo.

fit3=list()
lambda=c()
Meantime=c()
i=0
while(i<=4){
fit3[[i+1]]=survreg(Surv(t, status)~1, dist="exponential",data=dados[dados$x==i,])
lambda[i+1]=1/exp(fit3[[i+1]]$coefficients[1])
Meantime[i+1]=1/lambda[i+1]
i=i+1
}
lambda

## [1] 0.003715847 0.008820736 0.010601025 0.050626503 0.202347228

Meantime

## [1] 269.1177 113.3692  94.3305  19.7525   4.9420

Criando a função de sobrevivência - valores preditos pelo modelo

tempos=seq(0,50,1)
S_t=matrix(nrow=length(tempos),ncol=length(lambda))
S_t[,1]=exp(-lambda[1]*tempos)
S_t[,2]=exp(-lambda[2]*tempos)
S_t[,3]=exp(-lambda[3]*tempos)
S_t[,4]=exp(-lambda[4]*tempos)
S_t[,5]=exp(-lambda[5]*tempos)
S_t=data.frame(cbind(tempos,S_t))
names(S_t)=c("tempos","x_0","x_1","x_2","x_3","x_4")
datatable(S_t,options = list(pageLength = 5),cap="Função de sobrevivência para vários valores de x")

Representação Gráfica

colors <- c("x=0"="blue", "x=1"="red", "x=2"="green" ,"x=3" = "orange", "x=4" = "pink")

ggplot(S_t,aes(x=tempos,y=x_0,color="x=0"))+
    geom_line() +
geom_line(aes(x=tempos,y=x_1,color="x=1")) + 
geom_line(aes(x=tempos,y=x_2,color="x=2"))+
geom_line(aes(x=tempos,y=x_3,color="x=3"))+
geom_line(aes(x=tempos,y=x_4,color="x=4"))+
labs(x="t (em meses)",
         y="S(t)",
         title="Gráfico da função de sobrevivência - modelo exponencial",
         color = "Número de Incompatibilidades (x)") +
    scale_color_manual(labels=c(expression(x==0),expression(x==1),expression(x==2),expression(x==3),expression(x==4)),values = colors)

Abordagem Bayesiana

Entrada dos dados no RJAGS

t=dados$t
is.na(t)=dados$status==0
is.censored=1-dados$status
t.cen=dados$t+dados$status
dados_jags=list(n=nrow(dados),t=dados$t,is.censored=is.censored,t.cen=t.cen,x=dados$x)
dados_jags=list(n=nrow(dados),t=dados$t,is.censored=is.censored,t.cen=t.cen,x=dados$x,new.n=50,new.t=rep(seq(1,50,5),5),new.x=c(rep(0,10),rep(1,10),rep(2,10),rep(3,10),rep(4,10)))
dados_jags

## $n
## [1] 148
## 
## $t
##   [1]  0.035  0.068  0.100  0.101  0.167  0.168  0.197  0.213  0.233  0.234
##  [11]  0.508  0.508  0.533  0.633  0.767  0.768  0.770  1.066  1.267  1.300
##  [21]  1.600  1.639  1.803  1.867  2.180  2.667  2.967  3.328  3.393  3.700
##  [31]  3.803  4.311  4.867  5.180  6.233  6.367  6.600  6.600  7.180  7.667
##  [41]  7.733  7.800  7.933  7.967  8.016  8.300  8.410  8.607  8.667  8.800
##  [51]  9.100  9.233 10.541 10.607 10.633 10.667 10.869 11.067 11.180 11.443
##  [61] 12.213 12.508 12.533 13.467 13.800 14.267 14.475 14.500 15.213 15.333
##  [71] 15.525 15.533 15.541 15.934 16.200 16.300 16.344 16.600 16.700 16.933
##  [81] 17.033 17.067 17.475 17.667 17.700 17.967 18.115 18.115 18.933 18.934
##  [91] 19.508 19.574 19.733 20.148 20.180 20.900 21.167 21.233 21.600 22.100
## [101] 22.148 22.180 22.180 22.267 22.300 22.500 22.533 22.867 23.738 24.082
## [111] 24.180 24.705 25.705 25.213 29.705 30.443 31.667 31.934 32.180 32.367
## [121] 32.672 32.705 33.148 33.567 33.770 33.869 34.836 34.869 34.934 35.738
## [131] 36.180 36.213 39.410 39.433 39.672 40.001 41.733 41.734 42.311 42.869
## [141] 43.180 43.279 43.902 44.267 44.475 44.900 45.148 46.451
## 
## $is.censored
##   [1] 0 0 0 0 1 0 1 0 1 0 0 1 0 1 0 0 1 0 1 0 0 1 1 0 0 0 1 1 0 0 1 1 1 0 1 1 1
##  [38] 1 0 1 0 1 1 1 0 0 1 1 0 1 1 0 1 1 1 0 1 0 1 1 0 0 1 1 1 1 1 1 1 1 1 1 1 1
##  [75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $t.cen
##   [1]  1.035  1.068  1.100  1.101  0.167  1.168  0.197  1.213  0.233  1.234
##  [11]  1.508  0.508  1.533  0.633  1.767  1.768  0.770  2.066  1.267  2.300
##  [21]  2.600  1.639  1.803  2.867  3.180  3.667  2.967  3.328  4.393  4.700
##  [31]  3.803  4.311  4.867  6.180  6.233  6.367  6.600  6.600  8.180  7.667
##  [41]  8.733  7.800  7.933  7.967  9.016  9.300  8.410  8.607  9.667  8.800
##  [51]  9.100 10.233 10.541 10.607 10.633 11.667 10.869 12.067 11.180 11.443
##  [61] 13.213 13.508 12.533 13.467 13.800 14.267 14.475 14.500 15.213 15.333
##  [71] 15.525 15.533 15.541 15.934 16.200 16.300 16.344 16.600 16.700 16.933
##  [81] 17.033 17.067 17.475 17.667 17.700 17.967 18.115 18.115 18.933 18.934
##  [91] 19.508 19.574 19.733 20.148 20.180 21.900 21.167 21.233 21.600 22.100
## [101] 22.148 22.180 22.180 22.267 22.300 22.500 22.533 22.867 23.738 24.082
## [111] 24.180 24.705 25.705 25.213 29.705 30.443 31.667 31.934 32.180 32.367
## [121] 32.672 32.705 33.148 33.567 33.770 33.869 34.836 34.869 34.934 35.738
## [131] 36.180 36.213 39.410 39.433 39.672 40.001 41.733 41.734 42.311 42.869
## [141] 43.180 43.279 43.902 44.267 44.475 44.900 45.148 46.451
## 
## $x
##   [1] 3 0 0 1 4 2 1 1 1 2 0 2 3 0 3 4 0 4 2 3 1 2 2 4 3 4 1 2 3 4 3 1 0 1 2 2 1
##  [38] 0 3 1 1 2 1 1 2 1 0 1 1 1 0 1 2 3 1 2 3 2 0 0 1 3 2 0 2 0 4 1 1 0 1 2 1 0
##  [75] 1 0 1 0 1 3 3 0 1 1 1 1 2 1 0 1 2 3 0 2 0 2 0 0 3 1 2 0 0 0 2 1 1 1 1 1 0
## [112] 0 2 1 3 1 0 2 1 0 1 2 1 1 1 2 0 1 2 0 1 1 1 0 0 0 2 0 2 0 0 1 2 2 1 1 1 0
## 
## $new.n
## [1] 50
## 
## $new.t
##  [1]  1  6 11 16 21 26 31 36 41 46  1  6 11 16 21 26 31 36 41 46  1  6 11 16 21
## [26] 26 31 36 41 46  1  6 11 16 21 26 31 36 41 46  1  6 11 16 21 26 31 36 41 46
## 
## $new.x
##  [1] 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3
## [39] 3 3 4 4 4 4 4 4 4 4 4 4

Valors iniciais - Não preciso por valores iniciais - o algoritmo gera aleatoriamente, se necessário! - Também pode-se utilizar os valores iniciais sugeridos no material, mas não foi utilizado.

#tinits1<-dados$t+5
#is.na(tinits1)<-dados$status==1
#tinits2<-tinits1+5
#renalinits<-list(list(t=tinits1),list(t=tinits2))

Sintaxe do modelo

cat("
    model
    {
      # priors
      for(j in 1:5)
      {
        # prior lambda for each group
        lambda[j] ~ dgamma(0.001, 0.001)
        mu[j] <- 1/lambda[j] # mean time to death
      }
      # likelihood
      for(i in 1:n)
      {
        is.censored[i] ~ dinterval(t[i], t.cen[i])
        t[i] ~ dexp(lambda[x[i]+1])
      }
      # predictions
      for(i in 1:new.n)
      {
        S[i] = exp(-lambda[new.x[i]+1]*new.t[i])
      }
   }
   ", file="modelo6.txt")

Reconhece o modelo no arquivo .txt

modelo<-jags.model("modelo6.txt",data=dados_jags,n.chains=2)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 296
##    Unobserved stochastic nodes: 5
##    Total graph size: 811
## 
## Initializing model

Começa a rodar o modelo (update)

note que ainda não inicializamos os parâmetros
desta forma, as amostras iniciais serão descartadas (chamadas de burn-in), as amostras iniciais são dependentes dos valores iniciais.

update(modelo,5000)

Coleta as amostras desejadas:

params=c("lambda","mu")
samples=coda.samples(modelo,params,n.iter=50000,thin=5)

Traços e densidades a posteriori

mcmc_trace(samples)

mcmc_dens(samples)

Diagnóstico de Gelman Rubin

Satisfaz \(1 < valores de R <1.2\)

gelman.diag(samples, confidence = 0.95, transform=FALSE)

## Potential scale reduction factors:
## 
##           Point est. Upper C.I.
## lambda[1]          1          1
## lambda[2]          1          1
## lambda[3]          1          1
## lambda[4]          1          1
## lambda[5]          1          1
## mu[1]              1          1
## mu[2]              1          1
## mu[3]              1          1
## mu[4]              1          1
## mu[5]              1          1
## 
## Multivariate psrf
## 
## 1

gelman.plot(samples)

Gráficos de autocorrelação

Satisfaz autocorrelações tendem a zero a medida que os “lags” ou distancias crescem.

autocorr.plot(samples[[1]],col=1)

autocorr.plot(samples[[2]],col=2)

Sumários a posteriori

Pelos sumários a posteriori nós verificamos que os valores de \(\lambda\) para o modelo exponencial diferem de acordo com a covariável X: número de incompatibilidades;
O mesmo ocorre para a quantidade tempo médio de falha, dada pelo inverso do parâmetro \(\lambda\).

summary(samples)

## 
## Iterations = 5005:55000
## Thinning interval = 5 
## Number of chains = 2 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##               Mean       SD  Naive SE Time-series SE
## lambda[1]  0.04826 0.007703 5.447e-05      5.447e-05
## lambda[2]  0.05298 0.007182 5.079e-05      5.078e-05
## lambda[3]  0.05659 0.009959 7.042e-05      7.111e-05
## lambda[4]  0.10138 0.025398 1.796e-04      1.839e-04
## lambda[5]  0.28339 0.106291 7.516e-04      7.516e-04
## mu[1]     21.26481 3.491636 2.469e-02      2.469e-02
## mu[2]     19.22931 2.653796 1.877e-02      1.876e-02
## mu[3]     18.23309 3.306898 2.338e-02      2.361e-02
## mu[4]     10.52162 2.808602 1.986e-02      1.986e-02
## mu[5]      4.11805 1.884599 1.333e-02      1.333e-02
## 
## 2. Quantiles for each variable:
## 
##               2.5%      25%      50%      75%    97.5%
## lambda[1]  0.03427  0.04289  0.04784  0.05321  0.06434
## lambda[2]  0.03983  0.04798  0.05260  0.05763  0.06789
## lambda[3]  0.03882  0.04964  0.05599  0.06290  0.07791
## lambda[4]  0.05811  0.08335  0.09910  0.11699  0.15728
## lambda[5]  0.11328  0.20657  0.27058  0.34603  0.52563
## mu[1]     15.54125 18.79251 20.90490 23.31599 29.17816
## mu[2]     14.72974 17.35170 19.01075 20.84132 25.10447
## mu[3]     12.83560 15.89865 17.86007 20.14597 25.76238
## mu[4]      6.35829  8.54753 10.09091 11.99790 17.20872
## mu[5]      1.90248  2.88995  3.69577  4.84095  8.82738

Pedindo novas quantidades de interesse

Função de sobrevivência (foi colocada na sintaxe do modelo);
Não é necessário verificar a convergência novamente;
pois a função de sobrevivência é resultante dos parâmetros já analisados!

params_novo=c("S")
samples_novo=coda.samples(modelo,params_novo,n.iter=50000,thin=5)
S=summary(samples_novo)
mean_S=S[["statistics"]][,1]

Representação Gráfica:

S_t_Baye=matrix(mean_S,nrow=10,ncol=5)
t=seq(1,50,5)
S_t_Baye=data.frame(cbind(t,S_t_Baye))
names(S_t_Baye)=c("tempos","x_0","x_1","x_2","x_3","x_4")
datatable(S_t_Baye)

colors <- c("x=0"="blue", "x=1"="red", "x=2"="green" ,"x=3" = "orange", "x=4" = "pink")

ggplot(S_t_Baye,aes(x=tempos,y=x_0,color="x=0"))+
    geom_line() +
geom_line(aes(x=tempos,y=x_1,color="x=1")) + 
geom_line(aes(x=tempos,y=x_2,color="x=2"))+
geom_line(aes(x=tempos,y=x_3,color="x=3"))+
geom_line(aes(x=tempos,y=x_4,color="x=4"))+
labs(x="t (em meses)",
         y="S(t)",
         title="Gráfico da função de sobrevivência - modelo exponencial",
         color = "Número de Incompatibilidades (x)") +
    scale_color_manual(labels=c(expression(x==0),expression(x==1),expression(x==2),expression(x==3),expression(x==4)),values = colors)

Conclusões

A metodologia Bayesiana trouxe resultados satisfatórios - com prioris não informativas
Houve convergência do algoritmo.

Pendências

Os gráficos são parecidos com o survreg (abordagem não Bayesiana), no entanto a função de sobrevivência estimada via abordagem Bayesiana possui resultados bem menores (tendendo para zero mais rápido)
Verificar porque não fornece os mesmos resultados!

Propostas futuras

Considerar:

outras distribuições para o tempo de falha (ou tempos de sobrevivência): Weibull, Lognormal, etc;
o modelo de riscos proporcionais de Cox;
o modelo AFT - tempo de falha acelerado na presença de covariáveis;
Considerar outras abordagens para a covariável: além da abordagem em grupos (fatores) - pode-se considerar variável quantitativa contínua;
outros tipos de censura;
Representações gráficas para comparar Abordagem Clássica x Bayesiana;
Representações gráficas para comparar KM com modelos paramétricos.
Modelos de contagem de Poisson: para eventos recorrentes - indivíduos que apresentam reincidência da doença.

Análise de sobrevivência - Enfoque bayesiano com pacote RJAGS

Lia Hanna Martins Morita

2020