薛芳靜 黃 灝 許碧云 張慧敏 陳炳為△

?

·計(jì)算機(jī)應(yīng)用·

時(shí)間序列數(shù)據(jù)中的干預(yù)分析模型及SAS實(shí)現(xiàn)

薛芳靜1黃 灝1許碧云2張慧敏1陳炳為1△

目的 利用干預(yù)分析模型從定量分析的角度來(lái)評(píng)估政策干預(yù)或突發(fā)事件對(duì)時(shí)間序列數(shù)據(jù)的具體影響。方法 本文以2010年1月至2013年5月某醫(yī)院ICU抗生素多重耐藥性的時(shí)間序列為例進(jìn)行干預(yù)性分析。利用SAS中的PROC ARIMA程序可定量分析政策法規(guī)對(duì)抗生素多重耐藥性的影響。結(jié)果 多重耐藥性的干預(yù)分析模型中趨勢(shì)估計(jì)參數(shù)為-0.00174，抗生素多重耐藥性呈下降趨勢(shì)，干預(yù)的MA系數(shù)為-0.662，干預(yù)變量系數(shù)為-0.117，說(shuō)明限制抗生素使用的政策法規(guī)降低了醫(yī)院抗生素使用量，使抗生素多重耐藥性下降11.04%。結(jié)論 干預(yù)分析模型可以較好地處理“突變性”的時(shí)間序列問(wèn)題，SAS中的PROC ARIMA程序?yàn)樵撃Ｐ吞峁┝艘环N極為方便的實(shí)現(xiàn)形式。

干預(yù)分析模型 時(shí)間序列 抗生素

時(shí)間序列預(yù)測(cè)方法是以慣性原理推測(cè)其未來(lái)的狀態(tài)，常見的有指數(shù)平滑法、ARIMA方法等[1]。但這些方法不能解決“突變性”的時(shí)間序列，如政策的實(shí)施或干預(yù)，這必將對(duì)事件發(fā)展過(guò)程產(chǎn)生影響，甚至改變事件發(fā)展結(jié)果。Box和Tiao于1975年提出干預(yù)分析模型(intervention analysis model)，較好地解決了這種“突變性”的時(shí)間序列。其被廣泛用于經(jīng)濟(jì)政策的變化或突發(fā)事件對(duì)經(jīng)濟(jì)環(huán)境、經(jīng)濟(jì)過(guò)程或結(jié)果的具體影響的定量分析[2-3]。干預(yù)分析模型是將干預(yù)事件和時(shí)間序列模型結(jié)合起來(lái)進(jìn)行研究的，它較好地解決了政策的變化或突發(fā)事件的影響及影響程度。如Bordet等在德國(guó)一家社區(qū)醫(yī)院實(shí)施抗生素管理計(jì)劃后利用干預(yù)分析模型得到該醫(yī)院頭孢菌素使用下降了33%，氟喹諾酮類使用下降了31%，證實(shí)在干預(yù)后頭孢菌素類和氟喹諾酮類藥物的使用密度趨勢(shì)有顯著的變化，抗生素總使用降低了11%[4]。

由于受到我國(guó)開始正式實(shí)施《抗菌藥物臨床應(yīng)用管理辦法》的影響，ICU的抗生素多重耐藥性在干預(yù)前后發(fā)生了變化。本文利用SAS中的ARIMA過(guò)程，應(yīng)用干預(yù)分析模型對(duì)某醫(yī)院ICU的抗生素多重耐藥性進(jìn)行分析，以期得到干預(yù)前后抗生素多重耐藥性的變化情況。

原理與SAS程序

1.干預(yù)變量

干預(yù)分析模型是在常規(guī)時(shí)間序列中加入干預(yù)變量，干預(yù)變量有兩種形式[3,5]：一種是持續(xù)性的干預(yù)變量，即在T時(shí)刻發(fā)生后一直有影響，這時(shí)可以用階躍函數(shù)表示：

第二種是短暫性的干預(yù)變量，即在事件發(fā)生時(shí)刻T有影響，用單位脈沖函數(shù)表示：

2.干預(yù)事件的形式

對(duì)于干預(yù)時(shí)間序列模型，干預(yù)事件按其影響有四種最基本的形式。其他的干預(yù)模型，不管是否受到多種干預(yù)，也不管這些干預(yù)影響有多復(fù)雜，都可用這四種形式組合表達(dá)[3,5]：

(1)干預(yù)事件的影響突然開始并持續(xù)

設(shè)干預(yù)事件在T時(shí)刻突然開始，干預(yù)對(duì)因變量的影響是固定且未知的。這時(shí)干預(yù)模型可寫為：

其中ω表示干預(yù)影響程度的未知參數(shù)，B為后移算子。當(dāng)Xt不平穩(wěn)時(shí)，通過(guò)差分使序列成為平穩(wěn)序列；當(dāng)Xt平穩(wěn)時(shí)，則不需差分。Bb是代表干預(yù)事件在b個(gè)時(shí)期(b=1,2,…)滯后才產(chǎn)生影響。

(2)干預(yù)事件的影響逐漸增大且持續(xù)

干預(yù)事件在T時(shí)刻突然發(fā)生，但其影響不是立刻全部產(chǎn)生，而是隨著時(shí)間逐漸增大，這時(shí)模型可用下式表達(dá)：

(3)干預(yù)事件的影響突然開始且暫時(shí)性

這類干預(yù)現(xiàn)象可以用數(shù)學(xué)模型描述如下：

當(dāng)δ=0時(shí)，干預(yù)的影響只存在一個(gè)時(shí)期，當(dāng)δ=1時(shí)，干預(yù)的影響將長(zhǎng)期存在。

(4)干預(yù)事件的影響逐漸增大后衰減

干預(yù)事件在T時(shí)刻突然發(fā)生，其影響逐漸增加，在某個(gè)時(shí)刻到達(dá)峰值，然后又逐漸減弱甚至消失，模型可表達(dá)為：

3.干預(yù)分析模型的構(gòu)造與干預(yù)效應(yīng)的識(shí)別

(1)模型的構(gòu)造

干預(yù)模型主要是為了度量干預(yù)效應(yīng)，剔除了干預(yù)影響后的序列可看作普通的序列，進(jìn)而可以用一般的時(shí)間序列分析方法進(jìn)行擬合。本文以Box-Jenkins的ARIMA模型進(jìn)行干預(yù)模型的構(gòu)造。

(2)干預(yù)效應(yīng)的識(shí)別

根據(jù)序列的具體情況和干預(yù)變量的性質(zhì)進(jìn)行識(shí)別。利用干預(yù)變量產(chǎn)生影響之前或干預(yù)影響之后的凈化數(shù)據(jù)計(jì)算出自相關(guān)函數(shù)與偏自相關(guān)函數(shù)。首先識(shí)別ARIMA模型中的p和q，然后估計(jì)出θB、φB中的參數(shù)。

(3)模型的判斷

模型的擬合采用AIC統(tǒng)計(jì)量進(jìn)行評(píng)價(jià)，依據(jù)AIC最小原則，選出最優(yōu)ARIMA模型。加入干預(yù)變量后通過(guò)比較模型參數(shù)的t值和P值評(píng)價(jià)最終模型。

4.SAS軟件的實(shí)現(xiàn)

本文利用SAS9.3的ARIMA過(guò)程步來(lái)實(shí)現(xiàn)，具體程序如下：

(1)ARIMA模型的識(shí)別

建立數(shù)據(jù)集

data a;

input y @@;

format date monyy.;

date=INTNX(′month′,′01JAN10′d,_N_-1);

x=(date>=′01JUN11′d);

proc gplot data=a;

symbol1 i=spline v=dot c=red;

plot y*date=1;

在DATA步利用INTNX函數(shù)生成時(shí)間序列，并生成階梯函數(shù)X，利用GPLOT作出序列圖。再利用ARIMA過(guò)程的identify語(yǔ)句對(duì)模型進(jìn)行初步識(shí)別，觀察樣本自相關(guān)系數(shù)、偏相關(guān)系數(shù)及是否具有周期性等問(wèn)題。

proc arima data=a;

identify var=y(1);

estimate q=(1)method=uls plot;

forecast lead=12 out=b id=date interval=month id=date;

run;

在ARIMA(p,d,q)模型擬合中，identify語(yǔ)句后的y(1)代表進(jìn)行一階差分，estimateq=(1)考慮模型中q=1的模型，即此時(shí)模型為ARIMA(0,1,1)，同時(shí)分別考慮ARIMA(1,1,0)、ARIMA(1,1,1)，從三種模型中選出最優(yōu)模型。Forecast語(yǔ)句是按月份進(jìn)行12個(gè)月的預(yù)測(cè)。完成上述三個(gè)模型結(jié)果，選取AIC最小的模型進(jìn)行白噪聲檢驗(yàn)，同時(shí)可利用gplot程序作殘差序列圖，觀察殘差的特性。

(2)干預(yù)變量

proc arima data=a;

identify var=y(1)cross=(x(1));

estimate p=(1)q=(1)input=(x)method=cls;

forecast lead=12 out=a id=date interval=month;

run;

在identify語(yǔ)句中利用cross選項(xiàng)指定協(xié)變量x，在estimate利用input語(yǔ)句指定ARIMA模型中加入干預(yù)x，input后的變量一定要利用identify語(yǔ)句中的cross選項(xiàng)事先指定。

實(shí)例分析

據(jù)WHO資料顯示，我國(guó)門診感冒患者抗生率使用率為75%、外科手術(shù)高達(dá)95%、住院患者80%，其中使用廣譜抗生素和聯(lián)合應(yīng)用的患者占58%，遠(yuǎn)高于國(guó)際上規(guī)定的抗生素使用率不超過(guò)30%的標(biāo)準(zhǔn)，而且我國(guó)每年有8萬(wàn)多人直接或間接死于濫用抗生素，由此造成的機(jī)體損傷及病原菌耐藥性更是無(wú)法估量[6]。濫用抗生素將會(huì)使人類為對(duì)抗疾病付出越來(lái)越高的代價(jià)。2012年8月1日，我國(guó)正式開始實(shí)施《抗菌藥物臨床應(yīng)用管理辦法》，對(duì)抗生素臨床使用做了嚴(yán)格規(guī)定，此辦法的實(shí)施必然對(duì)我國(guó)抗生素的臨床使用量產(chǎn)生影響。本文利用干預(yù)分析模型對(duì)2010年1月至2013年5月某醫(yī)院ICU的抗生素多重耐藥性進(jìn)行干預(yù)前后間的差別分析，干預(yù)后的數(shù)據(jù)從2011年7月到2013年5月，數(shù)據(jù)見表1。使用上述SAS程序進(jìn)行分析，首先繪制序列圖(圖1)，觀察到該醫(yī)院抗生素的多重耐藥性從實(shí)施干預(yù)時(shí)開始減少。

隨之采用自相關(guān)圖檢驗(yàn)法來(lái)檢驗(yàn)平穩(wěn)性，圖2是樣本自相關(guān)系數(shù)散點(diǎn)圖，自相關(guān)系數(shù)緩慢趨于零，表明抗生素多重耐藥性序列是非平穩(wěn)的。因?yàn)樾蛄惺欠瞧椒€(wěn)的，采用一階差分進(jìn)行平穩(wěn)化處理。再次考察自相關(guān)系數(shù)圖，處理后的數(shù)據(jù)滿足平穩(wěn)性，同時(shí)通過(guò)白噪聲檢驗(yàn)，可以進(jìn)行ARIMA建模。用ARIMA(1,1,0)、ARIMA(1,1,1)、ARIMA(0,1,1)三種模型來(lái)進(jìn)行擬合。

表1 2010年1月至2013年5月某醫(yī)院ICU的抗生素多重耐藥性

圖1 某醫(yī)院2010年1月至2013年5月ICU的抗生素多重耐藥性時(shí)間序列圖

圖2 樣本自相關(guān)系數(shù)散點(diǎn)圖

模型ARIMA(1,1,0)ARIMA(1,1,1)ARIMA(0,1,1)AIC-125.632-125.029-126.977

表2分別列出了模型的擬合情況。根據(jù)AIC最小準(zhǔn)則，ARIMA(0,1,1)為最優(yōu)模型。加入干預(yù)變量進(jìn)行分析，結(jié)果見表3。

表3 加入干預(yù)變量后的ARIMA模型參數(shù)檢驗(yàn)結(jié)果

多重耐藥性的干預(yù)分析模型方程為Yt=-0.00174-0.117It+(1-0.662B)at，μ=-0.00174表明多重耐藥性呈下降趨勢(shì)，其下降量為1-e-0.00174=1.70%；干預(yù)的MA模型系數(shù)為-0.662，P<0.0001；干預(yù)變量的系數(shù)為-0.117，P=0.0046，說(shuō)明限制抗生素使用的政策法規(guī)降低了醫(yī)院抗生素多重耐藥性，進(jìn)而使抗生素多重耐藥性下降了1-e-0.117=11.04%。

討 論

延誤 ICU抗生素合理用藥(對(duì)病原體使用敏感藥物)與住院死亡率息息相關(guān)，抗生素用藥史是延誤用藥的最重要危險(xiǎn)因素之一，且抗生素用藥時(shí)間延長(zhǎng)也是危重癥患者出現(xiàn)耐藥菌感染的最重要誘發(fā)因素。因此，ICU臨床醫(yī)生必須權(quán)衡兩個(gè)問(wèn)題：一是抗生素用藥需滿足現(xiàn)在所救治患者的需求，二是兼顧 ICU其他患者可能暴露于該抗生素治療患者誘發(fā)的耐藥菌感染[7]。在ICU，臨床醫(yī)生常會(huì)遇到的棘手問(wèn)題是使用抗生素治療危重癥嚴(yán)重感染患者的同時(shí)，又要盡量減少耐藥菌的出現(xiàn)及其傳播。

干預(yù)分析模型是時(shí)間序列分析模型的一種推廣,如果不存在干預(yù)影響,這兩種模型就沒有什么差別。干預(yù)分析模型可以對(duì)具有“突變性”的時(shí)間序列進(jìn)行模擬化處理，定量分析政策干預(yù)或突發(fā)事件對(duì)事件過(guò)程結(jié)果的具體影響。干預(yù)分析模型與其他模型相比，有兩個(gè)優(yōu)點(diǎn)[8]：(1)干預(yù)事件影響存在，那么對(duì)時(shí)間序列分析和預(yù)測(cè)不論是用ARIMA 模型還是MARIMA模型(多變量時(shí)間序列的ARIMA模型)其預(yù)測(cè)精度都比干預(yù)分析模型差，也就是干預(yù)分析模型更接近實(shí)際，分析的結(jié)果更為客觀。(2)干預(yù)分析模型比回歸分析模型中加入虛擬量要靈活得多，適應(yīng)能力強(qiáng)一些，因?yàn)楦深A(yù)分析模型能對(duì)時(shí)間序列的動(dòng)態(tài)特性進(jìn)行合理的描述,甚至可以對(duì)未來(lái)的影響做出主觀估計(jì)，把先驗(yàn)信息反映到模型中。

干預(yù)分析模型從其出現(xiàn)時(shí)起，主要應(yīng)用于經(jīng)濟(jì)政策的變化或突發(fā)事件對(duì)經(jīng)濟(jì)環(huán)境、經(jīng)濟(jì)過(guò)程或結(jié)果的具體影響的定量分析，較少被引入醫(yī)學(xué)領(lǐng)域的研究。SAS中的PROC ARIMA程序?yàn)樵撃Ｐ吞峁┝艘环N極為方便的實(shí)現(xiàn)形式。

[1]王濤,苑新海,朱宗龍.ARIMA時(shí)間序列在乙肝發(fā)病預(yù)測(cè)中的應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2014,31(4):646-647,650.

[2]Box G,Tiao GC.Intervention Analysis with Applications to Economic and Environmental Problems.Journal of the American Statistical Association,1975,70(349):70-79.

[3]徐國(guó)祥.統(tǒng)計(jì)預(yù)測(cè)和決策.第2版.2005:上海財(cái)經(jīng)大學(xué)出版社,396.

[4]Borde JP,Litterst S,Ruhnke M,et al.Implementing an intensified antibiotic stewardship programme targeting cephalosporin and fluoroquinolone use in a 200-bed community hospital in Germany.Infection,2015,43(1):45-50.

[5]馮文權(quán),蔡基棟.干預(yù)分析模型及其應(yīng)用.武漢大學(xué)學(xué)報(bào)(自然科學(xué)版),1994,33(6):29-35.

[6]胡燕,白繼庚,胡先明,等.我國(guó)抗生素濫用現(xiàn)狀、原因及對(duì)策探討.中國(guó)社會(huì)醫(yī)學(xué)雜志,2013,30(2):128-130.

[7]Kollef MH,Micek ST.Rational Use of Antibiotics in the ICUBalancing Stewardship and Clinical Outcomes.JAMA,2014，321(14):1403-1404.

[8]馮文權(quán).政策干預(yù)或突發(fā)事件影響的統(tǒng)計(jì)分析模型.統(tǒng)計(jì)與決策,1990，(6):8-11.

(責(zé)任編輯：劉 壯)

1.東南大學(xué)流行病與衛(wèi)生統(tǒng)計(jì)學(xué)系(210009)

2.南京大學(xué)醫(yī)學(xué)院附屬鼓樓醫(yī)院

△通信作者：陳炳為，E-mail:drchenbw@126.com