Curso R - UFMG - Análise de sobrevivência
Ricardo Solar - UFMG (ricardosolar@ufmg.br)
Colaboração do Prof. Og De Souza - UFV (https://www.isoptera.ufv.br)
November, 2023
© Ricardo Solar/UFMG - compartilhamento e utilização não-comercial livres. Não reproduzir sem autorização prévia
1 Contextualizando a análise de sobrevivência
Algumas vezes nos interessa medir o tempo decorrido até que um determinado evento ocorra, por exemplo:
- expostos a um determinado tratamento, por quanto tempo sobrevivem os ratos?
- quanto tempo leva até uma semente germinar?
- quanto tempo é necessário para que um armadilha capture o primeiro inseto?
- dentre dois tipos de isca, qual deles é localizado primeiro?
- por quantas horas posso garantir o funcionamento de uma máquina?
Vejam que, em todos esses casos, eu espero que ao final do meu experimento o evento de fato possa ocorrer. Ou seja, ao final, independente do tratamento, eu posso ter 100% do evento em ambos os tratamentos. Por isso, a variável que me interessa é quanto tempo se decorreu até o evento, e não a porcentagem de ocorrência do mesmo. Isso representa uma mudança de paradigma nas análises que normalmente fazemos, pois mesmo que os níveis do tratamento acumulem a mesma quantidade de observações ao final do tempo do experimento, ainda poderemos distinguir diferenças entre um e outro.
1.1 A variável tempo
Determinar estatisticamente o momento em que algum evento ocorre é um problema muito mais comum e necessário do que pode parecer à primeira vista. Poderemos calcular, por exemplo, o tempo médio gasto para que um inseto chegue a uma isca. Mais comum ainda é medirmos a quantidade de tempo decorrida do início do experimento até a morte da maioria dos indivíduos de uma coorte, ou seja, o tempo médio de sobrevivência destes indivíduos.
O cálculo consiste em construirmos uma curva descendente que descreve a percentagem de indivíduos da coorte que ainda estão vivos a cada novo incremento do tempo, desde o início do experimento. Esta curva, que naturalmente começa no tempo zero, quando 100% estão vivos, e pode cair até o percentual de 0% de vivos, depentendo do tempo do experimento. Observe que vivo/morto aqui representa uma mudança de estado, ou seja, o mesmo raciocínio vale para atacado/não atacado, germinada/não germinada, funcionando/estragado … e por aí vai!
Pois bem, estamos convencidos que o tempo é uma variável crucial, contudo ela é traiçoeira. Frequentemente precisamos entender como se processa a dinâmica de um dado fenômeno e somos tentados a tomar medidas de uma mesma unidade experimental em vários momentos ao longo do tempo. Por exemplo, para saber a taxa de crescimento de uma planta, medimos esta planta diversas vezes ao longo do tempo, para produzir um gráfico de regressão onde tempo é a variável explicativa (x) e tamanho da planta é a variável resposta (y). Embora isto seja correto para se estimar a rapidez com que aquele indivíduo aumenta de tamanho, isto não pode ser extrapolado para os demais indivíduos daquela espécie, pois aí teríamos uma pseudo-repetição temporal. Nesse caso, medir várias vezes o mesmo indivíduo fere o princípio da independência amostral.
1.2 Tipos de curva de sobrevivência
O maior problema associado à variável tempo, entretanto, parece estar associado à avaliação da sobrevivência, pois há três padrões básicos:
- [Tipo I:] a maioria da população sobrevive até idades avançadas. Dito de outra forma: a probabilidade de morrer aumenta com a idade.
- [Tipo II:] a percentagem da população que morre é constante ao longo do tempo. Isto é, a probabilidade de morrer não se altera ao longo da idade
- [Tipo III:] a maioria da população morre em idade jovem, ou seja, a probabilidade de morrer é maior na infância do que em idades mais avançadas
Fonte da imagem: https://www.britannica.com/science/survivorship-curve
Como estamos falando de `probabilidade’, você já pode imaginar que cada um destes três padrões geraria diferentes funções de distribuição de frequência e, com isso, o cálculo de significância exigiria diferentes distribuições de erros. Isso seria um problemão, mas aqui entram as distribuições de sobrevivência, especialmente a Weibull. sta distribuição de probabilidades assume formatos que servem para os três tipos de sobrevivência descritos acima e, com isso, elimina a necessidade de testarmos várias distribuições como fizemos no caso de GLM clássico. Note que esta não é a única análise de sobrevência existente, nem a única análise de sobrevivência paramétrica. Há vários outros métodos, paramétricos ou não, mas aqui nos concentraremos na distribuição de Weibull por ser um método paramétrico (por definição mais conservador que os métodos não-paramétricos) muito comum na literatura biológica.
1.3 Variável censora
Um outro problema em experimentos que avaliam sobrevivência é o fato de algums indivíduos sobreviverem mais tempo do que o tempo do experimento! Isto é, eventualmente você precisa dar um basta nas observações e sempre tem um ou outro indivíduo que ainda está vivo. Como calcular o tempo médio que aquele grupo sobreviveu? Se eu deixar de considerar este indivíduo já estarei alterando a média daquele tratamento a priori. Se eu considerar que este indivíduo morreu naquele momento, também não estarei usando a informação correta e poderei privilegiar um tratamento em relação ao outro. O problema é que não sabemos quando o indivíduo morrerá. A dúvida clássica é: se eu esperasse mais um dia, será que todos os indivíduos morreriam?
Neste caso, a solução é usar uma variável censora, que é uma variável que informa ao programa quando foi o último momento em que aquele indivíduo foi visto vivo pela última vez. Isso informa que eu sei que o indivíduo estava vivo naquele momento, mesmo que eu não saiba quando ele de fato vai morrer. O algoritmo considerará este valor ao calcular o tempo médio de sobrevivência até aquele momento, mas o desconsiderará ao calcular médias maiores.
A variável censora também resolve aquele caso clássico do indivíduo que some ou que é adulterado durante o experimento. Isto é comum no caso de insetos presos em gaiolas ou iscas deixadas no campo. Eu sei que observei o indivíduo no tempo anterior, e ele estava lá “saudável”, porém na observação de agora ele não está lá mais, ou está com alguma adulteração que não deve-se ao experimento. O inseto pode ter fugido da gaiola e em seguinda morrido, pode ter havido canibalismo ou o danado pode estar vivo até agora!
O uso da variável censora é muito simples: se o indivíduo estiver vivo no momento em que encerramos o experimento, obviamente este é o momento em que ele foi visto vivo para última vez. Neste caso anotamos o valor de tempo que corresponde a este momento e usamos o valor 0 (zero) para a variável censora para indicar que o indivíduo ainda não tinha morrido. O mesmo se aplica no caso de indivíduos que desapareceram no meio do experimento, devo anotar o _último momento que o indivíduo foi visto vivo, e atribuir a ele o valor 0. Já para o caso de um indivíduo que eu vi efetivamente morto num determinado tempo, anoto aquele momento, nesse caso com o valor 1 (um).
2 Exemplo 1
Vamos usar os dados de Crawley (1997) que aplicou 3 tipos de inseticidas em baratas, objetivando verificar qual inseticida é mais eficiente. O problema é que as baratas variavam de tamanho e, por isso, teriam diferenças na sua resistência ao inseticida. Por isso, a análise inclui uma co-variável (peso da barata) para extrair este efeito.
2.1 Procedimentos iniciais
Este procedimento não é absolutamente necessário, mas pode ser bastante útil para limpar o histórico de análises e deixar o programa como o Environment vazio.
rm(list=ls()) ##remove todos os elementos existentes no Environment2.2 Carregando os pacotes que usaremos
library(survival)
library(ggplot2)
library(survminer)## Loading required package: ggpubr##
## Attaching package: 'survminer'## The following object is masked from 'package:survival':
##
## myeloma2.3 Carregando os dados
dados <- read.table("baratas.txt", h=T, stringsAsFactors = T) ## Warning in file(file, "rt"): cannot open file 'baratas.txt': No such file or
## directory## Error in file(file, "rt"): cannot open the connectionsummary(dados)## tempo censor peso grupo
## Min. : 1.00 Min. :0.0000 Min. : 0.055 g1:50
## 1st Qu.: 1.00 1st Qu.:1.0000 1st Qu.: 2.459 g2:50
## Median : 7.00 Median :1.0000 Median : 6.316 g3:50
## Mean :15.17 Mean :0.8667 Mean : 9.390
## 3rd Qu.:21.00 3rd Qu.:1.0000 3rd Qu.:11.955
## Max. :50.00 Max. :1.0000 Max. :42.0902.4 Criando o modelo
mod1 <- survreg(Surv(tempo,censor)~grupo*peso, data=dados)
anova(mod1)Pudemos ver que a interação e o peso não foram significativos, desta forma podemos começar a simplificar nosso modelo. Para simplificar o modelo, vamos retirar primeiro o fator mais complexo, em nosso caso a interação grupo:peso. Então nosso modelo fica assim:
mod2 <- survreg(Surv(tempo,censor)~grupo+peso, data=dados)
anova(mod1, mod2)anova(mod2)Com esta nova análise, podemos ver que a variável peso também não é significativa, portanto podemos retirá-la. Então vamos simplificar o modelo mais uma vez. Desta vez vamos retirar a variável peso.
mod3 <- survreg(Surv(tempo,censor)~grupo, data=dados)
anova(mod2, mod3)anova(mod3)Nosso modelo continua a ser significativo, como não temos mais nenhuma variável para retirar, este é nosso modelo final. Mas nossa simplificação ainda não é completa, nova variável grupo possui três níveis (três grupos) e sua significância apenas nos diz que existe pelo menos um nível (um grupo) diferente dos demais. Portanto, temos que verificar se nós temos dois ou três níveis (dois ou três grupos). Devemos analisar os contrastes. Quem juntar primeiro? Seguiremos igual o que fizemos antes, e vamos juntar aqueles níveis que têm menor diferença absoluta entre si.
summary(mod3)##
## Call:
## survreg(formula = Surv(tempo, censor) ~ grupo, data = dados)
## Value Std. Error z p
## (Intercept) 3.4593 0.2283 15.15 < 2e-16
## grupog2 -0.8222 0.3097 -2.65 0.0079
## grupog3 -1.5403 0.3016 -5.11 3.3e-07
## Log(scale) 0.3145 0.0705 4.46 8.1e-06
##
## Scale= 1.37
##
## Weibull distribution
## Loglik(model)= -470.5 Loglik(intercept only)= -483.3
## Chisq= 25.63 on 2 degrees of freedom, p= 2.7e-06
## Number of Newton-Raphson Iterations: 5
## n= 150Pelo summary do modelo, vemos que g1 e g2 são os mais próximos. Vamos então criar uma nova variável grupo.12 recategorizando os grupos.
grupo2 <- dados$grupo
levels(grupo2)## [1] "g1" "g2" "g3"levels(grupo2)[1] <- "g1.g2"
levels(grupo2)[2] <- "g1.g2"
levels(grupo2)## [1] "g1.g2" "g3"E devo testar se os modelos têm diferença significativa:
mod4 <- survreg(Surv(tempo,censor)~grupo2, data=dados)
anova(mod3,mod4)Os modelos são distintos, portanto ao juntar, eu pioro a qualidade da informação. Conclusão, g1 e g2 são diferentes. Faço o mesmo para g2 e g3.
grupo3 <- dados$grupo
levels(grupo3)## [1] "g1" "g2" "g3"levels(grupo3)[2] <- "g2.g3"
levels(grupo3)[3] <- "g2.g3"
levels(grupo3)## [1] "g1" "g2.g3"E devo testar se os modelos têm diferença significativa:
mod5 <- survreg(Surv(tempo,censor)~grupo3, data=dados)
anova(mod3,mod5)Os modelos também são distintos, portanto ao juntar, eu pioro a qualidade da informação. Conclusão, g1, g2 e g3 são diferentes.
2.5 Construindo o gráfico e obtendo as informações necessárias
Para se construir o gráfico de sobrevivência é necessário dois parâmetros. O primeiro parâmetro é o , que no R deve ser obtido dividindo-se 1/ scale . Esse valor de scale pode ser encontrado no summary(modelo), que no nosso exemplo é 1.37, portanto o \(\alpha = \frac{1}{1.37}\) = 0.729927.
O segundo parâmetro é o tempo médio para morte, denotado por \(\mu\) para cada tratamento. O tempo médio para morte é dado pelo antilog do parâmetro estimado, parâmetro este mostrado com o comando summary(modelo) ou coef(modelo), para isto basta fazer exponencial do valor de intercepto para cada grupo. Veja abaixo:
alfa <- 1/(summary(mod3)$scale)
summary(mod3)##
## Call:
## survreg(formula = Surv(tempo, censor) ~ grupo, data = dados)
## Value Std. Error z p
## (Intercept) 3.4593 0.2283 15.15 < 2e-16
## grupog2 -0.8222 0.3097 -2.65 0.0079
## grupog3 -1.5403 0.3016 -5.11 3.3e-07
## Log(scale) 0.3145 0.0705 4.46 8.1e-06
##
## Scale= 1.37
##
## Weibull distribution
## Loglik(model)= -470.5 Loglik(intercept only)= -483.3
## Chisq= 25.63 on 2 degrees of freedom, p= 2.7e-06
## Number of Newton-Raphson Iterations: 5
## n= 150coef(mod3)## (Intercept) grupog2 grupog3
## 3.4593448 -0.8222432 -1.5403091(mug1 <- exp(coef(mod3)[1]))## (Intercept)
## 31.79614(mug2 <- exp(coef(mod3)[1]+coef(mod3)[2]))## (Intercept)
## 13.97265(mug3 <- exp(coef(mod3)[1]+coef(mod3)[3]))## (Intercept)
## 6.814384Tendo em mãos estes valores, para fazer o gráfico basta usar a função da equação de sobrevivência:
\(exp(-\mu^{-\alpha}*t^{\alpha})\)
lembrando que t é a nossa variável x, \(\alpha\) foi o valor que calculamos a partir de scale e \(\mu\) é um valor por curva. Desta forma podemos construir a formula e mandar construir o gráfico em qualquer programa, no R pode-se usar a função curve.
2.6 Obtendo os pontos de ajuste do gráfico
Temos o resultado do modelo Survreg, agora para obter o ajuste, utilizaremos a função Survfit:
m.fit3 = survfit(Surv(tempo,censor)~grupo, data=dados)
summary(m.fit3)## Call: survfit(formula = Surv(tempo, censor) ~ grupo, data = dados)
##
## grupo=g1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 50 8 0.84 0.0518 0.744 0.948
## 2 42 1 0.82 0.0543 0.720 0.934
## 3 41 5 0.72 0.0635 0.606 0.856
## 4 36 1 0.70 0.0648 0.584 0.839
## 6 35 1 0.68 0.0660 0.562 0.822
## 7 34 2 0.64 0.0679 0.520 0.788
## 9 32 1 0.62 0.0686 0.499 0.770
## 10 31 1 0.60 0.0693 0.478 0.752
## 11 30 1 0.58 0.0698 0.458 0.734
## 12 29 1 0.56 0.0702 0.438 0.716
## 13 28 2 0.52 0.0707 0.398 0.679
## 17 26 1 0.50 0.0707 0.379 0.660
## 18 25 1 0.48 0.0707 0.360 0.641
## 20 24 1 0.46 0.0705 0.341 0.621
## 21 23 1 0.44 0.0702 0.322 0.602
## 22 22 1 0.42 0.0698 0.303 0.582
## 24 21 1 0.40 0.0693 0.285 0.562
## 26 20 1 0.38 0.0686 0.267 0.541
## 34 19 1 0.36 0.0679 0.249 0.521
## 36 18 2 0.32 0.0660 0.214 0.479
## 37 16 1 0.30 0.0648 0.196 0.458
## 46 15 1 0.28 0.0635 0.180 0.437
##
## grupo=g2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 50 12 0.76 0.0604 0.6504 0.888
## 2 38 7 0.62 0.0686 0.4991 0.770
## 3 31 2 0.58 0.0698 0.4581 0.734
## 4 29 3 0.52 0.0707 0.3984 0.679
## 5 26 5 0.42 0.0698 0.3032 0.582
## 7 21 1 0.40 0.0693 0.2849 0.562
## 8 20 2 0.36 0.0679 0.2488 0.521
## 11 18 1 0.34 0.0670 0.2311 0.500
## 12 17 1 0.32 0.0660 0.2136 0.479
## 17 16 2 0.28 0.0635 0.1795 0.437
## 19 14 1 0.26 0.0620 0.1629 0.415
## 21 13 2 0.22 0.0586 0.1305 0.371
## 38 11 1 0.20 0.0566 0.1149 0.348
## 40 10 1 0.18 0.0543 0.0996 0.325
## 41 9 1 0.16 0.0518 0.0848 0.302
## 46 8 2 0.12 0.0460 0.0566 0.254
##
## grupo=g3
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 50 19 0.62 0.0686 0.49906 0.770
## 2 31 3 0.56 0.0702 0.43801 0.716
## 3 28 1 0.54 0.0705 0.41811 0.697
## 5 27 4 0.46 0.0705 0.34067 0.621
## 6 23 2 0.42 0.0698 0.30324 0.582
## 7 21 3 0.36 0.0679 0.24877 0.521
## 9 18 2 0.32 0.0660 0.21363 0.479
## 10 16 1 0.30 0.0648 0.19644 0.458
## 11 15 4 0.22 0.0586 0.13054 0.371
## 13 11 1 0.20 0.0566 0.11489 0.348
## 14 10 1 0.18 0.0543 0.09962 0.325
## 17 9 1 0.16 0.0518 0.08478 0.302
## 20 8 3 0.10 0.0424 0.04354 0.230
## 21 5 1 0.08 0.0384 0.03125 0.205
## 23 4 1 0.06 0.0336 0.02003 0.180
## 24 3 1 0.04 0.0277 0.01029 0.156
## 30 2 1 0.02 0.0198 0.00287 0.139
## 46 1 1 0.00 NaN NA NAplot(m.fit3,
mark=c(1,2,3),
las=1,
col=c(2,4,7),
ylab= "Proporcao de sobreviventes",
xlab="Tempo (dias)",
bty="l")
# legenda
legend("bottomleft",legend=levels(dados$grupo),pch=1:3, col=c(2,4,7),bty="n", lty=1)
Esse é o gráfico mais simples que podemos ter para a sobrevivência. Mas podemos incrementar um pouco, colocando nossas curvas, para as quais obtivemos os parâmetros anteriormente:
curve(exp((-mug1^-alfa)*(x^alfa)), from=0,to=50, ylim=0:1, las=1, lwd=2, col=2, xlab="Tempo", ylab="Prop. Sobreviventes")
curve(exp((-mug2^-alfa)*(x^alfa)), from=0,to=50, add=T, lty=2, lwd=2, col=4)
curve(exp((-mug3^-alfa)*(x^alfa)), from=0,to=50, add=T, lty=3, lwd=2, col=7)
legend("bottomleft",legend=levels(dados$grupo),pch=1:3, col=c(2,4,7),bty="n", lty=1)
#Para adicionar os pontos, devo primeiro extraí-los do m.fit3:
s <- summary(m.fit3)
dfw <- data.frame(surv=s$surv, Grupos=s$strata, time=s$time)
cores <- c(2,4,7)
points(dfw$surv~dfw$time,
col=cores[as.numeric(dfw$Grupos)],
pch=as.numeric(dfw$Grupos),
cex=1.3)
Podemos fazer também a curva no pacote ggsurvplot:
ggsurvplot(
m.fit3, # survfit object with calculated statistics.
pval = TRUE, # show p-value of log-rank test.
conf.int = TRUE, # show confidence intervals for
# point estimaes of survival curves.
xlab = "Time in days", # customize X axis label.
break.time.by = 10, # break X axis in time intervals by 10.
#ggtheme = theme_light(), # customize plot and risk table with a theme.
risk.table = F, # print risk table. Allowed values for risk.table are:
#TRUE, FALSE, 'absolute', 'percentage', 'nrisk_cumcensor', 'nrisk_cumevents'
risk.table.y.text.col = F,# colour risk table text annotations.
risk.table.y.text = F,# show bars instead of names in text annotations
# in legend of risk table.
ncensor.plot = F, # plot the number of censored subjects at time t
surv.median.line = "h", # add the median survival pointer.
legend.title = "Grupos", # adiciona o nome da lengeda correto
legend.labs = c("G1", "G2", "G3"), # change legend labels.
palette = c("#E7B800", "#2E9FDF", "#2E6CDF") # custom color palettes.
)
2.7 Como obter o tempo médio para morte e o tempo para morte de 50% dos indivíduos?
2.7.1 Tempo médio para morte
O cálculo do tempo medio de de ocorrência do evento para cada grupo, ou seja, o tempo medio para morte é apresentado em John E. Pinder III, James G. Wiener and Michael H. Smith. 1978. The Weibull Distribution: A New Method of Summarizing Survivorship Data. Ecology, 59(1), 175-179
A fórmula é dada por:
\(\mu*\Gamma(1+\frac{1}{\alpha})\)
onde:
- \(\mu\) é o valor que calculamos para cada intercepto
- \(\alpha\) é o valor que calculamos dividindo 1 por scale
- \(\Gamma\) é uma função definida por Abramowitz and Stegun 1972 (https://en.wikipedia.org/wiki/Abramowitz_and_Stegun provides links to the full text which is in public domain. Chapter 6: Gamma and Related Functions.)
(TPM1 <- mug1*gamma(1+(1/alfa)))## (Intercept)
## 38.72844(TPM2 <- mug2*gamma(1+(1/alfa)))## (Intercept)
## 17.01901(TPM3 <- mug3*gamma(1+(1/alfa)))## (Intercept)
## 8.300082.7.2 Tempo para morte de 50% dos indivíduos
Nesse caso, basta eu utilizar a própria fórmula abaixo, com y = 0.5:
\(({\frac{\ln(y)}{-\mu^{-\alpha}}})^{\frac{1}{\alpha}}\)
(tmedio1 <- (log(0.5)/((-mug1^-alfa)))^(1/alfa))## (Intercept)
## 19.24749(tmedio2 <- (log(0.5)/((-mug2^-alfa)))^(1/alfa))## (Intercept)
## 8.45821(tmedio3 <- (log(0.5)/((-mug3^-alfa)))^(1/alfa))## (Intercept)
## 4.125023