金晶亮

（南通大學(xué)，江蘇南通226007）

0 引 言

在復(fù)雜系統(tǒng)中，具有很多不確定性的因素，建立數(shù)學(xué)模型時經(jīng)常需使用概率統(tǒng)計模型。一方面，隨著回歸分析理論的不斷發(fā)展，隨機(jī)效應(yīng)模型已成為目前重要的研究課題；另一方面，數(shù)據(jù)一般含有刪失或不精密的特點(diǎn)。刪失分為“右刪失”和“左刪失”，個體的確切壽命不知道，只知道壽命大于L，則稱該個體的壽命在L是右刪失的，并說L是右刪失數(shù)據(jù)［1］。由于刪失的引入，情況大為復(fù)雜。普通的統(tǒng)計學(xué)未討論這些，它只是討論每個數(shù)據(jù)都是完全數(shù)據(jù)的情形。生存分析的一大特點(diǎn)，就是討論含刪失數(shù)據(jù)的情形，因而發(fā)展出許多新的統(tǒng)計方法。綜合上述兩方面，顯然對既有隨機(jī)效應(yīng)又有刪失數(shù)據(jù)的研究也是十分必要的。基于實(shí)際問題和理論研究的需要，本文主要研究既有隨機(jī)效應(yīng)又有刪失數(shù)據(jù)的模型，側(cè)重于統(tǒng)計診斷技術(shù)的研究。

1 含右刪失數(shù)據(jù)的非線性隨機(jī)效應(yīng)模型

假設(shè)Y是一個（n+m）維的響應(yīng)向量，u是q維的隨機(jī)效應(yīng)因子，假設(shè)u服從正態(tài)分布N（0，σ2Iq），Y1|u，Y2|u，…Yn+m|u，相互獨(dú)立，且服從正態(tài)分布 N（μi，σ2），不失一般性，考慮最后 m 個生命時間數(shù)據(jù)由于試驗(yàn)的終止卻未壽終而刪失了。Yi|u在隨機(jī)效應(yīng)下的均值為：μi=（fxi，β）+u（i=1，2，…，n+m），其中β是p×1未知固定效應(yīng)向量，令φ（y）=，聯(lián)合似然函數(shù)方程由下式給出：L｝，對數(shù)似然函數(shù)可以寫成：l（β，σ2，u）=-uTu，此處方程對β求導(dǎo)，可得），其中，wi具有右刪失數(shù)據(jù)的原隨機(jī)效應(yīng)模型Y=（fX，β）+Zu+ε化為一般隨機(jī)效應(yīng)模型：

故式1可以改寫為W=f（β）+e，e～N（0，σ2）。由此可看出，W關(guān)于（β，σ2）的對數(shù)似然為式2前3項(xiàng)，因此lp（β）作為正態(tài)的非線性隨機(jī)效應(yīng)模型的邊緣似然是精確的（即Op（n-1）=0），Laplace的展開對于正態(tài)非線性隨機(jī)效應(yīng)模型是精確的。

1.1 模型的參數(shù)估計

介紹了含右刪失數(shù)據(jù)的非線性隨機(jī)效應(yīng)模型，討論參數(shù)β的估計方法及有關(guān)性質(zhì)。

引理1對于正態(tài)非線性隨機(jī)效應(yīng)模型1，似然函數(shù)lp（β）關(guān)于參數(shù)β的Score函數(shù)與觀察信息陣分別是］［G］，其中為方程?l（W，u；β）/?u=0的解，=（β），G=?2f（β）/?β?βT-?2W/?β?βT為 n× p× p立體陣，Ω=方括號乘積［·］［·］［2］。

證明式2對β求導(dǎo)得，-ip（β）T）-1代入即得結(jié)論。

1.2 算法

根據(jù)引理1，β的參數(shù)估計可采用通常的迭代算法，具體的算法如下：

（1）首先給出參數(shù)β的初值。這可采用不帶隨機(jī)效應(yīng)的一般指數(shù)族非線性回歸模型的常用算法［3，4］，得到參數(shù)β的估計并取為初值β0，并取u0=ZT（Y-μ）|；（2）根據(jù)迭代計算出的βi，ui，給出Wi，；（3）給定參數(shù)β，解方程?（lW，u；β）/?u=0，即采用下列迭代公式：ui+1=（；（4）對于給定的u，解方程?l（W，u；β）/?β=0，求解參數(shù)β，這時可采用Gauss-Newton迭代法［5］。由于：ip（β），以φDTΩ-1D代替-¨lp（β），則Gauss-Newton迭代公式為若|βi+1-βi|與|ui+1-ui|均小于給定的精度，則停止迭代，否則重復(fù)上述步驟2-4，直至達(dá)到給定的精度。

2 含右刪失數(shù)據(jù)的線性隨機(jī)效應(yīng)模型

2.1 模型的參數(shù)估計

假設(shè)Y是一個（n+m）維的響應(yīng)向量，u是q維的隨機(jī)效應(yīng)因子，假設(shè) u服從正態(tài)分布N（0，σ2uIq），Y1|u，Y2|u，…，Yn+m|u相互獨(dú)立，且服從正態(tài)分布N（μi，σ2ε），不失一般性，考慮最后m個生命時間數(shù)據(jù)由于試驗(yàn)的終止卻未壽終而刪失了。Yi|u在隨機(jī)效應(yīng)下的均值為：μi=xTiβ+ziTu（i=1，2，…，n+m），其中β是p×1未知固定效應(yīng)向量令φ（y聯(lián)合似然函數(shù)為：L=｝，記λ0=對數(shù)似然函數(shù)為：

2.2 實(shí)例分析

某產(chǎn)品的使用壽命與溫度T的高低有關(guān)，數(shù)據(jù)來自加速壽命試驗(yàn)如表1所示（數(shù)據(jù)來源：文獻(xiàn)1）。根據(jù)同類產(chǎn)品的知識，壽命時間t與絕對溫度T之間應(yīng)有關(guān)系式

表1 加速壽命試驗(yàn)數(shù)據(jù)（欄目中“1”表示壽終，“0”表示未壽終）

由于試驗(yàn)數(shù)據(jù)是重復(fù)測量的，可以考慮隨機(jī)效應(yīng)模型：Yij=μi+αj+eij，i=1，2，3，4，j=1，2，3，4，其中μi是固定效應(yīng)，α1，α2，α3，α4為隨機(jī)效應(yīng)，假定αj，eij都互不相關(guān)，且E（αj）=0 Var（αj）=，Var（eij）=。利用前面介紹的算法，可以得到估計值根據(jù)前面討論的統(tǒng)計診斷量進(jìn)行計算，得到如下結(jié)果如表2所示：顯示10號，12號，16號是強(qiáng)影響點(diǎn)。

表2 加速壽命試驗(yàn)診斷統(tǒng)計量（min）

［1］ 彭非，王偉.生存分析［M］.北京：中國人民大學(xué)出版社，2004：197.

［2］ 王松桂.線性模型理論［M］.北京：科學(xué)出版社，2004：17-43.

［3］ Wei BC.Exponential Family Nonlinear Models［M］.Singapore：Springer-Verlag，1998：215-228.

［4］ Wei B C.On confidence regions of embedded models in regular parametric families（A geometric approach）［J］.Austral：J Statist，1994，36（2）：327-328.

［5］ 韋博成，魯國斌，史建清，等.統(tǒng)計診斷引論［M］.南京：東南大學(xué)出版社，1991：23-56.