復旦大學公共衛(wèi)生學院生物統(tǒng)計學教研室和公共衛(wèi)生安全教育部重點實驗室(200032)

黃綠斕　趙耐青　秦國友△

?

·論著·

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

復旦大學公共衛(wèi)生學院生物統(tǒng)計學教研室和公共衛(wèi)生安全教育部重點實驗室(200032)

黃綠斕趙耐青秦國友△

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

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

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

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

(1)

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

考慮到藥物的效果可能受時間影響，模型引入交互項：

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

(2)

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

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

(3)

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

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

(4)

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

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

方　　法

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

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

隨機模擬研究

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

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

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

其中g1(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。

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

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

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

gij(rk))2

gij(rk))2

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

結　　果

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

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

表1　各方法的IMSE均值和標準差(IMSE)

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

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

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

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

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

實例分析

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

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

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

討　　論

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

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

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

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

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

[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)測連續(xù)非隨機缺失數(shù)據(jù)變系數(shù)模型及其應用.中國衛(wèi)生統(tǒng)計,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ù)回歸模型及模擬分析.測繪通報,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.

(責任編輯：鄧妍)

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

國家自然科學基金(11371100)

秦國友，Email:gyqin@fudan.edu.cn