李拂曉, 肖燕婷, 陳占壽

(1.西安理工大學理學院, 陜西 西安 710054;2.青海師范大學數(shù)學與統(tǒng)計學院,青海 西寧 810008)

1.引言

分類變量在現(xiàn)實生活中十分常見, 如每個學生的期末考試成績按分數(shù)段可分為“優(yōu)秀”,“良好”, “及格”和“不及格”四類; 每個人的血壓值可分為低血壓, 正常和高血壓三類.按照時間順序記錄分類變量值, 即構成分類時間序列.關于分類時間序列的建模, 研究者提出了多種模型, 如Markov鏈模型, 廣義線性模型, 整數(shù)自回歸模型和自回歸滑動平均模型等.廣義線性模型可以把響應變量的條件期望表示為過去觀測值和協(xié)變量的函數(shù), 且相比Markov鏈模型, 既不需要假設馬氏性, 也不需要假設平穩(wěn)性, 是分類時間序列建模比較常用的模型[1].累積Logistic回歸模型(Cumulative logistic regression model)是一類重要的廣義線性模型, 在生物醫(yī)學, 社會科學和遺傳學等領域有著廣泛的應用.

在數(shù)據(jù)分析中，數(shù)據(jù)所服從的模型由于受到某些因素的影響, 可能在某個時刻或位置發(fā)生改變, 即產(chǎn)生變點.過去的幾十年中，有許多學者研究了模型的結構變點檢驗和估計問題.[2-9]變點檢驗一般分為回顧性檢驗(Retrospective test)和序貫檢驗(Sequential test), 前者是分析已觀測的歷史數(shù)據(jù), 檢驗數(shù)據(jù)中是否存在變點; 后者是連續(xù)檢驗數(shù)據(jù)中是否存在變點, 即新觀測的數(shù)據(jù)是否仍符合原有的模型,也稱為在線監(jiān)測[10].針對回顧性變點檢驗,Antoch[11]基于整體最大值型(Overall maximum type)統(tǒng)計量研究了廣義線性模型的結構變點檢驗問題.Hudecova[12]利用得分型統(tǒng)計量研究了二分類自回歸模型(Binary autoregressive model)的結構變點檢驗問題.Fokianos[13]等基于偏似然得分過程構造統(tǒng)計量, 檢驗二分類Logistic回歸模型的系數(shù)是否存在變點, Gombay[14]等將其推廣到多分類Logistic回歸模型.WANG等[15]在高維同質性檢驗的基:上提出了一種新的方法, 研究了具有大量類別的多分類數(shù)據(jù)的變點估計問題.關于變點的序貫檢驗, XIA等[16]分別基于殘差的累計和(Cumulative sum, CUSUM)和滑動和(Moving sum, MOCUSM)構造統(tǒng)計量, 在線監(jiān)測廣義線性模型的結構是否存在變點.Hohle[17]利用CUSUM方法在線監(jiān)測分類時間序列所服從的Logistic回歸模型是否存在結構變化.LI等[18]基于偏最大似然得分過程構造統(tǒng)計量, 研究了多分類Logistic回歸模型的變點序貫檢驗問題.

本文基于偏似然得分過程構造序貫檢驗統(tǒng)計量, 在線監(jiān)測累積Logistic回歸模型的結構是否存在變點.原假設下推導出統(tǒng)計量的漸近分布, 備擇假設下證明其一致性.與文[17]提出的CUSUM變點監(jiān)測方法相比, 本文所提出的方法只需要在原假設下估計出模型參數(shù), 備擇假設下不需要重新估計.模擬試驗表明原假設下統(tǒng)計量的檢驗水平大部分都接近于顯著性水平0.05, 在備擇假設下雖然監(jiān)測到變點的平均運行長度較長, 但監(jiān)測到變點的誤報率較低, 且檢驗勢也較高.最后通過一組實際數(shù)據(jù)說明本文方法的有效性.

2.模型及變點檢驗問題

變點的序貫檢驗統(tǒng)計量一般是由監(jiān)測統(tǒng)計量Γ(m,k)和邊界函數(shù)g(m,k)構成.定義停時

Γ(m,k)和g(m,k)在原假設下滿足

其中α ∈(0,1).

3.變點的序貫檢驗

首先利用偏最大似然方法對θ進行估計[20], 偏似然函數(shù)為

則對數(shù)偏似然函數(shù)為

記

其中

假設1和假設2表明偏對數(shù)似然函數(shù)的二階導數(shù)是關于θ的連續(xù)函數(shù)[21].假設2保證了Ut(θ)是非奇異的.假設3表明

知了在樹上放肆地叫著，夕陽的余暉照過來，我看到西天堆起大山一樣的烏云，一重連一重。一注一注的陽光從烏云的后面射向天空，給綿綿云山綴上了一帶炫目的金邊。我強撐著爬起來，漫無目的朝前走著。身后不知么事時候，冒出幾個細伢兒跟著，他們齊聲唱道——

滿足正定性的概率是1.[21]若假設4成立, 則當m →∞時

其中G(θ)為正定矩陣.[21]

令Sm(θ)=0, 求解可得θ的偏最大似然估計量基于構造變點序貫檢驗統(tǒng)計量

g(s)是邊界函數(shù).

下面兩個定理給出統(tǒng)計量Stat在原假設下的極限分布以及備擇假設下的一致性.

定理3.1如果假設1-4和原假設H0成立, 則

(i) 若g(s)=c, 則

由Nκm」, 類似于文[13]中定理1和定理3的證明可得

則定理的結論成立.(ii)的證明與(i)類似.

根據(jù)假設1-4,

因δ0, 當m →∞時

則(i)的結論成立, 當m+1＜k ＜k*時, 證明是類似的.(ii)的結論可由(i)直接推出.

4.模擬試驗

本節(jié)利用Monte Carlo模擬研究變點序貫檢驗統(tǒng)計量Stat的有限樣本性質.在推導統(tǒng)計量的漸近分布時, 采用了兩種不同的邊界函數(shù)g(s)=c和g(s)=cg1(s)=cs0.8, 相應的臨界值可分別通過下式獲得：

其中α是顯著性水平.數(shù)據(jù)生成過程如下：

其中{Yt}是類別數(shù)為3的累積logistic回歸模型,α1=-0.5,α2=0.2,(φ1,φ2,φ3)T=(2,-0.5,1)T.假設歷史樣本量m= 100,200,500, 監(jiān)測樣本量N1=N -m= (κ-1)m,κ= 3,5,7.顯著性水平取α=0.05, 模擬試驗重復的次數(shù)為2500 次.

令θ= (α1,α2,φ1,φ2,φ3)T, 在現(xiàn)實應用中可能只關注部分參數(shù)是否發(fā)生變化.本文假設只監(jiān)測參數(shù)α1和φ1是否發(fā)生變化, 其余參數(shù)作為冗余參數(shù).表1給出了原假設下統(tǒng)計量Stat在兩種邊界函數(shù)下的檢驗水平.結果表明兩種情形下的檢驗水平大部分都接近于顯著性水平0.05, 少部分存在一些扭曲.

表1 兩種邊界函數(shù)下統(tǒng)計量的檢驗水平

考慮下述形式的備擇假設:

HA:α1在m+k*處從-0.5變?yōu)?2,φ1在m+k*處從2變?yōu)?,

其中k*= 0.01N,0.05N,0.1N,0.2N,0.3N.表2-5描述了統(tǒng)計量Stat的檢驗勢和監(jiān)測到變點運行長度(Run length)的四分之一分位數(shù)(Q1), 中位數(shù)(Q2), 均值(Average run length, ARL), 四分之三分位數(shù)(Q3), 最大值(Max) 以及誤報率(Pτ).運行長度是變點發(fā)生位置與監(jiān)測到變點的位置之間的樣本量, 其均值稱為平均運行長度(ARL).檢驗勢是監(jiān)測到變點發(fā)生在m+k*之后,N之前的次數(shù)所占的比率.ARL≥ˉm+k*表示監(jiān)測到變點位置發(fā)生在m+k*之后時, 其運行長度的平均值.誤報率是指監(jiān)測到變點發(fā)生在m+k*之前的次數(shù)所占的比率.做為對比, 考慮了文[17]中CUSUM變點檢測方法, 其臨界值可通過模擬平均運行長度的方式獲取.假設受控狀態(tài)下平均運行長度為370, 則相應的臨界值為2.5.CUSUM統(tǒng)計量的檢驗勢, 運行長度的數(shù)字特征和誤報率見表6.

表2 在HA下, 當m=100,N =400, g(s)=c時, Stat的檢驗勢, 運行長度的數(shù)字特征和誤報率

表3 在HA下, 當m=100,N =400, g(s)=cs0.8時, Stat的檢驗勢, 運行長度的數(shù)字特征和誤報率

表4 在HA下, 當m=200,N =800, g(s)=c時, Stat的檢驗勢, 運行長度的數(shù)字特征和誤報率

表5 在HA下, 當m=200,N =800, g(s)=cs0.8時, Stat的檢驗勢, 運行長度的數(shù)字特征和誤報率

從表2-5的結果可以看出, Stat的檢驗勢隨著樣本量的增加而增大, 但當變點發(fā)生的位置離監(jiān)測起始時刻越遠, 其檢驗勢越低.邊界函數(shù)不同, 其對應的檢驗勢也有差異.當m=100,N=400時, 邊界函數(shù)取g(s)=c的檢驗勢高于g(s)=cs0.8時的檢驗勢.但當樣本量增大,m=200,N=800時, 兩種邊界函數(shù)下的檢驗勢都趨于1.

運行長度的數(shù)字特征Q1,Q2, ARL,Q3, Max隨著樣本量的增大而增大.變點發(fā)生的位置也會對運行長度產(chǎn)生影響.當邊界函數(shù)取g(s) =c時, 運行長度受變點位置的影響較小.例如當m=100,N=400, 不同變點位置下的平均運行長度是230左右;m=200,N=800, 不同變點位置下的平均運行長度是293左右.但當g(s) =cs0.8時, 不同變點位置對運行長度的影響較大.例如m=100,N=400,k*=0.01N,0.3N時, 平均運行長度分別是156.66和259.13, 差異較大.相比g(s) =c, 邊界函數(shù)取g(s) =cs0.8時, 若變點的位置離監(jiān)測起始時刻較近時, ARL較短, 但當變點位置離監(jiān)測起始時刻較遠時, ARL較長.由于誤報率在大多數(shù)情形下都為0, 變點后的平均運行長度(ARL≥ˉm+k*) 與平均運行長度(ARL)幾乎都是相等的.

表6 在HA下, 當m=100,N =400, CUSUM變點檢測方法的檢驗勢, 運行長度的數(shù)字特征和誤報率

由表6可以看出, 當利用CUSUM方法監(jiān)測變點時, 其檢驗勢, 運行長度和誤報率受變點位置的影響較大.例如m= 100,N= 400, 當變點發(fā)生在監(jiān)測起始時刻時, 檢驗勢接近于1, 平均運行長度為11.5, 誤報率為0; 但當變點的位置離監(jiān)測起始時刻較遠時, 檢驗勢迅速下降, 誤報率迅速增加, 此時監(jiān)測到的變點大都發(fā)生在m+k*之前, 所以其平均運行長度為負值.

總的來說, 統(tǒng)計量Stat在兩種邊界函數(shù)下的平均運行長度較長, 但誤報率很小, 大多數(shù)情況下都為0.CUSUM方法在變點發(fā)生在監(jiān)測起始時刻時, ARL較短, 誤報率較低.但當變點發(fā)生在離起始時刻較遠時, 誤報率較高.

5.實例分析

圖1是一組新生兒睡眠狀態(tài)的觀測序列, 每隔30秒觀測一次, 共1000個觀測值.每個時刻的睡眠狀態(tài)一般可分為四類: 安靜睡眠(Quiet sleep), 不確定睡眠(Indeterminate sleep), 主動睡眠(Active sleep), 醒(Awake), 分別記為1, 2, 3, 4.根據(jù)新生兒的睡眠規(guī)律, 四個睡眠狀態(tài)有一定的次序, 即“4”＜“1”＜“2”＜“3”, 則觀測序列可看作有序分類時間序列.文[20]利用累積logistic回歸模型對其進行建模

其中α1=-14.722,α2=-10.389,α3=-4.078,β1= 18.663,β2= 12.173,β3= 7.566.采用文[13]的回顧性變點檢驗(Retrospective change-point detection)方法檢測上述模型是否存在結構變點, 結果發(fā)現(xiàn)α2在第596個觀測值處存在變化.在變點前后對模型參數(shù)重新進行估計, 得到調整模型.對比原模型和調整后模型的AIC值, 發(fā)現(xiàn)調整后的模型AIC值減小, 即調整后的模型更加合理.用本文所提出的變點序貫檢驗監(jiān)測變點, 以前200個數(shù)據(jù)為歷史數(shù)據(jù), 從第201個數(shù)據(jù)開始監(jiān)測, 發(fā)現(xiàn)統(tǒng)計量在第671個數(shù)據(jù)處停止監(jiān)測, 監(jiān)測延遲為175.雖然有一定的監(jiān)測延遲, 但也能有效監(jiān)測到變點.

圖1 新生兒睡眠狀態(tài)的觀測序列

6.結論

累積Logistic回歸模型是分類時間序列建模的一類重要模型.針對累積Logistic回歸模型結構變點的序貫檢驗問題, 本文基于偏似然得分過程構造統(tǒng)計量, 在線監(jiān)測模型的結構是否存在變點.在原假設下推導出統(tǒng)計量的漸近分布, 備擇假設下證明其一致性.模擬試驗采用了兩種形式的邊界函數(shù),結果表明原假設下統(tǒng)計量的檢驗水平大部分都接近于顯著性水平0.05, 少部分存在一些扭曲.在備擇假設下雖然監(jiān)測到變點的平均運行長度較長, 但監(jiān)測到變點的誤報率較低, 且檢驗勢也較高.CUSUM變點檢測方法在變點發(fā)生在監(jiān)測起始時刻時, ARL較短, 誤報率較低.但當變點發(fā)生在離起始時刻較遠時, 誤報率較高.最后通過一組新生兒睡眠狀態(tài)的觀測數(shù)據(jù)說明本文方法的有效性.