褪黑素時(shí)間序列中缺失值的填補(bǔ)方法研究

錢淑雯，李海燕，楊學(xué)智，李 慧，張尚尚，蘆 煜，王京平，張治霞，徐 靜

(北京中醫(yī)藥大學(xué)，北京 100029)

論 著

褪黑素時(shí)間序列中缺失值的填補(bǔ)方法研究

錢淑雯，李海燕，楊學(xué)智，李 慧，張尚尚，蘆 煜，王京平，張治霞，徐 靜

(北京中醫(yī)藥大學(xué)，北京 100029)

目的 模擬同一時(shí)間點(diǎn)數(shù)據(jù)完全缺失和部分缺失2種情況，通過填補(bǔ)值和實(shí)際值的對(duì)比，比較各填補(bǔ)方法對(duì)褪黑素(MT)時(shí)間序列的填補(bǔ)效果。方法 同一時(shí)間點(diǎn)完全缺失時(shí)，比較實(shí)際值與SPSS 5種填補(bǔ)方法填補(bǔ)結(jié)果；部分缺失時(shí)，除完全缺失的填補(bǔ)方法外，增加擬合時(shí)間序列模型填充。結(jié)果 完全缺失時(shí)，臨近點(diǎn)的中位數(shù)和線性插值法的填補(bǔ)結(jié)果和兩因素析因設(shè)計(jì)資料方差分析結(jié)果更接近于實(shí)際。實(shí)際值波動(dòng)幅度較小的時(shí)候，插值法擬合效果好；在實(shí)際值波動(dòng)較大時(shí)，臨近中位數(shù)擬合效果好。部分缺失時(shí)，擬合模型填充效果好。結(jié)論 完全缺失時(shí)，如排除缺失值大幅波動(dòng)，可以運(yùn)用臨近中位數(shù)和插值法對(duì)缺失值進(jìn)行填充。在臨近值波動(dòng)幅度較小時(shí)，選用插值法填充值；在臨近值波動(dòng)幅度較大時(shí)，選用臨近中位數(shù)填充值。部分缺失時(shí)，選用時(shí)間序列擬合模型填充。

時(shí)間序列;缺失值;填補(bǔ)方法

在中醫(yī)藥實(shí)驗(yàn)研究中，由于各種意外情況的發(fā)生，實(shí)驗(yàn)數(shù)據(jù)不可避免存在缺失的情況。在基于“天人相應(yīng)”的時(shí)間序列數(shù)據(jù)分析時(shí)，缺失值的出現(xiàn)對(duì)時(shí)間序列統(tǒng)計(jì)分析極為不利。時(shí)間序列是某一個(gè)指標(biāo)在不同時(shí)間上的不同數(shù)值按照時(shí)間的先后順序排列而成的數(shù)列。缺失值的出現(xiàn)不僅影響各時(shí)間點(diǎn)的連接，對(duì)時(shí)間序列數(shù)據(jù)，序列太短會(huì)有損參數(shù)估計(jì)的穩(wěn)健性，取較長序列則可以保證擬合模型的可靠性。而實(shí)際情形是，序列涵蓋的歷史值越長，序列中越有可能含有缺失數(shù)據(jù)，即使加強(qiáng)質(zhì)量控制，也不能完全避免。褪黑素(MT)是“天人相應(yīng)”研究中最常用的節(jié)律性指標(biāo)，有強(qiáng)大的抗氧化能力,調(diào)節(jié)生物鐘、保護(hù)小腸黏膜屏障、調(diào)節(jié)腸道細(xì)菌等廣泛的生理作用[1]。MT分泌值正是這樣一組典型的時(shí)間序列數(shù)據(jù)[2]。筆者以MT完整時(shí)間序列數(shù)為依據(jù)，創(chuàng)新地模擬同一時(shí)間點(diǎn)數(shù)據(jù)完全缺失和部分缺失2種情況，采用實(shí)際數(shù)據(jù)對(duì)比填補(bǔ)值的辦法，結(jié)合兩因素析因設(shè)計(jì)資料方差分析和時(shí)間序列圖，比較各填補(bǔ)方法對(duì)MT時(shí)間序列的填補(bǔ)效果。

1 實(shí)驗(yàn)資料

1.1 動(dòng)物及分組 192只健康SPF級(jí)昆明種(KM)雌性小鼠，體質(zhì)量14～20 g，由北京維通利華實(shí)驗(yàn)動(dòng)物技術(shù)有限公司提供。將192只雌性KM小鼠隨機(jī)分為4組，每組48只：①正常明-暗周期組(Ⅰ組)，在自然光條件下用普通飼料喂養(yǎng)；②正常明-暗周期+藥物組(Ⅱ組)，在自然光照條件下以含0.25%補(bǔ)腎方藥(杞蓉益精顆粒)的飼料喂養(yǎng)；③異常明-暗周期組(Ⅲ組)，使用電子定時(shí)器，每天18：00—24：00用日光燈光照6 h，普通飼料喂養(yǎng)；④異常明-暗周期+藥物組(Ⅳ組)，光照條件同Ⅲ組，以含0.25%補(bǔ)腎方藥(杞蓉益精顆粒)的飼料喂養(yǎng)。

1.2 取材方法 為保證造模和取材的有效性，各組實(shí)驗(yàn)動(dòng)物喂養(yǎng)30 d，造模時(shí)間相對(duì)較長。取材前4 h禁食，相近的時(shí)間點(diǎn)取材小鼠取自不同的飼養(yǎng)籠，盡量讓小鼠處于各飼養(yǎng)狀態(tài)，減少對(duì)指標(biāo)的干擾。為細(xì)化MT的節(jié)律變化，充分反映變化趨勢(shì)，取材于喂養(yǎng)結(jié)束后分別于次日2:00，4:00，6:00，8:00，10:00，12:00，14:00，16:00，18:00，20:00，22:00，24:00 12個(gè)時(shí)間點(diǎn)從小鼠眼眶后靜脈叢取血，每個(gè)時(shí)間段(24 h)取12只測(cè)量，連續(xù)取4個(gè)時(shí)間段。所取血樣經(jīng)3 500 r/min離心15～20 min分離血清，所取血清置于帶蓋的錐形小管內(nèi)于-70 ℃保存，用酶聯(lián)免疫法測(cè)定其血清MT含量。

1.3 實(shí)驗(yàn)方法

1.3.1 模擬缺失方案 文獻(xiàn)[3]研究表明，缺失比例在10%～30%時(shí)，各填補(bǔ)方法的填補(bǔ)結(jié)果較穩(wěn)定，因此模擬缺失值確定為每個(gè)24 h周期內(nèi)2個(gè)，缺失比例為16.7%，最接近穩(wěn)定比例的平均值。運(yùn)用EXCEL產(chǎn)生模擬缺失值的位置，以1～12的整數(shù)產(chǎn)生隨機(jī)數(shù)。

1.3.2 完全缺失填補(bǔ)方案 如果在取材4個(gè)時(shí)間段中，每個(gè)24 h周期的同一時(shí)間點(diǎn)取材的數(shù)據(jù)都缺失時(shí)，鑒于SPSS軟件醫(yī)學(xué)研究中應(yīng)用的普遍性，直接運(yùn)用SPSS 18.0軟件進(jìn)行缺失值填補(bǔ)。鑒于該時(shí)間點(diǎn)的數(shù)據(jù)完全缺失，擬合的時(shí)間序列周期模型無該時(shí)間點(diǎn)的實(shí)驗(yàn)數(shù)據(jù)，因此不采用擬合的時(shí)間序列周期模型填補(bǔ)缺失值。比較各方法的填補(bǔ)值和實(shí)際值，結(jié)合兩因素析因設(shè)計(jì)資料的方差分析和時(shí)間序列圖，得出最優(yōu)填補(bǔ)方法，并驗(yàn)證方法的合理性。①序列均值為取整列數(shù)據(jù)的均值；②臨近點(diǎn)均值為去該缺失值臨近的幾個(gè)點(diǎn)的均值，具體幾個(gè)點(diǎn)由附近點(diǎn)的跨度來設(shè)定；③臨近點(diǎn)的中位數(shù)為取該缺失值鄰近的幾個(gè)點(diǎn)的中位數(shù)，具體幾個(gè)點(diǎn)由附近點(diǎn)的跨度來設(shè)定；④線性插值法應(yīng)用線性插值法填補(bǔ)缺失值，用該列數(shù)據(jù)缺失值前一個(gè)數(shù)據(jù)和后一個(gè)數(shù)據(jù)建立插值直線，然后用缺失點(diǎn)在線性插值函數(shù)的函數(shù)值填充該缺失值；⑤點(diǎn)處的線性趨勢(shì)法應(yīng)用缺失值所在的整個(gè)序列建立線性回歸方程，然后用該回歸方程在缺失點(diǎn)的預(yù)測(cè)值填充缺失值。

1.3.3 部分缺失填補(bǔ)方案 在取材4個(gè)時(shí)間段中，每個(gè)24 h周期的同一時(shí)間點(diǎn)取材的數(shù)據(jù)部分缺失時(shí)，除完全缺失填補(bǔ)方法外，根據(jù)MT周期分泌規(guī)律，以現(xiàn)有4組的實(shí)驗(yàn)數(shù)據(jù)擬合的時(shí)間序列周期模型[4]，以周期模型值填充缺失值。

2 結(jié) 果

2.1 模擬缺失值位置 EXCEL隨機(jī)產(chǎn)生1～12之間的整數(shù)，結(jié)果為3和9，確定24 h周期的第3個(gè)、第9個(gè)時(shí)間點(diǎn)為模擬缺失值位置，即4:00和16:00的實(shí)驗(yàn)數(shù)據(jù)缺失。

2.2 完全缺失各方法填補(bǔ)結(jié)果

2.2.1 與實(shí)際數(shù)據(jù)對(duì)比 根據(jù)采集到的各實(shí)驗(yàn)組數(shù)據(jù)的情況，采用數(shù)據(jù)較全的Ⅳ組MT 4個(gè)24 h周期的時(shí)間段的時(shí)間序列數(shù)據(jù)，模擬在4個(gè)時(shí)間段的4:00和16:00數(shù)據(jù)都缺失。從4：00和16：00缺失數(shù)據(jù)的填補(bǔ)方法結(jié)果可以看出，臨近點(diǎn)的中位數(shù)和線性插值法的填補(bǔ)值更接近于實(shí)際值。這兩個(gè)時(shí)間點(diǎn)的實(shí)際值和SPSS 5種填補(bǔ)方法的填補(bǔ)結(jié)果見表1和表2。

表1 4：00 MT缺失數(shù)據(jù)的填補(bǔ)方法結(jié)果比較

表2 16：00 MT缺失數(shù)據(jù)的填補(bǔ)方法結(jié)果比較

2.2.2 填補(bǔ)結(jié)果的兩因素析因設(shè)計(jì)資料方差分析 對(duì)Ⅳ組MT 4個(gè)時(shí)間段的數(shù)據(jù)進(jìn)行兩因素析因設(shè)計(jì)資料的方差分析，因素為額外光照和藥物，額外光照有2個(gè)水平，“0”無額外光照，“1”表示有額外光照；藥物有2個(gè)水平， “0”表示無藥物，“1”表示有藥物。分析在4:00和16:00，額外光照作用、藥物作用和光照和藥物交互作用對(duì)小鼠MT分泌是否有顯著影響。見表3及表4。實(shí)際數(shù)據(jù)的方差分析結(jié)果顯示額外光照作用對(duì)該時(shí)間點(diǎn)MT分泌無影響顯著；藥物作用對(duì)該時(shí)間點(diǎn)MT分泌有顯著影響；光照與藥物交互作用對(duì)該時(shí)間點(diǎn)MT分泌有顯著影響。臨近點(diǎn)的中位數(shù)和線性插值法的P值和顯著性結(jié)果更接近于實(shí)際，但是存在與實(shí)際結(jié)果不一致的地方。因此，利用Ⅰ～Ⅲ組中較全的數(shù)據(jù)，進(jìn)一步分析臨近點(diǎn)的中位數(shù)和線性插值法適用的填補(bǔ)情形。見表3和表4。

2.2.3 各組填補(bǔ)前后的時(shí)間序列圖 利用Ⅰ～Ⅲ組中較全的數(shù)據(jù)，進(jìn)一步分析臨近點(diǎn)的中位數(shù)和線性插值法適用的填補(bǔ)情形。采用時(shí)間序列圖，分析臨近點(diǎn)的中位數(shù)和線性插值法的適用范圍，對(duì)比驗(yàn)證臨近點(diǎn)的中位數(shù)和線性插值法的合理性。采用Ⅰ～Ⅲ組中一個(gè)周期內(nèi)的各組較全的數(shù)據(jù)，模擬4:00和16:00的缺失值。結(jié)果顯示，在實(shí)際值波動(dòng)幅度較小的時(shí)候，插值法擬合效果好；在實(shí)際值波動(dòng)時(shí)，臨近中位數(shù)擬合效果好。缺失值填補(bǔ)方法不適用于臨近點(diǎn)大幅度增大和減小的情況。見圖1～3。

表3 4：00實(shí)際及缺失值填充結(jié)果方差分析

表4 16：00實(shí)際及缺失值填充結(jié)果方差分析

圖1 Ⅰ組實(shí)際及缺失值填充的時(shí)間序列圖

圖2 Ⅱ組實(shí)際及缺失值填充的時(shí)間序列圖

2.3 部分缺失各方法填補(bǔ)結(jié)果 由于各實(shí)驗(yàn)組第1個(gè)24 h時(shí)間段內(nèi)MT時(shí)間序列數(shù)據(jù)較完整，模擬4組的第一時(shí)間段4:00和16:00缺失實(shí)驗(yàn)數(shù)據(jù)，擬合的時(shí)間序列周期模型，以周期模型值填充各組缺失值。從填補(bǔ)值來看，周期模型值填充效果較臨近中位數(shù)和插值法好。見表5。

圖3 Ⅲ組實(shí)際及缺失值填充的時(shí)間序列圖

組別4:00實(shí)際模型插值臨近中位16:00實(shí)際模型插值臨近中位Ⅰ組0.3150.3400.3400.3900.4670.4500.7700.620Ⅱ組0.6210.7500.7700.6801.4431.0400.5700.540Ⅲ組0.6650.6500.6900.7000.4670.4500.7700.620Ⅳ組0.6170.5900.7900.7900.3320.3900.3300.380

3 討 論

缺失按其產(chǎn)生機(jī)制可分為3類[5]:①M(fèi)CAR完全隨機(jī)缺失。此種缺失機(jī)制下，數(shù)據(jù)的缺失不依賴于要研究的因素，是隨機(jī)的。②MAR隨機(jī)缺失。隨機(jī)是指條件隨機(jī)，即在給定觀察到的數(shù)據(jù)的情況下，缺失是一種隨機(jī)現(xiàn)象，數(shù)據(jù)缺失的可能性只依賴于觀察到的數(shù)據(jù)。③MNAR非隨機(jī)缺失。此種情況下，數(shù)據(jù)缺失的可能性依賴于未觀察到的變量。實(shí)驗(yàn)中時(shí)間序列數(shù)據(jù)缺失主要屬于①②2種情形。

當(dāng)數(shù)據(jù)存在缺失時(shí)，采用何種處理方法取決于數(shù)據(jù)缺失的特征，如若方法不適合數(shù)據(jù)，可能使得結(jié)果不可靠?，F(xiàn)有研究中所采用的方法大致可分為以下幾種：①完全集分析[6]，存在缺失的觀測(cè)將被剔除，僅有完整數(shù)據(jù)被用于分析。②單值填補(bǔ)，這種方法為每個(gè)缺失填補(bǔ)一個(gè)值，填補(bǔ)值可以是基于完整數(shù)據(jù)所計(jì)算的均值，也可能是根據(jù)完整數(shù)據(jù)擬合回歸模型再基于回歸模型的預(yù)測(cè)值做填補(bǔ)。③最大似然法，存在缺失數(shù)據(jù)時(shí)，似然法利用所觀察到的信息得到參數(shù)的似然函數(shù)，進(jìn)行參數(shù)估計(jì)。④多重填補(bǔ)[7]，這種方法將缺失視為一種不確定性，在特定假設(shè)下，認(rèn)為缺失的數(shù)據(jù)服從某種預(yù)測(cè)的概率分布;在填補(bǔ)時(shí)，每個(gè)缺失將被m個(gè)(m>l)從這個(gè)分布中所抽得的樣本值所替代，這m個(gè)值之間的變異代表缺失所引入的不確定性。由于每個(gè)缺失被m個(gè)填補(bǔ)值所代替，所以最終會(huì)得到m個(gè)完整的數(shù)據(jù)集，填補(bǔ)之后的數(shù)據(jù)分析將基于這m個(gè)完整的數(shù)據(jù)集;每一個(gè)填補(bǔ)后得到的完整集都可以采用完整數(shù)據(jù)的分析方法進(jìn)行分析。

在時(shí)間序列分析時(shí)，如果對(duì)缺失值采用簡單剔除的處理方式，那么對(duì)試驗(yàn)結(jié)果的解釋會(huì)產(chǎn)生影響數(shù)據(jù)的連續(xù)性、峰值的缺失、數(shù)據(jù)信息利用不充分等問題。對(duì)于時(shí)間序列數(shù)據(jù)，單個(gè)時(shí)間點(diǎn)數(shù)據(jù)缺失是較為常見的，單值填補(bǔ)是較為實(shí)用的。

研究結(jié)果顯示，同一時(shí)間點(diǎn)取材的數(shù)據(jù)都缺失時(shí)，結(jié)合相關(guān)文獻(xiàn)，能排除數(shù)據(jù)大幅變化，運(yùn)用SPSS填補(bǔ)缺失值的方法有其意義。在排除時(shí)間序列缺失數(shù)據(jù)臨近點(diǎn)大幅度增大和減小的情況下，可以運(yùn)用臨近中位數(shù)和插值法對(duì)缺失值進(jìn)行填充。在臨近值波動(dòng)幅度較小的時(shí)候，選用插值法填充值；在臨近值波動(dòng)幅度較大的時(shí)候，選用臨近中位數(shù)填充值。同一時(shí)間點(diǎn)取材的數(shù)據(jù)部分缺失時(shí)，建立時(shí)間序列模型擬合有很大意義，且同一時(shí)間點(diǎn)取材的數(shù)據(jù)部分缺失越少，擬合模型填充效果越好。但是本研究選用的填補(bǔ)方法相對(duì)簡單，如果對(duì)缺失值填補(bǔ)要求很高，需選用數(shù)據(jù)挖掘中更精確的辦法[8]。

目前尚未就缺失值處理方法達(dá)成廣泛一致，盡管如此，在處理缺失值時(shí)，結(jié)合具體情形加以考慮。在相關(guān)文獻(xiàn)中，缺失值的主要填補(bǔ)方法還有隨機(jī)回歸填補(bǔ)法、趨勢(shì)得分法(PS)、馬爾科夫鏈蒙特卡羅法(MCMC)和三次樣條插值法等。

作為出現(xiàn)缺失值的后果，數(shù)據(jù)分析困難會(huì)隨之發(fā)生，而且隨著缺失值數(shù)量的增加，這種情況愈加嚴(yán)重。因此，盡可能避免缺失值的出現(xiàn)就顯得極為重要，通過有效的設(shè)計(jì)最小化這種可能性，強(qiáng)化對(duì)數(shù)據(jù)的收集。如果是開創(chuàng)性的研究，相關(guān)文獻(xiàn)基礎(chǔ)較少時(shí)，應(yīng)在執(zhí)行實(shí)驗(yàn)過程中最低限度盡量減少缺失值。收集與研究目標(biāo)絕對(duì)必要的數(shù)據(jù)，防止施加一個(gè)研究人員不必要的負(fù)擔(dān)，提高收集的數(shù)據(jù)質(zhì)量。適當(dāng)?shù)臉颖玖抗烙?jì)，考慮實(shí)驗(yàn)中可能存在樣本的丟失。設(shè)置預(yù)實(shí)驗(yàn)或參考相關(guān)文獻(xiàn)，確定實(shí)驗(yàn)操作中可能發(fā)生的意外情況。提高實(shí)驗(yàn)人員的操作技能，規(guī)范實(shí)驗(yàn)操作，提高取材成功率。確定適合取材方法，保證指標(biāo)科學(xué)性，提高取材成功率。防止樣本的處理、保存和指標(biāo)測(cè)量中可能發(fā)生的失誤，提前進(jìn)行實(shí)驗(yàn)培訓(xùn)。

[1] 高卉,阮明鳳,龍浩成,等. 褪黑素對(duì)梗阻性黃疸大鼠小腸黏膜屏障保護(hù)作用的實(shí)驗(yàn)研究[J]. 現(xiàn)代中西醫(yī)結(jié)合雜志,2007,16(8):1024-1025

[2] 徐靜，馮前進(jìn)，牛欣. 時(shí)間序列插值法在天人相應(yīng)生物信號(hào)轉(zhuǎn)導(dǎo)數(shù)據(jù)分析中的應(yīng)用[J]. 中華中醫(yī)藥雜志，2012，29(4):895-899

[3] 高海威. 中醫(yī)藥臨床研究中缺失數(shù)據(jù)處理方法探討[D]. 廣州中醫(yī)藥大學(xué)，2012

[4] 張熙，李濟(jì)賓，張晉昕. 含有周期性的時(shí)間序列中隨機(jī)型缺失數(shù)據(jù)的填補(bǔ)方法[J]. 中國衛(wèi)生統(tǒng)計(jì)，2012，29(4):475-477

[5] Roderick JA Little,Donald B Rubin. 缺失數(shù)據(jù)統(tǒng)計(jì)分析[M]. 孫澤山，譯. 北京:中國統(tǒng)計(jì)出版社,2004：4-10

[6] Deznjssje S,LaValley MP,Horton NJ,et al. Bias due to missing exposure Data using complete-case analysis in the Proportional hazards regression model[J]. Statisticsin Medieine，2003,22(4):545-557

[7] Schafer JL，Graham JW. Missing data：our view of the state of the art[J]. Psychol Methods，2002，7(2)：147-177

[8] 陳光平. 基于時(shí)間序列數(shù)據(jù)特性的缺失值估計(jì)算法[J]. 計(jì)算機(jī)工程與應(yīng)用，2012，48(12):135-138

Study on the method to fill missing values of melatonin time series

Qian Shuwen, Li Haiyan, Yang Xuezhi, Li Hui, Zhang Shangshang, Lu Yu, Wang Jingping, Zhang Zhixia, Xu Jing

(Beijing University of TCM, Beijing 100029, China)

Objective It is to simulate two cases of data completely missing and partially missing at the same time, and to compare the fill effect of melatonin (MT) time series of every method by comparing the fill value and the actual value. Methods In the case of data completely missing at the same time, five SPSS fill methods were used to fill the simulated missing values. Filling results were compared with the actual value. In the case of data partially missing at the same time, in addition to the completely missing filling methods, the time series model fitting was added to fill. Results In the case of data completely missing at the same time, median fill method of approaching points and linear interpolation results were closer to the actual value. Two-factor factorial design analysis of variance results were better. When the actual value had small fluctuations, the linear interpolation result was closer to the actual value. When the actual value had small fluctuations, the linear interpolation result was closer to the actual value. When the actual value had the larger fluctuations, the median fill method of approaching points result was closer to the actual value. In the case of data partially missing at the same time, the time series model fitting had the better result. Conclusion In the case of data completely missing, excluding the very large increase and decrease, median fill method of approaching points and linear interpolation can be used to fill missing values. When the actual value has small fluctuations, choose linear interpolation. When the actual value has the larger fluctuations, choose median fill method of approaching points. In the case of data partially missing, choose time series model to fill.

times series; missing value; fill method

錢淑雯，女，碩士研究生在讀，研究方向?yàn)樘烊讼鄳?yīng)時(shí)間序列。

徐靜，E-mail：xuj@bucm.edu.cn

北京中醫(yī)藥大學(xué)自主課題(2013-JYBZZ-JS-041)；北京中醫(yī)藥大學(xué)自主課題(2014-JYBZZ-XS-031)；北京中醫(yī)藥大學(xué)可獲取的人體診斷信息關(guān)鍵技術(shù)創(chuàng)新團(tuán)隊(duì)(2011-CXTD-05)

10.3969/j.issn.1008-8849.2014.35.001

R-332

A

1008-8849(2014)35-3877-04

2014-05-30