Cara,

Primeiro um comentário a respeito do IC. Li umas paradas aqui que me levou a seguinte conclusão. O objeto y$std.chaz é identico ao objeto y$std.err. Então o cumhaz no y nada mais é que -log(y$surv). Dessa forma, pra achar o IC do cumhaz basta fazer o -log(y$lower) e -log(y$upper), independente do método utilizado, sendo que como tem o "-" na frente, o upper vira lower e vice versa. Está batendo aqui com os meus dados nas primeiras casas decimais. 
cumhaz.summ <- function(y, digits = 3){
y$cumhaz.upper  <- -log(y$lower)
y$cumhaz.lower  <- -log(y$upper)
      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)

A segunda parada é que essa função só foi testada para um preditor que estratifica a curva. No exemplo acima o "~ x". Possivelmente daria errado se não houver preditor (e dá errado mesmo no meu teste aqui retorna uma lista vazia pq o y$strata retorna NULL) ou mais de um preditor. Este deve ser o teu caso já que o erro dá a entender que o comprimento da condição está retornando objetos de tamanhos diferentes. Mas se colocar um exemplo reprodutível, talvez eu consiga te ajudar.  

Abraço forte,

Pedro Brasil
Pedro Brasil


Em qua., 15 de dez. de 2021 às 21:52, sznelwar--- por (R-br) <r-br@listas.c3sl.ufpr.br> escreveu:
Não consegui rodar esta função:
> 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)
Error in data.frame(Time = y$time[cond], `N risk` = y$n.risk[cond], `N events` = y$n.event[cond],  : 
  arguments imply differing number of rows: 3, 10
Além disso: Warning messages:
1: In y$cumhaz + 1.96 * y$std.chaz :
  comprimento do objeto maior não é múltiplo do comprimento do objeto menor
2: In y$cumhaz - 1.96 * y$std.chaz :
  comprimento do objeto maior não é múltiplo do comprimento do objeto menor
> 
 
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.
_______________________________________________
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.