R語言實例解析腫瘤治療藥物經濟學評價的分區生存模型

劉家一，李 偉

1.中國人民大學公共管理學院（北京 100872）

2.北京醫院麻醉科（北京 100005）

分區生存模型（partitioned survival model，PSM ） 是國外腫瘤治療藥物經濟學（pharmacoeconomics，PE）評價中應用較多的一種決策模型，用于腫瘤藥物的報銷決定[1-2]。PSM已在英國國家卓越健康與護理研究所（National Institute for Health and Care Excellence，NICE） 技術評估（Technology Appraisal，TA）計劃和加拿大藥物與衛生技術局（Canadian Agency for Drugs and Technologies in Health，CADTH）評價中廣泛使用，是目前NICE 評估晚期或轉移性癌癥干預措施的最常用方法[1,3]。與國內研究者常用的狀態轉移模型相比，PSM 具有簡潔、直觀的特點，它直接利用腫瘤治療藥物臨床報告中的無進展生存期（progression free survival，PFS）、總生存期（overall survival，OS）數據，避免在數據有限的情況下計算并驗證復雜的狀態間轉移概率[4-8]。

在PSM 應用中，盡管TreeAge 等軟件在建模方面具有一定優勢，但其固有設置要求定義狀態間轉移概率，限制了其在PSM 中的應用[9]。相較之下，R 語言在藥物成本效用分析中靈活性更強，專用程序包可直接應用PSM，同時提供了豐富的數據處理與可視化功能。國際研究者已在不斷嘗試和探索新的模型結構和變量選擇方法過程中產生了不同功能的專用程序包，如健康經濟模型分析的heemod[10-11]、提供靈活生存分析方法的flexsurv[12]、集成了健康經濟模型模擬分析的hesim[13]、標準生存數據分析的survival[14]、結合了生存分析與健康經濟評估的survHE[15-16]、執行貝葉斯成本效益分析的BCEA[17]等。利用R 語言進行腫瘤治療藥物的PSM 分析，能夠快速、直觀地復現結果，在缺乏患者個體數據（individual patient data，IPD）的情況下，研究者可直接提取臨床試驗報告中公開的PFS、OS 曲線數據[13,18]，以重新構建IPD 并完成整個PE 評價，但這要求研究者具有數理統計、健康經濟分析與編程的綜合知識背景，國內采用R 語言進行PSM 分析的應用報道不多，且缺乏具體實施案例。本研究以帕博利珠單抗治療非小細胞肺癌（non-small cell lung cancer，NSCLC）的 KEYNOTE189 臨床試驗（簡稱KN189）[19-20]為例，實現從數字化提取PFS、OS 曲線數據，然后由R 語言實現PE 分析與評價的全過程，以期為研究人員利用R 語言進行PE評價提供參考。

1 臨床試驗數據提取

基于帕博利珠單抗Pembrolizumab KN189臨床試驗數據[19-20]（圖1、圖2），使用GetDataW 2.26.0.20 軟件分別對其OS、PFS 曲線數據進行提取，同時提取number at risk 數據，一并保存為xlsx 文件。針對提取后數據進行整理，使后一個時間點的生存率小于或等于前一個時間點數值。

圖1 KN189臨床試驗OS曲線數據Figure 1.OS curve data of KN189

圖2 KN189臨床試驗PFS曲線數據Figure 2.PFS curve data of KN189

2 重構IPD

2.1 R語言環境設置

打開Rstudio 2023.09.1+494 軟件，設定工作路徑，安裝并下載軟件包：survial（v 3.5-7）、survHE（v 2.0.1）、IPDfromKM（v 0.1.10）、readbitmap（v 0.1.5）、survminer（v 0.4.9）、flexsurv（v 2.2.2）、data.table（v 1.15.2）、hesim （v 0.5.4）、openxlsx（v 4.25.2）。

2.2 定義數據集

基于KN189 試驗方案，將對照組命名為cemo，帕博利珠組定義為comb，OS 曲線命名為surv，PFS 曲線命名為pfs。分別通過read.xlsx()函數讀取第1 步中提取的xlsx 文件，整理并保存為如下數據集：survcomb、survcemo、nriskcomb、nriskcemo、pfscomb、pfscemo、nsirkcomb，表1、表 2分別展示了survcomb 與nriskcemo 數據集的數據結構。此處僅示例OS 曲線的數據重構過程，PFS 曲線與OS 曲線的數據提取過程類似。

表1 KN189試驗comb組OS曲線數據Table 1.OS curve data of the comb group in KN189

表2 KN189試驗cemo組OS曲線風險人數數據Table 2.OS curve data of the cemo group in KN189

2.3 預處理與重構IPD數據

2.3.1 預處理IPD對象格式

使用preprocess()函數預處理為IPD 對象適當格式，R 語言代碼如下：

pre_comb＜- preprocess (dat=survcomb,trisk=nriskcomb$time, nrisk=nriskcomb$nrisk,totalpts=410)

totalpts 為指定試驗的初始人數；

pre_cemo＜- preprocess (dat=survcemo,risk=nriskcemo$time, nrisk=nriskcemo$nrisk,totalpts=206)

206 為cemo 組試驗的初始人數。

2.3.2 重構IPD數據

使用getIPD 函數重建comb 組OS 曲線IPD數據，運行代碼如下：

est_comb＜ getIPD (prep=pre_comb, armID=1,tot.event=NULL)

summary (est_comb)

查看est_comb 數據集可知，采用IPDfromKM包的Kolmogorov_Smirnov(K-S) Test，得到D 值為0.053，p_value 為0.99，預測值est_comb 與read _in的臨床數據pre_comb 數據無顯著性差異，可進行后續分析，代碼如下：

est_cemo＜-getIPD (prep=pre_cemo, armID=2,tot.event=NULL)

同理獲得cemo 組OS 曲線IPD 數據，armID定義為2。

2.3.3 對比重構IPD數據與KN189的HR值

重構IPD 數據與KN189 報告的風險比（hazard ratio，HR）值對比分析的R 語言運行代碼如下：

est_IPD_comb_cemo＜-rbind (est_comb$IPD,est_cemo$IPD)

提取comb、cemo 組的重構IPD 為新數據集；

hr_est_IPD＜-coxph (Surv(time,status)～treat,data=est_IPD_comb_cemo)

計算重構IPD數據中comb、cemo組的HR值；

exp (coef (hr_est_IPD)

運行并提取HR 值為0.57，與KN189 報告中數值0.56 相近。

2.4 繪制重構IPD數據與原read_in曲線對比檢驗

使用plot 繪制comb 組與cemo 組重構IPD 數據與read_in 的臨床數據比較圖，圖3 顯示兩組數據重合度較高。若重構IPD 與read_in 曲線對比圖可見明顯差異，則需從第1 步開始重復以上開始步驟，并檢查原始數據，具體代碼如下：

圖3 cemo組重構IPD數據OS曲線與read-in KN189數據OS曲線對比Figure 3.Comparison of OS curves of reconstructed IPD data and read-in KN189 data in cemo group

plot(est_comb)

Estimated 為重構IPD 數據，Reported 為read_in KN189 臨床數據；

plot(est_cemo)

如圖3 所示，圖3-A 為基于重構IPD 數據的OS 曲線與數字化讀取的KN189 臨床試驗的OS 曲線的視覺對比圖。

將兩組IPD 數據分別加上treat 列并分別賦值為1、2，然后將comb 與cemo 組的重構IPD 數據整合為數據框（表3），具體代碼如下：

表3 OS曲線重構IPD數據Table 3.Reconstructed IPD data of OS curve

IPD_comb_cemo_os

surcombIPD＜- transform (est_comb$IPD,treat=1)

將治療方案comb 組定義為1；

surcemoIPD＜- transform (est_cemo$IPD,treat=2)

將治療方案cemo 組定義為2；

IPD_comb_cemo_os＜- rbind (surcombIPD,surcemoIPD)

重復以上2.2 至2.4 步驟，得到OS 曲線的重構數據框IPD_comb_cemo_os。

3 生存曲線參數分布擬合與評估

定義待擬合模型名稱的向量，同時考慮指數、威布爾、伽瑪、對數正態、對數邏輯及廣義伽瑪多個模型，運行代碼如下：

mods1＜-c (＂exponential＂, ＂weibull＂, ＂gamma＂,＂lognormal＂, ＂loglogistic＂, ＂gengamma＂)

定義模型公式，采用Surv()函數，以數據框中分類變量treat 作為協變量分析，具體代碼如下：

formula＜- Surv (time, status)～ as.factor (treat)

執行批量生存模型分析，采用最大似然估計（maximum likelihood estimates，MLE），使 用flexsurv 包fit.models()函數創建survHE 對象KN189case_mle，其中存儲了所考慮的多個參數模型的擬合結果，代碼如下：

KN189case_mle＜-fit.models (formula=formula,data=IPD_comb_cemo_os, distr=modsl, method=＂mle＂)

KN189case_mle$model.fitting

顯示批量擬合的結果AIC、BIC、DIC 值。

基于AIC、BIC 較小者模型擬合佳的原則，本例結果顯示AIC、BIC 值最小的擬合為loglogistic，即loglogistic 模型擬合度最優。同樣步驟進行MLE 批量擬合IPD_comb_cemo_pfs，結果也顯示loglogistic 擬合模型的AIC、BIC 值最小，擬合度最優。

4 構建生存、效用、費用模型

4.1 定義模型

調用R 軟件hesim 包的hesim_data 函數[13]，用以標準化后續生存效用模型的數據準備，代碼如下：

hesim_dat＜- hesim_data (strategies =data.table (strategy_id=1 : 2, strategy_name=c (＂comb＂,＂cemo＂))，

patients=data.table (patient_id=1, grp_id=1),states =data.table (state_id=1 : 2, state_name=c(＂Stable＂, ＂Progression＂)))

分別定義分區生存模型中治療方案strategies、患者類型patients、試驗疾病狀態states，死亡狀態不用列出。

labs＜- get_label (hesim_dat)

4.2 建立生存模型

將IPD_comb_cemo_os、IPD_comb_cemo_pfs數據合并，并整理為數據框surv_est_data，整理后的數據結構見表4。

表4 comb與cemo組OS及PFS曲線的重構患者數據Table 4.Reconstructed patient data of OS and PFS curves in comb group and cemo group

采用flessurvreg()函數，基于重構IPD 數據集surv_est_data，分別擬合OS 及PFS 生存模型，代碼如下：

flex_fit_os＜- flexsurvreg (Surv (os_time, os_status) ～strategy_name, data=surv_est_data,dist=“llogis”)

選擇步驟3 中pfs 及os 曲線最優擬合模型llogis：

flex_fit_pfs＜- flexsurvreg (Surv(pfs_time, pfs_status)～strategy_name, data=susurv_est_data,dist=＂llogis＂)

綜合flex_fit_os 及flex_fit_pfs 模型為psfit_wei，代碼如下：

psfit_wei＜- flexsurvreg_list (flex_fit_pfs, flex_fit_os)

包含了3 個狀態分別的生存擬合結果對象。

建立生存模型survmods，調用hesim 程序包的create_PsmCurves()函數。在其參數設定中，明確boostrap 為確定生成結果不確定性的方法，n為boostrap 的樣本數，運行代碼如下：

n_samples＜- 300

設定生存模型中將要模擬的樣本數量；

surv_input_data＜- expand(hesim_dat,by=c(＂strategies＂, ＂patients＂）

從hesim_dat 中指定參數進行模擬；

survmods＜- create_PsmCurves (psfit_wei, n=n_samples, input_data=surv_input_data,

uncertainty=＂bootstrap＂, est_data=＂surv_est_data＂)

input_data: hesim 數據框架，包含了需要生成生存函數的所有患者和策略組合形式；n 為用于bootstrap 的樣本數；uncertainty 用于確定生成結果的不確定性的方法，本例使用bootstrap 方法；est_data: hesim 數據集，類似于IPD 生存狀態與時間數據及其他如年齡等所有額外信息。

4.3 建立效用模型

疾病狀態的效用值采用文獻數據[21]，即pfs=0.815、pd=0.321，其效用值參數的變化范圍：DSA 在+15%范圍，PSA 為beta 分布。據此定義狀態效用值數據框（表5），運行代碼如下：

表5 不同狀態效用值Table 5.Utility values for different states

utility_tbl＜- stateval_tbl (utility_data,dist=＂beta＂)

utilitymod＜- create_StateVals (utility_tbl, n=n_samples, hesim_data=hesim_dat)

4.4 建立費用模型

根據KN189 臨床設計方案，針對使用藥品費用（包括疾病進展后使用二線治療方案藥物、不同方案治療比例）、檢查檢驗的周期計劃及費用、不良反應發生率及相應治療費用、臨終關懷費用等，采用藥智網及北京市醫療費用收費標準進行綜合計算，得出comb、cemo 組各自總體費用的均值（表6），運行代碼如下：

表6 不同治療方案費用Table 6.The cost of different plans

cost_strategy_id

cost_tbl＜- stateval_tbl (cost_strategy_id,dist=＂fixed＂)

costmod＜- create_StateVals (cost_tbl, n=n_samples, hesim_data=hesim_dat)

5 實例化PSM并計算質量調整生命年

5.1 定義模型

生存模型survmods 與后續的效用模型、費用模型匯總進行PSM 的實例化，并進行時間外推，分析代碼如下：

psm＜- Psm$new (survival_models=survmods,utility_model=utilitymod, cost_models=list(medical=costmod))

times＜- seq (0, 1000, by=50)

定義生存模型實例化外推時間的范圍為20年，與前數據統一以周為單位，將時間范圍帶入psm 進行生存預測；

psm$sim_survival (t=times)

將定義的外推時間輸入psm 進行生存預測。

5.2 繪制外推生存曲線圖

基于以上所定義的外推時間范圍times 內，將實例化PSM 中的生存概率繪制生存曲線，結果見圖4，分析代碼如下：

圖4 外推時間內的PFS和OS曲線Figure 4.PFS and OS curves over time

labs$curve＜- c(＂PSF_Curve＂=1, ＂OS_Curve＂=2)

autoplot (psm$survival_, lables=labs, ci=TRUE,ci_style=＂ribbon＂,

scale_x_continuous (breaks=seq (o, max (times), 50))

psm$sim_stateprobs()

整理狀態概率數據框，用于后續視圖表現；

stateprobs＜-psm$stateprobs_[sample==2&patient_id==1]

從psm 對象中任選一個子集，選擇sample=2且patient_id=1，將其狀態概率定義給到stateprobs；

stateprobs [, state: =factor (state_id, state,levels=rev (unique (state_id)))]

在stateprobs 數據框中增加state 列，將state_id 轉換為因子類型的變量賦值給state 列，并指定其級別為唯一值且逆序排列；同理增加strategy 列如下：

stateprobs [, strategy: =factor(strategy_id,labels=c (hesim_dat$strategies$srategy_name))]

基于整理的stareprobs 繪制狀態概率曲線，見圖5，分析代碼如下：

圖5 外推時間內的狀態曲線Figure 5.State curves over time

ggplot (stateprobs [strategy_id % in% c (1 : 2)]

指定數據來源stateprobs 中strategy_id 列數據繪圖；

aes(x=t, y=prob, fill=state, group=state))+

定義x、y 軸分，根據state 列不同值著色分組；

geom_area (alpha-.65), facet_wrap(～strategy)

將圖形分面，每個分面表示不同的策略(strategy)；

abs (x=＂Periods＂, y=＂Proportion in state＂)

設定x、y 軸的文字標識；

scale_fill_manual (name=＂Health state＂

設置圖例的標題；

values=c(＂gray＂, ＂red＂, ＂blue＂)

設置填充圖形的顏色映射；

lables=c (＂Death＂, rev (hesim_dat$states$stat e_name)))+

標識每個狀態的文字內容；

scale_x_continuous (breaks=seq(0, max(times), 50)+

設置x 軸的刻度標簽；

guides (fill= guide_legend (reverse=T, nrow=2,byrow=TRUE))+

調整圖例以逆序的方式且排2 行；

theme (legend.position=＂bottom＂)

設置圖例的位置在底部。

5.3 模擬成本與質量調整生命年

基于5.2 中的狀態概率進行數值積分模擬貼現成本與質量調整生命年（quality-adjusted life years，QALYs），選擇5%進行貼現，分析代碼如下：

psm$sim_costs (dr=0.05)

運行后結果見圖6。

圖6 基于PSM以5%的折現率模擬的貼現成本Figure 6.Discounted cost simulated based on PSM at a 5% discount rate

psm$sims_qalys (dr=0.05)

運行后結果見圖7。

圖7 基于PSM以5%的折現率模擬的貼現效用Figure 7.Discounted utilities simulated based on PSM at a 5% discount rate

6 視圖輔助決策的成本效用分析

6.1 創建模擬數據

調用hesim 包的cea()和cea_pw()函數分別執行單一或成對的兩種治療策略的成本效益分析，其運行結果對象中包含有每個PSA 樣本的平均成本與QALYs。先用summarize()方法創建一個具有平均值的hesim::ce 對象，并將其賦值給變量ce_ sim，該對象包括了在概率敏感分析PSA 中使用的模擬結果的摘要統計信息，分析代碼如下：

ce_sim＜- psm$summarize()

wtp＜- seq (from=0, to=20 000, by=2 000)

定義支付意愿閾值wtp 的范圍，以此進行CUA 分析；

cea_out＜- cea (ce_sim, dr_qalys=.05,dr_ costs=.05, k=wtp)

將cea()函數分析結果賦值給變量cea_out；

cea_pw_out＜- cea_pw (ce_sim, comparator=1,dr_qalys=.05, dr_costs=.05, k=wtp)

cea_pw()函數計算人群加權增量成本-效果比（populattion-weighted incremental costeffectiveness ratio，PWICER），comparator=2 即將strategy_id=2 的cemo 組為基準策略分析偏好權重，K 為支付意愿閾值。

6.2 繪制成本效果可接受曲線

從cea_out 中提取strategy_id=1 的comb 組成本效果數據，并繪制單組的成本效果可接受曲線，分析代碼如下：

strategy_1_data＜- cea_out$mce [strategy_id==1]

cea_ac_plot＜- ggplot (data=strategy_1_data,aes(x=strategy_1_data$k, y=strategy_1_data$prob,group=1))+

geom_line (linewidth=1), xlab (＂k＂)+ylab(＂prob＂)+

ggtitle (＂Cost-Effectiveness Acceptability Curve＂), theme_minimal()

提取strategy_id=2 的cemo 組數據，將兩組數據合并后，同時繪制cemo 與comb 組的成本效果可接受曲線：

strategy_2_data＜- cea_out$mce [strategy_id==2]

combined_data＜- rbind (strategy_1_data,strategy_2_data)

colors＜- c (＂blue＂, ＂red＂)

labels＜- c (＂Strategy 1＂, ＂Strategy 2＂)

cea_ac_plot＜- ggplot (data=combined_data, aes(x=k, y=prob, color=factor (strategy_id))) +

geom_line (linewidth=1), xlab (＂Willingness to pay＂) , ylab (＂Probability Acceptable＂)+

scale_color_manual (values=colors,lables=lables),

guides (color=guide_legend (title=＂Strategy＂))+

ggtitle (＂Cost-Effectiveness Acceptability Curve＂), theme_minimal()

print(cea_ac_plot)

由圖8 可知，當支付意愿值超過4 800 元時，comb 組具有成本效益的可能性開始大于cemo 組。

圖8 comb組與cemo組成本效果可接受曲線Figure 8.Acceptable cost-effectiveness curves of the comb group and the cemo group

6.3 繪制增量成本-效益比值散點圖

繪制增量成本-效益比值（incremental costeffectiveness ratio，ICER）散點圖，設定支付意愿值上下限并在圖中添加97.5%和2.5%可能性的正比例函數線條，運行代碼如下：

strategy_1＜- cea_pw_out$delta [strategy_id==1, ]

提取cea_pw_out 數據集中comb 組的數據；

startegy_1$icer＜- strategy_1$ic/strategy_1$ie

計算comb 組的ICER；

wtp_2＜- quantile (strategy_1$icer, probs=0.975)

wtp_3＜- quantile (strategy_1$icer, probs=0.025)

cea_sensitivity_plot＜- ggplot (data=cea_pw_out$delta,

aes (x=ie, y=ic, color=strategy_id, size=3,alpha=0.8))+

geom_point (color=＂red＂, size=3, alpha=0.8)+

geom_abline (slope=wtp_2, intercept=0,color=＂blue＂, linetype=＂dashed＂)+

geom_abline (slope=wtp_3, intercept=0,color=＂red＂, linetype=＂dashed＂)+

geom_text (data=data.frame (x=0.2, y=3 000,lable=paste0 (＂97.5% ICER=＂, round (wtp_2,2))) ,

geom_text (data=data.frame (x=1.8, y=3 000,lable=paste0 (＂2.5%ICER=＂, round (wtp_3,2))),

xlab(＂Incremental Effectiveness＂),ylab(＂Incremental Cost＂)+

ggtitle (＂Sensitivity Analysis＂), xlim (0, 1.25),ylim (1650000, 1800000)+

theme_minimal()

print (cea_sensitivity_plot)

由圖9 可知，當意愿支付閾值為14 312 時，comb 組干預方案具有成本-效益的可能性為97.5%；當意愿支付閾值為3 123 時，comb 組干預方案具有成本-效益的可能性僅為2.5%。

圖9 comb組Bootstrap自助法抽樣300次后得到的ICER散點圖Figure 9.Scatter plot of ICER obtained after 300 iterations of Bootstrap resampling in the comb group

7 討論

本研究展示了R 語言在成本效用分析中的高效性，具體表現在：一是代碼透明性，分析過程完全通過代碼記錄，確保模型運行的透明性；二是計算可重復性，僅需重新運行代碼即可重復整個計算過程并更新結果；三是錯誤識別與調整便捷性，模型設計及參數設定的錯誤或假設調整易于發現與更正；四是結果輸出直接性，PE 評價結果可直接通過R 代碼生成；五是敏感分析，可通過擴大分析的范圍評估關鍵參數和假設變化對模型輸出的影響，從而捕捉不確定性、識別關鍵參數及優化模型；六是結果輸出直接性，PE 評價結果可直接通過R 代碼生成。此外，還具有模型的靈活性，hesim 程序包[13]為健康經濟學模擬提供了一個靈活的框架，除PSM 外，還提供馬爾科夫模型（Markov models）、狀態轉換模型（state- transition models，STM）、隊列狀態轉換模型（cohort state-transition models）和個體狀態轉換模型（individual state-transition models）。上述模型均支持對成本和效用的模擬分析，適用于多種健康狀態和時間點的數據。

PSM 作為決策支持工具，在本研究中展現了其優勢與局限性。通過GetDataW 軟件提取KN189 試驗的PFS 和OS 數據，為生存分析提供了基礎，避免了獲取臨床試驗IPD 的需求。然而，PSM 的應用依賴于臨床試驗文獻對生存曲線的詳細報告，若需分析特定亞組或多狀態患者的時間相關數據，而文獻未提供必要信息時，分析可能受限。此外，PSM 能準確模擬疾病事件，但在隨訪時間超出臨床試驗報告時，需通過參數法外推PFS 和OS 曲線，這一過程對于晚期和轉移性腫瘤藥物的成本效用分析至關重要，但需要注意PSM 方法本身限制了敏感性分析的范圍，當長期的總生存期外推存在不確定性時，局限性就會顯現。因此，建議將PSM 與STM 一起使用，以支持對其推斷的評估[4]。研究者需仔細檢查外推擬合結果，利用個體臨床事件的試驗數據、外部數據和專家意見綜合評估PSM 和STM 的生存預測可信度，避免出現不合理的統計現象[13]，后續進一步通過實例進行研究和說明。采用R 語言進行腫瘤治療藥物的成本效用分析要求研究者具備臨床知識、數理統計能力及編程經驗，研究者需在臨床方案設計理解、參數設定、模型假設及驗證等方面謹慎設定與檢查，以確保分析結果的真實性和可靠性。

本研究聚焦于使用hesim 建立模型及PSM 輔助決策的視圖功能，未能詳盡闡述R 軟件及相關程序包安裝、PFS 和OS 曲線的數字化提取、參數分布的擬合與評估、生存分析的完整性及費用模型的數據來源和計算方法等內容，感興趣的讀者可參考文獻[10-17]，以獲取更多關于上述內容的詳細信息。本研究為了突出PSM 模型建立，采用簡單數值演示費用及效用模型，若需進行復雜的費用及效用模型分析，如疾病進展后多線治療方案的比例、不良反應發生率及治療費用等，則需單獨建立模型，然后并入第5 步的實例化PSM，并進行概率敏感性分析。