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

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@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@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@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

Engenheiro Agrônomo - Data Scientist 
CREA : 051984991-4
Técnico em Segurança do Trabalho 
Nº: 0012669/BA
Tel: +55 73 99151-9565


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.rabak@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.