Exemplo

O exemplo a ser utilizado será:

• Tempos de reincidência em meses de pacientes com câncer de bexiga submetidos a um procedimento cirúrgico feito por laser.

Para tanto, iremos considerar a distribuição Log-Logístico, comparado com o Kaplan-Meier e, para o segundo método de comparação, utilizando-se linearização, utilizaremos a Gamma-Generalizada.

Modelo Log-Logístico

name = 'Log-Logistico'
name_r = 'llogis'

library(flexsurv)

## Warning: package 'flexsurv' was built under R version 4.1.3

## Carregando pacotes exigidos: survival

library(survival)
library(survminer)

## Warning: package 'survminer' was built under R version 4.1.3

## Carregando pacotes exigidos: ggplot2

## Warning: package 'ggplot2' was built under R version 4.1.3

## Carregando pacotes exigidos: ggpubr

## Warning: package 'ggpubr' was built under R version 4.1.3

## 
## Attaching package: 'survminer'

## The following object is masked from 'package:survival':
## 
##     myeloma

library(lubridate)

## Warning: package 'lubridate' was built under R version 4.1.3

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

tempos<-c(3,5,6,7,8,9,10,10,12,15,15,18,19,20,22,25,28,30,40,45)
cens<-c(1,1,1,1,1,1,1,0,1,1,0,1,1,1,1,1,1,1,1,0)

Ajustando o modelo

Abaixo vamos ajustar o modelo utilizando-se da função flexsurvreg

name = 'Log-Logistico'
name_r = 'llogis'

fs2 <- flexsurvreg(Surv(tempos,cens)~1,
                   dist = name_r)

Vamos obter os valores da estimativa por verossimilhança do modelo:

shape = fs2$coefficients[1] # Parametro de shape
scale = fs2$coefficients[2] # Parametro de scale

# Funcao de sobrevivencia
S <- function(t, shape, scale){
  (1 + ((t/scale)**(shape)))**(-1)
}

#Obtendo as estimativas das fun??es de sobreviv?ncia
ekm<-survfit(Surv(tempos,cens)~1) #estimador de kaplan-meier
time<-ekm$time #tempos
st<-ekm$surv #sobreviv?ncia kaplan meier
stg <- S(t=time, shape=shape, scale=scale)

knitr::kable(cbind(time,st,stg))

Método 01

par(mfrow=c(1,2))
plot(stg,st,pch=16,ylim=range(c(min(st),max(st))), xlim=range(c(0,1)), ylab = "S(t): Kaplan-Meier",
     xlab=paste("S(t):", name))
lines(c(0,1), c(0,1), type="l", lty=1)
grid()

plot(ekm, conf.int=F, xlab="Tempos", ylab="S(t)")
lines(c(0,time),c(1,stg), lty=2)
legend("bottomleft",lty=c(1,2),c("Kaplan-Meier", name),bty="n",cex=0.8)
grid()

Método 02

Aqui utilizaremos uma transformação para linear a relação entre o tempo e a função de sobrevivencia.

• Demonstração:

Para tanto, sabemos que a função de sobrevivencia da gamma generalizada é dada por:

\[ S(t) = \frac{1}{1 + \bigg(\frac{t}{\alpha}\bigg)^{\beta}} \rightarrow \]

também podemos escrever da seguinte forma:

\[ \rightarrow 1 + \bigg(\frac{t}{\alpha}\bigg)^{\beta} = \frac{1}{S(t)} \rightarrow \]

isolando \(\bigg(\frac{t}{\beta}\bigg)^{\beta}\), temos

\[ \rightarrow \bigg(\frac{t}{\beta}\bigg)^{\beta} = \frac{1}{S(t)} - 1 = \frac{1 - S(t)}{S(t)} \rightarrow \]

\[ \rightarrow \log\bigg(\frac{1 - S(t)}{S(t)}\bigg) = \beta\cdot \log\bigg(\frac{t}{\alpha}\bigg) \rightarrow \] \[ \rightarrow \log\bigg(\frac{1 - S(t)}{S(t)}\bigg) = \beta\cdot [\log(t) - \log(\alpha)] \rightarrow \]

\[ \rightarrow \log(t) = \frac{1}{\beta}\cdot \log\bigg(\frac{ 1 - S(t)}{S(t)}\bigg) + \log(\alpha) \] seja \(k = \log(\alpha)\), \(\tau = \frac{1}{\beta}\), \(X = \log\bigg(\frac{1-S(t)}{S(t)}\bigg)\) e \(Y = \log(t)\), então notamos a seguinte expressão:

\[ Y = \tau\cdot X + k \]

• Computacionalmente:

lstg= (1/scale)*log((1-st)/st) + log(shape) # s(t): kaplan meier

par(mfrow=c(1,1))
plot(log(time), lstg,pch=16,
     xlab=expression(log(t)),
     ylab=expression((1/beta)*log((1-S(t))/S(t)) + log(alpha)),
     main=name)
abline(lm(lstg ~ log(time)))
grid()

Teste de RV

Aplicaremos um teste de razão de verossimilhanças.

fitg <- flexsurvreg(formula = Surv(tempos, cens) ~ 1,dist=name_r)
fitgg <- flexsurvreg(formula = Surv(tempos, cens) ~ 1,dist="gengamma")

#Obten??o das log-verossimilhan?as
loglikLLogis<-fitg$loglik[1]  #log-verossimilhan?a modelo log-logistico
loglikG<-fitgg$loglik[1]   #log-verossimilhan?a modelo generalizado 

#Obten??o das estat?sticas de teste

RVe<-2*(loglikG-loglikLLogis) #Estat?stica da raz?o de verossimilhan?as para o modelo exponencial

#Obten??o dos p-valores dos testes            
pRV<-1-pchisq(RVe,2) #Log-logistico

rbind(pRV,
      'LogLik.LLogis'=fitg$loglik[1],
      'LogLik.GG'=fitgg$loglik[1])

##                     [,1]
## pRV             0.711922
## LogLik.LLogis -66.030530
## LogLik.GG     -65.690743

Não rejeitamos a hipótese nula, ou seja, pelo teste de razão de verossimilhanças, e em suma pela analise gráfica, podemos ver que o modelo Log-Logístico não se ajustou bem para esse conjunto de dados.

References