讀者朋友你們好啊分苇，如果你關(guān)注過之前幾期公眾號(hào)（全哥的學(xué)習(xí)生涯）的內(nèi)容医寿，你會(huì)知道我著重分享了如何利用R作meta分析的方法掏颊。而由于我本人研究方向的問題乌叶，今天這一期的內(nèi)容將以文末參考文獻(xiàn)為例准浴，向大家分享展示如何利用heemod程序包構(gòu)建Markov模型進(jìn)行成本效果分析乐横，以期為衛(wèi)生經(jīng)濟(jì)評(píng)價(jià)的讀者伙伴們提供參考與幫助。需要本文完整代碼的請(qǐng)?jiān)诠娞?hào)（全哥的學(xué)習(xí)生涯）內(nèi)留言“衛(wèi)生經(jīng)濟(jì)學(xué)評(píng)價(jià)”獲取罐农。
heemod 程序包主要開發(fā)用于衛(wèi)生經(jīng)濟(jì)學(xué)評(píng)價(jià)宰睡，包括Markov 模型的構(gòu)建、運(yùn)行和敏感性分析等气筋；與Excel 和TreeAge 等衛(wèi)生經(jīng)濟(jì)學(xué)評(píng)價(jià)的常用軟件相比拆内，基于R 進(jìn)行衛(wèi)生經(jīng)濟(jì)學(xué)研究具有開源免費(fèi)、靈活性高和圖形功能強(qiáng)大等優(yōu)勢(shì)宠默。其效果圖展示如下：
1. heemod 包的安裝與加載
代碼如下：
install.packages(“heemod”)
library(“heemod”)
2. Markov 模型的構(gòu)建
2.1 確定評(píng)價(jià)方案
根據(jù)實(shí)際情況明確需要進(jìn)行比較的治療方案（兩種或兩種以上）麸恍。本文進(jìn)行比較的是兩種晚期乳腺癌（advanced breast cancer，ABC）的治療方案搀矫，拉帕替尼聯(lián)合卡培他濱方案和卡培他濱單藥方案抹沪。
2.2 Markov狀態(tài)的設(shè)定
Markov狀態(tài)主要依據(jù)具體的疾病和比較的方案進(jìn)行設(shè)定，在一個(gè)周期中瓤球，模型中的患者只能處于唯一的健康狀態(tài)中权她，腫瘤疾病中常用的是含有疾病穩(wěn)定蝴罪、疾病進(jìn)展和死亡三種狀態(tài)的Markov模型，也可以根據(jù)實(shí)際情況構(gòu)建含有更多健康狀態(tài)（state）的Markov模型，包含健康狀態(tài)越多，模型越復(fù)雜疮装。本例假設(shè)模型含有疾病穩(wěn)定刷袍、治療響應(yīng)滚局、疾病進(jìn)展和死亡4 種健康狀態(tài)。
2.3 轉(zhuǎn)移概率的設(shè)定
靜態(tài)Markov模型中省骂，各個(gè)健康狀態(tài)間的轉(zhuǎn)移概率不隨時(shí)間改變而變化涨缚，此種情況下可直接定義各狀態(tài)下的轉(zhuǎn)移概率并賦值。本文研究中假設(shè)轉(zhuǎn)移概率不隨時(shí)間改變，各狀態(tài)間的轉(zhuǎn)移概率源于關(guān)于拉帕替尼的兩項(xiàng)臨床試驗(yàn)惭嚣。轉(zhuǎn)移概率的定義代碼如下：
mod_par <- define_parameters(
pSTR_comb = 0.0620,
pSTS_comb = 0.7440,
pSTP_comb = 0.1400,
pSTD_comb = 0.0540,
pSTR_mono = 0.0430,
pSTS_mono = 0.6800,
pSTP_mono = 0.2150,
pSTD_mono =0.0620,
pRTR_comb = 0.8150,
pRTP_comb = 0.1450,
pRTD_comb = 0.0400,
pRTR_mono = 0.8150,
pRTP_mono = 0.1050,
pRTD_mono = 0.0800,
pPTP_comb = 0.9070,
pPTD_comb = 0.0930,
pPTP_mono = 0.9040,
pPTD_mono = 0.0960)
上述代碼中以S表示疾病穩(wěn)定谋国，P表示疾病進(jìn)展弯洗，R表示治療響應(yīng)，D表示死亡狀態(tài)；以comb表示拉帕替尼和卡培他濱聯(lián)合給藥方案，mono表示卡培他濱單藥方案沪摄；pSTP_comb表示聯(lián)合給藥組從疾病穩(wěn)定狀態(tài)轉(zhuǎn)移至疾病進(jìn)展?fàn)顟B(tài)的概率哄陶，余則依次類推。
2.4 成本和健康效用值的定義
成本可根據(jù)實(shí)際情況進(jìn)行區(qū)分，如可分為藥物治療成本、住院成本等。在本例中，250 mg拉帕替尼的價(jià)格為23美元（Clapa），卡培他濱價(jià)格為1.5美元/150 mg(uS)遗座，分別定義了疾病穩(wěn)定（0.70）、治療響應(yīng)（0.84）和疾病進(jìn)展（0.50）三種狀態(tài)下的效用值。此外，因?yàn)榭ㄅ嗨麨I按照體表面積給藥，所以也設(shè)定了模擬患者的基線體表面積（body surface area代箭，BSA）杜漠，具體代碼如下：
mod_par <- modify(mod_par,cLapa = 23,cCape = 1.5,uS = 0.70,uR = 0.84,uP = 0.50,BSA = 1.72)
2.5 轉(zhuǎn)移概率矩陣
動(dòng)態(tài)Markov模型中，各狀態(tài)間的轉(zhuǎn)移概率具有時(shí)間依賴性，通常需要進(jìn)行參數(shù)生存分析以便于獲得不同時(shí)間點(diǎn)的轉(zhuǎn)移概率矩陣击碗。當(dāng)可以獲取原始生存資料時(shí)稍途，參數(shù)定義完畢后械拍，即可通過define_transition()函數(shù)定義轉(zhuǎn)移概率矩陣，具體代碼如下：
mat_comb <- define_transition(state_names = c("S", "R","P", "D"),pSTS_comb,pSTR_comb, pSTP_comb, pSTD_comb,0, pRTR_comb,pRTP_comb,pRTD_comb,0,0,pPTP_comb,pPTD_comb,0, 0, 0, 1)
mat_mono <- define_transition(state_names = c("S", "R","P", "D",pSTS_mono, pSTR_mono, pSTP_mono, pSTD_mono,0, pRTR_mono,pRTP_mono, pRTD_mono,0, 0,pPTP_mono,pPTD_mono,0, 0, 0, 1)
每一種治療方案均有各自的轉(zhuǎn)移概率矩陣坷虑，此處mat_comb代表聯(lián)合給藥方案的轉(zhuǎn)移概率矩陣猖吴，mat_mono代表單藥方案的轉(zhuǎn)移概率矩陣摔刁，通過plot()函數(shù)可繪制狀態(tài)轉(zhuǎn)移圖（State Transition Diagram），命令如下：
plot(mat_comb)
plot(mat_mono)
結(jié)果如圖1海蔽，圖2所示党窜。
2.6 定義狀態(tài)值
定義狀態(tài)值在define_state()函數(shù)中進(jìn)行，主要定義各健康狀態(tài)下的成本和收益（效用值），如定義死亡狀態(tài)的效用值為0径簿，本例模型含有4種健康狀態(tài)，因此需要分別對(duì)4種狀態(tài)下的成本和收益進(jìn)行定義因惭。具體代碼如下：
S<- define_state(cost = dispatch_strategy(comb = 5*cLapa*6*7 + BSA*2000*cCape/150*4*7,mono = BSA*2500*cCape/150*4*7),qaly = uS)
R<- define_state(cost = dispatch_strategy(comb = 5*cLapa*6*7 + BSA*2000*cCape/150*4*7,mono = BSA*2500*cCape/150*4*7),qaly = uR)
P<- define_state(cost = dispatch_strategy(comb = 3631,mono = 3823),qaly = uP)
D<- define_state(cost = 0,qaly = 0)
以上代碼中岳锁，歡迎關(guān)注公眾號(hào)：全哥的學(xué)習(xí)生涯，獲取更多數(shù)據(jù)分析知識(shí)蹦魔。死亡狀態(tài)的成本和效用值設(shè)置為0激率。dispatch_strategy()函數(shù)區(qū)分同一狀態(tài)下不同治療方案的成本或效果。各個(gè)健康狀態(tài)下的成本是按照Markov周期長度計(jì)算的勿决，如原始研究根據(jù)臨床試驗(yàn)的隨訪時(shí)間將周期長度設(shè)為1.5月乒躺，計(jì)算的成本即為1.5月內(nèi)經(jīng)不同方案治療的患者在各狀態(tài)下的成本。
2.7 定義治療方案
完成以上的步驟以后低缩，即可對(duì)各治療方案進(jìn)行定義嘉冒，本例中進(jìn)行比較的方案為拉帕替尼聯(lián)合卡培他濱方案（strat_comb）和卡培他濱單藥方案（strat_mono）兩種，代碼如下：
strat_comb<- define_strategy(transition = mat_comb,S =S,R = R,P = P,D = D)
strat_mono<- define_strategy(transition = mat_mono,S = S,R = R,P = P,D =D)
3 運(yùn)行模型與結(jié)果展示
3.1 模型的運(yùn)行
完成參數(shù)咆繁、轉(zhuǎn)移概率矩陣讳推、狀態(tài)值和方案的設(shè)定后，Markov模型即構(gòu)建完成玩般，接下來就是運(yùn)行模型獲得分析結(jié)果银觅，模型運(yùn)行命令如下：
res_mod<-run_model(parameters = mod_par,comb = strat_comb,mono = strat_mono,cycles = 20,cost = cost,effect = qaly)
在run_model()函數(shù)中，可通過指定cycles參數(shù)確定循環(huán)的周期數(shù)（cycles）
3.2 模型結(jié)果的展示
模型運(yùn)行完成后壤短，可通過summary()函數(shù)查看模型運(yùn)行結(jié)果设拟，包括成本、增量成本久脯、效果纳胧、增量效果、增量成本-效果比（incremental cost-effectiveness ratio帘撰，ICER）和成本效果邊界等跑慕，結(jié)果如圖3所示。
使用plot()函數(shù)摧找，并指定type參數(shù)核行，即可生成隊(duì)列人數(shù)轉(zhuǎn)移圖和成本效果圖，命令如下：
plot(res_mod)
plot(res_mod, type = "ce")
結(jié)果分別如圖4和圖5所示.
圖4顯示了各健康狀態(tài)的人數(shù)隨模型周期進(jìn)行的變化情況蹬耘，模型初期芝雪，所有患者（1000人）均處于疾病穩(wěn)定狀態(tài)。圖5中默認(rèn)以成本最低的方案（單藥方案）為參照方案综苔，連線為成本效果邊界惩系，連接兩方案直線的斜率即為ICER位岔。
4.敏感性分析
4.1 單因素敏感性分析
因模型參數(shù)常存在不確定性，歡迎關(guān)注公眾號(hào)：全哥的學(xué)習(xí)生涯堡牡，獲取更多數(shù)據(jù)分析知識(shí)抒抬。所以通過敏感性分析來考察不確定因素對(duì)結(jié)論的影響，敏感性分析包括單因素敏感性分析和概率敏感性分析晤柄。在單因素敏感性分析中擦剑，單一變量的取值在一定范圍內(nèi)（如95%置信區(qū)間或在基線值上下浮動(dòng)10%等）變化，其他的變量固定為基線值芥颈，以觀察單一變量值對(duì)模型結(jié)論的影響惠勒，結(jié)果通常以龍卷風(fēng)圖（tornado diagram）的形式展示，在函數(shù)define_parameters()中直接定義的參數(shù)均可以進(jìn)行敏感性分析浇借，以下為單因素敏感性分析和繪制龍卷風(fēng)圖得代碼捉撮，結(jié)果如圖6所示。
define_dsa<- define_dsa(cLapa, 18.4, 27.6,cCape, 1.2, 1.8,BSA, 1.5, 1.9, uS, 0.63, 0.77,uR, 0.765, 0.924,uP, 0.45, 0.55)
result_dsa<- run_dsa(res_mod, dsa = def_dsa)
plot(result_dsa, result = "icer", type = "difference")
其中妇垢，cLapa為25mg拉帕替尼的成本巾遭，uS為疾病穩(wěn)定狀態(tài)健康效用值，uR為治療響應(yīng)狀態(tài)健康效用值闯估，uP為疾病進(jìn)展?fàn)顟B(tài)健康效用值灼舍，cCape為150mg卡培他濱成本,BSA為體表面積。
從龍卷風(fēng)圖中可見涨薪，拉帕替尼的成本對(duì)ICER的影響較大骑素，其次為疾病穩(wěn)定狀態(tài)、治療響應(yīng)狀態(tài)和疾病進(jìn)展?fàn)顟B(tài)的效用值刚夺，卡培他濱的價(jià)格和BSA對(duì)結(jié)果的影響較小献丑。
4.2 概率敏感性分析
概率敏感性分析中假定各變量均服從某一分布，如正態(tài)分布侠姑、對(duì)數(shù)正態(tài)分布创橄、gamma分布和beta分布等等，通過從概率分布中隨機(jī)取樣莽红，重復(fù)多次妥畏，觀察不確定性因素對(duì)模型的影響情況，概率敏感性分析和相關(guān)圖形繪制的代碼如下：
define_psa<- define_psa(cLapa ~gamma(mean = 23, sd = 5),
cCape~gamma(mean = 1.5, sd = 0.5),
BSA ~normal(mean = 1.75, sd = 0.8),
uS ~normal(mean = 0.7, sd = 0.07),
uR ~normal(mean = 0.84, sd = 0.08),
uP ~normal(mean = 0.5, sd = 0.05))
result_psa<- run_psa(res_mod, psa = def_psa, N = 1000)
plot(result_psa, type = "ce")
plot(result_psa, type = "ac")
可通過修改參數(shù)N調(diào)整模擬次數(shù)安吁，默認(rèn)為1000次醉蚁，當(dāng)指定type參數(shù)為“ce”時(shí)，繪制增量效果為橫坐標(biāo)鬼店，增量成本為縱坐標(biāo)繪制散點(diǎn)圖（圖7）网棍，當(dāng)type=“ac”時(shí)，繪制成本效果可接受曲線（圖8）妇智。
最后滥玷，如果屏幕前的你對(duì)R語言學(xué)習(xí)還有什么問題和看法捌锭，或者有什么建議，歡迎在公眾號(hào)（全哥的學(xué)習(xí)生涯）內(nèi)給我留言罗捎，或者直接添加我的個(gè)人微信（公眾號(hào)內(nèi)菜單欄”與我聯(lián)系—聯(lián)系方式“可獲得）
文中數(shù)據(jù)來源：
LE Q A, HAY J W. Cost-effectiveness analysis of lapatinib in HER-2-positive advanced breast cancer[J]. Cancer,2010,115(3):489-498.
