class: title-slide, left, bottom # 통계실습 II. 생존분석 ---- ## **가천대 길병원 고급통계교육** ### 방태모 ### 2022-06-21 --- # .green[목차] .left-column[ ] .right-column[ - .content-box-green[기초개념] - .content-box-green[생존 함수 추정] - Kaplan-Meier method - 실습 - .content-box-green[생존 함수 비교] - Log rank test - 실습 - .content-box-green[위험인자 분석] - Cox proportional hazard model - Time dependent Cox proportional hazard model - 실습 <!-- 위험인자 분석 또한 생존함수 추정방법에 해당 --> <!-- 위험함수 h(t)와 생존함수 S(t) 사이에는 관계식이 있고 하나를 알면, 하나를 유도할수 있음 h(t) = f(t)/S(t) --> ] --- class: inverse, middle, center # 기초개념 --- ## .green[생존 분석] - 생존자료(survival data)를 분석하는 통계적 방법 ## .green[생존 자료] - 관심 사건의 발생 여부(이진형 변수)와 생존시간(연속형 변수)을 관측한 자료 - 생존자료를 사건-시간 자료(time-to-event data)라 표현하기도 함 - 실제 생존 시간(true survival time, .red[T]) = 추적 시작 시점 ~ 사건 발생 지점 <!-- 관심사건은 사망, 질병 발생/재발/호전, 금연, 퇴사 등이 될 수 있음--> <!-- 단위: days, weeks, months, years --> - 생존 자료의 가장 큰 특징: 중도 절단(censored) - 중도 절단 시간 (censoring time, .green[C]) - 연구기간 내 사건 미발생 시 마지막 추적 시점까지의 시간 - 관측 생존 시간 (observed survival time) = min(.red[T], .green[C]) - 연구기간 내 생존 시간 또는 중도 절단 시간까지의 시간 <!-- 즉, 관측 생존 시간 `\(\neq\)` 실제 생존 시간 --> --- ## .green[생존 자료의 특징] .pull-left[ - 생존 시간은 정규분포를 따르지 않음 - 항상 양의 값 - 매우 치우친 분포 - 중도 절단 - Type I censoring - Type II censoring - Type III censoring - loss to follow up - drop out <!-- 사전에 중도절단시간 고정 e.g. 연구종료시점 정해놓은 경우 --> <!-- 사전에 관측 사망자수 설정 --> <!-- 임의중도절단 --> - 이진형 변수와 연속형 변수를 다루어야 함 - .green[생존 분석] 필요 ] .pull-right[ <div class="figure"> <img src="./fig/example.png" alt="fig 1. Survival data (Jim Gruman)" width="750" height="450" /> <p class="caption">fig 1. Survival data (Jim Gruman)</p> </div> ] <!-- right censored - 사건 발생을 관측하지 못한 경우 (연구 종료) --> <!-- left censored - 추적 시작 시점을 관측하지 못한 경우 (아주 드문 형태) --> <!-- interval censored - 연구 기간 중 일부를 관측하지 못한 경우 (아주 드문 형태) --> --- ## .green[생존 분석의 목적] .left-column[ ] .right-column[ - 추정(estimation) - 생존시간의 분포에 관한 정보 - .green[Kaplan-Meier Method] - 그룹간 비교(group comparison) - 실험군과 대조군이 있는 경우 두 처리군의 생존 분포 비교 - .green[Log-rank test] - 위험인자 분석(risk-factor analysis) - 유의한 예후인자(prognostic factor) 파악 - .green[Cox ph model] - .green[Time dependent Cox ph model] ] --- class: inverse, middle, center # 생존 함수 추정 ## Kaplan-Meier method --- ## .green[생존 함수 S(t)] .pull-left[ <!-- 생존시간을 T라 했을 때 --> `\begin{align} S(t) = Pr(T > t) \end{align}` - 관심사건이 `\(t\)` 시점까지 일어나지 않을 확률 - 관심사건: 사망 - 환자가 `\(t\)`시간 이상 생존할 확률 - `\(S(t)\)`: 생존율, 생존 곡선 - `\(F(t) = 1 - S(t)\)`: 누적 발생률 <!-- 생존율(survival rate), 생존곡선(surviaval curve) --> <!-- 누적 발생률(cumulative incidence rate) --> - e.g. `\(S(10)\)` = 0.1 ( `\(t\)` in years) - 10년 생존율 = `\(S(10)\)` = 10% - 10년 누적 사망률 = `\(F(10)\)` = 90% ] .pull-right[ <div class="figure"> <img src="./fig/survival_function.png" alt="fig 2. Survival function (Wikipedia)" width="450" height="450" /> <p class="caption">fig 2. Survival function (Wikipedia)</p> </div> ] --- ## .green[Kaplan-Meier method] .pull-left[ - 생존함수를 추정하는 대표적인 비모수적 방법론: <!-- 생존 시간의 분포를 특정 분포로 가정하지 않음 --> <!-- 또 다른 비모수적방법론: 생명표방법. 환자 수 상당히 많을때 유용. 개개의 생존시간을 고려하는 대신 생존시간을 몇 개의 구간으로 나누어 생존자료를 요약하는 방법. 보험통계에 많이 쓰임. --> <!-- 모수적 방법론: 지수분포, 곰페르츠분포, 와이블분포. 산업계에서 선호 --> <!-- 생존시간의 분포형태가 알려져 있을 경우에는 비모수적 방법보다 모수적 방법 이용하는 것이 더 효율적 --> `\begin{align} \hat{S}(t_j) = \prod_{i=1}^{j}(1- \frac{d_i}{n_i})^{\delta_i} \end{align}` - `\(t_j\)`: `\(j\)`-번째 사망 시간 - `\(d_i\)`: `\(t_i\)` 시점 사망자 수 - `\(n_i\)`: `\(t_i\)` 시점 바로 직전 생존자 수 - `\(\delta_i\)`: 중도절단여부(중도절단 = 0, 아니면 1) - 즉, 사건 발생 시점에서만 생존율 변화 ] .pull-right[ <div class="figure"> <img src="./fig/km_method.png" alt="fig 3. KM method (Jim Gruman)" width="550" height="300" /> <p class="caption">fig 3. KM method (Jim Gruman)</p> </div> ] --- ## .green[실습 - Kaplan-Meier method] .scroll-output[ ### .black[패키지 설치] ```r install.packages("dplyr") install.packages("survival") install.packages("survminer") ``` ### .black[패키지 로딩] ```r library(dplyr) library(survival) library(survminer) ``` ### .black[데이터 소개] - 228명 진행성 폐암 환자들의 생존 시간(`time`), 사망 여부(`status`) 데이터 - 1022일간 추적 연구(follow-up study) ```r glimpse(lung) ``` ``` ## Rows: 228 ## Columns: 10 ## $ inst <dbl> 3, 3, 3, 5, 1, 12, 7, 11, 1, 7, 6, 16, 11, 21, 12, 1, 22, 16… ## $ time <dbl> 306, 455, 1010, 210, 883, 1022, 310, 361, 218, 166, 170, 654… ## $ status <dbl> 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, … ## $ age <dbl> 74, 68, 56, 57, 60, 74, 68, 71, 53, 61, 57, 68, 68, 60, 57, … ## $ sex <dbl> 1, 1, 1, 1, 1, 1, 2, 2, 1, 1, 1, 2, 2, 1, 1, 1, 1, 1, 2, 1, … ## $ ph.ecog <dbl> 1, 0, 0, 1, 0, 1, 2, 2, 1, 2, 1, 2, 1, NA, 1, 1, 1, 2, 2, 1,… ## $ ph.karno <dbl> 90, 90, 90, 90, 100, 50, 70, 60, 70, 70, 80, 70, 90, 60, 80,… ## $ pat.karno <dbl> 100, 90, 90, 60, 90, 80, 60, 80, 80, 70, 80, 70, 90, 70, 70,… ## $ meal.cal <dbl> 1175, 1225, NA, 1150, NA, 513, 384, 538, 825, 271, 1025, NA,… ## $ wt.loss <dbl> NA, 15, 15, 11, 0, 0, 10, 1, 16, 34, 27, 23, 5, 32, 60, 15, … ``` ] --- ## .green[실습 - Kaplan-Meier method] .scroll-output[ ### .black[생존 자료 객체 만들기] ```r surv_lung <- Surv(lung$time, lung$status)[1:10] head(surv_lung) ``` ``` ## [1] 306 455 1010+ 210 883 1022+ ``` ### .black[Kaplan-Meier method을 이용한 생존 곡선 추정] ```r method_km <- survfit(Surv(time, status) ~ 1, data = lung) glimpse(method_km) ``` ``` ## List of 16 ## $ n : int 228 ## $ time : num [1:186] 5 11 12 13 15 26 30 31 53 54 ... ## $ n.risk : num [1:186] 228 227 224 223 221 220 219 218 217 215 ... ## $ n.event : num [1:186] 1 3 1 2 1 1 1 1 2 1 ... ## $ n.censor : num [1:186] 0 0 0 0 0 0 0 0 0 0 ... ## $ surv : num [1:186] 0.996 0.982 0.978 0.969 0.965 ... ## $ std.err : num [1:186] 0.0044 0.00885 0.00992 0.01179 0.01263 ... ## $ cumhaz : num [1:186] 0.00439 0.0176 0.02207 0.03103 0.03556 ... ## $ std.chaz : num [1:186] 0.00439 0.0088 0.00987 0.01173 0.01257 ... ## $ type : chr "right" ## $ logse : logi TRUE ## $ conf.int : num 0.95 ## $ conf.type: chr "log" ## $ lower : num [1:186] 0.987 0.966 0.959 0.947 0.941 ... ## $ upper : num [1:186] 1 1 0.997 0.992 0.989 ... ## $ call : language survfit(formula = Surv(time, status) ~ 1, data = lung) ## - attr(*, "class")= chr "survfit" ``` - `time`: unique한 생존 기간 및 중도 절단 기간 - `n.event`: 해당 기간의 관심사건 발생 수 - `surv`: 각 `time`의 생존율 ] --- ## .green[실습 - Kaplan-Meier method] ### .black[Kaplan-Meier plot 시각화] - {survminer} 패키지의 `ggsurvplot()` 이용 - 🔗 [cheatsheet](https://rpkgs.datanovia.com/survminer/survminer_cheatsheet.pdf) ```r ggsurvplot(method_km, xlab = "Days", ylab = "Overall survival probability") ``` <img src="2_survival_files/figure-html/unnamed-chunk-9-1.png" width="504" style="display: block; margin: auto;" /> <!-- 그래프의 tick 마크는 censored를 나타냄 --> --- ## .green[실습 - Kaplan-Meier method] ### .black[1년 생존율 추정하기] .scroll-box-20[ - 추정된 KM curve에 의하면 말기 암환자의 추정된 1년 생존율은 41%에 해당: ```r summary(method_km, 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 ``` - censoring 고려없이 해당 시점의 사망자 수만 고려하여 생존율을 계산하면 생존율이 과대추정됨: `\begin{align} (1-\frac{121}{228}) \times 100 = 47 \% \end{align}` - 해당 시점에는 42명의 censored된 환자가 존재하고, 이를 반영하여 KM method에 의해 더 정확하게 추정된 생존율은 41%에 해당 ] --- ## .green[실습 - Kaplan-Meier method] ### .black[1년 생존율 추정하기 - 시각화] - Number at risk: 해당 시점의 생존자 수 <img src="2_survival_files/figure-html/unnamed-chunk-11-1.png" width="504" style="display: block; margin: auto;" /> - 총 환자 수 228명 - 사망 121명 - 중도절단 42명 = 65명 생존 (1년 생존자 수) --- ## .green[실습 - Kaplan-Meier method] ### .black[1년 생존율 추정하기 - 시각화 코드] .scroll-box-20[ ```r plot_main <- ggsurvplot( data = lung, fit = method_km, xlab = "Months", legend = "none", xscale = 30.4, break.x.by = 182.4, risk.table = TRUE, risk.table.y.text = FALSE) plot1 <- plot_main plot1$plot <- plot1$plot + geom_segment(x = 365.25, xend = 365.25, y = -0.05, yend = 0.4092416, size = 1.5) + geom_segment(x = 365.25, xend = -40, y = 0.4092416, yend = 0.4092416, size = 1.5, arrow = arrow(length = unit(0.2, "inches"))) ``` ] --- ## .green[실습 - Kaplan-Meier method] .pull-left[ ### .black[중위 생존시간 추정] ```r method_km ``` ``` ## Call: survfit(formula = Surv(time, status) ~ 1, data = lung) ## ## n events median 0.95LCL 0.95UCL ## [1,] 228 165 310 285 363 ``` - 말기 암환자 중위 생존시간: 310일 ] .pull-right[ ### .black[중위 생존시간 추정 - 시각화] <img src="2_survival_files/figure-html/unnamed-chunk-14-1.png" width="504" style="display: block; margin: auto;" /> ] --- ## .green[실습 - Kaplan-Meier method] ### .black[중위생존시간 - 시각화 코드] ```r plot2 <- plot_main plot2$plot <- plot2$plot + geom_segment(x = -45, xend = 310, y = 0.5, yend = 0.5, size = 1.5) + geom_segment(x = 310, xend = 310, y = 0.5, yend = -0.03, size = 1.5, arrow = arrow(length = unit(0.2, "inches"))) plot2 ``` --- class: inverse, middle, center # 생존 함수 비교 ## Log-rank test --- ## .green[Log-rank test] .pull-left[ - 그룹 간 생존 곡선을 비교하는 대표적인 비모수 검정 방법 - 두 그룹 간 생존 함수에 차이가 있는지 살펴보고자 함 - 가설 - `\(H_0: S_1(t) = S_2(t)\)` - `\(H_1: S_1(t) \neq S_2(t)\)` ] .pull-right[ <div class="figure"> <img src="./fig/log-rank.png" alt="fig 4. Log-rank test" width="624" height="500" /> <p class="caption">fig 4. Log-rank test</p> </div> ] --- ## .green[실습 - Log-rank test] - 앞서 이용한 패키지 {survival}, {survminer} 그대로 이용 - KM method를 통해 추정한 말기 암환자들의 생존 곡선을 성별 그룹에따라 비교 ```r method_km2 <- survfit(Surv(time, status) ~ sex, data = lung) method_km2 ``` ``` ## Call: survfit(formula = Surv(time, status) ~ sex, data = lung) ## ## n events median 0.95LCL 0.95UCL ## sex=1 138 112 270 212 310 ## sex=2 90 53 426 348 550 ``` - 중위생존시간은 여성이 훨씬 김 --- ## .green[실습 - Log-rank test] ### .black[시각화] .pull-left[ ```r ggsurvplot(method_km2, pval = TRUE, conf.int = TRUE, risk.table = TRUE, # Add risk table risk.table.col = "strata", # Change risk table color by groups linetype = "strata", # Change line type by groups surv.median.line = "hv", # Specify median survival ggtheme = theme_bw(), # Change ggplot2 theme palette = c("#E7B800", "#2E9FDF")) ``` - 두 점선은 각 그룹의 중위생존시간을 의미 - `\(p < 0.01\)` ] .pull-right[ <img src="2_survival_files/figure-html/unnamed-chunk-19-1.png" width="504" style="display: block; margin: auto;" /> ] <!-- https://ticket.interpark.com/Ticket/Goods/TPBridge.asp?GoodsCode=22006468 --> --- class: inverse, middle, center # 위험인자 분석 ## Cox proportional hazard model ## Time dependent Cox proportional hazard model --- ## .green[위험함수] h(t) `\begin{align} h(t) = \lim_{dt\to0} P(t < T \leq t+dt) \end{align}` - 환자가 `\(t\)` 시점까지 생존했다가 `\(t\)` 시점 바로 직후에 사망하게 되는 순간 위험률 ## .green[Cox ph model] - 생존 자료에 대한 대표적인 준모수적 회귀 분석 방법론 - Cox regression 이라고도 부름 - 비례 위험 가정 (proportinal hazard assumption) - 각 설명변수의 효과는 시간에 관계없이 일정 `\begin{align} h(t) = h_0(t){\rm{exp}}(\beta_1x_1 + \beta_2x_2 + \cdots + \beta_kx_k) \end{align}` `\begin{align} h_0(t) = \rm{exp}(\beta_0): \rm{baseline \ hazard} \end{align}` <!-- 생존시간을 T라하고 이에 영향을 주는 `\(k\)`개 공변량이 있을 때 Cox 비례위험모형을 식 (5)와 같이 표현됨 --> <!-- 식 (6)은 공변량의 값들이 모두 0일때의 위험함수 값인 기저위험함수(baseline hazard function)에 해당 --> <!-- 그래서 이게 왜 준모수 모형이냐? 식 (5)에 자연로그를 취하고 logh0(t)를 알파(t)로 놓앗다고 생각했을 때, 알파(t)를 알파제로로 놓으면 지수분포 모형, 알파(t)가 알파제로 + 알파*t이면 곰페르츠 분포, 알파(t)가 알파제로 + 알파*로그t이면 와이블분포가 됨. 여기서 물론 분포의 서택이 필요하지 않고 알파(t)의 형태가 어떤 형태이든 무관함. 이에 따라 Cox ph 모형은 모수, 비모수도 아닌 준모수적 추정 방법으로 불림 --> <!-- 준모수적 추정방법으로 완전한 확률모형이 아니기 때문에 모수(회귀계수)들의 추정에 있어 편우도(partial likehood)를 이용해야함. 편우도 함수는 사실상 자료의 수치가 아닌 순위(rank)에만 이존하는 비모수적방법으로 , 기저분포에 의존하지 않는다는 장점을 가지고 있으면서도 완전한 우도함수에 비해 효율성은 별로 떨어지지 않는 것으로 알려져 있음 --> --- ## .green[Cox ph model] <!-- 연속형의 경우 분자 x+1, 분모 x --> - 위험비 HR (Hazard ratio) - 다른 설명변수들의 값이 고정되어 있을 때, `\(x_k\)`의 효과 - 시간에 대해 상수이므로 `\(x_k\)`에 대해서만 비례하여(proportional) 증가 - 즉, `\(x_k\)`가 한 단위 증가 시 사건 발생 위험은 `\({\rm{exp}}(\beta_k)\)`배가 됨 `\begin{align} {\rm{HR}}_k = \frac{h(t|X_k = 1)}{h(t|X_k = 0)} = \frac{h_0(t){\rm{exp}}(\beta_1x_1 + \cdots + \beta_kx_k)}{h_0(t){\rm{exp}}(\beta_1x_1 + \cdots \beta_{k-1}x_{k-1})} = {\rm{exp}}(\beta_k) \end{align}` - Cox ph model은 생존 자료에서 다른 변수들의 효과를 보정한, treatment 효과를 볼 수 있는 가장 대표적인 통계 모형에 해당 - Crude HR: Univariable Cox ph model에서 구한 HR - Adjusted HR: Multivariable Cox ph model에서 구한 HR --- ## .green[비례위험가정 검토] .left-column[ ] .right-column[ - Schoenfeld test - Schoenfeld residuals와 시간 `\(t\)` 간에 어떤 패턴이 존재하는가? - 시간을 설명변수, 잔차를 종속변수로 두고 단순회귀분석을를 수행하는 개념 - 설명변수별로 검토 - 귀무가설( `\(H_0\)` ): 잔차와 시간 간 패턴이 존재하지 않음 - 즉, 귀무가설을 기각할 수 없다면 비례위험가정을 만족함 <!-- 이 검정방법외에 로그-로그 생존함수를 통해 그래프로 검토하는 방법도 있음 --> ] --- ## .green[실습 - Cox ph model] ### .black[모형 적합] .scroll-box-20[ - {survival} 패키지의 `coxph()`를 이용해 univariable cox ph model 적합 - 폐암 환자들의 생존 시간이 성별에 따라 영향을 미치는가? ```r mod_cox <- coxph(Surv(time, status) ~ sex, data = lung) mod_cox ``` ``` ## Call: ## coxph(formula = Surv(time, status) ~ sex, data = lung) ## ## coef exp(coef) se(coef) z p ## sex -0.5310 0.5880 0.1672 -3.176 0.00149 ## ## Likelihood ratio test=10.63 on 1 df, p=0.001111 ## n= 228, number of events= 165 ``` - 적합 결과 tidy하게 출력하기: ```r install.packages("broom") broom::tidy(mod_cox, exp = TRUE) ``` |term | estimate| std.error| statistic| p.value| |:----|---------:|---------:|---------:|---------:| |sex | 0.5880028| 0.1671786| -3.176385| 0.0014912| - {gtsummary}를 이용하여 또 다른 형식으로도 출력이 가능함: ```r install.packages("gtsummary") gtsummary::tbl_regression(mod_cox, exp = TRUE) ```
Characteristic
HR
1
95% CI
1
p-value
sex
0.59
0.42, 0.82
0.001
1
HR = Hazard Ratio, CI = Confidence Interval
] --- ## .green[실습 - Cox ph model] ### .black[비례위험가정 검토] .pull-left[ - `cox.zph()` 함수를 이용하면 비례위험가정 검토를 손쉽게 수행할 수 있음 - `sex`의 경우 비례위험가정을 만족한다고 볼 수 있음 ```r mod_cox |> cox.zph() |> ggcoxzph() ``` ] .pull-right[ <img src="2_survival_files/figure-html/unnamed-chunk-26-1.png" width="504" style="display: block; margin: auto;" /> ] --- ## .green[실습 - Cox ph model] ### .black[누적 위험 시각화] .pull-left[ ```r mod_hr <- survfit(Surv(time, status) ~ sex, data = lung) ggsurvplot(data = lung, fit = mod_hr, xlab = "Months", xscale = 30.4, break.x.by = 182.4, fun = "cumhaz", legend.title = "", legend.labs = c("Male", "Female"), legend = "bottom", risk.table = TRUE, risk.table.y.text = FALSE) ``` ] .pull-right[ <img src="2_survival_files/figure-html/unnamed-chunk-28-1.png" width="504" style="display: block; margin: auto;" /> ] --- ## .green[Time dependent Cox ph model] - Cox ph model의 비례위험가정 위배시 고려해야하는 모형 - 편의상 고정된 공변량, 시간에 따라 변하는 공변량을 각각 1개씩 갖는다고 가정 `\begin{align} h(t) = h_0(t){\rm{exp}}(\beta_1x_1 + \beta_2x_2(t)) \end{align}` - 시간 `\(t\)`에서 위험이 `\(x_1\)`의 값과 시간 `\(t\)`에서의 `\(x_2\)` 값에 의존한다는 것을 의미 <!-- 이런 모형을 고려해야하는 상황은 언제일까? 예를 들어 환자의 혈압이나 백혈구 수 등과 같이 중요한 요인인데 연구 도중 변화하는 공변량들이 있을 때가 그런한 경우임. 이분된 변수도 시간에 따라 변할 수 있음. 예를 들자면 심장이식 후보에 들어있는 환자들에 대해 생존분석을 수행할 때가 있음. 이 환자들은 심장이식이 결정된 순간부터 생존시간을 계산하게 되지만 결정된 수순간부터 실제로 이식되는 시점까지 간격이 생기게됨. 어떤 이는 심장이식수술을 하기 전에 사망하기도 할 것. 즉, 심장이식이 생존시간에 영향을 미치는지 알기 위해 공변량으로 '이식상태'를 생각할 수 있음. 환자가 시점 t0에서 이식수술을 받았다면 0(t<t0), 아니면 1(t>=t0)인 공변량을 고려해 분석해야할 것. 이를 바탕으로 필요한 데이터는 생존시간, 중도절단여부, 이식이 이루어질때까지 기다린 시간이 될 것임. 심장이식을 받지 못햇다면 기다린 시간은 결측이 될 것. --> - Time dependent hazard ratio `\begin{align} {\rm{HR}}(t) = \frac{h(t|X_2=x_2+1)}{h(t|X_2=x_2)} = \frac{{\rm{exp}}(\beta_1(x_1+1)+\beta_2((x_2+1)t)}{{\rm{exp}}(\beta_1x_1+\beta_2(x_2t))} \newline = {\rm{exp}}(\beta_1 + \beta_2t) \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \end{align}` - `\(x\)`가 한단위 증가시 HR은 `\({\rm{exp}}(\beta_1+\beta_2t)\)`배가 됨 - 즉, Time dependent Cox ph model의 위험비는 `\(x\)` 뿐만이 아닌 시간에도 depend --- ## .green[실습 - Time dependent Cox ph model] ### .black[데이터 불러오기] .scroll-box-20[ - {SemiCompRisks} 패키지의 137명의 골수 이식 환자 데이터 `BMT` 이용 ```r install.packages("SemiCompRisks") data(BMT, package = "SemiCompRisks") BMT2 <- BMT |> select(T1, delta1, TA, deltaA) glimpse(BMT2) ``` ``` ## Rows: 137 ## Columns: 4 ## $ T1 <int> 2081, 1602, 1496, 1462, 1433, 1377, 1330, 996, 226, 1199, 1111,… ## $ delta1 <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, … ## $ TA <int> 67, 1602, 1496, 70, 1433, 1377, 1330, 72, 226, 1199, 1111, 38, … ## $ deltaA <int> 1, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 0, 1, 0, 1, 0, … ``` - `T1`: 사망 또는 last follow-up까지 기간 (단위: 일) - `delta1`: 사망 여부(1: 사망, 0: 생존) - `TA`: aGVHD(acture graft-versus-host disease) 발생될 때까지의 기간 (단위: 일) - `deltaA`: aGVHD 발생 여부(1: 발생, 0:발생X) - 분석 목적: aGVHD 발병이 골수이식환자의 생존시간에 영향을 미치는가? ] --- ## .green[실습 - Time dependent Cox ph model] ### .black[전처리] .scroll-box-20[ - aGVHD 여부를 Time depentent covariate으로 모형 적합을 수행할 예정 - R에서 Time dependent Cox ph의 적합에는 특별한 형태의 데이터를 요구함 - 1 행 번호 생성 - 2 `tmerge()` 이용 특별한 형태의 데이터 생성 ```r install.packages("tibble") bmt2 <- BMT2 |> tibble::rowid_to_column("my_id") bmt_time <- tmerge( data1 = bmt2 |> select(my_id, T1, delta1), data2 = bmt2, id = my_id, death = event(T1, delta1), agvhd = tdc(TA) ) |> as_tibble() bmt_time ``` ``` ## # A tibble: 163 × 7 ## my_id T1 delta1 tstart tstop death agvhd ## <int> <int> <dbl> <dbl> <int> <dbl> <int> ## 1 1 2081 0 0 67 0 0 ## 2 1 2081 0 67 2081 0 1 ## 3 2 1602 0 0 1602 0 0 ## 4 3 1496 0 0 1496 0 0 ## 5 4 1462 0 0 70 0 0 ## 6 4 1462 0 70 1462 0 1 ## 7 5 1433 0 0 1433 0 0 ## 8 6 1377 0 0 1377 0 0 ## 9 7 1330 0 0 1330 0 0 ## 10 8 996 0 0 72 0 0 ## # … with 153 more rows ``` - `tmerge()`: 각 환자의 covariate 값에 따른 multiple time interval 계산을 통해 long type 데이터 생성 - `event()`: `tmerge()`에 의해 새롭게 주어진 데이터에 맞게 관심 사건의 발생 여부 생성 - `tdc()`: `tmerge()`에 의해 새롭게 주어진 데이터에 맞게 time dependent covariate 지표 생성 ] --- ## .green[실습 - Time dependent Cox ph model] ### .black[전처리] .pull-left[ - 원 데이터 ```r head(bmt2, 5) ``` ``` ## my_id T1 delta1 TA deltaA ## 1 1 2081 0 67 1 ## 2 2 1602 0 1602 0 ## 3 3 1496 0 1496 0 ## 4 4 1462 0 70 1 ## 5 5 1433 0 1433 0 ``` ] .pull-right[ - `tmerge()`에 의해 변환된 데이터 ```r head(bmt_time, 5) ``` ``` ## # A tibble: 5 × 7 ## my_id T1 delta1 tstart tstop death agvhd ## <int> <int> <dbl> <dbl> <int> <dbl> <int> ## 1 1 2081 0 0 67 0 0 ## 2 1 2081 0 67 2081 0 1 ## 3 2 1602 0 0 1602 0 0 ## 4 3 1496 0 0 1496 0 0 ## 5 4 1462 0 0 70 0 0 ``` ] --- ## .green[실습 - Time dependent Cox ph model] ### .black[모형 적합] - Cox ph model과 동일한 함수 `coxph()`로 적합 가능 ```r coxph( Surv(time = tstart, time2 = tstop, event = death) ~ agvhd, data = bmt_time ) %>% gtsummary::tbl_regression(exp = TRUE) ```
Characteristic
HR
1
95% CI
1
p-value
agvhd
1.40
0.81, 2.43
0.2
1
HR = Hazard Ratio, CI = Confidence Interval
--- ## .green[책 추천] - 🔗[생명과학연구를 위한 통계적방법](http://www.kyobobook.co.kr/product/detailViewKor.laf?barcode=9788973385348) - 🔗[The Epidemiologist R Handbook](https://epirhandbook.com/en/index.html) --- class: inverse, middle, center # References --- ## .green[References] .scroll-output[ [1] 이재원, 박미라, and 유한나. 생명과학연구를 위한 통계적방법. 자유아카데미, 2006. http://www.yes24.com/Product/Goods/3444417. [2] Batra, Neale. The Epidemiologist R Handbook, 2021. https://epirhandbook.com/en/index.html. [3] Gruman, Jim. “Survival Analysis.” Personal blog. myTidyTuesday, October 12, 2021. https://opus1993.github.io/myTidyTuesday/Survival.html. [4] Harrison, Ewen, and Riinu Pius. R for Health Data Science. CRC Press, 2021. https://argoshare.is.ed.ac.uk/healthyr_book/. [5] Mohammed, Shariq. “Introduction to Survival Analysis Using R.” Workshop on Computational Biostatistics and Survival Analysis, December 23, 2019. https://shariq-mohammed.github.io/files/cbsa2019/1-intro-to-survival.html#1_setup. [6] Zabor, Emily. “Survival Analysis in R.” Personal blog. emilyzabor, March 2019. https://www.emilyzabor.com/tutorials/survival_analysis_in_r_tutorial.html#Part_1:_Introduction_to_Survival_Analysis. ] --- class: inverse # Thanks! .pull-right[.pull-down[ <a href="mailto:favorite@kakao.com"> .white[<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg"> <path d="M440 6.5L24 246.4c-34.4 19.9-31.1 70.8 5.7 85.9L144 379.6V464c0 46.4 59.2 65.5 86.6 28.6l43.8-59.1 111.9 46.2c5.9 2.4 12.1 3.6 18.3 3.6 8.2 0 16.3-2.1 23.6-6.2 12.8-7.2 21.6-20 23.9-34.5l59.4-387.2c6.1-40.1-36.9-68.8-71.5-48.9zM192 464v-64.6l36.6 15.1L192 464zm212.6-28.7l-153.8-63.5L391 169.5c10.7-15.5-9.5-33.5-23.7-21.2L155.8 332.6 48 288 464 48l-59.4 387.3z"></path></svg> favorite@kakao.com] </a> <a href="https://github.com/be-favorite"> .white[<svg viewBox="0 0 496 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg"> <path d="M165.9 397.4c0 2-2.3 3.6-5.2 3.6-3.3.3-5.6-1.3-5.6-3.6 0-2 2.3-3.6 5.2-3.6 3-.3 5.6 1.3 5.6 3.6zm-31.1-4.5c-.7 2 1.3 4.3 4.3 4.9 2.6 1 5.6 0 6.2-2s-1.3-4.3-4.3-5.2c-2.6-.7-5.5.3-6.2 2.3zm44.2-1.7c-2.9.7-4.9 2.6-4.6 4.9.3 2 2.9 3.3 5.9 2.6 2.9-.7 4.9-2.6 4.6-4.6-.3-1.9-3-3.2-5.9-2.9zM244.8 8C106.1 8 0 113.3 0 252c0 110.9 69.8 205.8 169.5 239.2 12.8 2.3 17.3-5.6 17.3-12.1 0-6.2-.3-40.4-.3-61.4 0 0-70 15-84.7-29.8 0 0-11.4-29.1-27.8-36.6 0 0-22.9-15.7 1.6-15.4 0 0 24.9 2 38.6 25.8 21.9 38.6 58.6 27.5 72.9 20.9 2.3-16 8.8-27.1 16-33.7-55.9-6.2-112.3-14.3-112.3-110.5 0-27.5 7.6-41.3 23.6-58.9-2.6-6.5-11.1-33.3 2.6-67.9 20.9-6.5 69 27 69 27 20-5.6 41.5-8.5 62.8-8.5s42.8 2.9 62.8 8.5c0 0 48.1-33.6 69-27 13.7 34.7 5.2 61.4 2.6 67.9 16 17.7 25.8 31.5 25.8 58.9 0 96.5-58.9 104.2-114.8 110.5 9.2 7.9 17 22.9 17 46.4 0 33.7-.3 75.4-.3 83.6 0 6.5 4.6 14.4 17.3 12.1C428.2 457.8 496 362.9 496 252 496 113.3 383.5 8 244.8 8zM97.2 352.9c-1.3 1-1 3.3.7 5.2 1.6 1.6 3.9 2.3 5.2 1 1.3-1 1-3.3-.7-5.2-1.6-1.6-3.9-2.3-5.2-1zm-10.8-8.1c-.7 1.3.3 2.9 2.3 3.9 1.6 1 3.6.7 4.3-.7.7-1.3-.3-2.9-2.3-3.9-2-.6-3.6-.3-4.3.7zm32.4 35.6c-1.6 1.3-1 4.3 1.3 6.2 2.3 2.3 5.2 2.6 6.5 1 1.3-1.3.7-4.3-1.3-6.2-2.2-2.3-5.2-2.6-6.5-1zm-11.4-14.7c-1.6 1-1.6 3.6 0 5.9 1.6 2.3 4.3 3.3 5.6 2.3 1.6-1.3 1.6-3.9 0-6.2-1.4-2.3-4-3.3-5.6-2z"></path></svg> @be-favorite] </a> <a href="https://twitter.com/TaemoBang"> .white[<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg"> <path d="M459.37 151.716c.325 4.548.325 9.097.325 13.645 0 138.72-105.583 298.558-298.558 298.558-59.452 0-114.68-17.219-161.137-47.106 8.447.974 16.568 1.299 25.34 1.299 49.055 0 94.213-16.568 130.274-44.832-46.132-.975-84.792-31.188-98.112-72.772 6.498.974 12.995 1.624 19.818 1.624 9.421 0 18.843-1.3 27.614-3.573-48.081-9.747-84.143-51.98-84.143-102.985v-1.299c13.969 7.797 30.214 12.67 47.431 13.319-28.264-18.843-46.781-51.005-46.781-87.391 0-19.492 5.197-37.36 14.294-52.954 51.655 63.675 129.3 105.258 216.365 109.807-1.624-7.797-2.599-15.918-2.599-24.04 0-57.828 46.782-104.934 104.934-104.934 30.213 0 57.502 12.67 76.67 33.137 23.715-4.548 46.456-13.32 66.599-25.34-7.798 24.366-24.366 44.833-46.132 57.827 21.117-2.273 41.584-8.122 60.426-16.243-14.292 20.791-32.161 39.308-52.628 54.253z"></path></svg> @TaemoBang] </a> <a href="https://github.com/be-favorite/Presentation_archive"> .white[<svg viewBox="0 0 512 512" style="height:1em;position:relative;display:inline-block;top:.1em;" xmlns="http://www.w3.org/2000/svg"> <path d="M326.612 185.391c59.747 59.809 58.927 155.698.36 214.59-.11.12-.24.25-.36.37l-67.2 67.2c-59.27 59.27-155.699 59.262-214.96 0-59.27-59.26-59.27-155.7 0-214.96l37.106-37.106c9.84-9.84 26.786-3.3 27.294 10.606.648 17.722 3.826 35.527 9.69 52.721 1.986 5.822.567 12.262-3.783 16.612l-13.087 13.087c-28.026 28.026-28.905 73.66-1.155 101.96 28.024 28.579 74.086 28.749 102.325.51l67.2-67.19c28.191-28.191 28.073-73.757 0-101.83-3.701-3.694-7.429-6.564-10.341-8.569a16.037 16.037 0 0 1-6.947-12.606c-.396-10.567 3.348-21.456 11.698-29.806l21.054-21.055c5.521-5.521 14.182-6.199 20.584-1.731a152.482 152.482 0 0 1 20.522 17.197zM467.547 44.449c-59.261-59.262-155.69-59.27-214.96 0l-67.2 67.2c-.12.12-.25.25-.36.37-58.566 58.892-59.387 154.781.36 214.59a152.454 152.454 0 0 0 20.521 17.196c6.402 4.468 15.064 3.789 20.584-1.731l21.054-21.055c8.35-8.35 12.094-19.239 11.698-29.806a16.037 16.037 0 0 0-6.947-12.606c-2.912-2.005-6.64-4.875-10.341-8.569-28.073-28.073-28.191-73.639 0-101.83l67.2-67.19c28.239-28.239 74.3-28.069 102.325.51 27.75 28.3 26.872 73.934-1.155 101.96l-13.087 13.087c-4.35 4.35-5.769 10.79-3.783 16.612 5.864 17.194 9.042 34.999 9.69 52.721.509 13.906 17.454 20.446 27.294 10.606l37.106-37.106c59.271-59.259 59.271-155.699.001-214.959z"></path></svg> Presentation archive] </a> <br><br><br> ]]