趙 志 周 倩 張晉昕△

·綜述·

時間序列分析方法及其進(jìn)展*

趙 志1周 倩2張晉昕1△

在醫(yī)學(xué)科研工作中，按某種(相等或不相等的)時間間隔對客觀事物進(jìn)行動態(tài)觀察，由于隨機(jī)因素的影響，各次觀察的指標(biāo)X1,X2,X3,…,Xi，…都是隨機(jī)變量，這種按時間順序排列的一系列隨機(jī)變量(或其觀測值)稱為時間序列。

例如流行病學(xué)家會關(guān)注在某地區(qū)觀察到的流感樣病例數(shù)、登革熱病例數(shù)隨時間的變化，時間序列模型可以考慮時間對病例數(shù)的影響，亦可納入不同的流行病學(xué)影響因素來預(yù)測病例數(shù)，以及探討疫情季節(jié)性特征。臨床中連續(xù)時間內(nèi)多次血壓的測量能夠評估藥物治療高血壓的效果；人腦的功能性核磁共振造影(fMRI)中能夠采集一系列時間序列，據(jù)以研究不同實驗條件下大腦對于刺激的反應(yīng)[1]。生物學(xué)家也會對基因表達(dá)譜中隱含的一些重要模式規(guī)律感興趣，以期與相應(yīng)的表觀或疾病關(guān)聯(lián)[2]。時間序列分析方法是處理這些問題的有力工具。本文將對時間序列分析的經(jīng)典模型及最新研究進(jìn)展進(jìn)行綜述。

時間序列的基本特征

1.平穩(wěn)性

時間序列是一種特殊的隨機(jī)過程，根據(jù)觀測記錄對隨機(jī)過程的結(jié)構(gòu)進(jìn)行統(tǒng)計推斷時，通常需要作出某些假設(shè)，其中最重要的假設(shè)是平穩(wěn)性，即決定過程特性的統(tǒng)計規(guī)律不隨著時間的變化而變化。從一定意義上講，過程位于統(tǒng)計的平衡點(diǎn)上。特別地，對一切時滯k和時點(diǎn)t1,t2,…,tn，都有Xt1,Xt2,…,Xtn與Xt1-k,Xt2-k,…,Xtn-k的聯(lián)合分布相同，則稱過程{Xt}為嚴(yán)平穩(wěn)過程[3](strictly stationary process)。

判斷一個隨機(jī)過程是否是嚴(yán)平穩(wěn)的，這在實際應(yīng)用中十分困難。因為它要求過程的各階矩都是時不變的(time-invariant)，即只與時間間隔有關(guān)，與初始時間無關(guān)，并且各階矩可以無窮大。而實際問題中，如果數(shù)據(jù)的方差無窮大，意味著變異無限大，那么就無法對數(shù)據(jù)進(jìn)行建模來預(yù)測或控制。因此，當(dāng)把條件適當(dāng)放寬，只限定存在常均值、有限且時不變的二階矩，便得到了寬平穩(wěn)過程(weakly stationary process)。若隨機(jī)過程的均值(一階距)和協(xié)方差(二階矩)存在，且滿足：

(1)E[X(t)]=μ,?t∈T;

(2)E[Xt+τ-μ][Xt-μ]=γ(τ),?t,t+τ∈T;

則稱{X(t),t∈T}為寬平穩(wěn)過程，均值μ為常數(shù)，γ(τ)是X(t)的協(xié)方差函數(shù)，它與時間t無關(guān)。

判斷時間序列是否平穩(wěn)是建模的前提，一般可通過繪制時序圖來判斷。若時序圖中的曲線圍繞一定均值呈現(xiàn)波動變化相似的情景，則可判定序列是平穩(wěn)的。更準(zhǔn)確地，可以通過單位根檢驗判斷平穩(wěn)性。

2.可逆性

可以假設(shè)隨機(jī)序列是隨機(jī)沖擊的線性組合，這種隨機(jī)序列就是用一般線性隨機(jī)模型來描述的[4]。例如對于隨機(jī)過程{Xt}有

Xt=εt-θεt-1

(1)

若將εt用現(xiàn)在和過去的Xt表示，則(1)式變?yōu)?/p>

εt=(1-θB)-1Xt

(2)

其中，B表示滯后算子，定義為BmXt=Xt-m。進(jìn)一步，若|θ|<1，(2)式可寫成無窮級數(shù)

(3)

這時就稱序列{Xt}是可逆的。

然而，如果|θ|>1，式(2)可寫成

(4)

即只能用過程的現(xiàn)在和將來值來刻畫當(dāng)前指標(biāo)值，這樣的模型不符合常理，被稱作不滿足可逆性條件。

時域分析方法的主流模型

常用的時間序列模型有很多，如移動平均模型、條件異方差模型、非線性時間序列模型等。這些模型是從時域(time domain)角度對序列的將來值用過去值建模，要求相鄰時間點(diǎn)序列的相關(guān)性能夠被過去值很好地刻畫。

1.ARIMA模型

ARIMA模型(autoregressive integrated moving average model)又稱博克斯-詹金斯(Box-Jenkins)模型[5]，是由美國統(tǒng)計學(xué)家G.E.P.Box和G.M.Jenkins于1970年首次系統(tǒng)提出的，該模型有三種基本模式：自回歸模型、移動平均模型與自回歸移動平均模型。當(dāng)序列為平穩(wěn)序列，其形式為ARMA模型：

Xt=φ1Xt-1+φ2Xt-2+…+φpXt-p+εt-θ1εt-1-

θ2εt-2-…-θqεt-q

(5)

當(dāng)序列為非平穩(wěn)序列，但通過d次差分可使序列平穩(wěn)時，采用的模型稱作自回歸求和移動平均模型，即ARIMA(p,d,q)模型。鑒于三種基本模式可視作ARIMA模型的特例，故ARIMA模型又常被作為這一族模型的總稱。對于判斷此模型是否平穩(wěn)，可通過DF(Dickey-Fuller法)或ADF(augmented Dickey-Fuller法)單位根檢驗，并經(jīng)過考察θ(B)的根符合“|θ|<1”判定序列滿足可逆性。

ARIMA模型能夠為傳染病或相關(guān)醫(yī)療衛(wèi)生設(shè)施方面的預(yù)測工作提供有效的指導(dǎo)。例如，2003年SARS爆發(fā)期間，Earnest和Chen等[6]利用ARIMA模型對新加坡的一家?？漆t(yī)院每日病床占用數(shù)進(jìn)行建模，所建立的ARMA(1,3)模型能夠較好地預(yù)測未來3天的病床需求數(shù)。

2.GARCH模型

經(jīng)典的ARIMA模型假定過程的當(dāng)前值依賴于過去值、當(dāng)前擾動和過去擾動，而擾動項εt為白噪聲。但現(xiàn)實中擾動不一定是理想的白噪聲，當(dāng)前的擾動可能會受到前一期甚至前幾期序列蘊(yùn)涵的信息影響，ARCH模型[7](autoregressive conditional heteroskedasticity model)由此而來。ARCH(q) 定義為：

Xt=φ1Xt-1+εt

這是對ARCH模型的推廣，稱之為GARCH(p,q)模型(generalizedautoregressiveconditionalheteroskedasticitymodel)。

條件異方差模型可看作是利用現(xiàn)在和過去的數(shù)據(jù)，基于AR模型對條件方差過程進(jìn)行建模，通過該模型可以預(yù)測序列未來的波動大小。

Moran和Solomon[8]使用GARCH模型對澳大利亞和新西蘭1995-2009年間重癥監(jiān)護(hù)室(ICU)成人月死亡數(shù)據(jù)進(jìn)行建模，從而對其實現(xiàn)有效的監(jiān)控。Modarres和Ouarda等[9]對加拿大蒙特利爾地區(qū)的氣候條件和髖部骨折率建立了多元GARCH模型，文章指出雪的深度、氣溫、氣壓以及白晝長度與髖部骨折率存在較大的關(guān)聯(lián)。

3.非線性時間序列

條件異方差模型在序列與擾動項之間建立了一種非線性的關(guān)系，同樣，也可以在序列之間建立非線性的函數(shù)關(guān)系，如

Xt=f1(Xt-1)+…+fp(Xt-p)+εt

(6)

其中，fi(Xt-i)可以是已知的非線性函數(shù)，i=1,2,…p，可取fi(Xt-i)=1/Xt-i。當(dāng)fi(Xt-i)為未知的函數(shù)時，模型(6)稱為可加自回歸模型(additiveautoregressivemodel)[10]，表示為Xt～AAP(p)，這類非線性關(guān)系未知的模型統(tǒng)稱為非參數(shù)時間序列模型。

另一種使用較為廣泛的非線性時間序列模型是門限自回歸模型(thresholdautoregressivemodel)。例如流感所引起的病死人數(shù)的增加速度往往比減少速度慢，并且當(dāng)病死人數(shù)達(dá)某一值時，過程會從一個方程(模式)突然切換到另一個方程。此時模型可設(shè)定為

其中ε1t和ε2t是白噪聲，d表示模式改變的時刻，稱之為延遲參數(shù)，c為臨界值。

頻域分析

頻域分析的思想認(rèn)為所有時間序列都可以看作是不同周期成分疊加的結(jié)果。周期成分廣泛存在于生物醫(yī)學(xué)的時間序列數(shù)據(jù)中，大到人群疾病流行強(qiáng)度的周期波動，小到細(xì)胞新陳代謝的周期生長，如生物醫(yī)學(xué)信號處理中的心電圖、腦電圖、醫(yī)院月度門診量等[11-12]。具有明顯周期成分的時間序列一般會在時域圖中顯示出周期性質(zhì)來，但更豐富的周期信息常常蘊(yùn)含于序列內(nèi)部，通過肉眼對時域圖的判讀難以發(fā)現(xiàn)并量化它，需要通過特定的方法將這種周期信息提取出來。早在 1929 年統(tǒng)計學(xué)家 R.A.Fisher就對時間序列周期性檢驗進(jìn)行過研究，他運(yùn)用傅立葉變換獲得時間序列周期圖并提出基于周期圖法的 Fisherg統(tǒng)計量用于檢測峰值[13]。傅立葉變換是頻域分析的基本工具。之后統(tǒng)計學(xué)家對長短不同、背景噪聲大小不同的周期性檢驗方法不斷予以改進(jìn)[14]。

定性時間序列

定性時間序列(categorical time series) 又稱分類時間序列，是指每個時間點(diǎn)觀測值的取值范圍為有限狀態(tài)空間C={C1，C2，…，Cn}的時間序列，其取值只能表示狀態(tài)或者類別。定性時間序列廣泛存在于各個領(lǐng)域，如生物醫(yī)學(xué)、行為學(xué)、流行病學(xué)、遺傳學(xué)等[15]。DNA序列、睡眠狀態(tài)監(jiān)測就是由有限狀態(tài)空間形成的定性時間序列[16-17]。

定性時間序列的分析方法早期主要集中在時域分析部分，常見的有馬爾科夫鏈分析法、連接函數(shù)法等[17]。隨著定性時間序列的周期性研究的應(yīng)用開展，頻域分析方法也開始出現(xiàn)并得到改進(jìn)。包括基于變換的功率譜分析法、基于序列啞元化的譜封分析法(spectral envelope analysis)[18]、基于非平穩(wěn)序列的小波分析法[19]等。其中，譜封分析法適用于各種類型時間序列的周期性分析，它是定性時間序列周期性分析中檢驗效果較好的一種方法。

多元時間序列

研究某一地區(qū)由于上呼吸道感染、肺炎、哮喘所引起的就診人數(shù)時間序列，如果出現(xiàn)異常的天氣條件，各類疾病的門診量時間序列就會同時發(fā)生改變，建模中不能忽視可能存在的關(guān)聯(lián)[20]。

關(guān)于一元時間序列的大部分基本理論和方法都可以推廣到多元時間序列，最常用的多元時間序列模型是向量自回歸模型(vector autoregressive model,VAR)?？紤]p階的多元時間序列Xt滿足VAR(p) 模型[21-22]，則

其中，Xt可以是上述例子中門診量的3維列向量，φ0是3維常數(shù)列向量，而φi是3×3維的系數(shù)矩陣，εt是3維的白噪聲過程。

多元情況下也存在向量ARMA模型，但其參數(shù)過多，往往存在模型的可識別性問題。針對非平穩(wěn)多元時間序列，一般需要檢驗序列的線性協(xié)整(cointegration)并建立誤差修正模型(error correction model)。另外，多元時間序列中還存在一種重要的模型結(jié)構(gòu)關(guān)系——格蘭杰因果關(guān)系(Granger causality)，它對于研究不同序列之間的繼起關(guān)系具有重要價值。

研究熱點(diǎn)與展望

時間序列模型的研究是為了不斷適應(yīng)實踐中遇到的新問題，建立更合理有效的模型。近年來在多元時間序列方面研究較多，該領(lǐng)域主要的發(fā)展動向是：(1)由一元非線性模型轉(zhuǎn)向多元非線性模型；(2)非平穩(wěn)的多元時間序列與小波分析等工具的結(jié)合；(3)高維多元時間序列的降維問題；(4)時間序列與統(tǒng)計過程控制相結(jié)合，可用于癥狀監(jiān)測，對于檢出疾病暴發(fā)等異常變動具有重要價值。

一般的研究對象往往處在紛繁復(fù)雜的系統(tǒng)中。例如傳染病發(fā)病人數(shù)的變化很可能與氣候因素、醫(yī)療衛(wèi)生條件等密切相關(guān)，這些序列各自都可能呈現(xiàn)非平穩(wěn)且非線性的特性，普通的VAR模型或非參數(shù)時間序列模型也不再適用。此時，建立起合理有效的多元非線性模型引起研究者的廣泛關(guān)注，主要研究角度是從非參數(shù)協(xié)整和半?yún)?shù)協(xié)整[23-26]對多元時間序列的非平穩(wěn)性進(jìn)行修正。它們對于非平穩(wěn)多元時間序列的推廣應(yīng)用具有重要意義。

針對非平穩(wěn)的多元時間序列，另一研究方向是從頻域角度出發(fā)與小波分析相結(jié)合[27]。小波分析不僅對非平穩(wěn)性具有很好的容忍性，還能夠同時結(jié)合時域分析與頻域分析的優(yōu)勢，使過程的信息得到更全面的利用。另外，還可以將小波變換與信息論相結(jié)合，利用小波熵來提取序列蘊(yùn)涵的信息[28-29]。

醫(yī)學(xué)研究中還存在一種具有很多變量或特征的時間序列資料，如fMRI影像數(shù)據(jù)中，每個像素就是一個變量或特征。對于高維的多元時間序列，為了便于統(tǒng)計分析以及實際意義的解釋，往往需要進(jìn)行降維處理，主成分分析[30]、因子模型[31-32]、LASSO[33]等降維方法在高維多元時間序列中使用較為廣泛。例如對于fMRI數(shù)據(jù)使用主成分分析和格蘭杰因果分析[34]，可以探討不同腦區(qū)之間的功能性聯(lián)系，為臨床診斷和治療提供參考。公共衛(wèi)生領(lǐng)域?qū)W者可以利用時間序列模型方法進(jìn)行疾病的長期趨勢預(yù)測與疫情監(jiān)測，歸納出疾病演變的規(guī)律；借助統(tǒng)計過程控制方法，提高監(jiān)測的準(zhǔn)確性[35]。目前的熱點(diǎn)是變點(diǎn)分析(change-point detection)，找出樣本序列的分布或特征突然發(fā)生變化的某個或某些時間點(diǎn)，來提示可能存在的疫情爆發(fā)并加以預(yù)警[36-38]；或評價公共衛(wèi)生干預(yù)手段是否可以起效，估計起效的遲滯時間長度。時間序列分析技術(shù)為前瞻性地了解疫情、主動地控制疾病流行提供了有效支持。

[1]Shumway RH,Stoffer DS.Time Series Analysis and Its Application With R Example.2nd Edition.北京:世界圖書出版社,2011.

[2]Vivian T,Liew WC,Yan H.Periodicity analysis of DNA microarray gene expression time series profiles in mouse segmentation clock data.Statistics and Its Interface,2010,3(3):413-418.

[3]Cryer JD,Chan KS.時間序列分析及應(yīng)用：R語言.北京:機(jī)械工業(yè)出版社,2011.

[4]Box GE,Jenkins GM,Reinsel GC.時間序列分析:預(yù)測與控制.北京:機(jī)械工業(yè)出版社,2011.

[5]方積乾,陸盈.現(xiàn)代醫(yī)學(xué)統(tǒng)計學(xué).北京:人民衛(wèi)生出版社,2002.

[6]Earnest A,Chen MI,Ng D,et al.Using autoregressive integrated moving average(ARIMA) models to predict and monitor the number of beds occupied during a SARS outbreak in a tertiary hospital in Singapore.BMC Health Serv Res,2005,5(36):1-8.

[7]Enders W.Applied Econometric Time Series.New Jersey:John Wiley & Sons,2004.

[8]Moran JL,Solomon PJ.Statistical process control of mortality series in the Australian and New Zealand Intensive Care Society(ANZICS) adult patient database:implications of the data generating process.BMC Medical Res Methodol,2013,13(66):1-12.

[9]Modarres R,Ouarda TBMJ,Vanasse A,et al.Modeling climate effects on hip fracture rate by the multivariate GARCH model in Montreal region,Canada.Int J Biometeorol,2014,58(5):921-930.

[10]Fan JQ,Yao QW.Nonlinear Time Series:Nonparametric and Parametric Methods.北京:科學(xué)出版社,2006.

[11]Stevenson NJ,O’Toole JM,Rankine LJ,et al.A nonparametric feature for neonatal EEG seizure detection based on a representation of pseudo-periodicity.Medical Eng & Phys,2011,34(4):437-446.

[12]周倩,張晉昕.時間序列周期性檢驗方法研究進(jìn)展.中國衛(wèi)生統(tǒng)計,2013,30(3):445-447.

[13]Fisher RA.Tests of significance in harmonic analysis.Proc Math Phys & Eng Sci,1929,125(796):54-59.

[14]Krafty RT,Xiong S,Stoffer DS,et al.Enveloping spectral surfaces:covariate dependent spectral analysis of categorical time series.J Time Series Anal,2011,33(5):797-806.

[15]Freeman WJ,Viana DPG.Relation of olfactory EEG to behavior:time series analysis.Behav Neurosci,1986,100(5):753-763.

[16]Stoffer DS,Tyler DE,Mcdougall AJ.Spectral analysis for categorical time series:scaling and the spectral envelope.Biometrika,1993,83(3):611-622.

[17]Mcgee M,Ensor K.Tests for harmonic components in the spectra of categorical time series.J Time Series Anal,1998,19(3):309-323.

[18]Stoffer DS,Tyler DE.Matching sequences:Cross-spectral analysis of categorical time series.Biometrika,1998,85(1):201-213.

[19]Rasheed F,Alshalalfa M,Alhajj R.Adaptive machine learning technique for periodicity detection in biological sequences.Int J Neural Syst,2009,19(1):11-24.

[20]Diggle P.Time Series:An Biostatistical Introduction.New York:Oxford University Press,1990.

[21]Tsay RS.Multivariate Time Series Analysis:With R and Financial Applications.New Jersey:John Wiley & Sons,2014.

[22]倪延延,張晉昕.向量自回歸模型擬合與預(yù)測效果評價.中國衛(wèi)生統(tǒng)計,2014,31(1):53-56.

[23]Gu JP,Liang ZW.Testing cointegration relationship in a semiparametric varying coefficient model.J Econometrics,2014,178(1):57-70.

[24]Boswijk HP,Lucas A.Semi-nonparametric cointegration testing.J Econometrics,2002,108(2):253-280.

[25]Wang QY,Phillips CB.Structural Nonparametric Cointegrating Regression.Econometrics ,2009,77(6):1901-1948.

[26]Kasparis I,Phillips CB.Dynamic misspecification in nonparametric cointegrating regression.J Econometrics,2012,168(2):270-284.

[27]Maharaj EA,Alonso AM.Discriminant analysis of multivariate time series:Application to diagnosis based on ECG signals.Comput Stat & Data Anal,2014,70:67-87.

[28]Chen JK,Li GQ.Tsallis Wavelet Entropy and Its Application in Power Signal Analysis.Entropy,2014,16(6):3009-3025.

[29]Fraiwan L,Lweesy K,Khasawneh N,et al.Classification of Sleep Stages Using Multi-wavelet Time Frequency Entropy and LDA.Methods Inf Med,2010,49(3):230-237.

[30]Li H.Asynchronism-based principal component analysis for time series data mining.Expert Syst Appl,2014,41(6):2842-2850.

[31]Lam C,Yao QW.Factor modelling for high-dimensional time series:inference for the number of factors.Annals Stat,2012,40(2):694-726.

[32]Pan JZ,Yao QW.Modelling multiple time series via common factors.Biometrika,2008,95(2):365-379.

[33]Hsu NJ,Hung HL,Chang YM.Subset selection for vector autoregressive processes using Lasso.Comput Stat & Data Anal,2008,52(7):3645-3657.

[34]Zhou ZY,Chen YH,Ding MZ,et al.Analyzing Brain Networks With PCA and Conditional Granger Causality.Hum Brain Mapp ,2009,30(7):2197-2206.

[35]詹思延,李立明,吳系科.流行病學(xué)進(jìn)展.第12卷.北京:人民衛(wèi)生出版社,2010.

[36]Chen MP,Shang N,Winston CA,et al.A Bayesian analysis of the 2009 decline in tuberculosis morbidity in the United States.Stat Med,2012,31(27):3278-3284.

[37]Kass-Hout TA,Xu ZH,McMurray P,et al.Application of change point analysis to daily influenza-like illness emergency department visits.J Am Med Inform Assoc,2012,19(6):1075-1081.

[38]Riffenburgh RH,Cummins KM.A simple and general change-point identifier.Stat Med,2006,25(6):1067-1077.

(責(zé)任編輯：郭海強(qiáng))

國家自然科學(xué)基金(30872182)

1.中山大學(xué)公共衛(wèi)生學(xué)院醫(yī)學(xué)統(tǒng)計與流行病學(xué)系(510080)

2.中山大學(xué)附屬第一醫(yī)院流行病學(xué)研究室

△通信作者：張晉昕，Email：zhjinx@mail.sysu.edu.cn