時間相依共變數之存活分析實作 慈濟醫學中心 / 研究部王仁宏 /Paul Wang paulwang@tzuchi.com.tw 分機 17651
課程目標 希望同仁在上完這堂課能學到 : 瞭解存活分析的基本概念及使用時機 評估模式適切性 & 解讀分析結果 如何處理 時間相依共變數 (time-varying) 學習利用統計軟體 (Stata) 完成常見的分析 2
Outline 存活分析資料的特質 完整資料 vs 設限資料 Kaplan-Meier Curve & Log-rank Test Cox Model 如何處理 時間相依共變數 Cox Model:Diagnostics 3
基本觀念 事件發生時間 (time-to-event) 的資料 醫學研究 存活分析 (Survival Analysis) 工程研究 可靠度分析 (Reliability Analysis) 癌症存活時間 手術後存活時間 疾病復發的時間 社會學研究 事件歷史分析 (Event History Analysis) 產品壽命時間 系統維修時間 經濟衰退結束時間 特定商品購買時間
PART I 存活分析資料的特質 5
存活分析資料的特質 存活分析是醫學統計中常見的一種分析方法, 分析研究對象追蹤直到事件 (ex: 死亡 ) 發生的時間 Ex: 癌症患者手術後的存活時間 資料特性 涵蓋兩種類型的資料 完整資料 (Complete Data) 個案發生了研究感興趣的事件 (ex: 死亡 ), 故有準確的時間 設限資料 (Censored Data) 不完整觀察值!! 個案直到可以追蹤到的時間點都未發生研究感興趣的事件 這樣的個案貢獻了至少活超過一段時間的資訊!! 6
存活分析資料蒐集示意圖 (1/2) 2000 2011 2012 1 2 3 4 5 A 療法 Cancer 手術 死亡 O O X X 存活 : X 1 2 3 4 5 : B 療法 O O O O X X X X 7
存活分析資料蒐集示意圖 (2/2) O O A 療法 O O 完整資料 O X 設限資料 B 療法 O O O X 8
常見名詞 (1/2) 失敗機率密度函數 f(t) 存活函數 /Survival function S(t) 不同時間點存活比例 失敗函數 /Failure function F(t) 累積死亡比例 危險函數 /Hazard function h(t) ) ( 1 ) ( ) ( ) ( t F dx x f t P T t S t ) ( 1 ) ( ) ( ) ( 0 t S dx x f t P T t F t dt t df t t t T t P t f t ) ( ) ( lim ) ( 0 ) ( ) ( ) ( lim ) ( 0 t S t f t t T t t T t P t h t 9
常見名詞 (2/2) f(t) S(t) Survival rate S(t) h(t) F(t) 10
範例資料 1:Example1.xlsx Excel 資料格式 比較不同手術療法存活率 ID Sex Group FU(Day) FU(Year) Death 1 1 0 25 0.07 1 2 0 0 890 2.44 1 3 0 0 560 1.53 1 4 1 0 302 0.83 1 5 1 0 48 0.13 1 6 0 1 1304 3.57 0 7 0 1 1701 4.66 1 8 0 1 300 0.82 1 9 1 1 400 1.10 1 10 0 1 1826 5.00 0 變項 0 1 Sex 女生 男生 Group A 療法 B 療法 Death 設限 完整 11
Excel to Stata Stata File Import Excel spreadsheet 1 選擇檔案!! 4 2 3 5 將第一列視為變數名稱!! <Stata 指令 > import excel "D:\Stata\Example1.xlsx", sheet("rawdata") firstrow 6 12
查看 Stata 資料檔 1 2 3 查看資料!! 13
Survival Setting(1/3) Statistics Survival analysis Setup and utilities Declare data to be survival-time data 1 將資料轉為存活分析所需資料格式!! 2 3 4 14
Survival Setting(2/3) 設定 Time variable 追蹤時間 Failure variable Failure 變項 Failure value Failure 數值 1 2 3 <Stata 指令 > stset FUYear, failure(death==1) 4 15
Survival Setting(3/3) Stata 資料檔自動產生下列欄位 : _st 存活分析是否將此筆資料納入 (0:No/1:Yes) _d 是否發生感興趣的事件 (0:No/1:Yes) _t 追蹤截止時間 _t0 追蹤起始時間 16
PART III KAPLAN-MEIER CURVE & LOG-RANK TEST 17
存活曲線 -- Kaplan-Meier 法 Kaplan-Meier 法是常用來估計存活曲線的方法, 此方法用每一個事件發生時間點及設限點來設定區間, 適用於一般臨床相關研究 估計追蹤過程不同時間點的存活率!! 18
The History of the Kaplan-Meier Estimator In a paper published in the Journal of the American Statistical Association in June 1958, Edward Kaplan and Paul Meier put forth a new, efficient method for estimating patient survival rates, taking into account the fact that some patients may have died during a research trial while others will survive beyond the end of the trial. The method, called the Kaplan-Meier estimator (also known as the product limit estimator), is based on a mathematical formula using information from those who have died and those who have survived to estimate the proportion of patients alive at any point during the trial. The estimator is plotted over time. The resulting curve is called the Kaplan-Meier curve, which is a series of horizontal steps of declining magnitude that, when a large enough sample is taken, approaches the true survival function for that population. 19
Selected Quotes in Honor of Meier s Work at His Death in 2011 The KM curve has become the standard tool used by medical researchers for determining the duration of survival in thousands of studies, ranging from cancer to AIDS to cardiovascular disease to diabetes, to name just a few -- New York Times (2011) The KM estimate was a very, very important advance. It seems so elementary now -- Washington Post (2011) Paul Meier s work and the KM analysis have been responsible for saving millions of lives --- The Significance Magazine (2011) The KM estimator is used ubiquitously in medical studies to estimate and depict the fraction of patients living for a certain amount of time after treatment. It has since been applied to data from clinical trials of therapies for every disease from cancer to cardiology to concussion -- Science 20 Life (2011)
Kaplan-Meier Survival Curve(1/4) Statistics Survival analysis Graphs Kaplan-Meier survivor function 1 3 4 2 21
Kaplan-Meier Survival Curve(2/4) Main 選擇圖形類型 設定 Grouping variables!! 其他細項設定 設定 Y axis X axis Titles Legend!! 1 2 <Stata 指令 > sts graph, by(group) ytitle(cumulative Survival Rate) xtitle(follow-up Time(years)) title(kaplan-meier Survival Curve) legend(order(1 "A 療法 " 2 "B 療法 "))
Y axis Kaplan-Meier Survival Curve(3/4) X axis Titles Legend 23
Kaplan-Meier Survival Curve(4/4) 24
Log-rank Test Kaplan-Meier 法僅能讓我們大致瞭解不同組別 (ex: 治療方式 ) 的存活曲線分佈, 若我們需要進一步比較組別間是否相同 or 是否有顯著差異, 則需要用到 Log-rank Test 一種統計檢定方法 ( 無母數檢定 )!! 25
Log-rank Test(1/2) Statistics Survival analysis Summary statistics, tests, and tables Test equality of survivor functions 1 3 2 4 26
Log-rank Test(2/2) Main 設定比較變項 (Group) 選擇統計檢定方法!! 1 2 <Stata 指令 > sts test Group, logrank 3 27
補充 : 存活時間估算 (1/2) Statistics Survival analysis Summary statistics, tests, and tables CIs for means and percentiles of survival time 1 3 2 28 4
補充 : 存活時間估算 (2/2) Main 設定比較變項 (Group) 選擇統計量 (ex: Median)!! 1 2 <Stata 指令 > stci, by(group) median 3 29
Log-rank Test + 存活時間估計 分析結果顯示 B 療法的存活時間 ( 中位數 ) 較 A 療法長!! 但其差異尚未達統計顯著 (p-value>0.05)!! 30
Kaplan-Meier Survival Curve + Log-rank Test Graph File Start Graph Editor Add Text!! 1 4 2 3 5 31
Kaplan-Meier Survival Curve + Log-rank Test 32
PART IV COX MODEL 33
Cox Proportional Hazard Model 當我們需要同時分析多個因子 (ex: Gender Age 手術 ) 對存活時間的影響性時, 則需要利用 Cox Regression Model 來進行比較!! 類似迴歸分析的模式 比較不同因素對不同時間點發生事件風險的影響!! 34
Cox Proportional Hazard Model(1/3) Statistics Survival analysis Regression models Cox proportional hazards model 1 3 4 2 35
Cox Proportional Hazard Model(2/3) Model 設定 Cox Model 所需考慮的變項 (ex: Sex Group)!! 1 36
Cox Proportional Hazard Model(3/3) 1 Step 1. 選擇變項類型 Type of variable Step 2. 選擇變數若是類別變項, 需設定比較基準 Step 3. 加入變項清單 varlist Step 4. 確認 OK 2 3 <Stata 指令 > stcox ib(first).sex ib(first).group 或 stcox i.sex i.group 4 37
Cox Proportional Hazard Model 男性相較於女性死亡的風險比為 11.5( 發生死亡風險較高 ), 且達統計顯著!! B 療法相較於 A 療法死亡的風險比為 0.248( 發生死亡風險較低 ), 但未達統計顯 38 著 (p-value=0.128 > 0.05)!!
Cox Model 的存活曲線圖形 (1/2) Statistics Survival analysis Regression models Plot survivor, hazard, cumulative hazard, or cumulative incidence function 1 3 4 2 39
Cox Model 的存活曲線圖形 (2/2) Main 選擇圖形類型 & 設定需要呈現 curve 的條件!! 1 2 [ 註 ] 針對未設定的變數, 軟體會自動以平均值代入模式!! <Stata 指令 > stcurve, survival at1( Group=0 ) at2( Group=1 ) 40
Cox Model 的存活曲線圖形 A 療法與 B 療法死亡風險無顯著差異!! (p=0.128 ) 男性相較於女性死亡的風險顯著較高!! (p=0.044) 41
PART IV COX MODEL 如何處理 時間相依共變數
範例資料 2:hip1.dta Hip-fracture data 研究新保護裝置是否能降低老人 hip fracture 資料筆數共 106 筆 48 位 60 歲女性, 隨機選取其中 28 位穿戴 資料包含 : 編號 (id) 起始時間 (time0) 終止時間 (time1) 是否發生 hip fracture(fracture) 是否穿戴保護裝置 (protect) 年齡 (age) 血液中鈣離子濃度 (calcium) 43
Survival Setting(hip1.dta hip2.dta) 設定 Time variable 起始 & 終止時間 ID variable 識別個案的欄位 Failure variable Failure 變項 <Stata 指令 > 先開啟 hip1.dta stset time1, id(id) time0(time0) failure(fracture) hip1.dta hip2.dta 增加的部分 44
Cox Proportional Hazard Model With time-varying covariate Calcium <Stata 指令 > stcox protect age calcium 穿戴新保護裝置相較於未使用者發生 hip fracture 的風險比為 0.11!! 發生 hip fracture 風險較低, 且達統計顯著!! 年齡大小 & 血液中鈣離子濃度與發生 hip fracture 的風險無顯著相關!! 45
Cox Model 的存活曲線圖形 <Stata 指令 > (1) 產生圖形 sts graph, by(protect) ytitle(cumulative Free of Hip Fracture Rate) xtitle(follow-up Time(months)) title(kaplan-meier Survival Curve) (2)Log-rank test sts test protect, logrank 穿戴新保護裝置相較於未使用者發生 hip fracture 的風險較低!! (i.e. 免於 hip fracture 機率較高, p<0.001*) 46
PART V COX MODEL DIAGNOSTICS 47
Cox Model:Diagnostics Cox Model 的基本假設是 Proportional Hazard 類似檢定迴歸模式是否合適的殘差分析 如果模式配適的好, 殘差 e i 應該服從特定分配!! 48
Cox Model:Diagnostics 方法一 殘差檢定法 Schoenfeld residuals 指令 :estat phtest 補充 :Cox-Snell residuals 整體模式評估 方法二 圖形法 適用於類別變項 指令 1:stphplot 指令 2:stcoxkm 49
EXAMPLE 1 比較不同手術療法存活率 50
方法一 殘差檢定法 (1/3) Statistics Survival analysis Regression models Test proportional-hazards assumption 1 3 4 2 51
方法一 殘差檢定法 (2/3) Main 選擇 Schoenfeld residuals(phtest) 方法!! 其他細項設定 選擇 Options!! 1 2 52
方法一 殘差檢定法 (3/3) Main 勾選 Test proportional-hazards assumption separately for each covariate!!( 針對每個變項都判斷 ) 1 2 <Stata 指令 > estat phtest, detail 3 53
殘差檢定結果 檢定結果顯示 global test p-value=0.921>0.05 未顯著違背 Proportional Hazard Assumption 54
補充 方法一 殘差檢定法 (3/3) Main 勾選 Plot smoothed, scaled Schoenfeld residuals versus time & 設定 Covariate 為 1.Sex 或 1.Group 1 可產生殘差圖!! 2 <Stata 指令 > estat phtest, plot(1.group) 3 55
補充 殘差圖 ( 針對個別變項 ) 若殘差圖形看起來散佈在 0 附近沒有特殊 pattern OK!! Sex 的殘差 <Stata 指令 > estat phtest, plot(sex) Group 的殘差 <Stata 指令 > estat phtest, plot(group) 56
補充 整體模式之殘差檢定 Cox-Snell 殘差對殘差估計的累積危險函數圖 藍色線看起來呈現 45 就代表模型配適效果還不錯!! <Stata 指令 > 藍色字體需對應修正 stcox i.sex i.group, efron predict cs, csnell stset cs, failure(death) sts generate H = na line H cs cs, sort xlab(0 1 to 4) ylab(0 1 to 4) 57
stphplot 方法二 圖形法 適用於類別變項 Statistics Survival analysis Regression models Graphically assess proportional-hazards assumption 1 3 4 2 58
stphplot 方法二 圖形法 適用於類別變項 Main 勾選 Fit separate Cox models & 設定 Independent variable 為 Sex & 設定 Adjust variables 為 Group!! 1 2 3 4 <Stata 指令 > stphplot, by(group) adjust(sex) zero 5 59
stphplot 圖形法判斷 若圖形看起來像兩條平行線 OK!! h( t x) h0 ( t)exp( x x) x S( t x) S exp( x ) 0 ( t) Y=aX+b ln[ ln{ S( t x)}] ln[ { S0( t)}] x x 60
stcoxkm 方法二 圖形法 適用於類別變項 Statistics Survival analysis Regression models Kaplan-Meier versus predicted survival 1 3 4 2 61
stcoxkm 方法二 圖形法 適用於類別變項 Main 勾選 Fit separate Cox models & 設定 Independent variable 為 Sex & 設定 Adjust variables 為 Group!! 1 2 <Stata 指令 > stcoxkm, by(sex) 3 62
stcoxkm 圖形法判斷 藉由 K-M Curve( 真實資料 ) 和 Cox Model( 假設成立 ) 的估計值的比較評估 proportional hazard 假設是否成立!! 63
EXAMPLE 2 HIP FRACTURE DATA 64
方法一 殘差檢定法 <Stata 指令 > estat phtest, detail 檢定結果顯示 global test p-value=0.958>0.05 未顯著違背 Proportional Hazard Assumption 65
補充 殘差圖 ( 針對個別變項 ) 若殘差圖形看起來散佈在 0 附近沒有特殊 pattern OK!! Age 的殘差 <Stata 指令 > estat phtest, plot(age) Protect 的殘差 <Stata 指令 > estat phtest, plot(protect) 66
補充 整體模式之殘差檢定 Cox-Snell 殘差對殘差估計的累積危險函數圖 藍色線看起來呈現 45 就代表模型配適效果還不錯!! 藍色字體需對應修正 <Stata 指令 > stcox protect age calcium, efron predict cs, csnell stset cs, failure(fracture) sts generate H = na line H cs cs, sort xlab(0 1 to 4) ylab(0 1 to 4) 67
stphplot 方法二 圖形法 適用於類別變項 若圖形看起來像兩條平行線 OK!! <Stata 指令 > stphplot, by(protect) adjust(age calcium) zero 68
stcoxkm 方法二 圖形法 適用於類別變項 藉由 K-M Curve( 真實資料 ) 和 Cox Model( 假設成立 ) 的估計值的比較評估 proportional hazard 假設是否成立!! <Stata 指令 > stcoxkm, by(protect) 69
Review 存活分析資料的特質 完整資料 vs 設限資料 Kaplan-Meier Curve & Log-rank Test 單變項個別分析 Cox Model 多變項相互校正分析 資料準備即可處理 時間相依共變數 Cox Model:Diagnostics
Reference 中國醫藥大學生物統計中心生統 e 報 http://biostatistics.cmu.edu.tw/ An Introduction to Survival Analysis Using Stata (3 rd ) 存活分析 ( 林建甫, 2008) 71
靜思語 : 知識要用心體會, 才能變成自己的智慧 感謝聆聽 Q & A
Example2 :Hip fracture data Stata 範例資料 : 研究目的 比較某種醫材是否能有效降低 hip fracture 的發生風險? use http://www.stata-press.com/data/cggm3/hip2 73
Stata 指令 Step 0.Import data from Stata website use http://www.stata-press.com/data/cggm3/hip2, clear Step 1.Cox Model stcox protect age calcium Step 2. 模式診斷 (1) 殘差法 estat phtest, detail estat phtest, plot(age) 或 estat phtest, plot(protect) (2) 圖形法 stphplot, by(protect) adjust(age) 74 stcoxkm, by(protect)
Example2 :Hip fracture data Stata 範例資料 : 研究目的 比較某種醫材是否能有效降低 hip fracture 的發生風險? use http://www.stata-press.com/data/cggm3/hip2 75