График кумулятивной заболеваемости и таблица в R

Контекст: у меня есть набор данных о пациентах, перенесших трансплантацию печени (некоторые из них находятся в группе лечения А, другие — в группе лечения В, в зависимости от того, какие иммунодепрессанты они принимают). Из-за трансплантации эти пациенты подвержены высокому риску развития донорских HBV-инфекций.

Что необходимо: Исследователя интересует время до начала инфекции (первый случай ВГВ) и процентное соотношение, в котором инфекция развивается с течением времени. Им также нужна кумулятивная заболеваемость ВГВ-инфекцией на исходном уровне и в каждый из периодов наблюдения после трансплантации (6 месяцев, 12 месяцев, 18 месяцев и 24 месяца) для группы A и группы B. Например, 6- данные за месяц будут представлять собой долю тех пациентов с 6-месячным последующим наблюдением, которые когда-либо болели ВГВ, данные за 12 месяцев будут представлять собой долю пациентов с 12-месячным последующим наблюдением, у которых когда-либо был ВГВ, и так далее.

Кумулятивная заболеваемость в этом конкретном случае просто означает 1 минус функция выживания, без учета каких-либо конкурирующих рисков. Анализируемая популяция не имеет смертельных случаев или потерь для последующего наблюдения.

Мои вопросы:

Есть ли способ добавить количество событий прямо под номером риска на графике?
Есть ли способ также получить кумулятивные показатели заболеваемости в каждый момент времени для каждой группы с std.err и 95% ДИ, подобно таблицам смертности, которые мы получаем, когда используем сводку (км) ниже? эти таблицы дожития дают мне вероятности выживания, поэтому я думаю, если мне нужна кумулятивная заболеваемость, я мог бы просто вручную сделать 1-вероятность выживания, но не знаю, как получить std.err и доверительные интервалы?

Ниже приведен тестовый набор данных, похожий на настоящий, и то, что я сделал до сих пор:

time<-c(1.5989,6.9433, 0.8890, 3.2691, 1.0514, 2.7625, 1.4319, 0.9681, 7.4416, 0.0268, 1.5168, 1.9647, 0.0657, 4.3571, 6.4490, 0.2198, 1.2028, 0.9555, 0.2601, 2.0096, 7.5156, 0.4463, 0.2355, 0.9391, 2.6996)

censor<-c(1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0)

group<-c(1, 2, 1, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2, 1, 2, 1)

df<-data.frame(ID, time, censor, group)
View(df)

km<-survfit(formula = Surv(time, censor) ~ group, data = df)
summary(km)

#cumulative incidence plot
plot(km, fun = function(x) 1-x)

#log rank test;
survdiff(Surv(time, censor) ~ group, data=df)
#plot survival curves for each treatment group
#cumulative incidence plot 
ggsurvplot(km,
           data = df,
           censor = T,
           risk.table = TRUE,
           legend.labs = c("group 1", "group 2"),
           xlim = c(0,10),
           ylim = c(0,1),
           pval = T,
           pval.method = T,
           pval.method.coord = c(2.5,0.5),
           pval.coord = c(4.2,0.5),
           xlab = "Months",
           ylab = "SURVIVAL PROBABILITY",
           linetype = c(1,2),
           legend.title = "",
           palette = c('red', 'blue'),
           fun = "event"
)

r survival

13.02.2023 06:49

Типы данных JavaScript

В JavaScript существует несколько типов данных, включая примитивные типы данных и ссылочные типы данных. Вот краткое объяснение различных типов данных...

Как сделать движок для футбольного матча? (простой вариант)

Футбол. Для многих людей, живущих на земле, эта игра - больше, чем просто спорт. И эти люди всегда мечтают стать футболистом или менеджером. Но, к...

Знайте свои исключения!

В Java исключение - это событие, возникающее во время выполнения программы, которое нарушает нормальный ход выполнения инструкций программы. Когда...

CSS Flex: что должен знать каждый разработчик

CSS Flex: что должен знать каждый разработчик Модуль flexbox, также известный как гибкий модуль разметки box, помогает эффективно проектировать и...

Введение в раздел "Заголовок" в HTML

Говорят, что лучшее о человеке можно увидеть только изнутри, и это относится и к веб-страницам HTML! Причина, по которой некоторые веб-страницы не...

Что такое CommonModule в Angular

Перейти к ответу Данный вопрос помечен как решенный

Ответы 1

Ответ принят как подходящий

Вы можете создать собственную таблицу, а затем удалить ggsurvplot таблицу по умолчанию и заменить ее своей пользовательской таблицей.

Как указано в Примечании о конкурирующих рисках при анализе данных о выживании статье в British Journal of Cancer о совокупной частоте смертности:

Это просто обратная сторона выживания. Другими словами, кумулятивная частота события в данный момент времени равна единице минус общая вероятность выживания в это время.

Таким образом, для кумулятивной частоты событий, как вы сказали, мы можем использовать 1-survival_probalility, поэтому их стандартная ошибка такая же, как и их доверительный интервал (хотя более низкий ci как верхний ci и верхний ci как нижний ci).

На пользовательском графике здесь я добавил количество рисков (первая строка), кумулятивное количество событий (вторая строка), вероятность выживания (третья строка) и ее доверительные интервалы (четвертая строка).

libs <- c("survminer", "survival", "tidyverse")
invisible(suppressMessages(sapply(libs,library,character.only=T)))

time <- c(1.5989,6.9433, 0.8890, 3.2691, 1.0514, 2.7625, 1.4319, 0.9681, 7.4416, 
          0.0268, 1.5168, 1.9647, 0.0657, 4.3571, 6.4490, 0.2198, 1.2028, 0.9555, 
          0.2601, 2.0096, 7.5156, 0.4463, 0.2355, 0.9391, 2.6996)

censor <- c(1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 0)

group <- c(1, 2, 1, 1, 2, 2, 1, 1, 1, 2, 1, 2, 1, 1, 2, 1, 2, 1, 1, 2, 2, 2, 1, 2, 1)

df <- data.frame(time, censor, group)

km <- survfit(formula = Surv(time, censor) ~ group, data = df)

#Draw cumulative incidence plot as usual
p1 <- ggsurvplot(km,
           data = df,
           censor = T,
           risk.table = T,
           legend.labs = c("group 1", "group 2"),
           xlim = c(0,7.5),
           ylim = c(0,1),
           pval = T,
           pval.method = T,
           pval.method.coord = c(2.5,0.5),
           pval.coord = c(4.2,0.5),
           xlab = "Months",
           ylab = "SURVIVAL PROBABILITY",
           linetype = c(1,2),
           legend.title = "",
           palette = c('red', 'blue'),
           fun = "event",
           cumevents = F, 
           tables.height = 0.4
) 

# extract data table of survival plot

p1_df <- p1$table$data
p1_df
#>     strata time n.risk pct.risk n.event cum.n.event n.censor cum.n.censor
#> 1  group 1    0     14      100       0           0        0            0
#> 2  group 1    2      4       29      10          10        0            0
#> 3  group 1    4      2       14       0          10        2            2
#> 4  group 1    6      1        7       0          10        1            3
#> 5  group 1    8      0        0       0          10        1            4
#> 6  group 2    0     11      100       0           0        0            0
#> 7  group 2    2      5       45       6           6        0            0
#> 8  group 2    4      3       27       0           6        2            2
#> 9  group 2    6      3       27       0           6        0            2
#> 10 group 2    8      0        0       0           6        3            5
#>    strata_size group llabels
#> 1           14     1      14
#> 2           14     1       4
#> 3           14     1       2
#> 4           14     1       1
#> 5           14     1       0
#> 6           11     2      11
#> 7           11     2       5
#> 8           11     2       3
#> 9           11     2       3
#> 10          11     2       0


# use summary of km survfit in time points of survival plot

sum_df <-summary(km, times = p1_df$time) 
sum_df
#> Call: survfit(formula = Surv(time, censor) ~ group, data = df)
#> 
#>                 group=1 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0     14       0    1.000   0.000        1.000        1.000
#>     0     14       0    1.000   0.000        1.000        1.000
#>     2      4      10    0.286   0.121        0.125        0.654
#>     2      4       0    0.286   0.121        0.125        0.654
#>     4      2       0    0.286   0.121        0.125        0.654
#>     4      2       0    0.286   0.121        0.125        0.654
#>     6      1       0    0.286   0.121        0.125        0.654
#>     6      1       0    0.286   0.121        0.125        0.654
#> 
#>                 group=2 
#>  time n.risk n.event survival std.err lower 95% CI upper 95% CI
#>     0     11       0    1.000    0.00        1.000        1.000
#>     0     11       0    1.000    0.00        1.000        1.000
#>     2      5       6    0.455    0.15        0.238        0.868
#>     2      5       0    0.455    0.15        0.238        0.868
#>     4      3       0    0.455    0.15        0.238        0.868
#>     4      3       0    0.455    0.15        0.238        0.868
#>     6      3       0    0.455    0.15        0.238        0.868
#>     6      3       0    0.455    0.15        0.238        0.868
str(sum_df)
#> List of 19
#>  $ n            : int [1:2] 14 11
#>  $ time         : num [1:16] 0 0 2 2 4 4 6 6 0 0 ...
#>  $ n.risk       : num [1:16] 14 14 4 4 2 2 1 1 11 11 ...
#>  $ n.event      : num [1:16] 0 0 10 0 0 0 0 0 0 0 ...
#>  $ n.censor     : num [1:16] 0 0 0 0 2 0 1 0 0 0 ...
#>  $ surv         : num [1:16] 1 1 0.286 0.286 0.286 ...
#>  $ std.err      : num [1:16] 0 0 0.121 0.121 0.121 ...
#>  $ cumhaz       : num [1:16] 0 0 1.17 1.17 1.17 ...
#>  $ std.chaz     : num [1:16] 0 0 0.39 0.39 0.39 ...
#>  $ strata       : Factor w/ 2 levels "group=1","group=2": 1 1 1 1 1 1 1 1 2 2 ...
#>  $ type         : chr "right"
#>  $ logse        : logi TRUE
#>  $ conf.int     : num 0.95
#>  $ conf.type    : chr "log"
#>  $ lower        : num [1:16] 1 1 0.125 0.125 0.125 ...
#>  $ upper        : num [1:16] 1 1 0.654 0.654 0.654 ...
#>  $ call         : language survfit(formula = Surv(time, censor) ~ group, data = df)
#>  $ table        : num [1:2, 1:9] 14 11 14 11 14 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:2] "group=1" "group=2"
#>   .. ..$ : chr [1:9] "records" "n.max" "n.start" "events" ...
#>  $ rmean.endtime: num [1:2] 7.52 7.52
#>  - attr(*, "class")= chr "summary.survfit"


# extract our columns of interest from `sum_df` dataframe

columns <- lapply(c(2:10,15:16), function(x) sum_df[x])
tbl <- do.call(data.frame, columns)
tbl
#>    time n.risk n.event n.censor      surv   std.err   cumhaz  std.chaz  strata
#> 1     0     14       0        0 1.0000000 0.0000000 0.000000 0.0000000 group=1
#> 2     0     14       0        0 1.0000000 0.0000000 0.000000 0.0000000 group=1
#> 3     2      4      10        0 0.2857143 0.1207363 1.168229 0.3903649 group=1
#> 4     2      4       0        0 0.2857143 0.1207363 1.168229 0.3903649 group=1
#> 5     4      2       0        2 0.2857143 0.1207363 1.168229 0.3903649 group=1
#> 6     4      2       0        0 0.2857143 0.1207363 1.168229 0.3903649 group=1
#> 7     6      1       0        1 0.2857143 0.1207363 1.168229 0.3903649 group=1
#> 8     6      1       0        0 0.2857143 0.1207363 1.168229 0.3903649 group=1
#> 9     0     11       0        0 1.0000000 0.0000000 0.000000 0.0000000 group=2
#> 10    0     11       0        0 1.0000000 0.0000000 0.000000 0.0000000 group=2
#> 11    2      5       6        0 0.4545455 0.1501314 0.736544 0.3072801 group=2
#> 12    2      5       0        0 0.4545455 0.1501314 0.736544 0.3072801 group=2
#> 13    4      3       0        2 0.4545455 0.1501314 0.736544 0.3072801 group=2
#> 14    4      3       0        0 0.4545455 0.1501314 0.736544 0.3072801 group=2
#> 15    6      3       0        0 0.4545455 0.1501314 0.736544 0.3072801 group=2
#> 16    6      3       0        0 0.4545455 0.1501314 0.736544 0.3072801 group=2
#>        lower     upper
#> 1  1.0000000 1.0000000
#> 2  1.0000000 1.0000000
#> 3  0.1248055 0.6540791
#> 4  0.1248055 0.6540791
#> 5  0.1248055 0.6540791
#> 6  0.1248055 0.6540791
#> 7  0.1248055 0.6540791
#> 8  0.1248055 0.6540791
#> 9  1.0000000 1.0000000
#> 10 1.0000000 1.0000000
#> 11 0.2379221 0.8684002
#> 12 0.2379221 0.8684002
#> 13 0.2379221 0.8684002
#> 14 0.2379221 0.8684002
#> 15 0.2379221 0.8684002
#> 16 0.2379221 0.8684002

# create custom table of our metrics and parameters of interests
custom_tbl <- data.frame(time=tbl$time, 
                         n.risk=tbl$n.risk, 
                         n.event=tbl$n.event,
                         cum.inc=1-tbl$surv, 
                         std.err=1-tbl$std.err,
                         ci=paste0(round((1-tbl$upper),2),"-",round((1-tbl$lower),2)),
                         strata=tbl$strata) %>% 
  mutate_if (is.numeric, round, digits=2)


# calculating cumulative incidence from n.event column, then remove duplicated rows and 
# replace them by last similar column

custom_tbl %<>% group_by(strata) %>% mutate(cum.n.event = cumsum(n.event)) %>% 
  group_by(time, strata) %>% 
  slice_tail(n=1)

custom_tbl
#> # A tibble: 8 × 8
#> # Groups:   time, strata [8]
#>    time n.risk n.event cum.inc std.err ci        strata  cum.n.event
#>   <dbl>  <dbl>   <dbl>   <dbl>   <dbl> <chr>     <fct>         <dbl>
#> 1     0     14       0    0       1    0-0       group=1           0
#> 2     0     11       0    0       1    0-0       group=2           0
#> 3     2      4       0    0.71    0.88 0.35-0.88 group=1          10
#> 4     2      5       0    0.55    0.85 0.13-0.76 group=2           6
#> 5     4      2       0    0.71    0.88 0.35-0.88 group=1          10
#> 6     4      3       0    0.55    0.85 0.13-0.76 group=2           6
#> 7     6      1       0    0.71    0.88 0.35-0.88 group=1          10
#> 8     6      3       0    0.55    0.85 0.13-0.76 group=2           6

# removing default risk table of our survival plot
p1$table <- NULL

# creating a custom ggplot by our custom table
p1$table <- ggplot(custom_tbl, aes(time, strata)) +
  geom_text(aes(label = n.risk), size = 3, vjust=-0.8) +
  geom_text(aes(label = cum.n.event), size = 3, vjust=1) + 
  geom_text(aes(label = cum.inc), size = 3, vjust=2.8) + 
  geom_text(aes(label = ci), size = 3, vjust=4) + 
  labs(x = NULL, y = NULL) + theme_classic()+
  coord_cartesian(xlim = c(0, 7.5))+
  scale_x_continuous(breaks=seq(0, 9, 2))+
  scale_y_discrete(
    labels=c("group 1", "group 2"))+
  theme(axis.text.y = element_text(colour=c("red", "blue"), size = 12),
        axis.text.x = element_text(size = 12))
#> Warning: Vectorized input to `element_text()` is not officially supported.
#> ℹ Results may be unexpected or may change in future versions of ggplot2.

p1

Created on 2023-02-13 with reprex v2.0.2

Большое спасибо, это именно то, что мне было нужно! Статья, которой вы поделились, также очень полезна для разъяснения некоторых вещей, которые меня немного смутили.

— 13.02.2023 21:06

Если этот ответ был вам полезен, подпишите его как принятый ответ и проголосуйте за него.

— 13.02.2023 21:21

13.02.2023 15:48