[R-br] Extraindo o "cumulative hazard" de objeto survfit

Cesar Rabak cesar.rabak em gmail.com
Sex Dez 10 13:32:17 -03 2021


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 <http://www.anovagro.com/>*
>> *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.
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20211210/d591bef1/attachment.htm>


Mais detalhes sobre a lista de discussão R-br