Analise de sobrevivência

Ja vimos como algumas vezes mesmo um valor p significativo não significa muita coisa.
Um tipo de modelo muito interessante e o que não tive muitas oportunidades de explorar são os de sobrevivência.

Para fazer um pequeno teste onde ele é melhor que uma análise de variância, vamos fazer um pequeno exemplo.

Uma aranha, vivendo sua vidinha numa flor qualquer corre o risco todo dia de ser predada. Podemos pensar nesse risco como uma chance.
Então podemos simular 100 dias da vida dessa aranha e ver se ela morre algum dia dado que ela tem 10% de chance de morrer por dia.

> rnos which(rnos<= 0.1) [1] 22 24 47 62 70 84 91

Então vemos aqui quais dias ela não passou na prova da sobrevivência, e quais dias ela morreu. Mas como todo mundo só morre uma vez, podemos dizer na verdade que essa aranha morreu no dia 22 de sua vida.

> which(rnos<= 0.1)[1] [1] 22

Certo e qual a relação do dia da morte com a chance de morte da aranha?
Bem se observamos varias aranhas, simulando 30 aranhas com o mesmo procedimento acima vemos que a chance de morte vai ser 1 dividido pela média em dias de vida

> death1 [1] 7 10 1 1 12 9 1 24 3 1 21 5 1 12 3 4 11 19 16 7 15 2 2 6 6 36 6 1 2 4 > 1/mean(death1) [1] 0.1209677

Usamos 10% na simulação, para 30 aranhas obtivemos 12% como estimativa grosseira de sobrevivência. Estimativa razoavelmente boa eu diria.
Agora poderíamos pensar em outra espécie com uma chance maior de morrer, ou ser predada, uma espécie com uma camuflagem ruim, o que aconteceria é que a média em dias seria menor correto? Os indivíduos tenderiam a morrer mais jovens. Vamos simular isso e ver o que acontece, então agora com 20% de chance de morrer.

> death2 [1] 1 1 1 2 8 5 2 3 20 1 4 1 5 9 6 3 1 3 2 3 2 5 11 7 3 8 2 5 2 2 > 1/mean(death2) [1] 0.234375

Certo vemos uma estimativa próxima do que foi usado na simulação, obtivemos 23%, usamos 20%, e vemos aranhas morrendo mais cedo.
Podemos perguntar então se em média aranhas da segunda espécie simulada aqui morrem mais cedo que da primeira e fazendo um modelo linear (teste t) temos.

> summary(model0) Call: lm(formula = death ~ factory) Residuals: Min 1Q Median 3Q Max -7.267 -3.267 -1.767 1.983 27.733 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8.267 1.188 6.956 3.48e-09 *** factory2 -4.000 1.681 -2.380 0.0206 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 6.51 on 58 degrees of freedom Multiple R-squared: 0.08896, Adjusted R-squared: 0.07326 F-statistic: 5.664 on 1 and 58 DF, p-value: 0.02063

Diferença significativa entre as espécies. Mas temos um problema, estamos quebrando fortemente uma premissa dos modelos lineares, a da homogeneidade das variâncias.
Olhem as variâncias das espécies:

> var(death1) [1] 68.75402 > var(death2) [1] 15.9954

Olhando o gráfico, vejam como a espécie um, que tinha 10% de chance de morrer tem uma variância muito maior que a aranha que tinha 20% de morrer. Isso não esta sendo considerado no modelo, o que o deixa muito ruim.

Uma solução é fazer um modelo geral linear com distribuição gamma.
Mas porque essa complicação? O que vamos ver sempre, é que quanto maior a chance de morrer, menor vai ser a média e a variação.
Não é isso que estamos vendo no gráfico ali em cima, dobramos a chance de morte, de 10% para 20% e a média e variância caíram vertiginosamente, podemos fazer mais simulações e ver se isso se confirma, caso não consigamos deduzir isso matematicamente, eu particularmente não consigo, mas deixando isso para depois, vamos ajustar o modelo geral linear e ver o que vamos obter.

> model1 summary(model1) Call: glm(formula = death ~ factory, family = Gamma) Deviance Residuals: Min 1Q Median 3Q Max -1.5705 -0.9383 -0.3327 0.3153 2.0701 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 0.12097 0.02144 5.642 5.26e-07 *** factory2 0.11341 0.04675 2.426 0.0184 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for Gamma family taken to be 0.9423881) Null deviance: 56.867 on 59 degrees of freedom Residual deviance: 50.422 on 58 degrees of freedom AIC: 337.43 Number of Fisher Scoring iterations: 6

As estimativas são nossas probabilidades, a da primeira espécie (tomado como base, intercepto) é 0.12, próximo do usado na simulação, e para a espécies dois temos uma diferença significativa (p de 0.0184), sendo essa maior que o intercepto 0.1231.
0.12 mais 0.11 da 0.23, próximo dos 20% que usamos na simulação, se pensarmos nos erros dessas medidas (Std. Error) temos que o valor da simulação esta ai dentro.

Se compararmos os resíduos, o que a parte determinística do modelo não consegue descrever, podemos ter uma noção de qual abordagem é mais apropriada.

Vemos como usando o modelo linear temos resíduos gigantes, além de uma tendência a sempre sobrar mais resíduos para baixo, enquanto o que queremos é mais o que vemos no resíduos gerados pelo glm gamma, pequenos e sem nenhuma tendência, esperamos errar igual tanto para mais quanto para menos.

Mas existe ainda uma terceira abordagem para esse problema, que é finalmente a analise de sobrevivência do titulo do post.

Então usamos uma função de sobrevivência para ajustar uma chance de sobrevivência para nossas arahas. Isso e muitas possibilidades estão disponíveis no pacote survival no R.

Então vamos la e primeiro ajustamos as curvas de sobrevivências para as duas espécies de aranhas que estamos trabalhando:

E temos uma escadinha de como a sobrevivência vai diminuindo ao longo do tempo. Lembre-se a cada dia você acumula mais chance de morrer, então o inverso, diminui sua chance de estar vivo a cada dia.

Bem, pelo gráfico já podemos ver como a diferença entre as aranhas, as do exemplo 2 (cor azul) tem uma mortalidade maior, logo sobrevivem menos, logo a curva cai bem mais rápido.
Podemos ajustar parâmetros, no caso essa chance de sobrevivência usando a função survreg() do pacote survival e temos:

> summary(model) Call: survreg(formula = Surv(death, event = rep(1, length(death))) ~ factory, dist = "exponential") Value Std. Error z p (Intercept) 2.112 0.183 11.57 5.91e-31 factory2 -0.661 0.258 -2.56 1.04e-02 Scale fixed at 1 Exponential distribution Loglik(model)= -166.9 Loglik(intercept only)= -170.1 Chisq= 6.45 on 1 degrees of freedom, p= 0.011 Number of Newton-Raphson Iterations: 4 n= 60

Vemos as estimativas, 2.112 para aranhas do exemplo 1 e as aranhas do exemplo 2 diferentes por -0.66. Bem esses valores precisam ser convertidos pela função logística de volta para probabilidade (chance) para podermos comparar com os valores da simulação, então convertendo…

> plogis(coefficients(model)[1]) (Intercept) 0.8920863 > plogis(coefficients(model)[1]+coefficients(model)[2]) (Intercept) 0.8101266

Hummmm 89% e 81%, ahhhh, mas aqui são estimativas da sobrevivência, e simulamos a chance de morrer, logo 1-sobrevivencia será a chance de mortalidade.

> 1-plogis(coefficients(model)[1]) (Intercept) 0.1079137 > 1-plogis(coefficients(model)[1]+coefficients(model)[2]) (Intercept) 0.1898734

Ae, 11% e 19%, valores muito próximos do usados na simulação.

Bem parece que analise de sobrevivência funciona muito bem recuperando os parâmetros que usamos para simular os dados.
Mas existem dois motivos que acho muito muito interessantes para usar analise de sobrevivências.

O primeiro, que está relacionado até com a motivação de criação desta técnica, é a de que para traçar a curva de sobrevivência, não precisamos acompanhar todos os indivíduos até a morte, pode parecer pouca coisa isso, mas se pensarmos por exemplo em pessoas e remédios, muitas vezes acompanhamos pessoas utilizando ou não um determinado tratamento, mas só podemos acompanhar essas pessoas durante o tratamento, então mesmo que não saibamos quando elas morrem, sabemos que ela sobreviveram o período de tempo que as acompanhamos, então temos ai informação a cerca da sobrevivência. Ou ainda temos que realizar um experimento num tempo determinado, então aquilo que não seria informação usando glm, aqui será informação, porque la esperamos todos os indivíduos morrerem.

Outra coisa é que em todos os casos, consideramos uma chance constante de mortalidade (ou sobrevivência se preferir), mas facilmente sabemos que isso pode ser diferente. O exemplo clássico é comparar pólipos com humanos com tartarugas, pólipos sim tem uma chance constante de morrer no mar, já que morrem quando são predados, nos humanos, muito diferente disso, quando nascemos, temos o cuidado de nossos pais, o que diminui muito a chance de morrer quando jovens, mas conforme envelhecemos, a partir de certa idade, doenças e muitos problemas vão aparecendo, diminuição da capacidade física logo podemos esperar um aumento gradativo da chance de mortalidade.
E o ultimo extremo seriam tartarugas, que quando nascem sozinhas na areia tem uma chance altíssima de mortalidade, mas esta vai diminuindo com o tempo, já que ela vai ficando grande e dai consegue se defender cada vez melhor. Isso seria uma complicação enorme nos modelos lineares e glms, mas podemos usando a regressão de sobrevivencia por exemplo, ajustar uma distribuição e modelo para essa alteração de sobrevivência, vejam que aqui usamos dist=”exponential”, ou seja a mortalidade é constante. Mas podemos ajustar vários modelos e comparar o ajuste por exemplo para ter uma visão mais clara de como a sobrevivência muda ao longo do tempo, se um tratamento demora a começar a fazer efeito, etc.

Bem é isso, uma primeira visita a possibilidade do uso de modelos de sobrevivência dentro da biologia. Segue o script, a maior parte dele esta no capitulo de analise de sobrevivência do livro do Crawley “The R Book” para quem quiser ler mais sobre o assunto.

##################################################
set.seed(3)
rnos<-runif(100)
which(rnos<= 0.1)
which(rnos<= 0.1)[1]
 
#exemplo 1
death1<-numeric(30)
for (i in 1:30) {
rnos<-runif(100)
death1[i]<- which(rnos<= 0.1)[1]
}
death1
 
1/mean(death1)
 
#exemplo 2
death2<-numeric(30)
for (i in 1:30) {
rnos<-runif(100)
death2[i]<- which(rnos<= 0.2)[1]
}
death2
1/mean(death2)
 
#juntando os dados
death<-c(death1,death2)
factory<-factor(c(rep(1,30),rep(2,30)))
 
plot(factory,death,ylab="Dia da Morte",frame.plot=F)
 
#Modelo linear
model0<-lm(death~factory)
summary(model0)
 
var(death1)
var(death2)
 
#Modelo gamma
model1<-glm(death~factory,Gamma)
summary(model1)
 
#grafico dos residuos
par(mfrow=c(2,1))
plot(resid(model0),ylim=c(-10,30),main="Modelo Linear (Distribuição normal)"
,ylab="Residuos",xlab="Aranhas")
abline(h=0,lwd=2,lty=2)
plot(resid(model1),main="Modelo Geral Linear (Distribuição gamma)",
ylim=c(-10,30),ylab="Residuos",xlab="Aranhas")
abline(h=0,lwd=2,lty=2)
 
#curva de sobrevivencia
library(survival)
plot(survfit(Surv(death,event=rep(1,length(death)))~factory),lwd=2,
col=c("red","blue"),frame.plot=F,xlab="Dias",ylab="Sobrevivencia")
legend("topright",legend=c("Exemplo 1","Exemplo 2"),lwd=2,col=c("red","blue"))
 
#Regressão de sobrevivencia
model<-survreg(Surv(death,event=rep(1,length(death)))~factory,dist="exponential")
summary(model)
 
#Probabilidades de sobrevivencias
plogis(coefficients(model)[1])
plogis(coefficients(model)[1]+coefficients(model)[2])
 
#Probabilidades de mortalidade
1-plogis(coefficients(model)[1])
1-plogis(coefficients(model)[1]+coefficients(model)[2])

Leave a Reply

Your email address will not be published. Required fields are marked *

* Copy This Password *

* Type Or Paste Password Here *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>