潘瑩麗，黃 河

（1.湖北大學a.數(shù)學與統(tǒng)計學學院；b.應用數(shù)學湖北省重點實驗室，武漢 430062；2.梧州學院 管理學院，廣西 梧州 543002)

0 引言

基于右刪失生存數(shù)據(jù)分析協(xié)變量與生存時間的相關性時，應用最廣泛的統(tǒng)計模型之一是加速失效時間模型（AFT模型）[1—5]，然而相關研究大多是基于存儲到磁盤的靜態(tài)數(shù)據(jù)進行的。由于現(xiàn)代很多功能強大的計算平臺的應用廣泛，有關生存分析中流數(shù)據(jù)的統(tǒng)計分析在大數(shù)據(jù)領域頻繁出現(xiàn)。生存分析中的流數(shù)據(jù)又稱為實時大數(shù)據(jù)，是一組連續(xù)的、海量的、快速的、連續(xù)到達的數(shù)據(jù)序列，其突出特點是“實時到達、速率多變、連續(xù)到達、次序獨立；規(guī)模宏大，有不可預測的極值；數(shù)據(jù)一旦被處理，除非經(jīng)過特別保存，否則不能取出來再次使用，或者再次提取數(shù)據(jù)的成本很高”。由于這些獨有的特征，因此依靠傳統(tǒng)的數(shù)據(jù)挖掘技術(shù)和方法很難實時獲取流數(shù)據(jù)中承載的重要信息和進行復雜的計算。基于流數(shù)據(jù)集的統(tǒng)計分析方法有很多，當優(yōu)化問題具有封閉形式的表達式時，Schifano等（2016）[6]提出可以直接通過對數(shù)據(jù)集的遞歸過程來更新估計量。封閉形式的表達式在實際應用中并不經(jīng)常存在。Robbins和Monro（1951）[7]提出一種梯度下降算法（SGD），然而SGD算法對學習速率非常敏感。為了解決這一問題，Toulis 等（2014）[8]提出了隱式SGD 算法和平均隱式SGD 算法。此外，一些學者還提出了在線二階算法來處理流數(shù)據(jù)[9]。為了保證估計量的一致性，基于SGD的方法和在線二階算法需要對數(shù)據(jù)集進行較強的約束假設。為了避免約束假設，在可再生廣義線性回歸的研究中，Luo和Song（2020）[10]采用構(gòu)造在線更新估計方程的方法，實現(xiàn)了對模型的在線統(tǒng)計推斷。在可再生分位數(shù)回歸研究中，Wang等（2022）[11]通過構(gòu)造在線更新似然函數(shù)，僅使用前期數(shù)據(jù)的統(tǒng)計量和當前批數(shù)據(jù)對提出的估計量實現(xiàn)了在線更新；在可再生參數(shù)、半?yún)?shù)和非參數(shù)模型研究中，Lin等（2020）[12]建立了針對流數(shù)據(jù)集的在線更新估計方法。然而，這些研究全部基于完全觀測的流數(shù)據(jù)。

對于流數(shù)據(jù)分析的架構(gòu)，Warren和Marz（2015）[13]提出了一種Lambda 架構(gòu)。Lambda 架構(gòu)是一個實時的大數(shù)據(jù)計算和存儲系統(tǒng)，由速度層、批處理層和服務層組成。通過使用增量算法在速度層捕獲瞬時的實時視圖，并且更新先前存儲的視圖，而批處理層存儲不斷增長的數(shù)據(jù)，并在新的數(shù)據(jù)批次到達時不斷重新計算批處理視圖，從而細化在速度層產(chǎn)生的無法保證估計準確性的結(jié)果。批處理層和速度層輸出的視圖被存儲在服務層中用于查詢。這種架構(gòu)非常靈活，但是卻完全忽略了實時統(tǒng)計推斷的需要。Luo和Song（2020）[10]提出添加一個新的推理層來擴展Lambda架構(gòu)中的速度層，并將這種新架構(gòu)命名為rho架構(gòu)，從而對流數(shù)據(jù)進行統(tǒng)計推斷?；贏FT模型，本文采用構(gòu)造在線更新估計方程的思想構(gòu)造可再生估計并引用rho架構(gòu)，提出一種新的基于AFT模型的可再生估計的流數(shù)據(jù)分析策略，這種策略僅使用當前批數(shù)據(jù)和歷史批數(shù)據(jù)的匯總統(tǒng)計來更新估計量，解決了流數(shù)據(jù)集量大無法被全部保存的問題。提出的流數(shù)據(jù)分析策略，在計算上是高效的，有利于流數(shù)據(jù)的存儲和算法的加速，也不會損失任何估計效率。

1 模型和可再生估計

1.1 帶流數(shù)據(jù)集的AFT模型

其中，β=( )β1，β2，…，βp為p維待估參數(shù)，ε是誤差項。將觀測時間y1，y2，…，ynk按從小到大的順序排列，得到y(tǒng)i的次序統(tǒng)計量y(1)≤y(2)≤…≤y(nk)。設δ(1)，δ(2)，…，δ(nk)為y(1)，y(2)，…，y(nk)的相應右刪失變量，x(1)，x(2)，…,x(nk)為相應的協(xié)變量。為了處理數(shù)據(jù)刪失帶來的信息不完整性，基于第k批數(shù)據(jù)，根據(jù)Stute 和Wang（1993）[14]的研究，構(gòu)造權(quán)重ω(1)，ω(2)，…，ω(nk)，其中ω(i)可以表示為：

1.2 可再生估計

由上述推導過程知，估計量β?2和β?2*是近似等價的，通過求解式（9），初始值β?1就可以更新為β?2。從數(shù)值上講，通過牛頓-拉夫遜迭代算法可以得到β?2的值，即β?2的第r+1次迭代為：

本文繪制了可再生過程rho結(jié)構(gòu)流程圖（見圖1）。

圖1 rho結(jié)構(gòu)流程圖

通過流程圖1，設計可再生估計更新算法（見算法1）。

算法1：基于AFT模型的可再生估計算法。

（1）輸入：流數(shù)據(jù)集D1，D2，…，Db，…。

（2）輸出：β?b和V?ar(β?b)，其中b=1，2，…。

（3）初始化：β?init，φ?0=0，J?0=0p×p。

（4）forb=1，2，…do。

（5）讀取流數(shù)據(jù)集Db。

（10）釋放數(shù)據(jù)集Db。

（11）end。

（12）返回β?b和V?ar(β?b)，forb=1，2，…。

算法解讀如下：

（1）第一行：所有流數(shù)據(jù)集都來自真實參數(shù)為β0的AFT模型。

（2）第二行：輸出數(shù)據(jù)包括每個批次b的β的可再生估計值和估計的漸近方差。

（3）第三行：為β0設置初始值，比如β?init=0。

（4）第四行：根據(jù)數(shù)據(jù)流運行在線更新過程。

（5）第六行：在推理層，計算負Hessian 矩陣Jb(Db;β?b-1)并與速度層通信。

（6）第七行：在速度層，運行更新算法，將β?b-1更新為β?b。

（7）第八行：在推理層，根據(jù)新的數(shù)據(jù)批次Db和加速層更新后的β?b，更新負Hessian矩陣J?b和估計量φ?b。

2 數(shù)值模擬

本文通過數(shù)值模擬評估帶流數(shù)據(jù)集的AFT 模型的可再生估計方法在有限樣本中的性能，并將其與基于整個數(shù)據(jù)集計算的Oracle 方法進行比較。采用如下四個指標來評估不同方法的性能：參數(shù)估計的偏差Bias（即參數(shù)估計值的均值與參數(shù)真實值的差）、估計值的樣本標準差SD、估計值的標準誤差SE和95%正態(tài)置信區(qū)間的經(jīng)驗覆蓋率CP。對于第k（k=1，2，…，b）批數(shù)據(jù)中的第i個個體，假設生存時間Ti由以下模型生成：

2.1 固定總樣本量Nb，變化每批樣本量nk

固定樣本量Nb=100000，每批樣本量nk分別取500、800、1000 和2000，表1 匯總了(β1，β2)=(-0.500，0.693)的模擬結(jié)果，表2 匯總了(β1，β2)=(0.500，-0.693)的模擬結(jié)果。

表1 固定Nb=100000，變化nk 時的模擬結(jié)果（( β1，β2 )=(-0.500，0.693)）

表2 固定Nb=100000，變化nk 時的模擬結(jié)果（β1，β2=(0.500，-0.693)）

2.2 固定每批樣本量nk，變化總樣本量Nb

固定每一批次的樣本量nk=100，總樣本量Nb分別取值1000、5000、8000和10000，表3匯總了(β1，β2)=(-0.500，0.693)的模擬結(jié)果，表4匯總了(β1，β2)=(0.500，-0.693)的模擬結(jié)果。

表3 固定nk=100，變化Nb 時的模擬結(jié)果（( β1，β2 )=(-0.500，0.693)）

表4 固定nk=100，變化Nb 時的模擬結(jié)果（( β1，β2 )=(0.500，-0.693)）

綜合比較表1至表4中的結(jié)果，可以得出以下結(jié)論：

（1）本文提出的可再生方法得到的估計值與Oracle方法可以相比較，偏差Bias接近于0，說明β1和β2的估計量是近似無偏的，而標準誤差SE和估計值的樣本標準差SD幾乎相同，說明估計的標準誤差為估計值的樣本標準差提供了較為良好的估計，此外覆蓋率CP也幾乎都在90%附近。

（2）固定總樣本量Nb=100000 保持不變，偏差Bias、估計值的樣本標準差SD 和標準誤差SE 隨著每一批次樣本量nk的增加幾乎保持不變，說明提出的可再生估計量的估計結(jié)果幾乎不受每一批次樣本量nk大小的影響。

（3）固定每一批次樣本量nk=100 保持不變，隨著批次從10 增加到100，總樣本量從1000 累積到了10000，可再生方法的標準差SD 和標準誤差SE 的值都在一定程度上有所減小，說明可再生方法的性能越來越好。此外，在這種情況下的模擬結(jié)果也證實了本文提出的更新方法對于數(shù)據(jù)批次是穩(wěn)健的。

3 實例

為了進一步證明本文提出的可再生估計方法的可操作性和有效性，將提出的方法應用于華盛頓大學Diana 等提供的單芽插條萌芽真實數(shù)據(jù)集中。該數(shù)據(jù)集的觀察期為2013—2016 年，在此期間，霞多麗（CH）和赤霞珠（CS）這兩個品種的芽從8月開始每周采樣一次，到12月第一周結(jié)束。在前三年中，每個樣品使用了18個插條，最后一年使用了30個插條。將插條放置在裝滿潮濕砂基材的鋁托盤上。將托盤置于生長室中，在24±2°C 的溫度下，光照15 小時。每兩天觀察一次插條，記錄萌芽的持續(xù)時間。該數(shù)據(jù)集包含總樣本量N=3288，將數(shù)據(jù)分成b=24 個批次來模擬一個動態(tài)過程，每個批次包含的樣本量為nk=137，這種模擬動態(tài)過程的分割策略在在線流數(shù)據(jù)的研究中被廣泛使用。

在分析中，考慮兩個潛在的協(xié)變量：品種（x1）和樣本的切割數(shù)量（x2）。本文最終關心的結(jié)果是萌芽持續(xù)的時間，其刪失率約為12.1%?？紤]如下AFT模型來評估上述兩個變量對萌芽持續(xù)時間T的影響：

log(T)=β1x1+β2x2

類似于模擬研究，將提出的可再生估計方法與Oracle方法進行比較。表5 匯總了參數(shù)估計（Est.），標準誤差（SE）和P 值（P-values）。從表5 可以看出，本文提出的方法得出協(xié)變量x1和x2對萌芽持續(xù)時間都具有顯著影響。對于協(xié)變量x2，Oracle方法和本文提出的方法均顯示對萌芽持續(xù)時間具有顯著影響。

表5 實證結(jié)果

4 結(jié)束語

本文基于帶有流數(shù)據(jù)集的AFT模型，通過構(gòu)造在線更新估計方程的方法，提出僅基于前期數(shù)據(jù)統(tǒng)計量和當前批數(shù)據(jù)的可再生估計，研究了一種基于AFT模型的流數(shù)據(jù)分析策略，有效解決了計算機存儲瓶頸的問題。模擬研究展示了所提出方法在有限樣本量下的優(yōu)良表現(xiàn)；實際數(shù)據(jù)分析展示了所提出方法的有效性與實用性。當每批右刪失流數(shù)據(jù)的樣本量非常大時，研究的主要預算和成本來自昂貴的重要協(xié)變量的測量，在有限的時間或預算下，觀測每批研究對象的昂貴協(xié)變量往往不可行。因此，基于AFT模型，未來繼續(xù)研究Case-Cohort 設計[15]、廣義的Case-Cohort設計[16]、Outcome-Dependent Sampling（ODS）設計[17]等有偏抽樣設計下右刪失流數(shù)據(jù)的建模、統(tǒng)計推斷、實時計算等將非常有意義。