劉家一，李 偉

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

2.北京醫(yī)院麻醉科（北京 100005）

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

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

1 臨床試驗(yàn)數(shù)據(jù)提取

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

圖1 KN189臨床試驗(yàn)OS曲線數(shù)據(jù)Figure 1.OS curve data of KN189

圖2 KN189臨床試驗(yàn)PFS曲線數(shù)據(jù)Figure 2.PFS curve data of KN189

2 重構(gòu)IPD

2.1 R語言環(huán)境設(shè)置

打開Rstudio 2023.09.1+494 軟件，設(shè)定工作路徑，安裝并下載軟件包：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 定義數(shù)據(jù)集

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

表1 KN189試驗(yàn)comb組OS曲線數(shù)據(jù)Table 1.OS curve data of the comb group in KN189

表2 KN189試驗(yàn)cemo組OS曲線風(fēng)險(xiǎn)人數(shù)數(shù)據(jù)Table 2.OS curve data of the cemo group in KN189

2.3 預(yù)處理與重構(gòu)IPD數(shù)據(jù)

2.3.1 預(yù)處理IPD對象格式

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

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

totalpts 為指定試驗(yàn)的初始人數(shù)；

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

206 為cemo 組試驗(yàn)的初始人數(shù)。

2.3.2 重構(gòu)IPD數(shù)據(jù)

使用getIPD 函數(shù)重建comb 組OS 曲線IPD數(shù)據(jù)，運(yùn)行代碼如下：

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

summary (est_comb)

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

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

同理獲得cemo 組OS 曲線IPD 數(shù)據(jù)，armID定義為2。

2.3.3 對比重構(gòu)IPD數(shù)據(jù)與KN189的HR值

重構(gòu)IPD 數(shù)據(jù)與KN189 報(bào)告的風(fēng)險(xiǎn)比（hazard ratio，HR）值對比分析的R 語言運(yùn)行代碼如下：

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

提取comb、cemo 組的重構(gòu)IPD 為新數(shù)據(jù)集；

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

計(jì)算重構(gòu)IPD數(shù)據(jù)中comb、cemo組的HR值；

exp (coef (hr_est_IPD)

運(yùn)行并提取HR 值為0.57，與KN189 報(bào)告中數(shù)值0.56 相近。

2.4 繪制重構(gòu)IPD數(shù)據(jù)與原read_in曲線對比檢驗(yàn)

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

圖3 cemo組重構(gòu)IPD數(shù)據(jù)OS曲線與read-in KN189數(shù)據(jù)OS曲線對比Figure 3.Comparison of OS curves of reconstructed IPD data and read-in KN189 data in cemo group

plot(est_comb)

Estimated 為重構(gòu)IPD 數(shù)據(jù)，Reported 為read_in KN189 臨床數(shù)據(jù)；

plot(est_cemo)

如圖3 所示，圖3-A 為基于重構(gòu)IPD 數(shù)據(jù)的OS 曲線與數(shù)字化讀取的KN189 臨床試驗(yàn)的OS 曲線的視覺對比圖。

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

表3 OS曲線重構(gòu)IPD數(shù)據(jù)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)

重復(fù)以上2.2 至2.4 步驟，得到OS 曲線的重構(gòu)數(shù)據(jù)框IPD_comb_cemo_os。

3 生存曲線參數(shù)分布擬合與評估

定義待擬合模型名稱的向量，同時(shí)考慮指數(shù)、威布爾、伽瑪、對數(shù)正態(tài)、對數(shù)邏輯及廣義伽瑪多個(gè)模型，運(yùn)行代碼如下：

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

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

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

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

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

KN189case_mle$model.fitting

顯示批量擬合的結(jié)果AIC、BIC、DIC 值。

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

4 構(gòu)建生存、效用、費(fèi)用模型

4.1 定義模型

調(diào)用R 軟件hesim 包的hesim_data 函數(shù)[13]，用以標(biāo)準(zhǔn)化后續(xù)生存效用模型的數(shù)據(jù)準(zhǔn)備，代碼如下：

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＂)))

分別定義分區(qū)生存模型中治療方案strategies、患者類型patients、試驗(yàn)疾病狀態(tài)states，死亡狀態(tài)不用列出。

labs＜- get_label (hesim_dat)

4.2 建立生存模型

將IPD_comb_cemo_os、IPD_comb_cemo_pfs數(shù)據(jù)合并，并整理為數(shù)據(jù)框surv_est_data，整理后的數(shù)據(jù)結(jié)構(gòu)見表4。

表4 comb與cemo組OS及PFS曲線的重構(gòu)患者數(shù)據(jù)Table 4.Reconstructed patient data of OS and PFS curves in comb group and cemo group

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

flex_fit_os＜- flexsurvreg (Surv (os_time, os_status) ～strategy_name, data=surv_est_data,dist=“l(fā)logis”)

選擇步驟3 中pfs 及os 曲線最優(yōu)擬合模型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 個(gè)狀態(tài)分別的生存擬合結(jié)果對象。

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

n_samples＜- 300

設(shè)定生存模型中將要模擬的樣本數(shù)量；

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

從hesim_dat 中指定參數(shù)進(jìn)行模擬；

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

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

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

4.3 建立效用模型

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

表5 不同狀態(tài)效用值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 建立費(fèi)用模型

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

表6 不同治療方案費(fèi)用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 實(shí)例化PSM并計(jì)算質(zhì)量調(diào)整生命年

5.1 定義模型

生存模型survmods 與后續(xù)的效用模型、費(fèi)用模型匯總進(jìn)行PSM 的實(shí)例化，并進(jìn)行時(shí)間外推，分析代碼如下：

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

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

定義生存模型實(shí)例化外推時(shí)間的范圍為20年，與前數(shù)據(jù)統(tǒng)一以周為單位，將時(shí)間范圍帶入psm 進(jìn)行生存預(yù)測；

psm$sim_survival (t=times)

將定義的外推時(shí)間輸入psm 進(jìn)行生存預(yù)測。

5.2 繪制外推生存曲線圖

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

圖4 外推時(shí)間內(nèi)的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()

整理狀態(tài)概率數(shù)據(jù)框，用于后續(xù)視圖表現(xiàn)；

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

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

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

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

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

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

圖5 外推時(shí)間內(nèi)的狀態(tài)曲線Figure 5.State curves over time

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

指定數(shù)據(jù)來源stateprobs 中strategy_id 列數(shù)據(jù)繪圖；

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

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

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

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

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

設(shè)定x、y 軸的文字標(biāo)識；

scale_fill_manual (name=＂Health state＂

設(shè)置圖例的標(biāo)題；

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

設(shè)置填充圖形的顏色映射；

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

標(biāo)識每個(gè)狀態(tài)的文字內(nèi)容；

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

設(shè)置x 軸的刻度標(biāo)簽；

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

調(diào)整圖例以逆序的方式且排2 行；

theme (legend.position=＂bottom＂)

設(shè)置圖例的位置在底部。

5.3 模擬成本與質(zhì)量調(diào)整生命年

基于5.2 中的狀態(tài)概率進(jìn)行數(shù)值積分模擬貼現(xiàn)成本與質(zhì)量調(diào)整生命年（quality-adjusted life years，QALYs），選擇5%進(jìn)行貼現(xiàn)，分析代碼如下：

psm$sim_costs (dr=0.05)

運(yùn)行后結(jié)果見圖6。

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

psm$sims_qalys (dr=0.05)

運(yùn)行后結(jié)果見圖7。

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

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

6.1 創(chuàng)建模擬數(shù)據(jù)

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

ce_sim＜- psm$summarize()

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

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

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

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

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

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

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

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

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 組數(shù)據(jù)，將兩組數(shù)據(jù)合并后，同時(shí)繪制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 可知，當(dāng)支付意愿值超過4 800 元時(shí)，comb 組具有成本效益的可能性開始大于cemo 組。

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

6.3 繪制增量成本-效益比值散點(diǎn)圖

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

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

提取cea_pw_out 數(shù)據(jù)集中comb 組的數(shù)據(jù)；

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

計(jì)算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 可知，當(dāng)意愿支付閾值為14 312 時(shí)，comb 組干預(yù)方案具有成本-效益的可能性為97.5%；當(dāng)意愿支付閾值為3 123 時(shí)，comb 組干預(yù)方案具有成本-效益的可能性僅為2.5%。

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

7 討論

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

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

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