[R-br] Extraindo o "cumulative hazard" de objeto survfit
Pedro Emmanuel Alvarenga Americano do Brasil
emmanuel.brasil em gmail.com
Qui Dez 16 15:47:19 -03 2021
Cara,
Primeiro um comentário a respeito do IC. Li umas paradas aqui que me levou a
seguinte conclusão. O objeto y$std.chaz é identico ao objeto y$std.err.
Então o cumhaz no y nada mais é que -log(y$surv). Dessa forma, pra achar o
IC do cumhaz basta fazer o -log(y$lower) e -log(y$upper), independente do
método utilizado, sendo que como tem o "-" na frente, o upper vira lower e
vice versa. Está batendo aqui com os meus dados nas primeiras casas
decimais.
cumhaz.summ <- function(y, digits = 3){
y$cumhaz.upper <- -log(y$lower)
y$cumhaz.lower <- -log(y$upper)
for(i in seq_along(levels(y$strata))){
cond <- y$strata == levels(y$strata)[i]
output[[i]] <- round(data.frame(Time = y$time[cond],
'N risk' = y$n.risk[cond],
'N events' = y$n.event[cond],
'Cum Hazard' = y$cumhaz[cond],
lower = y$cumhaz.lower[cond],
upper =
y$cumhaz.upper[cond]),digits)
}
names(output) <- levels(y$strata)
output
}
cumhaz.summ(y)
A segunda parada é que essa função só foi testada para um preditor que
estratifica a curva. No exemplo acima o "~ x". Possivelmente daria errado
se não houver preditor (e dá errado mesmo no meu teste aqui retorna uma
lista vazia pq o y$strata retorna NULL) ou mais de um preditor. Este deve
ser o teu caso já que o erro dá a entender que o comprimento da condição
está retornando objetos de tamanhos diferentes. Mas se colocar um exemplo
reprodutível, talvez eu consiga te ajudar.
Abraço forte,
Pedro Brasil
Pedro Brasil
Em qua., 15 de dez. de 2021 às 21:52, sznelwar--- por (R-br) <
r-br em listas.c3sl.ufpr.br> escreveu:
> Não consegui rodar esta função:
> > cumhaz.summ <- function(y, digits = 3){
> + cond <- y$strata == levels(y$strata)[1]
> + y$cumhaz.upper <- pmax(y$cumhaz + 1.96 * y$std.chaz,0)
> + y$cumhaz.lower <- pmax(y$cumhaz - 1.96 * y$std.chaz,0)
> + output <- list()
> + for(i in seq_along(levels(y$strata))){
> + cond <- y$strata == levels(y$strata)[i]
> + output[[i]] <- round(data.frame(Time = y$time[cond],
> + 'N risk' = y$n.risk[cond],
> + 'N events' = y$n.event[cond],
> + 'Cum Hazard' = y$cumhaz[cond],
> + lower = y$cumhaz.lower[cond],
> + upper = y$cumhaz.upper[cond]),digits)
> + }
> + names(output) <- levels(y$strata)
> + output
> + }
> > cumhaz.summ(y)
> Error in data.frame(Time = y$time[cond], `N risk` = y$n.risk[cond], `N
> events` = y$n.event[cond], :
> arguments imply differing number of rows: 3, 10
> Além disso: Warning messages:
> 1: In y$cumhaz + 1.96 * y$std.chaz :
> comprimento do objeto maior não é múltiplo do comprimento do objeto menor
> 2: In y$cumhaz - 1.96 * y$std.chaz :
> comprimento do objeto maior não é múltiplo do comprimento do objeto menor
> >
>
> Na verdade eu estou fuçando ainda como plot.survfit faz pq não esta
> escrito em qualquer lugar que eu tenha achado como essa função faz e a
> summary.survfit retorna o erro padrão. Mas como a duvida era como extrair o
> cumhaz... achei que poderia postar. Se voce souber como a função plot faz
> posta aí que eu acerto.
>
> Pedro Brasil
>
> Em sex., 10 de dez. de 2021 às 13:32, Cesar Rabak <cesar.rabak em gmail.com>
> escreveu:
>
>> Esse cálculo dos limites superiores e inferiores do risco acumulado me
>> parece estranho..
>>
>> A formulação é baseada em alguma referência?
>>
>>
>> On Fri, Dec 10, 2021 at 1:23 PM Pedro Emmanuel Alvarenga Americano do
>> Brasil por (R-br) <r-br em listas.c3sl.ufpr.br> wrote:
>>
>>> > fit <- survfit(Surv(time, status) ~ x, data = aml)
>>> > y <- summary(fit, times = c(14,28,35))
>>> > y
>>> Call: survfit(formula = Surv(time, status) ~ x, data = aml)
>>>
>>> x=Maintained
>>> time n.risk n.event survival std.err lower 95% CI upper 95% CI
>>> 14 8 2 0.818 0.116 0.619 1.000
>>> 28 6 2 0.614 0.153 0.377 0.999
>>> 35 3 2 0.368 0.163 0.155 0.875
>>>
>>> x=Nonmaintained
>>> time n.risk n.event survival std.err lower 95% CI upper 95% CI
>>> 14 7 5 0.583 0.142 0.3616 0.941
>>> 28 4 2 0.389 0.147 0.1854 0.816
>>> 35 2 2 0.194 0.122 0.0569 0.664
>>>
>>> > cumhaz.summ <- function(y, digits = 3){
>>> + cond <- y$strata == levels(y$strata)[1]
>>> + y$cumhaz.upper <- pmax(y$cumhaz + 1.96 * y$std.chaz,0)
>>> + y$cumhaz.lower <- pmax(y$cumhaz - 1.96 * y$std.chaz,0)
>>> + output <- list()
>>> + for(i in seq_along(levels(y$strata))){
>>> + cond <- y$strata == levels(y$strata)[i]
>>> + output[[i]] <- round(data.frame(Time = y$time[cond],
>>> + 'N risk' = y$n.risk[cond],
>>> + 'N events' = y$n.event[cond],
>>> + 'Cum Hazard' = y$cumhaz[cond],
>>> + lower = y$cumhaz.lower[cond],
>>> + upper =
>>> y$cumhaz.upper[cond]),digits)
>>> + }
>>> + names(output) <- levels(y$strata)
>>> + output
>>> + }
>>> > cumhaz.summ(y)
>>> $`x=Maintained`
>>> Time N.risk N.events Cum.Hazard lower upper
>>> 1 14 8 2 0.191 0.000 0.456
>>> 2 28 6 2 0.459 0.002 0.915
>>> 3 35 3 2 0.909 0.133 1.685
>>>
>>> $`x=Nonmaintained`
>>> Time N.risk N.events Cum.Hazard lower upper
>>> 1 14 7 5 0.492 0.056 0.928
>>> 2 28 4 2 0.858 0.187 1.530
>>> 3 35 2 2 1.442 0.385 2.499
>>>
>>> >
>>>
>>> Em qui., 9 de dez. de 2021 às 10:43, Cid Póvoas por (R-br) <
>>> r-br em listas.c3sl.ufpr.br> escreveu:
>>>
>>>> library(survival)
>>>> library(broom)
>>>> library(tidyverse)
>>>>
>>>> fit <- survfit(Surv(time, status) ~ x, data = aml)
>>>> df<-tidy(fit)
>>>> df$cumhaz <- fit$cumhaz
>>>> df %>% filter(time%in%c(9,8,28,33,48))
>>>>
>>>> fit1 <- summary(fit, times = c(14,28,35))
>>>> fit1
>>>>
>>>> tab <- data.frame(time = fit1$time,
>>>> n.risk = fit1$n.risk,
>>>> n.event = fit1$n.event,
>>>> survival = fit1$surv,
>>>> std.err = fit1$std.err,
>>>> `lower 95% CI` = fit1$lower,
>>>> `upper 95% CI` = fit1$upper,
>>>> cumhaz = fit1$cumhaz,
>>>> strata = fit1$strata)
>>>>
>>>> tab
>>>>
>>>> *Cid Edson Mendonça Póvoas*
>>>>
>>>> *AnovAgro*
>>>> *Engenheiro Agrônomo - **Data Scientist*
>>>> *CREA : 051984991-4*
>>>> *Técnico em Segurança do Trabalho *
>>>> *Nº: **0012669/BA*
>>>> *Tel: +55 73 99151-9565*
>>>> *Lattes : *http://lattes.cnpq.br/2303498368142537
>>>> *LinkedIn :* http://br.linkedin.com/in/cidedson/
>>>> *Whatsapp :* https://wa.me/5573991519565
>>>>
>>>> Em qui., 9 de dez. de 2021 às 09:26, Pedro Emmanuel Alvarenga Americano
>>>> do Brasil por (R-br) <r-br em listas.c3sl.ufpr.br> escreveu:
>>>>
>>>>> Ei Cesar,
>>>>>
>>>>> Ei sei que o cumhaz esta la. Eu so queria uma saida parecida com a do
>>>>> survival. Parece que vou ter que montar uma um esquema aqui pegar esses
>>>>> valores e montar uma tabela de saida.
>>>>>
>>>>> Valeu.
>>>>>
>>>>> Pedro Emmanuel Brasil
>>>>> (:)=
>>>>>
>>>>> Em qua, 8 de dez de 2021 20:42, Cesar Rabak <cesar.rabak em gmail.com>
>>>>> escreveu:
>>>>>
>>>>>> Isto não é suficiente?
>>>>>>
>>>>>> > fit$cumhaz
>>>>>> [1] 0.09090909 0.19090909 0.31590909 0.45876623 0.45876623 0.65876623
>>>>>> [7] 0.90876623 0.90876623 1.40876623 1.40876623 0.16666667 0.36666667
>>>>>> [13] 0.49166667 0.49166667 0.65833333 0.85833333 1.10833333 1.44166667
>>>>>> [19] 1.94166667 2.94166667
>>>>>> >
>>>>>> HTH
>>>>>> --
>>>>>> Cesar Rabak
>>>>>>
>>>>>> On Wed, Dec 8, 2021 at 5:12 PM Pedro Emmanuel Alvarenga Americano do
>>>>>> Brasil por (R-br) <r-br em listas.c3sl.ufpr.br> wrote:
>>>>>>
>>>>>>> Saudações amigos do R,
>>>>>>>
>>>>>>> Estou as voltas de estimar taxas de eventos e estou batendo cabeça.
>>>>>>> Antes de escrever uma função eu mesmo para fazer uma tabela com alguns
>>>>>>> valores gostaria de uma luz dos amigos de R.
>>>>>>>
>>>>>>> O banco aml está no pacote survival então bastaria carregar o pacote
>>>>>>> pra reproduzir o exemplo. Eu gostaria de ver em formato de tabela alguns
>>>>>>> momentos específicos da tabela de sobrevivência. Só que o summary.survfit
>>>>>>> só retorna a sobrevivência e não a taxa cumulativa.
>>>>>>>
>>>>>>> library("survival")
>>>>>>> fit <- survfit(Surv(time, status) ~ x, data = aml)
>>>>>>> > fit
>>>>>>> Call: survfit(formula = Surv(time, status) ~ x, data = aml)
>>>>>>>
>>>>>>> n events median 0.95LCL 0.95UCL
>>>>>>> x=Maintained 11 7 31 18 NA
>>>>>>> x=Nonmaintained 12 11 23 8 NA
>>>>>>> # Summary com todos os momentos de eventos
>>>>>>> > summary(fit)
>>>>>>> Call: survfit(formula = Surv(time, status) ~ x, data = aml)
>>>>>>>
>>>>>>> x=Maintained
>>>>>>> time n.risk n.event survival std.err lower 95% CI upper 95% CI
>>>>>>> 9 11 1 0.909 0.0867 0.7541 1.000
>>>>>>> 13 10 1 0.818 0.1163 0.6192 1.000
>>>>>>> 18 8 1 0.716 0.1397 0.4884 1.000
>>>>>>> 23 7 1 0.614 0.1526 0.3769 0.999
>>>>>>> 31 5 1 0.491 0.1642 0.2549 0.946
>>>>>>> 34 4 1 0.368 0.1627 0.1549 0.875
>>>>>>> 48 2 1 0.184 0.1535 0.0359 0.944
>>>>>>>
>>>>>>> x=Nonmaintained
>>>>>>> time n.risk n.event survival std.err lower 95% CI upper 95% CI
>>>>>>> 5 12 2 0.8333 0.1076 0.6470 1.000
>>>>>>> 8 10 2 0.6667 0.1361 0.4468 0.995
>>>>>>> 12 8 1 0.5833 0.1423 0.3616 0.941
>>>>>>> 23 6 1 0.4861 0.1481 0.2675 0.883
>>>>>>> 27 5 1 0.3889 0.1470 0.1854 0.816
>>>>>>> 30 4 1 0.2917 0.1387 0.1148 0.741
>>>>>>> 33 3 1 0.1944 0.1219 0.0569 0.664
>>>>>>> 43 2 1 0.0972 0.0919 0.0153 0.620
>>>>>>> 45 1 1 0.0000 NaN NA NA
>>>>>>>
>>>>>>> # Summary com os momentos desejados
>>>>>>> > summary(fit, times = c(14,28,35))
>>>>>>> Call: survfit(formula = Surv(time, status) ~ x, data = aml)
>>>>>>>
>>>>>>> x=Maintained
>>>>>>> time n.risk n.event survival std.err lower 95% CI upper 95% CI
>>>>>>> 14 8 2 0.818 0.116 0.619 1.000
>>>>>>> 28 6 2 0.614 0.153 0.377 0.999
>>>>>>> 35 3 2 0.368 0.163 0.155 0.875
>>>>>>>
>>>>>>> x=Nonmaintained
>>>>>>> time n.risk n.event survival std.err lower 95% CI upper 95% CI
>>>>>>> 14 7 5 0.583 0.142 0.3616 0.941
>>>>>>> 28 4 2 0.389 0.147 0.1854 0.816
>>>>>>> 35 2 2 0.194 0.122 0.0569 0.664
>>>>>>>
>>>>>>> > plot(fit)
>>>>>>> > plot(fit, cumhaz = T)
>>>>>>> >
>>>>>>> Reparem que há uma opção para o gráfico cumhaz = T, o que significa
>>>>>>> que o cumhaz está depositado no objeto, inclusive dentro do summary também.
>>>>>>> Tipo summary(fit)$cumhaz. Só que não há uma opção summary(fit, cumhaz = T)
>>>>>>> que retorne o cumhaz ao invés da sobrevivência. Alguém tem algum bizu pra
>>>>>>> fazer isso organizado, extrair a mesma tabela só que com o cumhaz, como summary(fit,
>>>>>>> times = c(14,28,35)) sem muito trabalho?
>>>>>>>
>>>>>>> Abraço forte,
>>>>>>> Pedro Brasil
>>>>>>> _______________________________________________
>>>>>>> R-br mailing list
>>>>>>> R-br em listas.c3sl.ufpr.br
>>>>>>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>>>>>>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e
>>>>>>> forneça código mínimo reproduzível.
>>>>>>
>>>>>> _______________________________________________
>>>>> R-br mailing list
>>>>> R-br em listas.c3sl.ufpr.br
>>>>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>>>>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>>>>> código mínimo reproduzível.
>>>>
>>>> _______________________________________________
>>>> R-br mailing list
>>>> R-br em listas.c3sl.ufpr.br
>>>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>>>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>>>> código mínimo reproduzível.
>>>
>>> _______________________________________________
>>> R-br mailing list
>>> R-br em listas.c3sl.ufpr.br
>>> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
>>> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
>>> código mínimo reproduzível.
>>
>> _______________________________________________
> R-br mailing list
> R-br em listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
> código mínimo reproduzível._______________________________________________
> R-br mailing list
> R-br em listas.c3sl.ufpr.br
> https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br
> Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça
> código mínimo reproduzível.
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20211216/cf9f7eae/attachment.htm>
Mais detalhes sobre a lista de discussão R-br