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

Pedro Emmanuel Alvarenga Americano do Brasil emmanuel.brasil em gmail.com
Sáb Dez 11 11:56:44 -03 2021


Cesar,
Não sei se entendi direito a tua colocação. Tanto as funções survfit como a
summary.survfit retornam o erro padrão dentro do objeto. Então, de qualquer
objeto survfit, se pedirmos um plot com cumhaz = T e conf.int = T o grafico
vem com a banda de confiança do hazard.

lsurv2 <- survfit(Surv(time, status) ~ x, aml, type='fleming')
plot(lsurv2, lty=2:3, fun="cumhaz", xlab="Months", ylab="Cumulative
Hazard", conf.int = T)

Funçando por achei o seguinte...
https://www.rdocumentation.org/packages/survival/versions/3.2-13/topics/survfit.formula

conf.type

One of "none", "plain", "log" (the default), "log-log", "logit" or "arcsin".
Only enough of the string to uniquely identify it is necessary. The first
option causes confidence intervals not to be generated. The second causes
the standard intervals curve +- k *se(curve), where k is determined from
conf.int. The log option calculates intervals based on the cumulative
hazard or log(survival). The log-log option bases the intervals on the log
hazard or log(-log(survival)), the logit option on
log(survival/(1-survival)) and arcsin on arcsin(survival).
O que significa que das varias formas de se estimar a banda, uma delas é a
"plain". Ainda não tenho certeza se o que ele chama de curve já é a cumhaz
ou se é S(t). Vou ter que fuçar melhor nas referencias.

Pedro Brasil


Em sex., 10 de dez. de 2021 às 20:38, Cesar Rabak <cesar.rabak em gmail.com>
escreveu:

> Pedro,
>
> A função plot.surfit não coloca os erros padrão da função cumulativa.
>
> Ademais, pegue os resultados do summary e veja se se o cálculo do
> intervalo do CI 95% "bate" com a conta 1,96 × d.p.
>
> HTH
>
> --
> Cesar Rabak
>
>
> On Fri, Dec 10, 2021 at 4:02 PM Pedro Emmanuel Alvarenga Americano do
> Brasil <emmanuel.brasil em gmail.com> wrote:
>
>> 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 <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/20211211/7e3904a4/attachment-0001.htm>


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