用的生存分析工作流)
生存分析跑過(guò)的人應(yīng)該都有體會(huì)數(shù)據(jù)清洗和變量編碼剛折騰完模型又要來(lái)回?fù)Q方向。先跑了一版 Cox統(tǒng)計(jì)師說(shuō)要看 KM 曲線和生存率老板問(wèn)能不能用 RMST 做敏感性分析隔壁組還想試一把貝葉斯生存模型做不確定性量化。每個(gè)方向都得改代碼、重跑、保存結(jié)果稍不注意就覆蓋了上一版的表格最后連哪份結(jié)果是哪個(gè)模型跑出來(lái)的都分不清。如果有人物化一點(diǎn)說(shuō)這就像代碼開發(fā)里的“多個(gè)功能分支同時(shí)開發(fā)”那分析場(chǎng)景其實(shí)也需要一套分支管理機(jī)制。這就是本文要展開的主題基于 IntelligenR 的分支分析讓 Cox、RMST、KM、BRM 四個(gè)分析方向在各自分支里獨(dú)立運(yùn)行、互不打架。本文適合三類讀者一是做臨床統(tǒng)計(jì)或生物統(tǒng)計(jì)的分析師經(jīng)常需要交叉驗(yàn)證多種生存分析模型二是用 R 做數(shù)據(jù)分析、希望把分析流程組織得更工程化的小伙伴三是剛接觸生存分析、想系統(tǒng)理解 KM、Cox、RMST、貝葉斯生存模型差異的初學(xué)者。讀完本文你不僅能理解分支分析的設(shè)計(jì)思路還能拿到一套可復(fù)用的 R 代碼模板直接用在自己的數(shù)據(jù)集上。1. 什么是 IntelligenR 分支分析為什么需要它1.1 傳統(tǒng)分析方式的痛點(diǎn)傳統(tǒng)做法里一個(gè)“完整”的生存分析流程往往長(zhǎng)這樣讀取數(shù)據(jù)、清洗、跑 KM、跑 Cox、畫圖、導(dǎo)出表格。改一個(gè)變量篩選條件就要從中間重新執(zhí)行換一個(gè)模型做敏感性分析又要復(fù)制一份腳本或者手動(dòng)備份結(jié)果。文件命名可能出現(xiàn)“final_v2”“final_v3_真的最終”這種狀況結(jié)果一多根本分不清。更麻煩的是不同分析方向之間還有依賴關(guān)系。KM 曲線可以幫你快速看數(shù)據(jù)分布Cox 模型能給出調(diào)整后的風(fēng)險(xiǎn)比RMST 在非比例風(fēng)險(xiǎn)場(chǎng)景下更穩(wěn)健貝葉斯模型則能給出參數(shù)的完整后驗(yàn)分布。這些模型彼此不是替代關(guān)系而是互補(bǔ)關(guān)系需要在同一份基礎(chǔ)數(shù)據(jù)上并行推進(jìn)再匯總判斷?？渴止し绞焦芾磉@種多方向分析效率低還容易出錯(cuò)。1.2 分支分析的核心思想關(guān)注過(guò) Git 工作流的朋友都知道“分支”的好處基于同一個(gè)主分支可以切出多個(gè) feature 分支各自開發(fā)互不阻塞最后再合并回主干。IntelligenR 里的分支分析思路與此高度相似只不過(guò)管理的對(duì)象從代碼變成了“分析任務(wù)”和“模型結(jié)果”。一套完整的 IntelligenR 分支分析通常包含三個(gè)層級(jí)基礎(chǔ)數(shù)據(jù)層所有分支共享同一份經(jīng)過(guò)校驗(yàn)的原始數(shù)據(jù)和分析基線。分析分支層每個(gè)分支對(duì)應(yīng)一個(gè)獨(dú)立分析方向比如 KM 分支、Cox 分支、RMST 分支、BRM 分支。結(jié)果匯總層分支完成后將關(guān)鍵結(jié)果匯總到統(tǒng)一的報(bào)告或?qū)Ρ缺碇泄Q策時(shí)取用。這樣做的好處是即時(shí)可見(jiàn)的每個(gè)分支獨(dú)立執(zhí)行模型 A 的改動(dòng)不會(huì)影響模型 B。同一數(shù)據(jù)可重復(fù)用于不同分析方向不需要復(fù)制多份數(shù)據(jù)集。分支結(jié)果可追溯哪個(gè)分支用了什么參數(shù)一目了然。遇到結(jié)果沖突時(shí)可以直接定位到具體分支排查不用從頭檢查整個(gè)流程。1.3 適用場(chǎng)景分支分析不是銀彈它特別適合以下幾類場(chǎng)景生存分析中的多模型對(duì)比KM、Cox、RMST、BRM 并行跑最后匯總。敏感性分析更換時(shí)間截點(diǎn)、變量集合、缺失值處理策略時(shí)開不同分支。團(tuán)隊(duì)協(xié)作統(tǒng)計(jì)分析員和數(shù)據(jù)科學(xué)家在同一份數(shù)據(jù)基礎(chǔ)上各自負(fù)責(zé)不同分支。分析流程復(fù)用同一套模板換一批數(shù)據(jù)自動(dòng)生成全套分支結(jié)果。如果你的分析任務(wù)只需要“一次性跑一個(gè)模型、得到一張表”那分支分析確實(shí)有點(diǎn)重但只要你的分析結(jié)果需要反復(fù)迭代、多方評(píng)審、結(jié)論要經(jīng)得起追問(wèn)分支分析的價(jià)值就非常明顯。2. 核心模型定位KM / Cox / RMST / BRM 分別解決什么問(wèn)題在進(jìn)入代碼之前先把四個(gè)模型的關(guān)系梳理清楚。這幾位出場(chǎng)頻率高但很多人對(duì)它們“什么時(shí)候用哪個(gè)”還是模糊的。2.1 KM 生存曲線先描述生存狀況Kaplan-Meier 估計(jì)是一種非參數(shù)方法用來(lái)估計(jì)生存函數(shù)。它不需要假設(shè)生存時(shí)間服從某個(gè)特定分布只會(huì)利用已知的事件和刪失數(shù)據(jù)計(jì)算出每個(gè)時(shí)間點(diǎn)的生存概率。輸出是一張階梯狀的生存曲線以及中位生存時(shí)間、某時(shí)間點(diǎn)生存率等描述性指標(biāo)。KM 主要用于探索和展示不同分組之間的生存曲線長(zhǎng)什么樣兩組曲線是否明顯分開它配合 log-rank 檢驗(yàn)可以給出組間差異的 p 值但它無(wú)法同時(shí)調(diào)整多個(gè)協(xié)變量也無(wú)法直接回答“某個(gè)因素每增加一個(gè)單位風(fēng)險(xiǎn)會(huì)變化多少”這種問(wèn)題。2.2 Cox 比例風(fēng)險(xiǎn)回歸量化風(fēng)險(xiǎn)因素Cox 回歸是生存分析里最常用的半?yún)?shù)模型。它不關(guān)心基線風(fēng)險(xiǎn)函數(shù)的具體形狀直接對(duì)協(xié)變量與風(fēng)險(xiǎn)之間的關(guān)系做回歸建模核心輸出是風(fēng)險(xiǎn)比 HRHazard Ratio。比如 HR 0.65意思是治療組的風(fēng)險(xiǎn)是對(duì)照組的 0.65 倍也就是風(fēng)險(xiǎn)降低了 35%。Cox 模型的優(yōu)勢(shì)在于可以同時(shí)納入多個(gè)協(xié)變量調(diào)整混雜因素。但它有一個(gè)重要前提就是比例風(fēng)險(xiǎn)假設(shè)PH 假設(shè)即各組的風(fēng)險(xiǎn)比在時(shí)間上保持恒定。如果 KM 曲線明顯交叉說(shuō)明不同時(shí)間區(qū)間內(nèi)風(fēng)險(xiǎn)優(yōu)勢(shì)在變化此時(shí) Cox 模型給出的單一 HR 可能會(huì)誤導(dǎo)結(jié)論。2.3 RMST非比例風(fēng)險(xiǎn)下的穩(wěn)健比較RMSTRestricted Mean Survival Time限制平均生存時(shí)間衡量的是在某個(gè)預(yù)先指定的時(shí)間窗口比如 5 年內(nèi)研究對(duì)象平均存活的時(shí)間。它不依賴比例風(fēng)險(xiǎn)假設(shè)即使生存曲線交叉也能給出一個(gè)直觀的組間差異估計(jì)比如“治療組在 5 年內(nèi)平均比對(duì)照組多活了 4.2 個(gè)月”。RMST 在臨床研究里經(jīng)常作為 Cox 模型的敏感性分析。如果 Cox 的 PH 假設(shè)不滿足或者你想給非統(tǒng)計(jì)背景的讀者講清楚“兩組到底差了多少”RMST 是更直觀的選擇。2.4 BRM 貝葉斯回歸模型把先驗(yàn)知識(shí)和不確定性帶進(jìn)分析BRM 在這里指貝葉斯回歸模型Bayesian Regression Model。在生存分析場(chǎng)景下就是對(duì)生存時(shí)間或風(fēng)險(xiǎn)函數(shù)建立一個(gè)貝葉斯回歸框架比如貝葉斯 Weibull 生存回歸、貝葉斯 Cox 模型等。貝葉斯方法的核心是“先驗(yàn)分布 數(shù)據(jù)似然 后驗(yàn)分布”。它的好處在于可以顯式納入先驗(yàn)知識(shí)比如參考?xì)v史研究認(rèn)為某藥物的 HR 在 0.6 到 0.8 之間。輸出不是單一點(diǎn)估計(jì)而是參數(shù)的后驗(yàn)分布能更自然地表達(dá)不確定性。在小樣本、事件數(shù)較少的情況下適當(dāng)使用信息先驗(yàn)可以降低模型不穩(wěn)定性。BRM 的代價(jià)是計(jì)算量大對(duì)使用者有更高的建模功底要求而且在團(tuán)隊(duì)協(xié)作里解釋成本也更高。因此它通常作為進(jìn)階驗(yàn)證手段而不是第一默認(rèn)模型。這四種模型其實(shí)是層層遞進(jìn)的KM 負(fù)責(zé)“看”Cox 負(fù)責(zé)“量”RMST 負(fù)責(zé)“穩(wěn)”BRM 負(fù)責(zé)“深”。分支分析恰好可以把它們組織在一個(gè)工作流里讓每個(gè)方向各司其職。3. 環(huán)境準(zhǔn)備與示例數(shù)據(jù)說(shuō)明3.1 軟件環(huán)境本文示例以 R 語(yǔ)言為主因?yàn)樯娣治鱿嚓P(guān)的生態(tài)最成熟。你需要準(zhǔn)備一個(gè)可用的 R 環(huán)境并安裝以下包install.packages(c(survival, survminer, survRM2, brms, dplyr))版本說(shuō)明survival是 R 自帶生態(tài)里的核心包survminer用于繪制 KM 曲線survRM2是跑 RMST 的經(jīng)典包brms用于貝葉斯生存回歸。具體版本不需要完全一致但建議保持較新的穩(wěn)定版文中代碼以常見(jiàn)版本的接口為準(zhǔn)如果你使用的版本有接口變化以官方文檔為準(zhǔn)。3.2 示例數(shù)據(jù)結(jié)構(gòu)為了演示分支分析我們先構(gòu)造一份模擬的臨床試驗(yàn)數(shù)據(jù)。假設(shè)我們想比較某種治療方案對(duì)患者生存時(shí)間的影響同時(shí)調(diào)整年齡和性別兩個(gè)協(xié)變量。數(shù)據(jù)字段如下字段含義id患者編號(hào)time生存時(shí)間單位月status事件狀態(tài)1 表示發(fā)生終點(diǎn)事件0 表示刪失group分組Treatment 治療組 / Control 對(duì)照組age年齡sex性別生成模擬數(shù)據(jù)的 R 代碼如下set.seed(2024) n - 300 group - factor(sample(c(Treatment, Control), n, replace TRUE)) age - rnorm(n, 60, 10) sex - factor(sample(c(Male, Female), n, replace TRUE)) true_hr - ifelse(group Treatment, 0.65, 1.0) scale_t - 50 / (true_hr^(1/2)) time - rweibull(n, shape 2, scale scale_t) cens - runif(n, 0, 80) status - ifelse(time cens, 1, 0) time - pmin(time, cens) data - data.frame(id 1:n, time time, status status, group group, age age, sex sex) write.csv(data, data/raw_data.csv, row.names FALSE)這段代碼的思路是先設(shè)定治療組真實(shí)風(fēng)險(xiǎn)比 HR 0.65然后基于 Weibull 分布生成生存時(shí)間再隨機(jī)生成一個(gè)刪失時(shí)間如果刪失時(shí)間早于生存時(shí)間則該樣本視為刪失。這樣做出來(lái)的模擬數(shù)據(jù)帶有真實(shí)效應(yīng)方便后面驗(yàn)證模型是否能識(shí)別出治療組的風(fēng)險(xiǎn)降低。3.3 IntelligenR 分支工作區(qū)結(jié)構(gòu)為了方便后續(xù)擴(kuò)展建議先建立一套清晰的分析目錄結(jié)構(gòu)intelligenr_branch_demo/ ├── data/ │ └── raw_data.csv ├── branches/ │ ├── km_analysis.R │ ├── cox_analysis.R │ ├── rmst_analysis.R │ └── brm_analysis.R ├── reports/ │ ├── km/ │ ├── cox/ │ ├── rmst/ │ └── brm/ ├── logs/ └── scripts/ └── run_all_branches.shdata放原始數(shù)據(jù)branches放各分支分析腳本reports按分支分別保存輸出結(jié)果logs記錄運(yùn)行日志scripts放批量執(zhí)行腳本。這樣的結(jié)構(gòu)本身就是“分支”思想在文件系統(tǒng)上的落地。4. 實(shí)戰(zhàn)在 IntelligenR 中建立四個(gè)分析分支下面進(jìn)入實(shí)操環(huán)節(jié)。我們將基于同一份模擬數(shù)據(jù)分別跑 KM、Cox、RMST、BRM 四個(gè)分支最后匯總結(jié)果。4.1 分支 AKM 曲線與 log-rank 檢驗(yàn)KM 分支是整個(gè)分析流程的地基。它的作用是描述性展示兩組生存曲線是否分開生存率隨時(shí)間怎么變化組間差異是否顯著。這是后續(xù)所有模型解釋的背景。文件路徑branches/km_analysis.Rlibrary(survival) library(survminer) data - read.csv(data/raw_data.csv) # 基于分組擬合 KM 生存曲線 fit_km - survfit(Surv(time, status) ~ group, data data) # 查看整體生存估計(jì) print(fit_km) # 輸出各時(shí)間點(diǎn)的生存率 summary(fit_km, times c(12, 24, 36, 48, 60)) # log-rank 檢驗(yàn)比較兩組生存曲線差異 logrank_test - survdiff(Surv(time, status) ~ group, data data) print(logrank_test) # 繪制 KM 曲線并疊加 p 值和風(fēng)險(xiǎn)表 pdf(reports/km/km_curve.pdf, width 8, height 6) ggsurvplot( fit_km, data data, pval TRUE, pval.method TRUE, risk.table TRUE, conf.int TRUE, xlab Time (months), ylab Overall Survival Probability ) dev.off()關(guān)鍵輸出會(huì)包括中位生存時(shí)間median survival timeTreatment 組假設(shè)為 39 個(gè)月Control 組假設(shè)為 28 個(gè)月12、24、36、48、60 個(gè)月的生存率。log-rank 檢驗(yàn) p 值如果小于 0.05說(shuō)明兩組生存分布有顯著差異。這里要提醒一句KM 曲線的 p 值只告訴我們“不同分組間生存分布有差異”并不等于“治療有效”因?yàn)樗鼪](méi)有調(diào)整年齡、性別等混雜因素。真正的因果性結(jié)論需要回到 Cox 或更嚴(yán)格的試驗(yàn)設(shè)計(jì)。4.2 分支 BCox 比例風(fēng)險(xiǎn)回歸Cox 分支回答的問(wèn)題是在同時(shí)調(diào)整年齡、性別的情況下治療組相對(duì)對(duì)照組風(fēng)險(xiǎn)降低了多少如果 HR 的可信區(qū)間不包含 1說(shuō)明效應(yīng)有統(tǒng)計(jì)學(xué)意義。文件路徑branches/cox_analysis.Rlibrary(survival) data - read.csv(data/raw_data.csv) # 擬合 Cox 回歸 fit_cox - coxph(Surv(time, status) ~ group age sex, data data) # 輸出模型結(jié)果 summary(fit_cox) # 提取 HR 和置信區(qū)間 hr_table - data.frame( variable names(coef(fit_cox)), HR exp(coef(fit_cox)), lower_95 exp(confint(fit_cox)[, 1]), upper_95 exp(confint(fit_cox)[, 2]), p_value summary(fit_cox)$coefficients[, 5] ) write.csv(hr_table, reports/cox/hr_table.csv, row.names FALSE) # 比例風(fēng)險(xiǎn)假設(shè)檢驗(yàn) ph_test - cox.zph(fit_cox) print(ph_test) pdf(reports/cox/ph_assumption.pdf, width 8, height 6) plot(ph_test) dev.off()在結(jié)果解讀上優(yōu)先看以下幾點(diǎn)groupTreatment對(duì)應(yīng)的 HR如果約為 0.65說(shuō)明治療組風(fēng)險(xiǎn)比對(duì)照組低約 35%。HR 的 95% 置信區(qū)間如果區(qū)間完全不包含 1說(shuō)明結(jié)果穩(wěn)定。cox.zph的 p 值如果某個(gè)變量的 p 值小于 0.05說(shuō)明該變量可能不滿足比例風(fēng)險(xiǎn)假設(shè)Cox 模型的結(jié)論需要謹(jǐn)慎此時(shí)應(yīng)去看 RMST 分支的結(jié)果。補(bǔ)充一點(diǎn)cox.zph全局檢驗(yàn)和單變量檢驗(yàn)都要看不要只看某一個(gè) p 值。如果 PH 假設(shè)不滿足常規(guī)做法是改用分層 Cox、時(shí)變系數(shù)模型或者直接用 RMST 作為主要分析。4.3 分支 CRMST 限制平均生存時(shí)間RMST 分支是 Cox 分支的重要補(bǔ)充尤其在生存曲線交叉或 PH 假設(shè)不成立時(shí)RMST 可以提供更穩(wěn)健的組間差異估計(jì)。文件路徑branches/rmst_analysis.Rlibrary(survRM2) data - read.csv(data/raw_data.csv) # 設(shè)定時(shí)間窗口 tau一般取臨床上有意義的時(shí)間點(diǎn) tau - 60 # 構(gòu)造治療的二元向量1 表示治療組0 表示對(duì)照組 arm - as.numeric(data$group Treatment) # RMST 分析 rmst_res - rmst2(time data$time, status data$status, arm arm, tau tau) # 輸出結(jié)果 print(rmst_res) # 手動(dòng)整理常用結(jié)果 rmst_table - data.frame( group c(Treatment, Control), rmst c(rmst_res$RMST.arm1$est[1], rmst_res$RMST.arm0$est[1]) ) write.csv(rmst_table, reports/rmst/rmst_summary.csv, row.names FALSE)rmst2的輸出核心是兩組 RMST 的差值RMST difference以及對(duì)應(yīng)的 95% 置信區(qū)間。如果差值為正且置信區(qū)間不包含 0說(shuō)明治療組在 60 個(gè)月內(nèi)的平均生存時(shí)間顯著長(zhǎng)于對(duì)照組。注意tau的設(shè)定要有臨床依據(jù)不能隨便取一個(gè)極大值。如果tau大于最大隨訪時(shí)間RMST 實(shí)際上就退化成整個(gè)隨訪期內(nèi)的平均生存時(shí)間可能會(huì)被右刪失影響不夠穩(wěn)健。通常建議把tau設(shè)為一個(gè)大多數(shù)患者都已經(jīng)出現(xiàn)終點(diǎn)事件的時(shí)點(diǎn)比如本研究中的 60 個(gè)月。4.4 分支 DBRM 貝葉斯生存回歸BRM 分支走的是貝葉斯路線。我們會(huì)用brms包擬合一個(gè)貝葉斯 Weibull 生存回歸把治療組、年齡、性別納入模型得到每個(gè)參數(shù)的后驗(yàn)分布。文件路徑branches/brm_analysis.Rlibrary(brms) data - read.csv(data/raw_data.csv) # brms 中 status1 表示事件需要轉(zhuǎn)換為刪失指示變量 data$censored - ifelse(data$status 1, 0, 1) # 貝葉斯 Weibull 生存回歸 fit_brm - brm( bf(time | cens(censored) ~ group age sex), data data, family weibull(), chains 2, cores 2, iter 2000, seed 2024 ) # 查看后驗(yàn)匯總 summary(fit_brm) # 提取關(guān)鍵參數(shù)的后驗(yàn)樣本 posterior_samples - as.data.frame(fit_brm) write.csv( summary(fit_brm)$fixed, reports/brm/fixed_effects_summary.csv, row.names TRUE ) # 繪制治療組效應(yīng)的后驗(yàn)分布 pdf(reports/brm/posterior_plot.pdf, width 8, height 6) bayesplot::mcmc_areas(fit_brm, pars b_groupTreatment) dev.off()在解釋貝葉斯結(jié)果時(shí)不要套用頻率學(xué)派的 p 值表達(dá)而是關(guān)注后驗(yàn)均值b_groupTreatment表示治療效果在 Weibull 尺度上的回歸系數(shù)負(fù)值意味著治療組的風(fēng)險(xiǎn)更低、生存時(shí)間更長(zhǎng)。90% 或 95% 可信區(qū)間如果區(qū)間不跨越 0說(shuō)明后驗(yàn)證據(jù)比較明確。樣本量較小時(shí)后驗(yàn)分布的尾部信息可以告訴你結(jié)果有多不確定。如果你對(duì)貝葉斯模型不熟悉建議先在經(jīng)典模型Cox上確認(rèn)結(jié)果方向再用 BRM 做穩(wěn)健性驗(yàn)證。貝葉斯模型的計(jì)算更慢鏈路更多適合作為“驗(yàn)證分支”而不是“起步分支”。這里需要說(shuō)明一點(diǎn)brms的刪失參數(shù)約定可能隨版本變化cens()的具體寫法請(qǐng)以你安裝版本的官方文檔為準(zhǔn)。如果brms安裝不便也可以先用rstanarm::stan_surv或者使用 JAGS 手寫 Weibull 生存模型思路完全相同。4.5 分支并行執(zhí)行與結(jié)果匯總四個(gè)分支腳本都準(zhǔn)備好后可以寫一個(gè)批量執(zhí)行腳本讓它們依次或并行運(yùn)行。并行執(zhí)行能省時(shí)間但要注意每個(gè)分支的日志和輸出目錄必須隔離避免互相覆蓋。文件路徑scripts/run_all_branches.sh#!/bin/bash set -e mkdir -p reports/km reports/cox reports/rmst reports/brm mkdir -p logs echo Run KM Branch Rscript branches/km_analysis.R logs/km.log 21 echo Run Cox Branch Rscript branches/cox_analysis.R logs/cox.log 21 echo Run RMST Branch Rscript branches/rmst_analysis.R logs/rmst.log 21 echo Run BRM Branch Rscript branches/brm_analysis.R logs/brm.log 21 echo All branches completed.如果你的機(jī)器是多核環(huán)境也可以使用parallel或 R 的future包做并行調(diào)度。但要注意貝葉斯分支本身已經(jīng)占用了多個(gè) CPU 核再和其他分支并行可能會(huì)造成資源競(jìng)爭(zhēng)反而變慢。結(jié)果匯總階段四個(gè)分支產(chǎn)出如下分支核心輸出文件位置KM生存曲線圖、中位生存時(shí)間、log-rank p 值reports/km/CoxHR 表、PH 假設(shè)檢驗(yàn)結(jié)果reports/cox/RMST兩組 RMST、RMST 差值與置信區(qū)間reports/rmst/BRM后驗(yàn)分布、可信區(qū)間、后驗(yàn)診斷圖reports/brm/5. 分支對(duì)比、合并與沖突處理多個(gè)分支跑完最關(guān)鍵的一步是“對(duì)比和匯總”。如果四個(gè)分支結(jié)論一致那結(jié)論就比較穩(wěn)固但如果分支之間出現(xiàn)了沖突怎么辦5.1 四類模型結(jié)果之間的互補(bǔ)性一個(gè)比較理想的結(jié)果應(yīng)該是這樣的KM 曲線顯示治療組生存曲線整體在對(duì)照組上方log-rank p 值 0.05。Cox 模型顯示治療組 HR 0.6595% 置信區(qū)間為 0.48 ~ 0.87PH 假設(shè)檢驗(yàn)不顯著。RMST 顯示治療組 60 個(gè)月 RMST 比對(duì)照組高 4.5 個(gè)月95% 置信區(qū)間不包含 0。BRM 后驗(yàn)分布顯示治療組效應(yīng)參數(shù)的后驗(yàn)均值遠(yuǎn)離 0可信區(qū)間不含 0。這種“四方結(jié)論一致”的情況最理想可以直接寫進(jìn)結(jié)論部分。5.2 結(jié)果沖突怎么辦如果 KM 曲線交叉、Cox 卻給出了顯著的 HR這就進(jìn)入了一個(gè)高頻沖突場(chǎng)景。此時(shí)應(yīng)該優(yōu)先懷疑 PH 假設(shè)是否滿足。KM 曲線交叉往往說(shuō)明不同時(shí)間段內(nèi)兩組的風(fēng)險(xiǎn)比在變化Cox 模型給出的“平均 HR”可能掩蓋了這種變化。此時(shí)應(yīng)該查看cox.zph的 p 值確認(rèn)是否違反 PH 假設(shè)。改看 RMST 分支用“某個(gè)時(shí)間窗內(nèi)的平均生存時(shí)間差異”做組間比較。對(duì) Cox 分支做分層分析或引入時(shí)間交互項(xiàng)而不是直接刪掉有問(wèn)題的變量。還有一種沖突是頻率學(xué)派和貝葉斯學(xué)派結(jié)論不一致Cox 給出的 p 值大于 0.05但 BRM 的后驗(yàn)可信區(qū)間不包含 0。這種情況通常和先驗(yàn)設(shè)置有關(guān)也可能是因?yàn)閮膳蓹z驗(yàn)邏輯本身不同。處理方式是把先驗(yàn)設(shè)置寫清楚并在報(bào)告中同時(shí)呈現(xiàn)兩套結(jié)果由讀者或評(píng)審專家判斷。5.3 分支合并與回退規(guī)則分支分析不是所有分支都必須合并。你可以根據(jù)分析目的制定簡(jiǎn)單的規(guī)則場(chǎng)景處理方式所有分支結(jié)論一致合并所有分支結(jié)論寫入最終報(bào)告Cox 和 RMST 沖突以 RMST 為主Cox 作為敏感性分析BRM 分支不收斂檢查先驗(yàn)和采樣參數(shù)必要時(shí)回退到經(jīng)典模型KM 曲線交叉放棄單一 HR 的解釋采用時(shí)變效應(yīng)或 RMST分支管理的核心原則是寧可保留多個(gè)分支的分歧結(jié)果也不要為了得到一個(gè)“漂亮結(jié)論”強(qiáng)行統(tǒng)一。6. 常見(jiàn)問(wèn)題與排查思路在實(shí)際跑 IntelligenR 分支分析時(shí)大家最常遇到的還是環(huán)境和模型細(xì)節(jié)問(wèn)題。下面整理一份排查清單。問(wèn)題現(xiàn)象常見(jiàn)原因解決思路中文標(biāo)簽在圖里顯示為方塊中文字體缺失在 R 里配置中文字體或把圖片內(nèi)文字改為英文cox.zph檢驗(yàn) p 值 0.05比例風(fēng)險(xiǎn)假設(shè)不成立改用分層 Cox、時(shí)變系數(shù)或直接看 RMST 分支RMST 的 tau 設(shè)太大最遠(yuǎn)隨訪時(shí)間前事件率不足調(diào)整 tau 為臨床有意義的截點(diǎn)如 48 或 60 個(gè)月貝葉斯模型采樣警告鏈長(zhǎng)不足或先驗(yàn)不合適增加 iter調(diào)整先驗(yàn)檢查 Rhat 是否接近 1分支之間結(jié)果不一致數(shù)據(jù)版本或隨機(jī)種子不一致確保所有分支使用同一份 raw_data.csv固定 set.seed事件數(shù)太少模型不穩(wěn)定樣本量不足或刪失率過(guò)高合并分組、減少協(xié)變量、考慮容錯(cuò)更強(qiáng)的貝葉斯模型并行執(zhí)行時(shí)內(nèi)存不足多個(gè) R 進(jìn)程同時(shí)吃內(nèi)存限制并行度BRM 分支單獨(dú)執(zhí)行這些問(wèn)題的共性是先定位是數(shù)據(jù)問(wèn)題、模型問(wèn)題還是環(huán)境問(wèn)題再對(duì)癥下藥不要一上來(lái)就改模型公式。7. 最佳實(shí)踐與工程化建議7.1 分支命名與分析記錄分支的名字要能體現(xiàn)分析方向、數(shù)據(jù)版本和參數(shù)變更。比如“cox_v2_adjust_age_sex”“rmst_tau60_sensitivity”。每次跑完分支在日志里記錄時(shí)間、數(shù)據(jù) hash、依賴包版本和運(yùn)行結(jié)果摘要方便后面追溯。7.2 數(shù)據(jù)不可變?cè)瓌t分支分析的前提是“基礎(chǔ)數(shù)據(jù)一致”。在實(shí)際操作中建議洗好一份raw_data.csv后就不允許再改動(dòng)。任何變量轉(zhuǎn)換、缺失值填充都放到分支腳本內(nèi)部做這樣分支之間不會(huì)互相污染。如果確實(shí)要修改基礎(chǔ)數(shù)據(jù)請(qǐng)生成一個(gè)新的數(shù)據(jù)版本而不是在原文件上改。7.3 隨機(jī)種子與可復(fù)現(xiàn)性所有涉及隨機(jī)數(shù)的環(huán)節(jié)包括模擬數(shù)據(jù)、樣本拆分、貝葉斯采樣都要顯式設(shè)置set.seed()這樣別人拿到你的代碼能復(fù)現(xiàn)出相同結(jié)果。貝葉斯模型還要額外記錄采樣器參數(shù)比如 chains、iter、warmup、adapt_delta這些都會(huì)影響結(jié)果穩(wěn)定性。7.4 輸出格式統(tǒng)一每個(gè)分支的輸出建議采用統(tǒng)一的命名和格式。比如所有表格都輸出 CSV所有圖片都輸出 PDF并且以模型名_內(nèi)容_日期命名。匯總時(shí)不僅給數(shù)字還要給出置信區(qū)間或可信區(qū)間避免只列點(diǎn)估計(jì)。7.5 合規(guī)與倫理提醒如果你的數(shù)據(jù)來(lái)自臨床研究或真實(shí)患者信息務(wù)必確認(rèn)分析方案已通過(guò)倫理審批數(shù)據(jù)使用已獲得授權(quán)并且分析過(guò)程中遵守最小權(quán)限和數(shù)據(jù)脫敏原則。分支分析只是技術(shù)手段它不能代替合規(guī)審查。任何時(shí)候都不應(yīng)該在未授權(quán)數(shù)據(jù)上開展分析也不應(yīng)該把內(nèi)部數(shù)據(jù)結(jié)果隨意對(duì)外發(fā)布。8. 總結(jié)與下一步學(xué)習(xí)路線本文圍繞 IntelligenR 分支分析展開把 KM、Cox、RMST、BRM 四種生存分析模型組織成了四條互不干擾的分析分支。前半部分幫你建立概念框架搞清四個(gè)模型各自的適用場(chǎng)景后半部分用一份模擬數(shù)據(jù)完整演示了從數(shù)據(jù)清洗、分支腳本編寫、批量執(zhí)行到結(jié)果匯總的全流程。下一步你可以從三個(gè)方向繼續(xù)深入如果對(duì)生存分析本身還不熟悉建議先吃透 KM 曲線和 Cox 模型的輸出解讀再嘗試 RMST 和貝葉斯模型。如果想把分支分析做得更自動(dòng)化可以學(xué)習(xí)targets、drake這類 R 工作流管理包把分支依賴關(guān)系寫得更規(guī)范。如果對(duì)貝葉斯生存模型感興趣可以系統(tǒng)學(xué)習(xí)brms的刪失數(shù)據(jù)建模、先驗(yàn)敏感性分析以及后驗(yàn)預(yù)測(cè)檢查。在真實(shí)項(xiàng)目中最優(yōu)先要關(guān)注的永遠(yuǎn)不是模型有多復(fù)雜而是數(shù)據(jù)質(zhì)量、模型假設(shè)和結(jié)果可復(fù)現(xiàn)性。分支分析給你的是一套“多條腿走路”的安全機(jī)制但最終能不能站住還是要看每條腿是否都踩在扎實(shí)的數(shù)據(jù)和嚴(yán)謹(jǐn)?shù)慕y(tǒng)計(jì)假設(shè)上。如果你在跑分支、匯總結(jié)果或解讀模型時(shí)遇到其他報(bào)錯(cuò)歡迎帶著具體報(bào)錯(cuò)和模型輸出繼續(xù)交流。 IntelligenR分支分析：打造高效可復(fù)用的生存分析工作流
生存分析跑過(guò)的人應(yīng)該都有體會(huì)數(shù)據(jù)清洗和變量編碼剛折騰完模型又要來(lái)回?fù)Q方向。先跑了一版 Cox統(tǒng)計(jì)師說(shuō)要看 KM 曲線和生存率老板問(wèn)能不能用 RMST 做敏感性分析隔壁組還想試一把貝葉斯生存模型做不確定性量化。每個(gè)方向都得改代碼、重跑、保存結(jié)果稍不注意就覆蓋了上一版的表格最后連哪份結(jié)果是哪個(gè)模型跑出來(lái)的都分不清。如果有人物化一點(diǎn)說(shuō)這就像代碼開發(fā)里的“多個(gè)功能分支同時(shí)開發(fā)”那分析場(chǎng)景其實(shí)也需要一套分支管理機(jī)制。這就是本文要展開的主題基于 IntelligenR 的分支分析讓 Cox、RMST、KM、BRM 四個(gè)分析方向在各自分支里獨(dú)立運(yùn)行、互不打架。本文適合三類讀者一是做臨床統(tǒng)計(jì)或生物統(tǒng)計(jì)的分析師經(jīng)常需要交叉驗(yàn)證多種生存分析模型二是用 R 做數(shù)據(jù)分析、希望把分析流程組織得更工程化的小伙伴三是剛接觸生存分析、想系統(tǒng)理解 KM、Cox、RMST、貝葉斯生存模型差異的初學(xué)者。讀完本文你不僅能理解分支分析的設(shè)計(jì)思路還能拿到一套可復(fù)用的 R 代碼模板直接用在自己的數(shù)據(jù)集上。1. 什么是 IntelligenR 分支分析為什么需要它1.1 傳統(tǒng)分析方式的痛點(diǎn)傳統(tǒng)做法里一個(gè)“完整”的生存分析流程往往長(zhǎng)這樣讀取數(shù)據(jù)、清洗、跑 KM、跑 Cox、畫圖、導(dǎo)出表格。改一個(gè)變量篩選條件就要從中間重新執(zhí)行換一個(gè)模型做敏感性分析又要復(fù)制一份腳本或者手動(dòng)備份結(jié)果。文件命名可能出現(xiàn)“final_v2”“final_v3_真的最終”這種狀況結(jié)果一多根本分不清。更麻煩的是不同分析方向之間還有依賴關(guān)系。KM 曲線可以幫你快速看數(shù)據(jù)分布Cox 模型能給出調(diào)整后的風(fēng)險(xiǎn)比RMST 在非比例風(fēng)險(xiǎn)場(chǎng)景下更穩(wěn)健貝葉斯模型則能給出參數(shù)的完整后驗(yàn)分布。這些模型彼此不是替代關(guān)系而是互補(bǔ)關(guān)系需要在同一份基礎(chǔ)數(shù)據(jù)上并行推進(jìn)再匯總判斷?？渴止し绞焦芾磉@種多方向分析效率低還容易出錯(cuò)。1.2 分支分析的核心思想關(guān)注過(guò) Git 工作流的朋友都知道“分支”的好處基于同一個(gè)主分支可以切出多個(gè) feature 分支各自開發(fā)互不阻塞最后再合并回主干。IntelligenR 里的分支分析思路與此高度相似只不過(guò)管理的對(duì)象從代碼變成了“分析任務(wù)”和“模型結(jié)果”。一套完整的 IntelligenR 分支分析通常包含三個(gè)層級(jí)基礎(chǔ)數(shù)據(jù)層所有分支共享同一份經(jīng)過(guò)校驗(yàn)的原始數(shù)據(jù)和分析基線。分析分支層每個(gè)分支對(duì)應(yīng)一個(gè)獨(dú)立分析方向比如 KM 分支、Cox 分支、RMST 分支、BRM 分支。結(jié)果匯總層分支完成后將關(guān)鍵結(jié)果匯總到統(tǒng)一的報(bào)告或?qū)Ρ缺碇泄Q策時(shí)取用。這樣做的好處是即時(shí)可見(jiàn)的每個(gè)分支獨(dú)立執(zhí)行模型 A 的改動(dòng)不會(huì)影響模型 B。同一數(shù)據(jù)可重復(fù)用于不同分析方向不需要復(fù)制多份數(shù)據(jù)集。分支結(jié)果可追溯哪個(gè)分支用了什么參數(shù)一目了然。遇到結(jié)果沖突時(shí)可以直接定位到具體分支排查不用從頭檢查整個(gè)流程。1.3 適用場(chǎng)景分支分析不是銀彈它特別適合以下幾類場(chǎng)景生存分析中的多模型對(duì)比KM、Cox、RMST、BRM 并行跑最后匯總。敏感性分析更換時(shí)間截點(diǎn)、變量集合、缺失值處理策略時(shí)開不同分支。團(tuán)隊(duì)協(xié)作統(tǒng)計(jì)分析員和數(shù)據(jù)科學(xué)家在同一份數(shù)據(jù)基礎(chǔ)上各自負(fù)責(zé)不同分支。分析流程復(fù)用同一套模板換一批數(shù)據(jù)自動(dòng)生成全套分支結(jié)果。如果你的分析任務(wù)只需要“一次性跑一個(gè)模型、得到一張表”那分支分析確實(shí)有點(diǎn)重但只要你的分析結(jié)果需要反復(fù)迭代、多方評(píng)審、結(jié)論要經(jīng)得起追問(wèn)分支分析的價(jià)值就非常明顯。2. 核心模型定位KM / Cox / RMST / BRM 分別解決什么問(wèn)題在進(jìn)入代碼之前先把四個(gè)模型的關(guān)系梳理清楚。這幾位出場(chǎng)頻率高但很多人對(duì)它們“什么時(shí)候用哪個(gè)”還是模糊的。2.1 KM 生存曲線先描述生存狀況Kaplan-Meier 估計(jì)是一種非參數(shù)方法用來(lái)估計(jì)生存函數(shù)。它不需要假設(shè)生存時(shí)間服從某個(gè)特定分布只會(huì)利用已知的事件和刪失數(shù)據(jù)計(jì)算出每個(gè)時(shí)間點(diǎn)的生存概率。輸出是一張階梯狀的生存曲線以及中位生存時(shí)間、某時(shí)間點(diǎn)生存率等描述性指標(biāo)。KM 主要用于探索和展示不同分組之間的生存曲線長(zhǎng)什么樣兩組曲線是否明顯分開它配合 log-rank 檢驗(yàn)可以給出組間差異的 p 值但它無(wú)法同時(shí)調(diào)整多個(gè)協(xié)變量也無(wú)法直接回答“某個(gè)因素每增加一個(gè)單位風(fēng)險(xiǎn)會(huì)變化多少”這種問(wèn)題。2.2 Cox 比例風(fēng)險(xiǎn)回歸量化風(fēng)險(xiǎn)因素Cox 回歸是生存分析里最常用的半?yún)?shù)模型。它不關(guān)心基線風(fēng)險(xiǎn)函數(shù)的具體形狀直接對(duì)協(xié)變量與風(fēng)險(xiǎn)之間的關(guān)系做回歸建模核心輸出是風(fēng)險(xiǎn)比 HRHazard Ratio。比如 HR 0.65意思是治療組的風(fēng)險(xiǎn)是對(duì)照組的 0.65 倍也就是風(fēng)險(xiǎn)降低了 35%。Cox 模型的優(yōu)勢(shì)在于可以同時(shí)納入多個(gè)協(xié)變量調(diào)整混雜因素。但它有一個(gè)重要前提就是比例風(fēng)險(xiǎn)假設(shè)PH 假設(shè)即各組的風(fēng)險(xiǎn)比在時(shí)間上保持恒定。如果 KM 曲線明顯交叉說(shuō)明不同時(shí)間區(qū)間內(nèi)風(fēng)險(xiǎn)優(yōu)勢(shì)在變化此時(shí) Cox 模型給出的單一 HR 可能會(huì)誤導(dǎo)結(jié)論。2.3 RMST非比例風(fēng)險(xiǎn)下的穩(wěn)健比較RMSTRestricted Mean Survival Time限制平均生存時(shí)間衡量的是在某個(gè)預(yù)先指定的時(shí)間窗口比如 5 年內(nèi)研究對(duì)象平均存活的時(shí)間。它不依賴比例風(fēng)險(xiǎn)假設(shè)即使生存曲線交叉也能給出一個(gè)直觀的組間差異估計(jì)比如“治療組在 5 年內(nèi)平均比對(duì)照組多活了 4.2 個(gè)月”。RMST 在臨床研究里經(jīng)常作為 Cox 模型的敏感性分析。如果 Cox 的 PH 假設(shè)不滿足或者你想給非統(tǒng)計(jì)背景的讀者講清楚“兩組到底差了多少”RMST 是更直觀的選擇。2.4 BRM 貝葉斯回歸模型把先驗(yàn)知識(shí)和不確定性帶進(jìn)分析BRM 在這里指貝葉斯回歸模型Bayesian Regression Model。在生存分析場(chǎng)景下就是對(duì)生存時(shí)間或風(fēng)險(xiǎn)函數(shù)建立一個(gè)貝葉斯回歸框架比如貝葉斯 Weibull 生存回歸、貝葉斯 Cox 模型等。貝葉斯方法的核心是“先驗(yàn)分布 數(shù)據(jù)似然 后驗(yàn)分布”。它的好處在于可以顯式納入先驗(yàn)知識(shí)比如參考?xì)v史研究認(rèn)為某藥物的 HR 在 0.6 到 0.8 之間。輸出不是單一點(diǎn)估計(jì)而是參數(shù)的后驗(yàn)分布能更自然地表達(dá)不確定性。在小樣本、事件數(shù)較少的情況下適當(dāng)使用信息先驗(yàn)可以降低模型不穩(wěn)定性。BRM 的代價(jià)是計(jì)算量大對(duì)使用者有更高的建模功底要求而且在團(tuán)隊(duì)協(xié)作里解釋成本也更高。因此它通常作為進(jìn)階驗(yàn)證手段而不是第一默認(rèn)模型。這四種模型其實(shí)是層層遞進(jìn)的KM 負(fù)責(zé)“看”Cox 負(fù)責(zé)“量”RMST 負(fù)責(zé)“穩(wěn)”BRM 負(fù)責(zé)“深”。分支分析恰好可以把它們組織在一個(gè)工作流里讓每個(gè)方向各司其職。3. 環(huán)境準(zhǔn)備與示例數(shù)據(jù)說(shuō)明3.1 軟件環(huán)境本文示例以 R 語(yǔ)言為主因?yàn)樯娣治鱿嚓P(guān)的生態(tài)最成熟。你需要準(zhǔn)備一個(gè)可用的 R 環(huán)境并安裝以下包install.packages(c(survival, survminer, survRM2, brms, dplyr))版本說(shuō)明survival是 R 自帶生態(tài)里的核心包survminer用于繪制 KM 曲線survRM2是跑 RMST 的經(jīng)典包brms用于貝葉斯生存回歸。具體版本不需要完全一致但建議保持較新的穩(wěn)定版文中代碼以常見(jiàn)版本的接口為準(zhǔn)如果你使用的版本有接口變化以官方文檔為準(zhǔn)。3.2 示例數(shù)據(jù)結(jié)構(gòu)為了演示分支分析我們先構(gòu)造一份模擬的臨床試驗(yàn)數(shù)據(jù)。假設(shè)我們想比較某種治療方案對(duì)患者生存時(shí)間的影響同時(shí)調(diào)整年齡和性別兩個(gè)協(xié)變量。數(shù)據(jù)字段如下字段含義id患者編號(hào)time生存時(shí)間單位月status事件狀態(tài)1 表示發(fā)生終點(diǎn)事件0 表示刪失group分組Treatment 治療組 / Control 對(duì)照組age年齡sex性別生成模擬數(shù)據(jù)的 R 代碼如下set.seed(2024) n - 300 group - factor(sample(c(Treatment, Control), n, replace TRUE)) age - rnorm(n, 60, 10) sex - factor(sample(c(Male, Female), n, replace TRUE)) true_hr - ifelse(group Treatment, 0.65, 1.0) scale_t - 50 / (true_hr^(1/2)) time - rweibull(n, shape 2, scale scale_t) cens - runif(n, 0, 80) status - ifelse(time cens, 1, 0) time - pmin(time, cens) data - data.frame(id 1:n, time time, status status, group group, age age, sex sex) write.csv(data, data/raw_data.csv, row.names FALSE)這段代碼的思路是先設(shè)定治療組真實(shí)風(fēng)險(xiǎn)比 HR 0.65然后基于 Weibull 分布生成生存時(shí)間再隨機(jī)生成一個(gè)刪失時(shí)間如果刪失時(shí)間早于生存時(shí)間則該樣本視為刪失。這樣做出來(lái)的模擬數(shù)據(jù)帶有真實(shí)效應(yīng)方便后面驗(yàn)證模型是否能識(shí)別出治療組的風(fēng)險(xiǎn)降低。3.3 IntelligenR 分支工作區(qū)結(jié)構(gòu)為了方便后續(xù)擴(kuò)展建議先建立一套清晰的分析目錄結(jié)構(gòu)intelligenr_branch_demo/ ├── data/ │ └── raw_data.csv ├── branches/ │ ├── km_analysis.R │ ├── cox_analysis.R │ ├── rmst_analysis.R │ └── brm_analysis.R ├── reports/ │ ├── km/ │ ├── cox/ │ ├── rmst/ │ └── brm/ ├── logs/ └── scripts/ └── run_all_branches.shdata放原始數(shù)據(jù)branches放各分支分析腳本reports按分支分別保存輸出結(jié)果logs記錄運(yùn)行日志scripts放批量執(zhí)行腳本。這樣的結(jié)構(gòu)本身就是“分支”思想在文件系統(tǒng)上的落地。4. 實(shí)戰(zhàn)在 IntelligenR 中建立四個(gè)分析分支下面進(jìn)入實(shí)操環(huán)節(jié)。我們將基于同一份模擬數(shù)據(jù)分別跑 KM、Cox、RMST、BRM 四個(gè)分支最后匯總結(jié)果。4.1 分支 AKM 曲線與 log-rank 檢驗(yàn)KM 分支是整個(gè)分析流程的地基。它的作用是描述性展示兩組生存曲線是否分開生存率隨時(shí)間怎么變化組間差異是否顯著。這是后續(xù)所有模型解釋的背景。文件路徑branches/km_analysis.Rlibrary(survival) library(survminer) data - read.csv(data/raw_data.csv) # 基于分組擬合 KM 生存曲線 fit_km - survfit(Surv(time, status) ~ group, data data) # 查看整體生存估計(jì) print(fit_km) # 輸出各時(shí)間點(diǎn)的生存率 summary(fit_km, times c(12, 24, 36, 48, 60)) # log-rank 檢驗(yàn)比較兩組生存曲線差異 logrank_test - survdiff(Surv(time, status) ~ group, data data) print(logrank_test) # 繪制 KM 曲線并疊加 p 值和風(fēng)險(xiǎn)表 pdf(reports/km/km_curve.pdf, width 8, height 6) ggsurvplot( fit_km, data data, pval TRUE, pval.method TRUE, risk.table TRUE, conf.int TRUE, xlab Time (months), ylab Overall Survival Probability ) dev.off()關(guān)鍵輸出會(huì)包括中位生存時(shí)間median survival timeTreatment 組假設(shè)為 39 個(gè)月Control 組假設(shè)為 28 個(gè)月12、24、36、48、60 個(gè)月的生存率。log-rank 檢驗(yàn) p 值如果小于 0.05說(shuō)明兩組生存分布有顯著差異。這里要提醒一句KM 曲線的 p 值只告訴我們“不同分組間生存分布有差異”并不等于“治療有效”因?yàn)樗鼪](méi)有調(diào)整年齡、性別等混雜因素。真正的因果性結(jié)論需要回到 Cox 或更嚴(yán)格的試驗(yàn)設(shè)計(jì)。4.2 分支 BCox 比例風(fēng)險(xiǎn)回歸Cox 分支回答的問(wèn)題是在同時(shí)調(diào)整年齡、性別的情況下治療組相對(duì)對(duì)照組風(fēng)險(xiǎn)降低了多少如果 HR 的可信區(qū)間不包含 1說(shuō)明效應(yīng)有統(tǒng)計(jì)學(xué)意義。文件路徑branches/cox_analysis.Rlibrary(survival) data - read.csv(data/raw_data.csv) # 擬合 Cox 回歸 fit_cox - coxph(Surv(time, status) ~ group age sex, data data) # 輸出模型結(jié)果 summary(fit_cox) # 提取 HR 和置信區(qū)間 hr_table - data.frame( variable names(coef(fit_cox)), HR exp(coef(fit_cox)), lower_95 exp(confint(fit_cox)[, 1]), upper_95 exp(confint(fit_cox)[, 2]), p_value summary(fit_cox)$coefficients[, 5] ) write.csv(hr_table, reports/cox/hr_table.csv, row.names FALSE) # 比例風(fēng)險(xiǎn)假設(shè)檢驗(yàn) ph_test - cox.zph(fit_cox) print(ph_test) pdf(reports/cox/ph_assumption.pdf, width 8, height 6) plot(ph_test) dev.off()在結(jié)果解讀上優(yōu)先看以下幾點(diǎn)groupTreatment對(duì)應(yīng)的 HR如果約為 0.65說(shuō)明治療組風(fēng)險(xiǎn)比對(duì)照組低約 35%。HR 的 95% 置信區(qū)間如果區(qū)間完全不包含 1說(shuō)明結(jié)果穩(wěn)定。cox.zph的 p 值如果某個(gè)變量的 p 值小于 0.05說(shuō)明該變量可能不滿足比例風(fēng)險(xiǎn)假設(shè)Cox 模型的結(jié)論需要謹(jǐn)慎此時(shí)應(yīng)去看 RMST 分支的結(jié)果。補(bǔ)充一點(diǎn)cox.zph全局檢驗(yàn)和單變量檢驗(yàn)都要看不要只看某一個(gè) p 值。如果 PH 假設(shè)不滿足常規(guī)做法是改用分層 Cox、時(shí)變系數(shù)模型或者直接用 RMST 作為主要分析。4.3 分支 CRMST 限制平均生存時(shí)間RMST 分支是 Cox 分支的重要補(bǔ)充尤其在生存曲線交叉或 PH 假設(shè)不成立時(shí)RMST 可以提供更穩(wěn)健的組間差異估計(jì)。文件路徑branches/rmst_analysis.Rlibrary(survRM2) data - read.csv(data/raw_data.csv) # 設(shè)定時(shí)間窗口 tau一般取臨床上有意義的時(shí)間點(diǎn) tau - 60 # 構(gòu)造治療的二元向量1 表示治療組0 表示對(duì)照組 arm - as.numeric(data$group Treatment) # RMST 分析 rmst_res - rmst2(time data$time, status data$status, arm arm, tau tau) # 輸出結(jié)果 print(rmst_res) # 手動(dòng)整理常用結(jié)果 rmst_table - data.frame( group c(Treatment, Control), rmst c(rmst_res$RMST.arm1$est[1], rmst_res$RMST.arm0$est[1]) ) write.csv(rmst_table, reports/rmst/rmst_summary.csv, row.names FALSE)rmst2的輸出核心是兩組 RMST 的差值RMST difference以及對(duì)應(yīng)的 95% 置信區(qū)間。如果差值為正且置信區(qū)間不包含 0說(shuō)明治療組在 60 個(gè)月內(nèi)的平均生存時(shí)間顯著長(zhǎng)于對(duì)照組。注意tau的設(shè)定要有臨床依據(jù)不能隨便取一個(gè)極大值。如果tau大于最大隨訪時(shí)間RMST 實(shí)際上就退化成整個(gè)隨訪期內(nèi)的平均生存時(shí)間可能會(huì)被右刪失影響不夠穩(wěn)健。通常建議把tau設(shè)為一個(gè)大多數(shù)患者都已經(jīng)出現(xiàn)終點(diǎn)事件的時(shí)點(diǎn)比如本研究中的 60 個(gè)月。4.4 分支 DBRM 貝葉斯生存回歸BRM 分支走的是貝葉斯路線。我們會(huì)用brms包擬合一個(gè)貝葉斯 Weibull 生存回歸把治療組、年齡、性別納入模型得到每個(gè)參數(shù)的后驗(yàn)分布。文件路徑branches/brm_analysis.Rlibrary(brms) data - read.csv(data/raw_data.csv) # brms 中 status1 表示事件需要轉(zhuǎn)換為刪失指示變量 data$censored - ifelse(data$status 1, 0, 1) # 貝葉斯 Weibull 生存回歸 fit_brm - brm( bf(time | cens(censored) ~ group age sex), data data, family weibull(), chains 2, cores 2, iter 2000, seed 2024 ) # 查看后驗(yàn)匯總 summary(fit_brm) # 提取關(guān)鍵參數(shù)的后驗(yàn)樣本 posterior_samples - as.data.frame(fit_brm) write.csv( summary(fit_brm)$fixed, reports/brm/fixed_effects_summary.csv, row.names TRUE ) # 繪制治療組效應(yīng)的后驗(yàn)分布 pdf(reports/brm/posterior_plot.pdf, width 8, height 6) bayesplot::mcmc_areas(fit_brm, pars b_groupTreatment) dev.off()在解釋貝葉斯結(jié)果時(shí)不要套用頻率學(xué)派的 p 值表達(dá)而是關(guān)注后驗(yàn)均值b_groupTreatment表示治療效果在 Weibull 尺度上的回歸系數(shù)負(fù)值意味著治療組的風(fēng)險(xiǎn)更低、生存時(shí)間更長(zhǎng)。90% 或 95% 可信區(qū)間如果區(qū)間不跨越 0說(shuō)明后驗(yàn)證據(jù)比較明確。樣本量較小時(shí)后驗(yàn)分布的尾部信息可以告訴你結(jié)果有多不確定。如果你對(duì)貝葉斯模型不熟悉建議先在經(jīng)典模型Cox上確認(rèn)結(jié)果方向再用 BRM 做穩(wěn)健性驗(yàn)證。貝葉斯模型的計(jì)算更慢鏈路更多適合作為“驗(yàn)證分支”而不是“起步分支”。這里需要說(shuō)明一點(diǎn)brms的刪失參數(shù)約定可能隨版本變化cens()的具體寫法請(qǐng)以你安裝版本的官方文檔為準(zhǔn)。如果brms安裝不便也可以先用rstanarm::stan_surv或者使用 JAGS 手寫 Weibull 生存模型思路完全相同。4.5 分支并行執(zhí)行與結(jié)果匯總四個(gè)分支腳本都準(zhǔn)備好后可以寫一個(gè)批量執(zhí)行腳本讓它們依次或并行運(yùn)行。并行執(zhí)行能省時(shí)間但要注意每個(gè)分支的日志和輸出目錄必須隔離避免互相覆蓋。文件路徑scripts/run_all_branches.sh#!/bin/bash set -e mkdir -p reports/km reports/cox reports/rmst reports/brm mkdir -p logs echo Run KM Branch Rscript branches/km_analysis.R logs/km.log 21 echo Run Cox Branch Rscript branches/cox_analysis.R logs/cox.log 21 echo Run RMST Branch Rscript branches/rmst_analysis.R logs/rmst.log 21 echo Run BRM Branch Rscript branches/brm_analysis.R logs/brm.log 21 echo All branches completed.如果你的機(jī)器是多核環(huán)境也可以使用parallel或 R 的future包做并行調(diào)度。但要注意貝葉斯分支本身已經(jīng)占用了多個(gè) CPU 核再和其他分支并行可能會(huì)造成資源競(jìng)爭(zhēng)反而變慢。結(jié)果匯總階段四個(gè)分支產(chǎn)出如下分支核心輸出文件位置KM生存曲線圖、中位生存時(shí)間、log-rank p 值reports/km/CoxHR 表、PH 假設(shè)檢驗(yàn)結(jié)果reports/cox/RMST兩組 RMST、RMST 差值與置信區(qū)間reports/rmst/BRM后驗(yàn)分布、可信區(qū)間、后驗(yàn)診斷圖reports/brm/5. 分支對(duì)比、合并與沖突處理多個(gè)分支跑完最關(guān)鍵的一步是“對(duì)比和匯總”。如果四個(gè)分支結(jié)論一致那結(jié)論就比較穩(wěn)固但如果分支之間出現(xiàn)了沖突怎么辦5.1 四類模型結(jié)果之間的互補(bǔ)性一個(gè)比較理想的結(jié)果應(yīng)該是這樣的KM 曲線顯示治療組生存曲線整體在對(duì)照組上方log-rank p 值 0.05。Cox 模型顯示治療組 HR 0.6595% 置信區(qū)間為 0.48 ~ 0.87PH 假設(shè)檢驗(yàn)不顯著。RMST 顯示治療組 60 個(gè)月 RMST 比對(duì)照組高 4.5 個(gè)月95% 置信區(qū)間不包含 0。BRM 后驗(yàn)分布顯示治療組效應(yīng)參數(shù)的后驗(yàn)均值遠(yuǎn)離 0可信區(qū)間不含 0。這種“四方結(jié)論一致”的情況最理想可以直接寫進(jìn)結(jié)論部分。5.2 結(jié)果沖突怎么辦如果 KM 曲線交叉、Cox 卻給出了顯著的 HR這就進(jìn)入了一個(gè)高頻沖突場(chǎng)景。此時(shí)應(yīng)該優(yōu)先懷疑 PH 假設(shè)是否滿足。KM 曲線交叉往往說(shuō)明不同時(shí)間段內(nèi)兩組的風(fēng)險(xiǎn)比在變化Cox 模型給出的“平均 HR”可能掩蓋了這種變化。此時(shí)應(yīng)該查看cox.zph的 p 值確認(rèn)是否違反 PH 假設(shè)。改看 RMST 分支用“某個(gè)時(shí)間窗內(nèi)的平均生存時(shí)間差異”做組間比較。對(duì) Cox 分支做分層分析或引入時(shí)間交互項(xiàng)而不是直接刪掉有問(wèn)題的變量。還有一種沖突是頻率學(xué)派和貝葉斯學(xué)派結(jié)論不一致Cox 給出的 p 值大于 0.05但 BRM 的后驗(yàn)可信區(qū)間不包含 0。這種情況通常和先驗(yàn)設(shè)置有關(guān)也可能是因?yàn)閮膳蓹z驗(yàn)邏輯本身不同。處理方式是把先驗(yàn)設(shè)置寫清楚并在報(bào)告中同時(shí)呈現(xiàn)兩套結(jié)果由讀者或評(píng)審專家判斷。5.3 分支合并與回退規(guī)則分支分析不是所有分支都必須合并。你可以根據(jù)分析目的制定簡(jiǎn)單的規(guī)則場(chǎng)景處理方式所有分支結(jié)論一致合并所有分支結(jié)論寫入最終報(bào)告Cox 和 RMST 沖突以 RMST 為主Cox 作為敏感性分析BRM 分支不收斂檢查先驗(yàn)和采樣參數(shù)必要時(shí)回退到經(jīng)典模型KM 曲線交叉放棄單一 HR 的解釋采用時(shí)變效應(yīng)或 RMST分支管理的核心原則是寧可保留多個(gè)分支的分歧結(jié)果也不要為了得到一個(gè)“漂亮結(jié)論”強(qiáng)行統(tǒng)一。6. 常見(jiàn)問(wèn)題與排查思路在實(shí)際跑 IntelligenR 分支分析時(shí)大家最常遇到的還是環(huán)境和模型細(xì)節(jié)問(wèn)題。下面整理一份排查清單。問(wèn)題現(xiàn)象常見(jiàn)原因解決思路中文標(biāo)簽在圖里顯示為方塊中文字體缺失在 R 里配置中文字體或把圖片內(nèi)文字改為英文cox.zph檢驗(yàn) p 值 0.05比例風(fēng)險(xiǎn)假設(shè)不成立改用分層 Cox、時(shí)變系數(shù)或直接看 RMST 分支RMST 的 tau 設(shè)太大最遠(yuǎn)隨訪時(shí)間前事件率不足調(diào)整 tau 為臨床有意義的截點(diǎn)如 48 或 60 個(gè)月貝葉斯模型采樣警告鏈長(zhǎng)不足或先驗(yàn)不合適增加 iter調(diào)整先驗(yàn)檢查 Rhat 是否接近 1分支之間結(jié)果不一致數(shù)據(jù)版本或隨機(jī)種子不一致確保所有分支使用同一份 raw_data.csv固定 set.seed事件數(shù)太少模型不穩(wěn)定樣本量不足或刪失率過(guò)高合并分組、減少協(xié)變量、考慮容錯(cuò)更強(qiáng)的貝葉斯模型并行執(zhí)行時(shí)內(nèi)存不足多個(gè) R 進(jìn)程同時(shí)吃內(nèi)存限制并行度BRM 分支單獨(dú)執(zhí)行這些問(wèn)題的共性是先定位是數(shù)據(jù)問(wèn)題、模型問(wèn)題還是環(huán)境問(wèn)題再對(duì)癥下藥不要一上來(lái)就改模型公式。7. 最佳實(shí)踐與工程化建議7.1 分支命名與分析記錄分支的名字要能體現(xiàn)分析方向、數(shù)據(jù)版本和參數(shù)變更。比如“cox_v2_adjust_age_sex”“rmst_tau60_sensitivity”。每次跑完分支在日志里記錄時(shí)間、數(shù)據(jù) hash、依賴包版本和運(yùn)行結(jié)果摘要方便后面追溯。7.2 數(shù)據(jù)不可變?cè)瓌t分支分析的前提是“基礎(chǔ)數(shù)據(jù)一致”。在實(shí)際操作中建議洗好一份raw_data.csv后就不允許再改動(dòng)。任何變量轉(zhuǎn)換、缺失值填充都放到分支腳本內(nèi)部做這樣分支之間不會(huì)互相污染。如果確實(shí)要修改基礎(chǔ)數(shù)據(jù)請(qǐng)生成一個(gè)新的數(shù)據(jù)版本而不是在原文件上改。7.3 隨機(jī)種子與可復(fù)現(xiàn)性所有涉及隨機(jī)數(shù)的環(huán)節(jié)包括模擬數(shù)據(jù)、樣本拆分、貝葉斯采樣都要顯式設(shè)置set.seed()這樣別人拿到你的代碼能復(fù)現(xiàn)出相同結(jié)果。貝葉斯模型還要額外記錄采樣器參數(shù)比如 chains、iter、warmup、adapt_delta這些都會(huì)影響結(jié)果穩(wěn)定性。7.4 輸出格式統(tǒng)一每個(gè)分支的輸出建議采用統(tǒng)一的命名和格式。比如所有表格都輸出 CSV所有圖片都輸出 PDF并且以模型名_內(nèi)容_日期命名。匯總時(shí)不僅給數(shù)字還要給出置信區(qū)間或可信區(qū)間避免只列點(diǎn)估計(jì)。7.5 合規(guī)與倫理提醒如果你的數(shù)據(jù)來(lái)自臨床研究或真實(shí)患者信息務(wù)必確認(rèn)分析方案已通過(guò)倫理審批數(shù)據(jù)使用已獲得授權(quán)并且分析過(guò)程中遵守最小權(quán)限和數(shù)據(jù)脫敏原則。分支分析只是技術(shù)手段它不能代替合規(guī)審查。任何時(shí)候都不應(yīng)該在未授權(quán)數(shù)據(jù)上開展分析也不應(yīng)該把內(nèi)部數(shù)據(jù)結(jié)果隨意對(duì)外發(fā)布。8. 總結(jié)與下一步學(xué)習(xí)路線本文圍繞 IntelligenR 分支分析展開把 KM、Cox、RMST、BRM 四種生存分析模型組織成了四條互不干擾的分析分支。前半部分幫你建立概念框架搞清四個(gè)模型各自的適用場(chǎng)景后半部分用一份模擬數(shù)據(jù)完整演示了從數(shù)據(jù)清洗、分支腳本編寫、批量執(zhí)行到結(jié)果匯總的全流程。下一步你可以從三個(gè)方向繼續(xù)深入如果對(duì)生存分析本身還不熟悉建議先吃透 KM 曲線和 Cox 模型的輸出解讀再嘗試 RMST 和貝葉斯模型。如果想把分支分析做得更自動(dòng)化可以學(xué)習(xí)targets、drake這類 R 工作流管理包把分支依賴關(guān)系寫得更規(guī)范。如果對(duì)貝葉斯生存模型感興趣可以系統(tǒng)學(xué)習(xí)brms的刪失數(shù)據(jù)建模、先驗(yàn)敏感性分析以及后驗(yàn)預(yù)測(cè)檢查。在真實(shí)項(xiàng)目中最優(yōu)先要關(guān)注的永遠(yuǎn)不是模型有多復(fù)雜而是數(shù)據(jù)質(zhì)量、模型假設(shè)和結(jié)果可復(fù)現(xiàn)性。分支分析給你的是一套“多條腿走路”的安全機(jī)制但最終能不能站住還是要看每條腿是否都踩在扎實(shí)的數(shù)據(jù)和嚴(yán)謹(jǐn)?shù)慕y(tǒng)計(jì)假設(shè)上。如果你在跑分支、匯總結(jié)果或解讀模型時(shí)遇到其他報(bào)錯(cuò)歡迎帶著具體報(bào)錯(cuò)和模型輸出繼續(xù)交流。