基于广义交互的光滑样条非参数回归分析

基于广义交互的光滑样条非参数回归分析

回归分析的模型通常分为两种类型：参数回归模型和非参数回归模型。如果表示反应变量与解释变量之间数量关系的回归函数属于由有限个参数所决定的一类函数时,即回归函数的形式已知,而其中的参数未知,则模型叫参数回归模型。例如:直线回归模型和多重回归模型等。如果回归函数只限制属于某一光滑函数类(如函数是连续且可导的,并有平方可积的二阶导数),即属于某个无穷维的函数集合,则模型叫做非参数回归模型。例如:样条函数回归和kernel回归等。在实际工作中,参数回归模型经常会受到限制而不适用。如:反应变量与解释变量间的具体依存关系不明确、反应变量的分布不易判定等。在这些情况下,参数回归模型难以进行拟合处理,而非参数回归模型则能进行有效的分析,这是因为非参数模型在模型假定方面具有适应性很强的特点,即它所需的假定要比参数模型弱得多,因而应用范围更广,适应数据变化的能力更强。本文立足粗糙度惩罚方法(roughnesspenaltyapproach)和三次样条函数(cubicsplinefunction)的有机结合,构造了惩罚平方和(penalizedsumofsquares),使光滑样条非参数回归分析综合考虑了曲线回归的两个方面:曲线拟合优度和拟合曲线的光滑度。gt算法所谓样条函数,从数学上说,就是按一定光滑性要求“对接”起来的分段多项式。具体地说,对所考察的区间[a,b]作一分划Δ∶a=t0＜t1＜⋯＜tn＜tn+1=b如果函数Sk(t)在分划Δ的每个子段[ti,ti+1](i=0,1,…,n)上都是k次多项式,且在每个结点(knot)ti(i=1,2,…,n)上具有直到k-1阶连续导数,则称分段k次多项式Sk(t)为k次样条函数。在实际应用中,最常用的样条函数为三次样条函数,显然,三次样条函数S3(t)可表示为:S3(t)=g(t)=di(t-ti)3+ci(t-ti)2+bi(t-ti)+ai,ti≤t≤ti+1,i=0,1,⋯,n(1)如果三次样条函数g(t)的二阶、三阶导数在a与b处为0,则g(t)叫做三次自然样条(naturalcubicspline,NCS),这些约束条件叫自然边界条件(naturalboundaryconditions),即d0=c0=dn=cn=0,因此在[a,t1]和[tn,b]上,g(t)为线性函数。设在区间[a,b]上有n对数据(ti,yi),i=1,2,…,n,其中a<t1<t2<…<tn<b(若有必要,可进行排序),则可对数据(ti,yi)拟合以下模型yi=g(ti)+εi(2)假设g(t)为具有分划Δ的三次自然样条函数,则可规定g=(g1,…,gn)′,γ=(γ2,…,γn-1)′,其中gi=g(ti),γi=g″(ti),i=1,2,…,n,由自然边界条件可知,γ1=γn=0。令Q为n×(n-2)阶矩阵,元素为qij,i=1,2,…,n;j=2,…,n-1。qj+1,j=h-1j,qj-1,j=h-1j-1,qjj=-h-1j-1-h-1j,其中hi=ti+1-ti,i=1,2,…,n-1。如果|i-j|≥2,则qij=0。又令R为(n-2)×(n-2)阶对称阵,元素为rij,i,j=2,3,…,n-1。rii=13(hi-1+hi),i=2,⋯‚n-1ri,i+1=ri+1,i=16hi,i=2,⋯‚n-2如果|i-j|≥2,则rij=0。求解待定参数g和γ时,只需使以下惩罚平方和达到最小值即可。S(g)=n∑i=1{yi-g(ti)}2+α∫ba{g″(t)}2dt(3)其中a>0为光滑参数。由于g和γ完全确定三次自然样条函数g(t)的首要条件为Q′g=Rγ,则∫ba{g″(t)}2dt=γ′Rγ=g′Κg(4)其中K=QR-1Q′令Y=(y1,…,yn)′,则最小化S(g)可求得g,g=(I+αK)-1Y,其中I为单位矩阵,对于待定参数γ的求解,则可由以下线性方程组解得(R+αQ′Q)γ=Q′Y(5)由于Q′g=Rγ,因此,g=Y-αQγ。得到g和γ后,就可按以下三种情况来构造三次自然样条函数g(t)。(1)gt计算令hi=ti+1-ti,vi=γi+1-γihi‚ui=γi+1+γi+γiti+1-γi+1tihi,则g(t)=ai+bit+cit2+dit3其中ai=ti+1gi-tigi+1hi+ti+1tiui6,bi=gi+1-gihi-(ti+1+ti)ui-ti+1tivi6‚ci=-(ti+1+ti)vi-ui6‚di=vi6。(2)非参数回归分析g(t)=a0+b0t其中a0=g1-b0t1‚b0=g2-g1t2-t1-(t2-t1)γ26(3)t≥tn。g(t)=an+bnt其中an=gn-bntn,bn=gn-gn-1tn-tn-1+(tn-tn-1)γn-16。非参数回归分析需要对光滑参数α进行选择,α越大,拟合曲线越光滑,反之亦然。对于光滑参数α的选择,可以直接指定某个正数,或者通过模式搜索法迭代求解以使得下面的广义交互有效(generalizedcross-validation)得分函数取最小值。GCV(α)=n×SS残差(EDF)2(6)其中残差平方和SS残差=n∑i=1[yi-g(ti)]2,EDF=tr{I-A(α)},I为单位阵,A(α)=I-αQ(R+αQ′Q)-1Q′。拟合优度R2为R2=1-SS残差n∑i=1(yi-ˉy)2(7)其中ˉy=1nn∑i=1yi,拟合曲线的均方差MSE为ΜSE=SS残差EDF(8)本文利用6.11版SAS软件的IML模块进行编程来实现以上分析过程。拟合误差分析为研究血糖水平与胰岛素水平的关系,对16名糖尿病人的血糖水平与胰岛素水平进行测定,测定结果见表1。由表1的数据绘一散点图(图1),从散点图可知,y关于t基本上呈下降趋势,由病理分析可知,胰岛素水平越高,血糖水平就越低,由于血糖水平的分布未知,故采用光滑样条非参数回归分析方法。经计算,可得到使(6)式达到最小值时模式搜索法自动给出的光滑参数α=0.9811,拟合优度R2=0.9213,误差的自由度EDF=9.7189,残差平方和SS残差=1612.58,均方差MSE=165.92。非参数回归给出y的估计值y^和拟合误差见表2。同时,表2也给出了Monotonic回归与最小二乘线性回归的拟合结果,便于相互比较。从拟合误差看,非参数回归优于其他二种回归,Monotonic回归与LS回归的残差平方和分别为2094.95和3467.11,均方差分别为174.58和247.65。虽然t与y的线性相关程度较高,相关系数r=-0.9115,R2=0.8309,线性回归方程为y^=311.6987-6.9439t,但由于样本的随机性,从表2和图1可以看出,个别数据偏离较大,从而使拟合误差增大。而采用Monotonic秩次回归,虽然在一定程度上弱化了数据离散性,改善了拟合效果(R2=0.8978),但它考虑的是变量的秩次,不是具体的原始数据,因而损失了部分信息,而且其回归曲线过于粗糙抖动,非参数回归克服了上述缺陷,不但弱化了数据的离散性,而且进一步改善了拟合效果,其回归曲线也较为平滑。而且对于任意给定的点t,它都能给出预测值y^,例如t=12.5μu/100ml,y^=224.50mg/100ml。曲线光滑度检验1.惩罚平方和S(g)综合考虑了曲线拟合的两个方面:拟合优度和拟合光滑度。从(3)式可以看出,S(g)不仅与拟合优度(由残差平方和∑i=1n{yi-g(ti)}2度量)有关,而且也与所拟合的曲线的粗糙度∫ab{g″(t)}2dt有关。对于给定的光滑参数α值,最小化S(g)可得到最佳兼顾光滑度和拟合优度的曲线,即可得到g(t)的最小惩罚二乘(penalizedleastsquares)估计g^(t)。2.光滑参数α决定了拟合曲线的光滑度,如果α较大,则S(g)中的主要影响项为粗糙度惩罚项α∫ab{g″(t)}2dt,因此g^(t)的曲率较小,当α→∞时,∫ab{g^″(t)}2dt→0‚g^(t)接近直线拟合,曲线很光滑;另一方面,若α较小,则S(g)中的主要影响项为残差平方和,g^(t)接近数据点,曲线凸凹不平,当α→0时,g^(t)逼近插值曲线,曲线很粗糙(凸凹不平)。3.光滑样条非参数回归分析不必事先对结点进行选择,克服了以往利用样条函数进行曲线拟合时所存在的缺点,即避免了结点选择的盲目性以及由此所造成的计算不便,既提高了拟合优度,又在一定程度上保证了拟合曲线的光滑度,使得拟合曲线较为美观。4.非参数回归分析应用范围广,适应性强,它所需的假定要比参数回归模型弱得多,适用于任意分布的资料,尤其当反应变量与解释变量间的函数形式不清楚时,参数模型难以进行拟合处理,而非参数法却能加以有效分析,它只需知道函数具有某些性质即可,如:函数可微且其二阶导数的平方可积。因此,作为曲线拟合的一个很有效方法,光滑样条非参