連續(xù)型變量界值點(diǎn)敏感度分析的SAS宏程序*

昆明醫(yī)科大學(xué)公共衛(wèi)生學(xué)院流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(650500) 肖媛媛 李曉梅 何利平 孟 瓊

【提 要】 目的 探討如何用SAS宏程序?qū)崿F(xiàn)連續(xù)型變量不同截?cái)嘀祵?duì)研究結(jié)果影響的敏感度分析。方法 采用自行編制的SAS宏程序，結(jié)合具體實(shí)例，演示連續(xù)型變量界值點(diǎn)敏感度分析的實(shí)現(xiàn)及結(jié)果的呈現(xiàn)。結(jié)果 該SAS宏程序可以按照研究者不同分析目的及設(shè)計(jì)精度進(jìn)行連續(xù)型變量不同截?cái)嘀档拿舾卸确治觥＝Y(jié)論 當(dāng)前研究所提供的SAS宏程序簡(jiǎn)潔靈活，可在涉及定量變量二分類處理的醫(yī)學(xué)資料分析中推廣應(yīng)用，對(duì)研究結(jié)果的穩(wěn)健性進(jìn)行更為全面的評(píng)價(jià)。

連續(xù)型變量的分類轉(zhuǎn)換在醫(yī)學(xué)資料的統(tǒng)計(jì)分析實(shí)踐中十分常見，尤其是根據(jù)某特定指標(biāo)的取值將研究對(duì)象進(jìn)行二分類(dichotomization)。常用的分類依據(jù)有兩種：一是采用公認(rèn)的醫(yī)學(xué)參考值(medical reference value)來界定；而大多無參考值或無固定參考值的指標(biāo)，則一般通過繪制受試者工作特征(receiver operating characteristic，ROC)曲線來確定最為適宜的截?cái)帱c(diǎn)(cut-off point)[1-2]。

敏感度分析(sensitivity analysis)是經(jīng)濟(jì)學(xué)領(lǐng)域的一種常用研究方法，概括而言，其研究特定模型自變量(輸入)的變動(dòng)對(duì)因變量(輸出)的影響[3]。由上不難看出，不論何種方法，對(duì)于連續(xù)型變量，用某一特定截?cái)帱c(diǎn)取值或參考值區(qū)間將其判定為截然不同的兩個(gè)分類顯然會(huì)有些武斷，從而導(dǎo)致對(duì)分類后分析結(jié)果可靠性的質(zhì)疑。因此，十分有必要分析當(dāng)界值點(diǎn)在其取值附近變動(dòng)時(shí)，分析結(jié)果是否穩(wěn)定，即對(duì)分析結(jié)果對(duì)于界值點(diǎn)取值變化進(jìn)行敏感度分析。

由于連續(xù)型變量在其定義域內(nèi)的任意區(qū)間均存在無窮多個(gè)取值，因此可行的敏感度分析策略一定是選擇某個(gè)或某些區(qū)間內(nèi)、一定精度下、有限但足夠多個(gè)取值，分別探討其對(duì)因變量的影響，這樣的分析意圖往往涉及次數(shù)龐大的迭代(iteration)，借助現(xiàn)有的編程式統(tǒng)計(jì)分析軟件(如SAS和R)方可便捷實(shí)現(xiàn)。但文獻(xiàn)檢索發(fā)現(xiàn)，目前尚無相關(guān)研究報(bào)道。本文將提出一種連續(xù)型變量界值點(diǎn)敏感度分析的SAS宏程序，并結(jié)合實(shí)例演示其分析結(jié)果，為推廣敏感度分析在醫(yī)學(xué)研究中的應(yīng)用提供便利，也為提升醫(yī)學(xué)研究分析結(jié)果的穩(wěn)健性(robustness)提供有力工具。

設(shè)計(jì)思路

SAS宏程序

根據(jù)不同研究目的，本文提供多重線性回歸模型(multiple linear regression model)、logistic回歸模型(logistic regression model)、Cox比例風(fēng)險(xiǎn)模型(Cox proportional hazards model)、動(dòng)態(tài)生存分析的Anderson-Gill計(jì)數(shù)過程模型(Anderson-Gill counting process model，AG模型)四種模型中連續(xù)型變量界值點(diǎn)敏感度分析的備選宏程序模塊，讀者也可依據(jù)SAS編程規(guī)則自行拓展至其他一般線性模型(generalized linear model)。

具體宏程序、備選模塊及注釋如下：

%macrosensitivity(start=，end=，interval=，)；

%do i=&start %to &end %by &interval；

data database2；

set database1；

if variable<=&i then group=0；

if variable>&i then group=1；

run；

/*注釋：variable為擬進(jìn)行敏感度分析的連續(xù)型變量，group為根據(jù)截?cái)帱c(diǎn)產(chǎn)生的二分類變量，database1為初始數(shù)據(jù)庫，database2為根據(jù)新的截?cái)帱c(diǎn)調(diào)整生成的數(shù)據(jù)庫*/

ods output ParameterEstimates=estimate1；/*注釋：estimate1為模型擬合結(jié)果輸出數(shù)據(jù)庫*/

/*以下提供Cox比例風(fēng)險(xiǎn)模型和動(dòng)態(tài)生存分析的AG模型兩種模型中連續(xù)型變量界值點(diǎn)敏感度分析的宏程序備選模塊，讀者可根據(jù)分析需要，結(jié)合Cox比例風(fēng)險(xiǎn)模型的宏程序?qū)崿F(xiàn)多重線性回歸分析及l(fā)ogistic回歸分析*/

/*Cox比例風(fēng)險(xiǎn)模型*/

proc phreg data=database2；

class var1 var2；

model(start end)*censor(0)=variable var1 var2 var3 var4 … varN；

run；

/*注釋：variable為擬進(jìn)行敏感度分析的連續(xù)型變量，var1～varN為Cox比例風(fēng)險(xiǎn)模型中的其他協(xié)變量*/

/*動(dòng)態(tài)生存分析的AG模型*/

proc phreg data=database2 covs(aggregate)covm；

class var1 var2；

model(start end)*censor(0)=variable var1 var2 var3 var4 … varN；

id subject；

run；

/*注釋：covs(aggregate)語句給出由Lin和Wei提出的穩(wěn)健的夾心方差估計(jì)[5]，model語句中，variable為擬進(jìn)行敏感度分析的連續(xù)型變量，var1～varN為Cox比例風(fēng)險(xiǎn)模型中的其他協(xié)變量，id語句為夾心方差的調(diào)整依據(jù)，subject為研究對(duì)象的唯一識(shí)別編碼*/

ods output close；

data estimate2；

set estimate1；

cutoff=&i；

logistic回歸模型，效應(yīng)值為OR

Cox比例風(fēng)險(xiǎn)模型，效應(yīng)值為HR

AG模型，效應(yīng)值為HR

where Parameter=“variable” and StdErrRatio ne.；

keep point lower upper cutoff；

run；

/*注釋：point，lower，upper分別是效應(yīng)值點(diǎn)估計(jì)值、95%可信區(qū)間下限和上限*/

%if &i=&start %then %do；

data estimate3；

set estimate2；

%end；

%else %do；

proc append base=estimate3 data=estimate2；run；

%end；

/*注釋：estimate3為最后輸出的效應(yīng)估計(jì)值及其95%可信區(qū)間匯總數(shù)據(jù)庫*/

dm ′clear output′；

dm ′clear log′；

dm ′odsresults；clear′；

%end；

%mend；

宏程序分析結(jié)果示例

對(duì)于本例，我們?cè)O(shè)置精度為1 IU/L。則選擇上述宏程序模塊3，提交如下宏程序執(zhí)行語句：

%simulation(start=70，end=120，interval=1)；

執(zhí)行后輸出的效應(yīng)值及其95%可信區(qū)間匯總結(jié)果estimate3內(nèi)容如表1所示。

表1 匯總結(jié)果數(shù)據(jù)庫內(nèi)容

可以選擇采用繪圖軟件如Sigmaplot等更好地直觀展示敏感度分析結(jié)果，所繪制圖形如下：

從圖1不難看出：隨著ALP截?cái)帱c(diǎn)在70到120間變化，對(duì)于ALP較高的早期胰腺癌病人來說，其總體死亡風(fēng)險(xiǎn)基本上均顯著高于ALP較低的病人?？梢哉J(rèn)為研究結(jié)果穩(wěn)健性較好。

圖1 不同ALP截?cái)帱c(diǎn)敏感度分析結(jié)果

討 論

廣泛意義上的敏感度分析涵蓋很多方面，除了本文研究?jī)?nèi)容所涉及的研究變量不同定義外，還包括探討不同統(tǒng)計(jì)分析方法、實(shí)驗(yàn)流程變動(dòng)、缺失值(missing data)及離群值(outliner)處理等對(duì)研究結(jié)果的影響[6]。由此不難看出，對(duì)于醫(yī)學(xué)資料的統(tǒng)計(jì)分析而言，敏感度分析是十分必要的。但與之形成鮮明對(duì)比的是，目前在醫(yī)學(xué)研究領(lǐng)域已發(fā)表的論文中，敏感度分析的使用比例卻非常低。在一項(xiàng)2013年的研究中，Thabane等對(duì)64篇已發(fā)表的英文論文進(jìn)行梳理，發(fā)現(xiàn)僅有13篇在統(tǒng)計(jì)分析中采用了任意形式的敏感度分析，占20.3%[7]。

針對(duì)敏感度分析在醫(yī)學(xué)研究中應(yīng)用不足、亟待推進(jìn)的現(xiàn)狀，本文根據(jù)醫(yī)學(xué)資料統(tǒng)計(jì)分析實(shí)踐中經(jīng)常涉及對(duì)連續(xù)型變量根據(jù)某一截?cái)嘀颠M(jìn)行二分類處理的問題，提出了一種如何根據(jù)截?cái)帱c(diǎn)的不同取值，對(duì)研究結(jié)果進(jìn)行敏感度分析的SAS宏程序。該宏程序簡(jiǎn)潔、靈活，能根據(jù)研究人員的不同精度要求及分析意圖執(zhí)行敏感度分析，為提升醫(yī)學(xué)研究統(tǒng)計(jì)分析結(jié)果的穩(wěn)健性提供有力支持。