[R-br] Extraindo o "cumulative hazard" de objeto survfit
Cesar Rabak
cesar.rabak em gmail.com
Sáb Dez 11 14:11:36 -03 2021
Pedro,
Em retrospecto, vejo que não fui muito claro no que escrevi...
O que quis dizer foi que (segundo sua primeira postagem), a chamada:
> plot(fit, cumhaz = T)
não inclui o intervalo de confiança do risco acumulado.
Secundariamente, o cálculo do desvio padrão para esse resultado estatístico
me parecia incorreto, e usando os dados que já haviam sido lançados na sua
1ª postagem pode-se ver que nem para cada ponto de risco (no tempo) a
fórmula valeria.
Com o avanço da sua pesquisa, estou seguro que verá as formas corretas de
calcular esses IC.
Por ora estou sem acesso à minhas referências em polpa de celulose e com
acesso bastante restrito à Web para encontrar algo mais assertivo.
HTH
--
Cesar Rabak
On Sat, Dec 11, 2021 at 11:56 AM Pedro Emmanuel Alvarenga Americano do
Brasil <emmanuel.brasil em gmail.com> wrote:
> 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/189a6cd0/attachment.htm>
Mais detalhes sobre a lista de discussão R-br