CÀI ĐẶT VÀ SỬ DỤNG PACKAGE SURVIVAL

Để sử dụng package survival trên RStudio, bạn có thể thực hiện các bước sau:

Bước 1: Cài đặt package survival

Bạn có thể cài đặt package survival bằng cách chạy lệnh sau trong RStudio:

install.packages("survival")

Nếu package đã được cài đặt trước đó, bạn có thể bỏ qua bước này.

Bước 2: Tải package survival

Sau khi cài đặt, bạn cần tải package survival bằng lệnh sau:

library(survival)

Bước 3: Sử dụng các hàm của package survival

Sau khi đã tải package survival, bạn có thể sử dụng các hàm của package để phân tích dữ liệu thời gian đến sự kiện. Các hàm này bao gồm:

survfit(): tính toán các ước lượng sống sót bằng phương pháp Kaplan-Meier hoặc phương pháp Nelson-Aalen.
coxph(): phân tích mô hình hồi quy Cox để tìm các yếu tố ảnh hưởng đến thời gian đến sự kiện.
survreg(): phân tích mô hình hồi quy Cox để tìm các yếu tố ảnh hưởng đến thang đo thời gian đến sự kiện không phải là dạng số nguyên (ví dụ: thời gian sống sót được đo bằng giờ, ngày, tháng,…).
cox.zph(): kiểm định giả thuyết về sự đồng nhất của các hệ số hồi quy Cox theo thời gian.

Gọi package

library(survival)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

Xử lý dữ liệu

cancer <- data(cancer)
head(lung, 6)

##   inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1    3  306      2  74   1       1       90       100     1175      NA
## 2    3  455      2  68   1       0       90        90     1225      15
## 3    3 1010      1  56   1       0       90        90       NA      15
## 4    5  210      2  57   1       1       90        60     1150      11
## 5    1  883      2  60   1       0      100        90       NA       0
## 6   12 1022      1  74   1       1       50        80      513       0

lung <- 
  lung %>% 
  mutate(
    status = recode(status, `1` = 0, `2` = 1)
  )
head(lung,6)

##   inst time status age sex ph.ecog ph.karno pat.karno meal.cal wt.loss
## 1    3  306      1  74   1       1       90       100     1175      NA
## 2    3  455      1  68   1       0       90        90     1225      15
## 3    3 1010      0  56   1       0       90        90       NA      15
## 4    5  210      1  57   1       1       90        60     1150      11
## 5    1  883      1  60   1       0      100        90       NA       0
## 6   12 1022      0  74   1       1       50        80      513       0

head(lung[, c("time", "status", "sex")])

##   time status sex
## 1  306      1   1
## 2  455      1   1
## 3 1010      0   1
## 4  210      1   1
## 5  883      1   1
## 6 1022      0   1

Tạo đối tượng

Surv(lung$time, lung$status)[1:10]

##  [1]  306   455  1010+  210   883  1022+  310   361   218   166

s1 <- survfit(Surv(time, status) ~ 1, data = lung)
s1

## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
## 
##        n events median 0.95LCL 0.95UCL
## [1,] 228    165    310     285     363

Tạo đường cong sinh tồn

library(ggsurvfit)

## Loading required package: ggplot2

survfit2(Surv(time, status) ~ 1, data = lung) %>% 
  ggsurvfit() +
  labs(
    x = "Days",
    y = "Overall survival probability"
  )

survfit2(Surv(time, status) ~ 1, data = lung) %>% 
  ggsurvfit() +
  labs(
    x = "Days",
    y = "Overall survival probability"
  ) + 
  add_confidence_interval()

survfit2(Surv(time, status) ~ 1, data = lung) %>% 
  ggsurvfit() +
  labs(
    x = "Days",
    y = "Overall survival probability"
    ) + 
  add_confidence_interval() +
  add_risktable()

1 năm

summary(survfit(Surv(time, status) ~ 1, data = lung), times = 365.25)

## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   365     65     121    0.409  0.0358        0.345        0.486

library(gtsummary)
survfit(Surv(time, status) ~ 1, data = lung) %>% 
  tbl_survfit(
    times = 365.25,
    label_header = "**1-year survival (95% CI)**"
  )

Characteristic	1-year survival (95% CI)
Overall	41% (34%, 49%)

Ước tính thời gian sống sót tb

survfit(Surv(time, status) ~ 1, data = lung)

## Call: survfit(formula = Surv(time, status) ~ 1, data = lung)
## 
##        n events median 0.95LCL 0.95UCL
## [1,] 228    165    310     285     363

lung %>% 
  filter(status == 1) %>% 
  summarize(median_surv = median(time))

##   median_surv
## 1         226

survfit(Surv(time, status) ~ 1, data = lung) %>% 
  tbl_survfit(
    probs = 0.5,
    label_header = "**Median survival (95% CI)**")

Characteristic	Median survival (95% CI)
Overall	310 (285, 363)

So sánh giữa các nhóm

survdiff(Surv(time, status) ~ sex, data = lung)

## Call:
## survdiff(formula = Surv(time, status) ~ sex, data = lung)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=1 138      112     91.6      4.55      10.3
## sex=2  90       53     73.4      5.68      10.3
## 
##  Chisq= 10.3  on 1 degrees of freedom, p= 0.001

Fit complex survival curves

require("survival")
fit2 <- survfit( Surv(time, status) ~ sex,
                data = lung )

library("survminer")

## Loading required package: ggpubr

## 
## Attaching package: 'survminer'

## The following object is masked _by_ '.GlobalEnv':
## 
##     myeloma

## The following object is masked from 'package:survival':
## 
##     myeloma

ggsurv <- ggsurvplot(fit2, fun = "event", conf.int = TRUE,
                     ggtheme = theme_bw())
ggsurv$plot

fit3 <- survfit( Surv(time, status) ~ age,
                data = lung )
ggsurv2 <- ggsurvplot(fit3, fun = "event", conf.int = TRUE,
                     ggtheme = theme_bw())
ggsurv2$plot

UntitledT7

Khánh An

2023-07-11