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.bra...@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.ra...@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.bra...@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.ra...@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@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@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@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.ra...@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@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@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@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@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@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@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.