變系數(shù)模型中穩(wěn)健估計(jì)方法的比較和應(yīng)用*

復(fù)旦大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計(jì)學(xué)教研室和公共衛(wèi)生安全教育部重點(diǎn)實(shí)驗(yàn)室(200032)

黃綠斕　趙耐青　秦國(guó)友△

?

·論著·

變系數(shù)模型中穩(wěn)健估計(jì)方法的比較和應(yīng)用*

復(fù)旦大學(xué)公共衛(wèi)生學(xué)院生物統(tǒng)計(jì)學(xué)教研室和公共衛(wèi)生安全教育部重點(diǎn)實(shí)驗(yàn)室(200032)

黃綠斕趙耐青秦國(guó)友△

【提要】目的在變系數(shù)模型中比較七種常見的穩(wěn)健估計(jì)方法與最小二乘法的表現(xiàn)，為變系數(shù)模型中估計(jì)方法的選擇提供依據(jù)。方法通過(guò)R軟件隨機(jī)模擬，以變系數(shù)模型產(chǎn)生數(shù)據(jù)并對(duì)其進(jìn)行污染，比較穩(wěn)健估計(jì)方法和最小二乘法估計(jì)結(jié)果的偏差、方差、均方誤差以及積分均方誤差的差異。結(jié)果當(dāng)數(shù)據(jù)存在擾動(dòng)時(shí)，尤其是存在X方向上的異常點(diǎn)時(shí)，M-Huber、最小絕對(duì)離差(least absolute deviation,LAD)估計(jì)、MM以及R這幾種穩(wěn)健方法的四項(xiàng)指標(biāo)幾乎都小于最小二乘法，其中，MM表現(xiàn)最好。而最小截?cái)嗥椒椒?least trimmed squares,LTS)、最小中位數(shù)平方法(least median of squares,LMS)以及S由于在R軟件中穩(wěn)定性較差，并不適用于變系數(shù)模型。結(jié)論在變系數(shù)模型中，當(dāng)有異常點(diǎn)存在時(shí)，采用MM估計(jì)能得到更加準(zhǔn)確的結(jié)果。

變系數(shù)模型穩(wěn)健異常點(diǎn)

通常，我們使用一般線性模型來(lái)研究變量之間的關(guān)系。比如，為比較A藥和B藥在療程為6個(gè)月中持續(xù)減肥的療效，我們將10個(gè)女性肥胖志愿者隨機(jī)分成2組(group=0為A藥組，group=1為B藥組)，分別考察這2組肥胖者在服藥前、服藥后3個(gè)月和6個(gè)月的體重變化(分別對(duì)應(yīng)t=0,1,2)。一般線性模型可表示為：

Δt=β0+β1t+β2group+ε

(1)

其中Δt表示個(gè)體在t時(shí)刻體重增量，即Δt=Yt-Y0。β2表示時(shí)間固定時(shí)，group變化一個(gè)單位，Δt平均變化β2個(gè)單位。

考慮到藥物的效果可能受時(shí)間影響，模型引入交互項(xiàng)：

Δt=β0+β1t+β2group+β3group·t+ε

(2)

此時(shí)，group的效應(yīng)為β2+β3t，對(duì)應(yīng)t=0,1,2。即時(shí)間固定時(shí)，group變化一個(gè)單位，Δt平均變化β2+β3t個(gè)單位。所以，藥物的效應(yīng)是關(guān)于時(shí)間的函數(shù)，隨時(shí)間線性變化。推廣到更一般的情況，模型可變?yōu)?/p>

Δt=g0(t)+group·g1(t)+ε

(3)

其中g(shù)0(t)，g1(t)為光滑函數(shù)，分別對(duì)應(yīng)公式(2)中的(β0+β1t)和(β2+β3t)。g1(t)表示固定基線Y0和時(shí)間t，group變化一個(gè)單位，Δt平均變化g1(t)個(gè)單位，g0(t)的解釋也類似。此模型實(shí)際上是變系數(shù)模型(varyingcoefficientmodel,VCM)[1]的特殊形式，不僅包含了上述所有的模型，而且更靈活、更容易解釋。而變系數(shù)模型更一般的形式為

Y=g0(r0)+X1g1(r1)+…+Xpgp(rp)+ε

(4)

其中Y為響應(yīng)變量，X1,X2,…,Xp以及r1,r2,…,rp均為協(xié)變量，gj(rj)(j=1,2…p)是未知的光滑函數(shù)，ε是隨機(jī)誤差且E(ε)=0,var(ε)=σ2。其中g(shù)j(rj)表示固定其他因素時(shí)Xj變化一個(gè)單位，Y平均變化gj(rj)個(gè)單位，這個(gè)平均變化量隨rj而發(fā)生改變。變系數(shù)模型是經(jīng)典線性模型的推廣，具有適應(yīng)性和解釋性強(qiáng)的特點(diǎn)，在經(jīng)濟(jì)金融、流行病學(xué)、環(huán)境科學(xué)以及生物醫(yī)學(xué)等領(lǐng)域也有著廣泛應(yīng)用[2-6]。

變系數(shù)模型中系數(shù)函數(shù)的估計(jì)通?？梢圆捎没诤说木植慷囗?xiàng)式估計(jì)以及樣條等方法[5]。目前，這些方法主要是建立在最小二乘法(ordinary least square,OLS)之上的，但是OLS方法對(duì)數(shù)據(jù)中異常點(diǎn)非常敏感，可能導(dǎo)致估計(jì)結(jié)果產(chǎn)生偏差，甚至得到完全錯(cuò)誤的結(jié)論[5]。因此，很多學(xué)者[7-11]提出了穩(wěn)健估計(jì)方法，這些方法對(duì)異常點(diǎn)有一定的抵抗能力。在線性模型下，穩(wěn)健估計(jì)方法的比較研究很多[12-14]，但是，在較線性模型復(fù)雜的變系數(shù)模型中穩(wěn)健估計(jì)方法研究并不多見。本文通過(guò)隨機(jī)模擬的方法比較在各種數(shù)據(jù)污染情況下變系數(shù)模型中幾種常見的穩(wěn)健方法與OLS的表現(xiàn)，為今后變系數(shù)模型中估計(jì)方法的選擇提供依據(jù)。

方　　法

本研究在自然立方樣條的基礎(chǔ)上，使用各種估計(jì)方法得到變系數(shù)模型中系數(shù)函數(shù)的估計(jì)。

自然樣條函數(shù)實(shí)際上是一分段多項(xiàng)式，首先對(duì)區(qū)間[a,b]進(jìn)行劃分:a=t1

隨機(jī)模擬研究

本節(jié)將通過(guò)隨機(jī)模擬比較各穩(wěn)健估計(jì)在變系數(shù)模型中的表現(xiàn)。

我們考慮如下變系數(shù)模型[18]：

Y=g1(r)+g2(r)X1+g3(r)X2+ε

其中g(shù)1(r)=exp(2r-1),g2(r)=8r(r-1),g3(r)=2sin(2πr)2,X1～N(0,12)，X2～B[1,0.6],r～U[0,1],ε～N(0,12)，由模型產(chǎn)生Y，從而建立未污染數(shù)據(jù)UC。

為了研究估計(jì)的穩(wěn)健性，我們通過(guò)隨機(jī)選取np個(gè)點(diǎn)替換成其他點(diǎn)的方式來(lái)對(duì)原始數(shù)據(jù)進(jìn)行污染，其中n表示樣本量，p表示污染比例。四種污染方式分別為：C1，y方向上的污染，對(duì)隨機(jī)選擇響應(yīng)變量Y的np個(gè)值乘以3產(chǎn)生異常點(diǎn)；C2，x方向上的污染，對(duì)隨機(jī)選擇協(xié)變量中連續(xù)變量X1的np個(gè)值加3來(lái)產(chǎn)生異常點(diǎn)，即均值漂移異常點(diǎn)；C3，誤差項(xiàng)的污染，通過(guò)替換誤差項(xiàng)來(lái)產(chǎn)生異常點(diǎn)，誤差項(xiàng)服從N(0，102)并與原誤差分布獨(dú)立；以及C4，誤差項(xiàng)的污染，誤差項(xiàng)服從自由度為3的t分布。其中樣本量為n=500，污染比例一般設(shè)為p=0.20，由于均值漂移異常點(diǎn)通常會(huì)對(duì)經(jīng)典的估計(jì)產(chǎn)生很大的影響，C2中p設(shè)為0.05。

在以上設(shè)定的每一參數(shù)組合下，隨機(jī)模擬Nsim=500次。以df=5的自然樣條為基礎(chǔ)，使用穩(wěn)健方法以及OLS對(duì)變系數(shù)模型進(jìn)行估計(jì)。為了衡量各方法對(duì)g(r)的估計(jì)精度，我們報(bào)告了500次模擬中積分均方誤差(integrated mean square error，IMSE)[18]的均值和標(biāo)準(zhǔn)差。其中每次模擬IMSE定義如下：

此外，我們還比較了各估計(jì)的平均絕對(duì)偏差、平均方差以及平均均方誤差，分別以ABIAS、AVAR、AMSE表示，并定義如下

gij(rk))2

gij(rk))2

其中i=1,…Nsim;j=1,2,3。{rk,k=1,…,ngrid}是r在[0,1]內(nèi)平均分布的柵格點(diǎn)，ngrid=200。

結(jié)　　果

本模擬考慮了未污染數(shù)據(jù)和擾動(dòng)數(shù)據(jù)，各方法500次模擬IMSE的均值和標(biāo)準(zhǔn)差結(jié)果見表1。首先，在C1～C3中MM估計(jì)表現(xiàn)最好，即IMSE的均值和標(biāo)準(zhǔn)差最小，其在UC和C4中表現(xiàn)也不錯(cuò)。其次，在UC中OLS的IMSE均值和標(biāo)準(zhǔn)差最小，但在擾動(dòng)數(shù)據(jù)C1～C4中，有了顯著的增大，尤其在C1、C2、C3中，遠(yuǎn)大于除LTS、LMS以及S估計(jì)以外的其他穩(wěn)健方法的結(jié)果。最后，LTS、LMS以及S估計(jì)即使在UC中的IMSE均值和標(biāo)準(zhǔn)差也遠(yuǎn)遠(yuǎn)大于其他方法，在擾動(dòng)數(shù)據(jù)中更甚，提示這三種穩(wěn)健方法可能不適合變系數(shù)模型的估計(jì)。

圖1表示通過(guò)穩(wěn)健方法以及OLS方法，在未污染數(shù)據(jù)UC和擾動(dòng)數(shù)據(jù)C1～C4中估計(jì)出的g(r)曲線。由于LTS、LMS以及S估計(jì)不穩(wěn)定，遠(yuǎn)遠(yuǎn)偏離真實(shí)的曲線，圖中并未畫出。從中我們可以發(fā)現(xiàn)，在這幾種方法中，真實(shí)的g(r)曲線與MM估計(jì)的曲線最接近，與OLS曲線相差最遠(yuǎn)。

表1　各方法的IMSE均值和標(biāo)準(zhǔn)差(IMSE)

*：僅取兩位小數(shù)，其中每一列的最小值用粗體表示。

圖1　穩(wěn)健估計(jì)方法以及OLS方法估計(jì)的g(r)曲線

各估計(jì)方法的ABIAS、AVAR以及AMSE見圖2，LTS、LMS以及S估計(jì)的結(jié)果由于不穩(wěn)定同樣并未給出。與表1的結(jié)果類似，這幾種方法在UC和C4下表現(xiàn)良好，有較小的ABIAS、AVAR以及AMSE。在擾動(dòng)數(shù)據(jù)C1～C3中，OLS的ABIAS、AVAR以及AMSE遠(yuǎn)遠(yuǎn)大于穩(wěn)健方法，MM的ABIAS、AVAR以及AMSE在幾乎所有的擾動(dòng)情形下都明顯小于或至少不大于其他方法。

圖2　各穩(wěn)健估計(jì)方法及OLS估計(jì)的ABIAS、AVAR以及AMSE

圖3　CD4數(shù)據(jù)通過(guò)變系數(shù)模型估計(jì)的gj(t)曲線(j=0,1,2,3)

實(shí)例分析

為了研究艾滋病自然史及其影響因素，多中心艾滋病研究[19]收集了283名感染HIV病毒的男同性戀患者1984-1991年隨訪情況。其中t表示患者艾滋病病毒診斷后的隨訪觀測(cè)時(shí)間(年)，Y是診斷后t時(shí)刻個(gè)體的CD4濃度，smoke表示該患者診斷前是否吸煙，age表示該患者經(jīng)中心化后診斷時(shí)的年齡，preCD4表示患者診斷時(shí)中心化的CD4濃度。為研究吸煙、年齡以及基線的CD4濃度對(duì)個(gè)體CD4濃度的影響，模型可設(shè)為

Y(t)=g0(t)+g1(t)smoke+g2(t)age+g3(t)preCD4+ε

圖3給出了該模型的估計(jì)結(jié)果。圖(a)可表示基線濃度為42的34歲不吸煙男同性戀患者CD4濃度隨時(shí)間不斷下降。圖(b)中吸煙的效應(yīng)在0附近波動(dòng)。圖(c)表明年齡的效應(yīng)幾乎是負(fù)向的，并且近似一條斜率為-0.08的直線。表明診斷時(shí)年齡越大CD4濃度越低，時(shí)間以及其他變量固定時(shí)，年齡增加一個(gè)單位個(gè)體t時(shí)刻CD4濃度平均下降0.08個(gè)單位。圖(d)表明基線CD4濃度越高，個(gè)體當(dāng)前CD4濃度越高，另外，基線的影響隨時(shí)間發(fā)生改變，前兩年急劇減小，之后趨于平緩。另外對(duì)gj(t)是否恒等于0進(jìn)行檢驗(yàn)，除g1(t)外P值均小于0.05，表明個(gè)體t時(shí)刻CD4濃度與年齡和基線情況有關(guān)，與吸煙無(wú)關(guān)。

討　　論

變系數(shù)模型實(shí)際上是更加一般的交互作用模型，交互作用呈非線性變化，在許多實(shí)際應(yīng)用中往往被忽略，在這類模型的理論研究中往往注重理論性質(zhì)而忽略了該模型具有很好的應(yīng)用價(jià)值和結(jié)果的詮釋，本文對(duì)模型的結(jié)果做了初步的詮釋。

本文模擬研究在變系數(shù)模型中，比較了幾種常見的穩(wěn)健方法和OLS在處理不同類型的異常點(diǎn)的結(jié)果差異，發(fā)現(xiàn)MM估計(jì)在各情形中綜合表現(xiàn)最好。

首先，不存在異常點(diǎn)時(shí)，使用自然樣條的OLS、MM、Huber-M、LAD以及R估計(jì)都能較準(zhǔn)確地得到變系數(shù)模型的估計(jì)，但LTS、LMS以及S估計(jì)效果較差。由于R中LTS、LMS以及S估計(jì)都是通過(guò)lqs函數(shù)來(lái)估計(jì)，這些方法很難得到精確的估計(jì)值并且計(jì)算量巨大。lqs函數(shù)用的是一種重抽樣的近似算法[20]，由于每次都是隨機(jī)抽樣，所以穩(wěn)定性相對(duì)較差。故而，在R軟件中，LTS、LMS以及S估計(jì)不適用于變系數(shù)模型。

其次，存在Y方向上的異常點(diǎn)時(shí)，基于OLS估計(jì)的結(jié)果準(zhǔn)確性以及穩(wěn)定性都不及MM、Huber-M、LAD和R估計(jì)這幾種穩(wěn)健方法，其中MM估計(jì)效果最好。X方向的異常點(diǎn)對(duì)所有的估計(jì)方法都有較大的影響，尤其是對(duì)OLS方法，其結(jié)果準(zhǔn)確性以及穩(wěn)定性都不及以上幾種穩(wěn)健方法。誤差項(xiàng)混合方差較大的正態(tài)分布時(shí)，對(duì)OLS影響較大。但誤差項(xiàng)混合t分布時(shí)對(duì)結(jié)果的影響不大，增加模擬加大污染比例至0.3以及調(diào)整自由度df=4,10等也幾乎沒有影響。

本研究的局限在于：(1)本研究只模擬了存在一個(gè)二分類變量和一個(gè)連續(xù)型變量、兩個(gè)協(xié)變量的情況，沒有對(duì)多分類以及多個(gè)協(xié)變量進(jìn)行模擬。(2)模型中由于沒有考慮兩個(gè)或以上的連續(xù)型協(xié)變量，并未考慮變量之間的相關(guān)性。(3)本研究考慮了4種類型的污染，但是實(shí)際數(shù)據(jù)往往更加復(fù)雜，一個(gè)數(shù)據(jù)中可能存在多種類型的污染。對(duì)于這些不足，我們將在以后做進(jìn)一步的研究。

[1]Hastie T,Tibshirani R.Varying-coefficient Models.Journal of the Royal Statistical Society.Series B(Methodological),1993,55(4):757-796.

[2]Fan J,Zhang W.Statistical methods with varying coefficient models.Stat Interface,2008,1(1):179-195.

[3]Fan J,Zhang W.Statistical estimation in varying coefficient models.Ann Stat,1999,27(5):1491-1518.

[4]Park B,Mammen E,Lee Y,et al.Varying Coefficient Regression Models:A Review and New Developments.Int Stat Rev,2015,83(1):36-64.

[5]Feng L,Zou C,Wang Z,et al.Robust spline-based variable selection in varying coefficient model.Metrika,2015,78(1):185-118.

[6]徐麗紅,劉志永,劉桂芬,等.縱向監(jiān)測(cè)連續(xù)非隨機(jī)缺失數(shù)據(jù)變系數(shù)模型及其應(yīng)用.中國(guó)衛(wèi)生統(tǒng)計(jì),2012,29(3):314-317.

[7]Yohai V,Zamar R.High Breakdown-Point Estimates of Regression by Means of the Minimization of an Efficient Scale.J Am Stat Assoc,1988,83(402):406-413.

[8]Rousseeuw P.Least median of squares regression.J Am Stat Assoc,1984,79(388):871-880.

[9]Rousseeuw P,Yohai V.Robust regression by means of S-estimators.Springer,1984.

[10]Jaeckel L.Estimating Regression Coefficients by Minimizing the Dispersion of the Residuals.The Annals of Mathematical Statistics,1972,43(5):1449-1458.

[11]Huber P.Robust estimation of a location parameter.The Annals of Mathematical Statistics,1964,35(1):173-101.

[12]Alma.Comparison of Robust Regression Methods in Linear Regression.Int J Contemp Math Sciences,2011,6(9):409-421.

[13]Anderson C,Schumacker R.A comparison of five robust regression methods with ordinary least squares regression:Relative efficiency,bias,and test of the null hypothesis.Understanding Statistics:Statistical Issues in Psychology,Education,and the Social Sciences,2003,2(2):179-103.

[14]Schumacker R,Monahan M,Mount R.A comparison of OLS and robust regression using S-PLUS.Multiple Linear Regression Viewpoints,2002,28(2).

[15]丁士俊,陶本藻.自然樣條非參數(shù)回歸模型及模擬分析.測(cè)繪通報(bào),2004-1-25(1):17-19.

[16]Hastie T,Tibshirani R,Friedman J.The Elements of Statistical Learning:Data Mining,Inference,and Prediction.Springer,2011.

[17]Rousseeuw P,Hubert M.Robust statistics for outlier detection.Wires Data Min Knowl,2011,1(1):73-79.

[18]Wang L,Kai B,Li R.Local Rank Inference for Varying Coefficient Models.J Am Stat Assoc,2009,104(488):1631-1645.

[19]Kaslow R,Ostrow D,Detels R,et al.The Multicenter AIDS Cohort Study:rationale,organization,and selected characteristics of the participants.Am J Epidemiol,1987,126(2):310-318.

[20]Rousseeuw P,Hubert M.Recent developments in PROGRESS.L1-Statistical Procedures and Related Topics IMS Lecture Notes,1997,3:201-214.

(責(zé)任編輯：鄧妍)

Comparison of Robust Methods for Varying Coefficient Model

Huang Lvlan,Zhao Naiqing,Qin Guoyou.

(Department of Biostatistics,School of Public Health and Key Laboratory of Public Health Safety,Fudan University(200032),Shanghai)

ObjectiveTo compare the performance of several common robust methods and Ordinary Least Square(OLS)in varying coefficient model.MethodsWe used R software to simulate uncontaminated data and contaminated data.Bias,variance,mean square error(MSE)and integrated mean square error(IMSE)were used for the evaluation indices to compare the performance of these robust methods and OLS.ResultsWhen outliers were present,especially occured in x-space,M-Huber,LAD(Least Absolute Deviation),MM and R performed much better than OLS with smaller Bias,variance,MSE and IMSE in almost all cases.Among them,MM performed best overall against a comprehensive set of outlier conditions.Furthermore,LTS(Least Trimmed Squares),LMS(Least Median of Squares)and S did not seem to apply in varying coefficient model for their instability in R software.ConclusionWhen outliers occured,MM resulted in more accurate results in varying coefficient model.

Varying coefficient model;Robustness;Outlier

國(guó)家自然科學(xué)基金(11371100)

秦國(guó)友，Email:gyqin@fudan.edu.cn