方園
安徽省第二人民醫(yī)院 (安徽合肥 230000)
時間序列為按照時間順序獲得的一系列觀測值的組合。隨著科學技術的不斷發(fā)展,對于時間序列的研究被應用于金融、氣象、農(nóng)業(yè)生產(chǎn)等領域[1-3]。時間序列數(shù)據(jù)挖掘工作的主要任務是獲取數(shù)據(jù)中蘊含的與時間相關的信息,掌握事物發(fā)展的規(guī)律。通過對時間歷史數(shù)據(jù)的挖掘,我們可以預測未來一段時間可能發(fā)生的數(shù)據(jù),并以此分析事物未來的發(fā)展趨勢[4]。
在放射治療中,醫(yī)用直線加速器是不可或缺的設備,其輸出劑量直接決定了患者治療的預期效果,因此,保證輸出劑量的穩(wěn)定性非常重要[5]。由于醫(yī)用直線加速器是一類結構復雜的高精度大型醫(yī)療設備,某些元器件的故障或老化均會導致加速器輸出劑量發(fā)生突發(fā)的或緩慢的變化,因此,對加速器輸出劑量的日常安全監(jiān)測成了放射治療質(zhì)量保證(quality assurance,QA)工作的重點[6-7]。很多學者利用統(tǒng)計學處理方法開展了大量關于醫(yī)用直線加速器劑量穩(wěn)定性的研究,且取得了較好的研究成果,為加速器劑量穩(wěn)定性的監(jiān)測工作提供了參考依據(jù)[8-10]。也有學者利用時間序列分析方法對醫(yī)用直線加速器質(zhì)量穩(wěn)定性進行了分析和預測,并證實了此方法在加速器劑量預測工作中的應用可行性[11]。鑒于此,本研究對醫(yī)用直線加速器質(zhì)控工作中的劑量監(jiān)測數(shù)據(jù)進行了研究,并利用時間序列挖掘知識進行了數(shù)據(jù)分析,以達到對加速器劑量偏移進行預測的目的。
本研究的數(shù)據(jù)來源于我院放射治療中心醫(yī)科達precise 直線加速器,選取3個調(diào)整周期內(nèi)共計66條加速器劑量實測數(shù)據(jù)作為實驗樣本。
首先,選取1個調(diào)整周期內(nèi)的前19條劑量監(jiān)測數(shù)據(jù)作為原始研究數(shù)據(jù),對其進行穩(wěn)定性分析;然后,利用時間序列研究中常用的差分自回歸移動平均(auto regressive integrated moving average,ARIMA)模型對其進行建模并預測后5次的測量結果;最后,將預測結果與另外兩個測量周期用同樣預測方法得出的預測結果進行對比分析。
ARIMA 模型于1976年由Box 和Jenkins 等學者提出,之后被廣泛應用于經(jīng)濟、金融和醫(yī)療等領域[12-14]。ARIMA(p,d,q)模型的通式為:
該模型分為前后兩部分,前部分?1yt-1+?2yt-2+...+?pyt-p是自回歸方程,后部分et-(θ1et-1+θ2et-2+...+θqet-q)是誤差移動方程,其中,yt是當前值,?1、?2...是自相關系數(shù),θ1、θ2...是誤差系數(shù),yt-1、yt-2...及et-1、et-2...是滯后項,et是誤差,p 是自回歸項數(shù),q 為移動平均項數(shù),d 為時間序列成為平穩(wěn)時所作的差分次數(shù)。
存在一定發(fā)展趨勢的時間序列均是非平穩(wěn)的,而ARIMA模型必須建立在平穩(wěn)的基礎上才有意義[15]。因此,時間序列建模首先應去除確定性因素,包括趨勢和季節(jié)性因素,然后檢測剩下數(shù)據(jù)的平穩(wěn)性,最終通過檢驗的數(shù)據(jù)才可進行ARIMA建模,建模的具體參數(shù)可以通過殘差自相關函數(shù)(autocorrelation function,ACF)和偏自相關函數(shù)(partial autocorrelation function,PACF)值來確定,R語言中有1個被命名為“auto.arima”的函數(shù),可以用來直接定參,然后用此參數(shù)進行模型的擬合和預測[16]。
auto.arima 函數(shù)是Hyndman-Khandakar 算法的1個變種,其結合了單位根檢驗、最小化赤池信息準則(Akaike information criterion,AIC)和極大似然估計(maximum likelihood estimate,MLE)等評價標準來獲得1個ARIMA 模型。
Hyndman-Khandakar 自動ARIMA 建模算法步驟如下:
步驟1:通過重復地進行單位根檢驗(Kwiatkowski Phillips Schmidt Shin test,KPSS)測試來確定差分階數(shù)d(0≤d≤2)。
步驟2:對數(shù)據(jù)差分d 次之后,通過最小化AIC來選擇最優(yōu)的p,q;算法通過stepwise search 而不是遍歷所有可能的p,q 組合來尋找最優(yōu)的p,q 組合。
我們通過對周期內(nèi)的數(shù)據(jù)進行觀察,發(fā)現(xiàn)其呈現(xiàn)逐步上升的趨勢且不存在季節(jié)性因素,劑量實測數(shù)據(jù)序列見圖1。
圖1 劑量實測
經(jīng)過一階差分和R 語言中的auto.arima 函數(shù),我們得到1個p,q 的最優(yōu)組合,即ARIMA(1,1,0)。劑量預測模型的擬合預測見圖2。
圖2 劑量預測模型的擬合預測
通過對預測模型的殘差ACF 和PACF 值的觀察,發(fā)現(xiàn)取值均在置信區(qū)間內(nèi),因此殘差檢驗可通過。預測模型的殘差ACF 和PACF 值見圖3;預測模型的參數(shù)估計見表1,由表可知,ARIMA(1,1,0)模型的參數(shù)在統(tǒng)計學上顯著。
表1 ARIMA 模型的參數(shù)
圖3 預測模型的殘差ACF 和PACF 值
經(jīng)過ARIMA 模型參數(shù)的選取、檢驗和確認,發(fā)現(xiàn)運用此模型獲得的在1個預測周期內(nèi)的預測值能夠較好地反映醫(yī)用直線加速器劑量變化的范圍和趨勢,預測結果見圖4和表2。表2是周期內(nèi)第20~24次的劑量預測值與實測值的對比,通過對比可知,預測值對實測值有較好的體現(xiàn),預測誤差為0.18%~0.80%。
表2 劑量預測數(shù)據(jù)與實際數(shù)據(jù)的對比
圖4 劑量預測值
在對數(shù)據(jù)結果進行總結時,我們發(fā)現(xiàn)預測是基于1個固定周期內(nèi)的數(shù)據(jù)進行的,為了明確ARIMA模型在其他周期內(nèi)的數(shù)據(jù)表現(xiàn),本研究就此增加了兩個測量周期的數(shù)據(jù)并用上述方法進行了預測和對比,發(fā)現(xiàn)ARIMA(1,1,0)模型對其他測量周期的數(shù)據(jù)預測也具有較好的表現(xiàn)。周期二、三的預測結果見圖5和表3。表3中周期二、三分別代表的是兩個測量周期第20~24次的實測值和預測值,由表可知,運用此模型兩組預測值與實測值的對比表現(xiàn)良好,誤差為-1.28%~0.50%。
表3 周期二、三劑量預測數(shù)據(jù)與實際數(shù)據(jù)的對比
本研究提出的利用時間序列進行醫(yī)用直線加速器劑量預測方法,具有較好的預測精度,對加速器劑量參數(shù)的及時調(diào)整起到了參考和提示作用。但是,ARIMA 模型也存在一定的局限性,僅適用于正常檢測和調(diào)整內(nèi)的數(shù)據(jù)預測,當醫(yī)用直線加速器出現(xiàn)維修和保養(yǎng),且維修和保養(yǎng)的項目影響到劑量輸出時,此模型將不再適用,檢測周期應以維修和保養(yǎng)后的時間開始計算。