對生存數據進行網狀Meta分析時往往依賴恒定比例風險假設，但當生存曲線存在相交時，這一假設便不成立。隨著免疫治療等創新藥物和療法出現，基于非恒定比例風險的網狀Meta分析在當前腫瘤藥物循證評價的重要性日益突出。分數多項式模型不依賴恒定比例風險假設，能靈活捕捉生存曲線特征，擬合效果較優于恒定比例風險模型。本研究展示了采用非恒定比例風險的分數多項式模型進行網狀Meta分析在R語言中的完整操作流程。
引用本文： 趙明燁, 邵泰航, 唐文熙. 基于分數多項式的非恒定比例風險網狀Meta分析：R語言實例研究. 中國循證醫學雜志, 2022, 22(7): 853-861. doi: 10.7507/1672-2531.202202001 復制
在循證醫學中，頭對頭的隨機對照試驗（randomized clinical trial，RCT）是評價藥物有效性和安全性的金標準。但實際上，不同創新療法之間通常缺乏頭對頭的RCT，這種情況下通常需要對現有的臨床證據進行間接比較（網狀Meta分析，network meta-analysis，NMA）[1]。目前絕大多數文獻進行生存數據NMA是利用治療方案與標準對照間的風險比值（hazards ratios，HR），基于恒定比例風險（proportional hazards，PH）進行分析[2-3]。然而，當干預措施間的生存曲線及其風險函數存在交叉時，PH假設便不再成立。近年來，免疫治療等創新療法的出現極大地改善了患者的生存獲益，相比傳統的化療，長期生存獲益更加突出，此時PH假設難以成立[4-6]。若此時在干預措施評價中強行假設PH，生存數據擬合會出現嚴重偏差。因此，當PH假設不成立時，迫切需要在NMA中引入非恒定比例風險假設（non-proportional hazards，Non-PH）和相應計算模型。
近年來，越來越多基于Non-PH假設的方法出現。Ouwens等[7]開發出不依賴PH假定的參數化生存曲線模型，該方法主要是將形狀參數和尺度參數進行再參數化的基礎上對風險值建模。Lu等[8]提出分段指數模型，假定生存曲線在不同階段服從PH假設，通過觀察對數累計風險圖來確定分隔時點。Petit等[9]應用限制性平均生存時間（restricted mean survival time，RMST）作為HR的代替指標進行NMA，將不同治療方案的曲線下面積之差作為最終效應量。Jansen等[10]提出分數多項式（fractional polynomials，FP）模型，可對風險值靈活建模，不依賴PH假設。以上4種方法都有一定的應用，且出現在了NICE指南中[11]。張天嵩等[12]提出，在有足夠的KM（Kaplan-Meier）曲線證據時，推薦采用FP方法進行建模。然而，國內還沒有相關文獻對該方法具體介紹。因此，本研究以治療非小細胞肺癌的4種藥物生存數據為例，納入7個研究，在廣義線性模型（generalized Linear model，GLM）的框架下，展示如何基于FP模型進行Non-PH的NMA，基本流程見圖1。
圖1
基于FP模型進行Non-PH的NMA流程圖
1 模型原理簡介
FP模型包括一階FP與二階FP，可表達成式1。其中，p1和p2可在（?2，?1，?0.5，0，0.5，1，2，3）中任取一值。因此，總計8種一階FP模型，36種二階FP模型。在一階FP模型中，當p=0時，FP模型與weibull模型一致；當p=1時，此時FP模型等同于Gompertz模型[13]。
 |
研究表明分組生存數據可采用二項似然分布［rjkt～Bin（pjkt，njkt）］，rjkt表示發生事件人數，njkt表示危險人數，pjkt表示t時刻事件發生率）及相應地loglog連接函數建模[14]。一般形式為：
 |
式2中，
jkt表示t時刻j研究的k治療方案的線性預測方程，ln（
tjkt）表示時間間隔造成的偏移量[12]。連續時間線性風險方程通過如下方式轉化：
 |
上式最后一步中，假設風險值在［
時段內恒定。在
趨近于0時，這一假設是成立的。此公式稍作變形便可得到Jansen研究中的風險方程[13]：
 |
式4中，對方程的左右兩邊取對數值，即可得到：
 |
由此可見，對于分組數據而言，任何線性對數風險方程都可以在GLM中表達。值得注意的是，“線性”并不意味著模型的風險函數隨時間呈線性變化，而是可依賴于任一變量，因此風險函數的形狀是多變的。一般地，在固定效應分散多項式的GLM模型中，
jkt可以寫成以下形式[14]：
 |
式6中，
為j研究的特異效應系數，
為治療方法k的特異效應系數。M為FP的階數，一般而言，M取1或2均可。假設A方案在j試驗中，B方案在w試驗中，t時刻A相對于B的HR可表示為式7。
 |
此外，GLM模型也可采用隨機效應模型進行擬合，原理與固定效應模型相似，只是在固定效應模型的基礎上加入了隨機效應量。一般可表達成式8：
 |
2 軟件安裝與程序加載
2.1 R語言安裝
截止筆者完稿時，R語言最新版本為4.1.2，讀者可在以下網站免費安裝：https://www.r-project.org/。
2.2 程序包加載
R語言中，通過install.packages函數可實現程序包的安裝，library函數可加載相關軟件包。在本研究中，需要的程序包包括：dplyr、ggplot2。通過以下代碼可實現程序的安裝及運行：
install.packages("dplyr")
install.packages("ggplot2")
3 數據預處理與加載
3.1 示例數據的選擇
本研究采用與Jansen[10]研究中一致的數據。共納入7個研究[15-21]、涉及4種治療非小細胞肺癌的治療方法，包括最佳支持治療（best support care，BSC）、多西他賽、吉非替尼和培美曲塞。
3.2 重構個體數據
缺失個體數據時，Guyot[22]的方法是目前最常用的重構方法，其根據KM曲線相關信息即可較好地重構患者個體數據。在R語言中可通過survHE程序包實現，已有文獻詳細講解了該方法[23]。由于篇幅限制，本文不再贅述。
3.3 數據格式整理
為方便在R語言中實現模型的運算，數據需整理成一定格式，見表1。由于數據過多，本研究只展示各項治療方案前6個月的生存數據。表1中，studyn表示研究方案的編號，“1～7”分別表示：“Chang2006”、“Cufer2006”、“Hanna2004”、“Kim2008”、“Lee2010”、“Maruyama2008”及“Shepherd2000”；trtn表示治療方法的編號，“1～4”分別表示：BSC、多西他賽，吉非替尼和培美曲塞；time為研究時限（月）；timeDelta為時間間隔（月），nevents和natrisk分別表示對應時間的發生事件數和危險人數。
表1
示例試驗數據整理
3.4 數據讀取與錄入
將表1保存為NMA.csv文件后，通過“read.csv”命令讀取該表格，使用“data.frame”生成數據矩陣，并命名為“km”；“ref”命令用來規定對照研究及治療方法。“studies”、“treatments”結合factor函數重命名各研究和治療方法。（注：代碼部分中“#”后內容是代碼釋義）
data<- read.csv(C:/Users/lenovo/Desktop/NMA.csv)#括號中為NMA.csv文件的存儲路徑
km <- data.frame(studyn = data$studyn,trtn = data$trtn,time = data$time.c,timeDelta = data$timeDelta,nevents = data$nevents,natrisk = data$narisk)
ref.study <- "Kim2008"
ref.trt <- "Docetaxel"
studies <- c("Chang2006", "Cufer2006", "Hanna2004", "Kim2008", "Lee2010", "Maruyama2008", "Shepherd2000")
treatments <- c("BSC", "Docetaxel", "Gefitinib", "Pemetrexed")
km$studyf <- factor(km$studyn, labels = studies)
km$trtf <- factor(km$trtn, labels = treatments)
km$studyf <- relevel(km$studyf, ref=ref.study)#將對照研究放在最前列
km$trtf <- relevel(km$trtf, ref=ref.trt)
km$trtn=as.numeric(km$trtf)#將因子型變量轉化為數值型
km$studyn=as.numeric(km$studyf)
4 模型構建
4.1 模型結構設置
首先，對全部的FP模型進行擬合。因此，需要生成一個涵蓋所有模型的列表，通過list函數來實現。function函數構造相應模型，具體代碼如下：
models <- list(
"First order FP,p1=0" = list(g1=function(x){log(x)},g2=function(x){0},f1=function(x){log(x)},f2=function(x){0}),
"First order FP,p1=1" = list(g1=function(x){x},g2=function(x){0},f1=function(x){x},f2=function(x){0}),
"Second order FP,p1=-2, p2=1" = list(g1=function(x){x^-2},g2=function(x){x},f1=function(x){x^-2},f2=function(x){x}),
"Second order FP,p1=-1, p2=2" = list(g1=function(x){x^-1},g2=function(x){x^2},f1=function(x){x^-1},f2=function(x){x^2}))
由于篇幅所限，本研究只展示部分模型代碼。讀者可通過改變功能函數中p1和p2自行添加。
4.2 模型擬合
將擬合好的數據保存于“km.new”數據集中，并引入FP模型相關參數。由function構造擬合所需要的功能函數，g1、g2、f1及f2已在模型篩選部分設定，g0和f0為常數項對應系數，默認為1。cblind函數用于構造包含“nevents”及“natrisk-nevents”的矩陣，在GLM框架下擬合相應模型，最后通過lapply函數將所有模型在設置好的結構（“fit.KM.NMA”）下進行擬合。
fit.KM.NMA<-function(bf){
km.new=km
km.new$g0=1
km.new$f0=1
km.new$g1=bf[[1]](km.new$time)
km.new$g2=bf[[2]](km.new$time)
km.new$f1=bf[[3]](km.new$time)
km.new$f2=bf[[4]](km.new$time)
f=cbind(nevents,natrisk-nevents)~trtf*f0+studyf*g0+trtf*f1+trtf*f2+studyf*g1+studyf*g2
glm(f,family=binomial(link=cloglog),data=km.new,offset = log(timeDelta)}#模型選擇二項分布，cloglog函數作為連接函數，時間間隔的對數作為偏移量
fits=lapply(models,fit.KM.NMA)
4.3 擬合度檢驗
在模型擬合過程中，R語言會給出相應的AIC值，使用lapply函數可生成全部擬合模型AIC值。data.frame函數進一步將AIC值整理成可視矩陣。
aics=lapply(fits, AIC)
data.frame(AIC=round(unlist(aics), 2))
根據赤池原則，AIC值越小表明模型擬合效果越佳。結果表明，二階FP（p1=?2，p2=1）總體擬合效果最好。此外，一階FP（p=?2）在所有一階FP模型中擬合最佳。FP模型擬合效果優于標準參數法（p=0或1）。
然而，NICE指南提出AIC最小并不意味外推性能佳。例如，在對成熟度不足的生存數據進行擬合時，AIC值小的二階FP模型可能存在對尾部數據過度擬合的問題，最終的外推效果尚不如一階FP模型。因此，建議考慮先采用多個模型，并通過觀察模型對應的生存曲線確定合適模型。本研究納入一階和二階FP模型各三個，見表2。
表2
納入模型基本信息
4.4 生存數據外推
通過function函數構造模型的結構，pred函數用于模型的預測。km.pred用于保存模型預測的數據及模型結構。具體代碼如下：
pred.KM.NMA<-function(bf){
trts=data.frame(trtf=unique(km$trtf))#剔除重復的治療方案
trts$studyf=sort(unique(km$studyf))[1]#確定對照研究
timehorizon=data.frame(time=2*(1:30))#確定模擬周期
km.pred=merge.data.frame(timehorizon,trts)#在數據矩陣中引入“timehorizon”及“trts”
km.pred$g0=1
km.pred$f0=1
km.pred$g1=bf[[1]](km.pred$time)
km.pred$g2=bf[[2]](km.pred$time)
km.pred$f1=bf[[3]](km.pred$time)
km.pred$f2=bf[[4]](km.pred$time)#生成模型結構
km.pred$timeDelta<-2#設置時間間隔
km.pred}
pred.KM.data=lapply(models,pred.KM.NMA)#批量計算全部模型預測結果
4.5 HR計算
由“predict.glm”命令生成風險預測值“pred”，“%>%”為管道函數，可將左側數據或表達式傳遞到右側。通過filter函數可以篩選出對照治療方案，進而計算其風險函數值“pred.doce”。最后由“mutate(lnhr = pred-pred.doce)”命令計算出各時點全部治療方案與對照方案相比的lnHR。具體代碼如下：
for(i in 1:length(models)){
pred.KM.data[[i]]$pred=predict.glm(fits[[i]], pred.KM.data[[i]])
d=pred.KM.data[[i]]
d1<-d %>% #d賦值于d1
filter(trtf == ref.trt)%>%#篩選出對照治療方案
mutate(pred.doce = pred)%>% #生成對照方案的風險值“pred.doce”作為新的一列
select(pred.doce, time) %>% #選擇“pred.doce”、“time”兩列
arrange(time) %>% #按時間排序“pred.doce”
left_join(d) %>%#生成數據從左邊加入d數據集
mutate(lnhr = pred-pred.doce) %>% #生成lnhr
mutate(hr=exp(lnhr))#%>%#將lnhr轉化為hr
d1$modelc=names(models)[[i]]
if(i==1)dpred<-d1#對照方案的“dpred”賦值于“d1”
if(i!=1)dpred<-rbind(dpred,d1)}#非對照方案，將“dpred”和“d1”合并輸出
dpred$Model<-factor(dpred$modelc,levels=names(models))#dpred中生成新的一列“Model”
生成HR后，還可利用ggplot2對其進行可視化分析，具體代碼如下：
f1<-ggplot() +
geom_line(data=dpred, aes(x=time, y=hr, group=trtf, colour=trtf), size=1)+#圖形相關參數設置
geom_hline(yintercept=1, lty=2)+#添加水平線
#geom_vline(xintercept=14, lty=2)+#添加垂直線
facet_wrap(~Model, nrow=3)+#設置數據排列格式
scale_color_discrete(name="Treatment")+#采用離散標尺給圖形上色
scale_y_log10(limits = c(0.1, 10), breaks = c(0.1, 0.25, 0.5, 1, 2, 4, 10))+#定義縱坐標的刻度及范圍
scale_x_continuous(breaks = c(seq(from=0, to=60, by = 12)))+#定義橫坐標的刻度及范圍
ylab("Hazard Ratio")+#設置縱坐標軸名稱
xlab("Time(months)")+#設置橫坐標軸名稱
theme(legend.position = "bottom")+#設置圖例位置
theme_bw()#使用默認主題
ggsave(plot = f1, filename = "hr.png", width = 12, height = 9)#命名并保存hr圖
運行上述代碼后，即可得到HR的可視化結果，見圖2。
圖2
隨時間變化的HR值
4.6 生存率計算
通過加總各時點風險值（haz），可獲取累計風險值（cumhaz）。由式9，可獲取各試點全部治療方案的生存數據[14]。本部分代碼同3.4基本一致，因此僅對部分代碼進行解釋。
 |
for(i in 1:length(models)){
pred.KM.data[[i]]$pred=predict.glm(fits[[i]], pred.KM.data[[i]])
d=pred.KM.data[[i]]
d1<-d %>%
dplyr::mutate(haz = exp(pred))%>%
dplyr::group_by(trtf) %>%
dplyr::arrange(time) %>%
dplyr::mutate(cumhaz = cumsum(haz)) %>%
dplyr::mutate(survProp = exp(-1*cumhaz))
zero <- unique(d1[,c("trtf", "studyf")])
zero$time=0#添加一行，t=0
zero$survProp=1#添加0時刻的生存率
d1=rbind(zero, d1)#合并“zero”與“d1”
d1$modelc=names(models)[[i]]
if(i==1)dpred<-d1
if(i!=1)dpred<-rbind(dpred, d1)}
dpred$Model<-factor(dpred$modelc, levels=names(models))
同樣地，利用ggplot2程序包可繪制各治療方案的生存曲線，見圖3。代碼如下：（代碼相關釋義與3.4一致）
圖3
各模型對應生存曲線
f2= ggplot() +
geom_line(data=dpred, aes(x=time, y=survProp, group=trtf, colour=trtf), size=1)+
scale_color_discrete(name="Treatment")+
expand_limits(y=c(0,1), x=c(0,60))+
facet_wrap(~Model, nrow=3)+
scale_x_continuous(breaks = c(seq(from=0, to=60, by = 12)))+
ylab("Proportion surviving")+
xlab("Time(months)")+
guides(color = guide_legend(ncol = 1))+
theme(legend.position = "bottom")+
theme_bw()
ggsave(plot = f2, filename = "paperSurv.png", width = 12, height = 9)
本研究中BSC生存數據相對較短，部分FP模型對其尾部數據過度擬合。通過圖2可知，2階FP模型（p1=1，p2=2或3）對于BSC的外推結果較差，其長期生存獲益優于積極治療，與臨床實際效果不符，需排除。
4.7 最優模型的確定和回歸系數生成
根據各模型AIC值及生存曲線圖，“p1=?2，p2=1”的二階FP模型總體AIC值最小，且生存曲線擬合較好。因此本研究與Jansen等研究一致，選擇該模型作為最優模型[13-14]。采用以下代碼合成最終的模型：
pred.KM.data=pred.KM.NMA(models[["Second order FP,p1=-2, p2=1"]])#“p1=-2,p2=1”的二階FP模型預測數據作為最終的pred.KM.data
最后，通過以下代碼可提取各治療方案相比于對照方案的回歸系數，即式6中的“
，其中，m取0，1，2。將回歸系數代入式7即可算出各時點各方案相比于對照方案的HR值。見表3。
表3
最優模型參數估計
timeh<-data.frame(time=1:60,a=1)#確定模擬時限
fit<-fits[["Second order FP,p1=-2, p2=1"]]#確定最終擬合模型
pred.data<-fit$data
pred.data$est1<-predict.glm(fit)#獲取比較系數均值
pred.data$est2=NA
pred.data$se<-predict.glm(fit,se.fit=T)$se.fit#獲取比較系數標準誤
coeff.data<- as.data.frame(coef(summary(fits[["Second order FP,p1=-2, p2=1"]])))[,1:2]#選取模型擬合回歸結果的前兩列，構建數據集“coeff.data”
names(coeff.data)<- c("est","std.err")#對前兩列重命名，分別為"est","std.err"
attach(coeff.data)#綁定數據集
coeff.data$conf.int.lower<- est-std.err*qnorm(0.975)#生成CI下限
coeff.data$conf.int.upper<- est+std.err*qnorm(0.975)#生成CI上限
coeff.data[,1:4]<- round(coeff.data[,1:4],3)#數值保留三位小數
coeff.data$pn<- rownames(coeff.data)#添加治療方案名稱
coeff.data$contrast<- paste(levels(km$trtf)[-1],"vs",ref.trt)#添加對比方案名稱
detach(coeff.data)#解綁數據集
coeff.data#生成回歸系數
5 討論
本研究采用治療非小細胞肺癌的4種藥物的生存曲線構建基于FP模型的非恒定比例風險NMA，對其原理及方法進行介紹，并采用實際案例展示了如何在R語言中構建模型、擬合生存曲線。本研究分析結果表明，相比于標準曲線，如岡茲分布、威布爾分布，FP模型的擬合效果更優。通過視覺檢查，基于FP模型重構的生存曲線擬合也能較好地貼合臨床實際。
隨著PD-（L）1等創新藥物的出現，腫瘤藥物生存曲線形狀更加多樣的同時，相比于傳統的化療往往有明顯的拖尾效應，因此PH模型不再適用。Non-PH模型中，RMST模型只能截取KM時限內的生存數據進行分析，無法外推的特點也限制了其應用；而參數化生存曲線模型及分段指數模型雖然一定程度上放松了PH假設，但當生存曲線的風險率不再單調變化時，這類模型便不再適用。
FP模型的優勢包括：① 相比于傳統NMA模型中，FP模型更能靈活捕捉生存曲線的特征，允許各治療方案生存曲線有不同形狀，不再依賴于PH假設，因此基于FP的NMA擬合效果較傳統NMA好。② FP模型可擬合到生存曲線的漸近線，相比于標準多項式或樣條模型，末端數據更加穩定。③ FP模型可用于預測KM時點外的生存數據，方便測算干預措施終生投入和產出，對于衛生決策意義重大[13]。同時，FP模型有一定的不足，容易受尾部稀疏數據影響，產生過度擬合。因此需觀察多個生存曲線或風險圖以確定最合適的模型參數。
鑒于FP模型的優越性，其應用越來越多。如，Wang等在對前列腺癌進行NMA時采用了一階FP模型[24]；Vickers等[25]采用FP模型評價非小細胞肺癌二線治療方案的療效。此外，NICE相關HTA評價中也多次用到了FP模型[10]。
由于篇幅有限，本研究只在頻率學框架中展示了基于FP模型在NMA中的應用，后續研究將進一步闡述如何在貝葉斯框架下實現模型的搭建。
聲明 所有作者均聲明無利益沖突。
在循證醫學中，頭對頭的隨機對照試驗（randomized clinical trial，RCT）是評價藥物有效性和安全性的金標準。但實際上，不同創新療法之間通常缺乏頭對頭的RCT，這種情況下通常需要對現有的臨床證據進行間接比較（網狀Meta分析，network meta-analysis，NMA）[1]。目前絕大多數文獻進行生存數據NMA是利用治療方案與標準對照間的風險比值（hazards ratios，HR），基于恒定比例風險（proportional hazards，PH）進行分析[2-3]。然而，當干預措施間的生存曲線及其風險函數存在交叉時，PH假設便不再成立。近年來，免疫治療等創新療法的出現極大地改善了患者的生存獲益，相比傳統的化療，長期生存獲益更加突出，此時PH假設難以成立[4-6]。若此時在干預措施評價中強行假設PH，生存數據擬合會出現嚴重偏差。因此，當PH假設不成立時，迫切需要在NMA中引入非恒定比例風險假設（non-proportional hazards，Non-PH）和相應計算模型。
近年來，越來越多基于Non-PH假設的方法出現。Ouwens等[7]開發出不依賴PH假定的參數化生存曲線模型，該方法主要是將形狀參數和尺度參數進行再參數化的基礎上對風險值建模。Lu等[8]提出分段指數模型，假定生存曲線在不同階段服從PH假設，通過觀察對數累計風險圖來確定分隔時點。Petit等[9]應用限制性平均生存時間（restricted mean survival time，RMST）作為HR的代替指標進行NMA，將不同治療方案的曲線下面積之差作為最終效應量。Jansen等[10]提出分數多項式（fractional polynomials，FP）模型，可對風險值靈活建模，不依賴PH假設。以上4種方法都有一定的應用，且出現在了NICE指南中[11]。張天嵩等[12]提出，在有足夠的KM（Kaplan-Meier）曲線證據時，推薦采用FP方法進行建模。然而，國內還沒有相關文獻對該方法具體介紹。因此，本研究以治療非小細胞肺癌的4種藥物生存數據為例，納入7個研究，在廣義線性模型（generalized Linear model，GLM）的框架下，展示如何基于FP模型進行Non-PH的NMA，基本流程見圖1。
圖1
基于FP模型進行Non-PH的NMA流程圖
1 模型原理簡介
FP模型包括一階FP與二階FP，可表達成式1。其中，p1和p2可在（?2，?1，?0.5，0，0.5，1，2，3）中任取一值。因此，總計8種一階FP模型，36種二階FP模型。在一階FP模型中，當p=0時，FP模型與weibull模型一致；當p=1時，此時FP模型等同于Gompertz模型[13]。
|
研究表明分組生存數據可采用二項似然分布［rjkt～Bin（pjkt，njkt）］，rjkt表示發生事件人數，njkt表示危險人數，pjkt表示t時刻事件發生率）及相應地loglog連接函數建模[14]。一般形式為：
|
式2中，jkt表示t時刻j研究的k治療方案的線性預測方程，ln（tjkt）表示時間間隔造成的偏移量[12]。連續時間線性風險方程通過如下方式轉化：
|
上式最后一步中，假設風險值在［時段內恒定。在趨近于0時，這一假設是成立的。此公式稍作變形便可得到Jansen研究中的風險方程[13]：
|
式4中，對方程的左右兩邊取對數值，即可得到：
|
由此可見，對于分組數據而言，任何線性對數風險方程都可以在GLM中表達。值得注意的是，“線性”并不意味著模型的風險函數隨時間呈線性變化，而是可依賴于任一變量，因此風險函數的形狀是多變的。一般地，在固定效應分散多項式的GLM模型中，jkt可以寫成以下形式[14]：
|
式6中，為j研究的特異效應系數，為治療方法k的特異效應系數。M為FP的階數，一般而言，M取1或2均可。假設A方案在j試驗中，B方案在w試驗中，t時刻A相對于B的HR可表示為式7。
|
此外，GLM模型也可采用隨機效應模型進行擬合，原理與固定效應模型相似，只是在固定效應模型的基礎上加入了隨機效應量。一般可表達成式8：
|
2 軟件安裝與程序加載
2.1 R語言安裝
截止筆者完稿時，R語言最新版本為4.1.2，讀者可在以下網站免費安裝：https://www.r-project.org/。
2.2 程序包加載
R語言中，通過install.packages函數可實現程序包的安裝，library函數可加載相關軟件包。在本研究中，需要的程序包包括：dplyr、ggplot2。通過以下代碼可實現程序的安裝及運行：
install.packages("dplyr")
install.packages("ggplot2")
3 數據預處理與加載
3.1 示例數據的選擇
本研究采用與Jansen[10]研究中一致的數據。共納入7個研究[15-21]、涉及4種治療非小細胞肺癌的治療方法，包括最佳支持治療（best support care，BSC）、多西他賽、吉非替尼和培美曲塞。
3.2 重構個體數據
缺失個體數據時，Guyot[22]的方法是目前最常用的重構方法，其根據KM曲線相關信息即可較好地重構患者個體數據。在R語言中可通過survHE程序包實現，已有文獻詳細講解了該方法[23]。由于篇幅限制，本文不再贅述。
3.3 數據格式整理
為方便在R語言中實現模型的運算，數據需整理成一定格式，見表1。由于數據過多，本研究只展示各項治療方案前6個月的生存數據。表1中，studyn表示研究方案的編號，“1～7”分別表示：“Chang2006”、“Cufer2006”、“Hanna2004”、“Kim2008”、“Lee2010”、“Maruyama2008”及“Shepherd2000”；trtn表示治療方法的編號，“1～4”分別表示：BSC、多西他賽，吉非替尼和培美曲塞；time為研究時限（月）；timeDelta為時間間隔（月），nevents和natrisk分別表示對應時間的發生事件數和危險人數。
表1
示例試驗數據整理
3.4 數據讀取與錄入
將表1保存為NMA.csv文件后，通過“read.csv”命令讀取該表格，使用“data.frame”生成數據矩陣，并命名為“km”；“ref”命令用來規定對照研究及治療方法。“studies”、“treatments”結合factor函數重命名各研究和治療方法。（注：代碼部分中“#”后內容是代碼釋義）
data<- read.csv(C:/Users/lenovo/Desktop/NMA.csv)#括號中為NMA.csv文件的存儲路徑
km <- data.frame(studyn = data$studyn,trtn = data$trtn,time = data$time.c,timeDelta = data$timeDelta,nevents = data$nevents,natrisk = data$narisk)
ref.study <- "Kim2008"
ref.trt <- "Docetaxel"
studies <- c("Chang2006", "Cufer2006", "Hanna2004", "Kim2008", "Lee2010", "Maruyama2008", "Shepherd2000")
treatments <- c("BSC", "Docetaxel", "Gefitinib", "Pemetrexed")
km$studyf <- factor(km$studyn, labels = studies)
km$trtf <- factor(km$trtn, labels = treatments)
km$studyf <- relevel(km$studyf, ref=ref.study)#將對照研究放在最前列
km$trtf <- relevel(km$trtf, ref=ref.trt)
km$trtn=as.numeric(km$trtf)#將因子型變量轉化為數值型
km$studyn=as.numeric(km$studyf)
4 模型構建
4.1 模型結構設置
首先，對全部的FP模型進行擬合。因此，需要生成一個涵蓋所有模型的列表，通過list函數來實現。function函數構造相應模型，具體代碼如下：
models <- list(
"First order FP,p1=0" = list(g1=function(x){log(x)},g2=function(x){0},f1=function(x){log(x)},f2=function(x){0}),
"First order FP,p1=1" = list(g1=function(x){x},g2=function(x){0},f1=function(x){x},f2=function(x){0}),
"Second order FP,p1=-2, p2=1" = list(g1=function(x){x^-2},g2=function(x){x},f1=function(x){x^-2},f2=function(x){x}),
"Second order FP,p1=-1, p2=2" = list(g1=function(x){x^-1},g2=function(x){x^2},f1=function(x){x^-1},f2=function(x){x^2}))
由于篇幅所限，本研究只展示部分模型代碼。讀者可通過改變功能函數中p1和p2自行添加。
4.2 模型擬合
將擬合好的數據保存于“km.new”數據集中，并引入FP模型相關參數。由function構造擬合所需要的功能函數，g1、g2、f1及f2已在模型篩選部分設定，g0和f0為常數項對應系數，默認為1。cblind函數用于構造包含“nevents”及“natrisk-nevents”的矩陣，在GLM框架下擬合相應模型，最后通過lapply函數將所有模型在設置好的結構（“fit.KM.NMA”）下進行擬合。
fit.KM.NMA<-function(bf){
km.new=km
km.new$g0=1
km.new$f0=1
km.new$g1=bf[[1]](km.new$time)
km.new$g2=bf[[2]](km.new$time)
km.new$f1=bf[[3]](km.new$time)
km.new$f2=bf[[4]](km.new$time)
f=cbind(nevents,natrisk-nevents)~trtf*f0+studyf*g0+trtf*f1+trtf*f2+studyf*g1+studyf*g2
glm(f,family=binomial(link=cloglog),data=km.new,offset = log(timeDelta)}#模型選擇二項分布，cloglog函數作為連接函數，時間間隔的對數作為偏移量
fits=lapply(models,fit.KM.NMA)
4.3 擬合度檢驗
在模型擬合過程中，R語言會給出相應的AIC值，使用lapply函數可生成全部擬合模型AIC值。data.frame函數進一步將AIC值整理成可視矩陣。
aics=lapply(fits, AIC)
data.frame(AIC=round(unlist(aics), 2))
根據赤池原則，AIC值越小表明模型擬合效果越佳。結果表明，二階FP（p1=?2，p2=1）總體擬合效果最好。此外，一階FP（p=?2）在所有一階FP模型中擬合最佳。FP模型擬合效果優于標準參數法（p=0或1）。
然而，NICE指南提出AIC最小并不意味外推性能佳。例如，在對成熟度不足的生存數據進行擬合時，AIC值小的二階FP模型可能存在對尾部數據過度擬合的問題，最終的外推效果尚不如一階FP模型。因此，建議考慮先采用多個模型，并通過觀察模型對應的生存曲線確定合適模型。本研究納入一階和二階FP模型各三個，見表2。
表2
納入模型基本信息
4.4 生存數據外推
通過function函數構造模型的結構，pred函數用于模型的預測。km.pred用于保存模型預測的數據及模型結構。具體代碼如下：
pred.KM.NMA<-function(bf){
trts=data.frame(trtf=unique(km$trtf))#剔除重復的治療方案
trts$studyf=sort(unique(km$studyf))[1]#確定對照研究
timehorizon=data.frame(time=2*(1:30))#確定模擬周期
km.pred=merge.data.frame(timehorizon,trts)#在數據矩陣中引入“timehorizon”及“trts”
km.pred$g0=1
km.pred$f0=1
km.pred$g1=bf[[1]](km.pred$time)
km.pred$g2=bf[[2]](km.pred$time)
km.pred$f1=bf[[3]](km.pred$time)
km.pred$f2=bf[[4]](km.pred$time)#生成模型結構
km.pred$timeDelta<-2#設置時間間隔
km.pred}
pred.KM.data=lapply(models,pred.KM.NMA)#批量計算全部模型預測結果
4.5 HR計算
由“predict.glm”命令生成風險預測值“pred”，“%>%”為管道函數，可將左側數據或表達式傳遞到右側。通過filter函數可以篩選出對照治療方案，進而計算其風險函數值“pred.doce”。最后由“mutate(lnhr = pred-pred.doce)”命令計算出各時點全部治療方案與對照方案相比的lnHR。具體代碼如下：
for(i in 1:length(models)){
pred.KM.data[[i]]$pred=predict.glm(fits[[i]], pred.KM.data[[i]])
d=pred.KM.data[[i]]
d1<-d %>% #d賦值于d1
filter(trtf == ref.trt)%>%#篩選出對照治療方案
mutate(pred.doce = pred)%>% #生成對照方案的風險值“pred.doce”作為新的一列
select(pred.doce, time) %>% #選擇“pred.doce”、“time”兩列
arrange(time) %>% #按時間排序“pred.doce”
left_join(d) %>%#生成數據從左邊加入d數據集
mutate(lnhr = pred-pred.doce) %>% #生成lnhr
mutate(hr=exp(lnhr))#%>%#將lnhr轉化為hr
d1$modelc=names(models)[[i]]
if(i==1)dpred<-d1#對照方案的“dpred”賦值于“d1”
if(i!=1)dpred<-rbind(dpred,d1)}#非對照方案，將“dpred”和“d1”合并輸出
dpred$Model<-factor(dpred$modelc,levels=names(models))#dpred中生成新的一列“Model”
生成HR后，還可利用ggplot2對其進行可視化分析，具體代碼如下：
f1<-ggplot() +
geom_line(data=dpred, aes(x=time, y=hr, group=trtf, colour=trtf), size=1)+#圖形相關參數設置
geom_hline(yintercept=1, lty=2)+#添加水平線
#geom_vline(xintercept=14, lty=2)+#添加垂直線
facet_wrap(~Model, nrow=3)+#設置數據排列格式
scale_color_discrete(name="Treatment")+#采用離散標尺給圖形上色
scale_y_log10(limits = c(0.1, 10), breaks = c(0.1, 0.25, 0.5, 1, 2, 4, 10))+#定義縱坐標的刻度及范圍
scale_x_continuous(breaks = c(seq(from=0, to=60, by = 12)))+#定義橫坐標的刻度及范圍
ylab("Hazard Ratio")+#設置縱坐標軸名稱
xlab("Time(months)")+#設置橫坐標軸名稱
theme(legend.position = "bottom")+#設置圖例位置
theme_bw()#使用默認主題
ggsave(plot = f1, filename = "hr.png", width = 12, height = 9)#命名并保存hr圖
運行上述代碼后，即可得到HR的可視化結果，見圖2。
圖2
隨時間變化的HR值
4.6 生存率計算
通過加總各時點風險值（haz），可獲取累計風險值（cumhaz）。由式9，可獲取各試點全部治療方案的生存數據[14]。本部分代碼同3.4基本一致，因此僅對部分代碼進行解釋。
|
for(i in 1:length(models)){
pred.KM.data[[i]]$pred=predict.glm(fits[[i]], pred.KM.data[[i]])
d=pred.KM.data[[i]]
d1<-d %>%
dplyr::mutate(haz = exp(pred))%>%
dplyr::group_by(trtf) %>%
dplyr::arrange(time) %>%
dplyr::mutate(cumhaz = cumsum(haz)) %>%
dplyr::mutate(survProp = exp(-1*cumhaz))
zero <- unique(d1[,c("trtf", "studyf")])
zero$time=0#添加一行，t=0
zero$survProp=1#添加0時刻的生存率
d1=rbind(zero, d1)#合并“zero”與“d1”
d1$modelc=names(models)[[i]]
if(i==1)dpred<-d1
if(i!=1)dpred<-rbind(dpred, d1)}
dpred$Model<-factor(dpred$modelc, levels=names(models))
同樣地，利用ggplot2程序包可繪制各治療方案的生存曲線，見圖3。代碼如下：（代碼相關釋義與3.4一致）
圖3
各模型對應生存曲線
f2= ggplot() +
geom_line(data=dpred, aes(x=time, y=survProp, group=trtf, colour=trtf), size=1)+
scale_color_discrete(name="Treatment")+
expand_limits(y=c(0,1), x=c(0,60))+
facet_wrap(~Model, nrow=3)+
scale_x_continuous(breaks = c(seq(from=0, to=60, by = 12)))+
ylab("Proportion surviving")+
xlab("Time(months)")+
guides(color = guide_legend(ncol = 1))+
theme(legend.position = "bottom")+
theme_bw()
ggsave(plot = f2, filename = "paperSurv.png", width = 12, height = 9)
本研究中BSC生存數據相對較短，部分FP模型對其尾部數據過度擬合。通過圖2可知，2階FP模型（p1=1，p2=2或3）對于BSC的外推結果較差，其長期生存獲益優于積極治療，與臨床實際效果不符，需排除。
4.7 最優模型的確定和回歸系數生成
根據各模型AIC值及生存曲線圖，“p1=?2，p2=1”的二階FP模型總體AIC值最小，且生存曲線擬合較好。因此本研究與Jansen等研究一致，選擇該模型作為最優模型[13-14]。采用以下代碼合成最終的模型：
pred.KM.data=pred.KM.NMA(models[["Second order FP,p1=-2, p2=1"]])#“p1=-2,p2=1”的二階FP模型預測數據作為最終的pred.KM.data
最后，通過以下代碼可提取各治療方案相比于對照方案的回歸系數，即式6中的“，其中，m取0，1，2。將回歸系數代入式7即可算出各時點各方案相比于對照方案的HR值。見表3。
表3
最優模型參數估計
timeh<-data.frame(time=1:60,a=1)#確定模擬時限
fit<-fits[["Second order FP,p1=-2, p2=1"]]#確定最終擬合模型
pred.data<-fit$data
pred.data$est1<-predict.glm(fit)#獲取比較系數均值
pred.data$est2=NA
pred.data$se<-predict.glm(fit,se.fit=T)$se.fit#獲取比較系數標準誤
coeff.data<- as.data.frame(coef(summary(fits[["Second order FP,p1=-2, p2=1"]])))[,1:2]#選取模型擬合回歸結果的前兩列，構建數據集“coeff.data”
names(coeff.data)<- c("est","std.err")#對前兩列重命名，分別為"est","std.err"
attach(coeff.data)#綁定數據集
coeff.data$conf.int.lower<- est-std.err*qnorm(0.975)#生成CI下限
coeff.data$conf.int.upper<- est+std.err*qnorm(0.975)#生成CI上限
coeff.data[,1:4]<- round(coeff.data[,1:4],3)#數值保留三位小數
coeff.data$pn<- rownames(coeff.data)#添加治療方案名稱
coeff.data$contrast<- paste(levels(km$trtf)[-1],"vs",ref.trt)#添加對比方案名稱
detach(coeff.data)#解綁數據集
coeff.data#生成回歸系數
5 討論
本研究采用治療非小細胞肺癌的4種藥物的生存曲線構建基于FP模型的非恒定比例風險NMA，對其原理及方法進行介紹，并采用實際案例展示了如何在R語言中構建模型、擬合生存曲線。本研究分析結果表明，相比于標準曲線，如岡茲分布、威布爾分布，FP模型的擬合效果更優。通過視覺檢查，基于FP模型重構的生存曲線擬合也能較好地貼合臨床實際。
隨著PD-（L）1等創新藥物的出現，腫瘤藥物生存曲線形狀更加多樣的同時，相比于傳統的化療往往有明顯的拖尾效應，因此PH模型不再適用。Non-PH模型中，RMST模型只能截取KM時限內的生存數據進行分析，無法外推的特點也限制了其應用；而參數化生存曲線模型及分段指數模型雖然一定程度上放松了PH假設，但當生存曲線的風險率不再單調變化時，這類模型便不再適用。
FP模型的優勢包括：① 相比于傳統NMA模型中，FP模型更能靈活捕捉生存曲線的特征，允許各治療方案生存曲線有不同形狀，不再依賴于PH假設，因此基于FP的NMA擬合效果較傳統NMA好。② FP模型可擬合到生存曲線的漸近線，相比于標準多項式或樣條模型，末端數據更加穩定。③ FP模型可用于預測KM時點外的生存數據，方便測算干預措施終生投入和產出，對于衛生決策意義重大[13]。同時，FP模型有一定的不足，容易受尾部稀疏數據影響，產生過度擬合。因此需觀察多個生存曲線或風險圖以確定最合適的模型參數。
鑒于FP模型的優越性，其應用越來越多。如，Wang等在對前列腺癌進行NMA時采用了一階FP模型[24]；Vickers等[25]采用FP模型評價非小細胞肺癌二線治療方案的療效。此外，NICE相關HTA評價中也多次用到了FP模型[10]。
由于篇幅有限，本研究只在頻率學框架中展示了基于FP模型在NMA中的應用，后續研究將進一步闡述如何在貝葉斯框架下實現模型的搭建。
聲明 所有作者均聲明無利益沖突。
