徐禮文，瞿開毅

(1.中國人民大學(xué)統(tǒng)計(jì)學(xué)院，北京100872；2.北方工業(yè)大學(xué)理學(xué)院，北京100144)

帶有等相關(guān)誤差結(jié)構(gòu)生長曲線模型的參數(shù)boo tstrap檢驗(yàn)

徐禮文1,2，瞿開毅2

(1.中國人民大學(xué)統(tǒng)計(jì)學(xué)院，北京100872；2.北方工業(yè)大學(xué)理學(xué)院，北京100144)

文章研究了具有等相關(guān)誤差結(jié)構(gòu)的生長曲線模型回歸系數(shù)的檢驗(yàn)問題，構(gòu)造了參數(shù)bootstrap(PB)檢驗(yàn)統(tǒng)計(jì)量和相應(yīng)的PB檢驗(yàn)，并與已有的廣義p值(GP)檢驗(yàn)進(jìn)行了比較。模擬研究表明，PB方法和GP方法在單處理組情形下的表現(xiàn)趨于一致，均能很好的控制第一類錯誤率；在多處理組情形下，GP方法在一些情形下不能很好地控制犯第一類錯誤的概率，而PB方法則在很好地保證檢驗(yàn)名義水平的前提下，同時也具有良好的勢表現(xiàn)。

生長曲線模型；重復(fù)觀測；bootstrap重抽樣；廣義p值

0　引言

生長曲線模型在生物學(xué)、醫(yī)藥學(xué)、社會經(jīng)濟(jì)學(xué)、心理學(xué)等領(lǐng)域具有非常廣泛的應(yīng)用，許多學(xué)者對該模型進(jìn)行了大量研究[1-4]。由于生長曲線模型用于建模重復(fù)觀測的數(shù)據(jù)，且模型協(xié)方差矩陣因包含方差分量通常未知，傳統(tǒng)檢驗(yàn)方法大多基于近似理論，難以得到精確檢驗(yàn)。Weerahandi和Berger[5]提出運(yùn)用廣義p值(GP)方法構(gòu)造了具有誤差獨(dú)立結(jié)構(gòu)的簡單生長曲線模型回歸系數(shù)的精確檢驗(yàn)。Lin和Lee[6]進(jìn)一步研究了等相關(guān)誤差結(jié)構(gòu)的生長曲線模型回歸系數(shù)的廣義p值檢驗(yàn)。其中文獻(xiàn)[5,6]所說精確性其實(shí)是指廣義p值的計(jì)算公式是有精確表達(dá)式的。但廣義p值和經(jīng)典p值的定義卻有很大差異，有時無法保證經(jīng)典p值檢驗(yàn)的優(yōu)良性。

大量的研究表明，bootstrap方法在檢驗(yàn)問題中具有維持名義水平的優(yōu)良性質(zhì)，在許多情形下要優(yōu)于GP方法[7-9]。而鮮有學(xué)者研究生長曲線模型回歸系數(shù)檢驗(yàn)問題中兩種方法的優(yōu)良性比較。因此本文針對回歸系數(shù)檢驗(yàn)問題構(gòu)造了參數(shù)bootstrap(PB)檢驗(yàn)，并與GP方法進(jìn)行比較研究。

1　單處理組生長曲線模型回歸系數(shù)的單邊檢驗(yàn)

單處理情形生長曲線模型的一般形式為:

其中，Yit表示因變量在第i個體第t個時刻的觀測，Xt是協(xié)變量向量，β為相應(yīng)的回歸系數(shù)，αi和εit分別表示個體隨機(jī)效應(yīng)和隨機(jī)誤差項(xiàng)。其矩陣表達(dá)式為：

該協(xié)方差矩陣具有等相關(guān)結(jié)構(gòu)。在該協(xié)方差陣結(jié)構(gòu)下，可以證明回歸系數(shù)的廣義最小二乘估計(jì)和普通最小二乘估計(jì)是相同的[10]，也是此模型下最大似然估計(jì)。模型(2)的最小二乘估計(jì)為

1.1GP檢驗(yàn)

(bk，ssw，ssb)表示的一組樣本觀測，TGP的樣本觀測值T0=bk-βk，與除檢驗(yàn)參數(shù)外的其他討厭參數(shù)無關(guān)。因此，廣義p值可表示為：

1.3PB-S檢驗(yàn)

2　多處理組生長曲線模型回歸系數(shù)的相等性檢驗(yàn)

上文討論了只包含單個處理組時生長曲線模型回歸系數(shù)的檢驗(yàn)問題，下面考慮包含多個處理組時生長曲線模型回歸系數(shù)的相等性檢驗(yàn)，即H0:β1=…=βm.考慮的模型一般形式為：

αij和εijt分別為個體隨機(jī)效應(yīng)和隨即干擾項(xiàng)。其矩陣形式為：

2.1廣義F(GF)檢驗(yàn)

對于模型(8)，稱協(xié)方差矩陣Σi不相等的情形為異方差。模型的殘差平方和為其中

且有：

eij～N(0，IT).令為原假設(shè)下的標(biāo)準(zhǔn)化殘差平方和，為備擇假設(shè)下的標(biāo)準(zhǔn)化殘差平方和。對于檢驗(yàn)H0:β1=…=βm的廣義檢驗(yàn)變量為：

其中，

2.2PB-F檢驗(yàn)

于是,PB-F檢驗(yàn)統(tǒng)計(jì)量為:

3　隨機(jī)模擬研究

為了比較參數(shù)bootstrap方法與廣義p值方法在上述不同情形下的優(yōu)良性，我們采用Monte Carlo方法進(jìn)行模擬研究。具體地，分別模擬兩種檢驗(yàn)方法對于檢驗(yàn)問題的第一類錯誤概率和勢函數(shù)。為了便于模擬，對模型作如下假定

3.1單處理組回歸系數(shù)的模擬檢驗(yàn)

考慮β2及線性組合β2+β3的單邊檢驗(yàn)。設(shè)定回歸系數(shù)的真值為β=(10，2，1)′，模擬模型在不同樣本量n和不同組合情形下的第一類錯誤率和勢函數(shù)的表現(xiàn)?，F(xiàn)給出單個β2的PB-S檢驗(yàn)的p值算法如下：

內(nèi)循環(huán)：l=1 to L；

(4)利用(6)式計(jì)算TBS.如果TBS＞T0，記countl=1;否則記為0；

結(jié)束內(nèi)循環(huán)；

檢驗(yàn)的勢可通過類似的算法得到?？紤]到該模擬涉及到內(nèi)外兩層循環(huán)，本文設(shè)定內(nèi)循環(huán)和外循環(huán)的次數(shù)為(L，M)=(3000，2000)。計(jì)算過程是在SAS軟件環(huán)境下進(jìn)行的。第一類錯誤概率和檢驗(yàn)勢的模擬結(jié)果分別見表1和表2。

表1　單處理組情形下回歸系數(shù)的第一類錯誤概率的估計(jì)值

從表1的第一類錯誤率估計(jì)中可以看出，兩種方法在各種情形下都能很好的控制第一類錯誤概率；而從表2可以看出，兩種方法的檢驗(yàn)勢表現(xiàn)基本一致。值得注意的是，在不同的組合情形下，兩種方法的檢驗(yàn)勢差別較大。

表2　單處理組情形下回歸系數(shù)的檢驗(yàn)勢的估計(jì)值

3.2多處理組的回歸系數(shù)的相等性檢驗(yàn)?zāi)M比較

對于含3個處理組的生長曲線模型，我們模擬比較GF和PB-FS三種方法在不同樣本組合和不同協(xié)方差矩陣情形下的第一類錯誤率和檢驗(yàn)勢表現(xiàn)。PB-FS方法的p值算法如下：

循環(huán)：for i=1 to L;

(3)根據(jù)(13)式計(jì)算TPB-F,若TPB-F＞T0,記counti=1;否則記為0；

(4)即為p值的估計(jì)值；

與單處理組類似，假定回歸系數(shù)向量的真值為β1=β2=β3=(10，2，1)′，設(shè)定內(nèi)循環(huán)和外循環(huán)的次數(shù)為(L，M)=(3000，2000).并且在檢驗(yàn)勢的模擬中，為了便于比較，保持回歸系數(shù)向量βi的第一和第三分量不變，僅改變第二分量。檢驗(yàn)的模擬結(jié)果如表3和表4。

從表3可看出，當(dāng)樣本量非單調(diào)、三個處理組隨機(jī)誤差項(xiàng)的方差相差較大時，GF方法的第一類錯誤率顯著大于0.05，尤其當(dāng)樣本量=(8，4，6)時，GF的第一類錯誤率達(dá)到0.08，表現(xiàn)出一定的隨意性；而PB方法在各種情形下均表現(xiàn)穩(wěn)定。

表3　多處理組回歸系數(shù)相等性檢驗(yàn)的第一類錯誤率估計(jì)

表4　多處理組回歸系數(shù)相等性檢驗(yàn)勢估計(jì)

從表4可看出GF方法在樣本量非單調(diào)、三個處理組隨機(jī)誤差項(xiàng)的方差相差較大時，GF的檢驗(yàn)勢顯著大于PF-F方法，其原因是此情形下具有較高的第一類錯誤率。在其他情形下，兩種方法的檢驗(yàn)勢趨于一致。綜上分析，PB-F方法表現(xiàn)較為穩(wěn)定，而GF方法在一些情形下表現(xiàn)過于自由。

4　結(jié)論

當(dāng)考慮多個處理組的生長曲線模型回歸系數(shù)相等性比較時，廣義p值檢驗(yàn)無法控制第一類錯誤的問題沒有受到足夠重視。本文從單處理組和多處理組兩方面分別比較了參數(shù)bootstrap方法(PB)和廣義p值方法(GP)在生長曲線模型回歸系數(shù)檢驗(yàn)的功效表現(xiàn)。模擬結(jié)果表明，兩種方法在單處理組情形下的表現(xiàn)是較為一致的，均能很好的控制第一類錯誤率；而在多處理組情形下，當(dāng)樣本量非單調(diào)、處理組隨機(jī)誤差項(xiàng)的方差相差較大時，GP方法無法保證名義水平，第一類錯誤率偏高，而PB方法在各種情形表現(xiàn)更加穩(wěn)健，且具有良好勢的表現(xiàn)。

[1]Rao CR.Comparison ofGrowth Curves[J].Biometrics,1958，(14).

[2]Kshirsagar A M,Smith W B.Growth Curves[M].New York,Dekker, 1995.

[3]Pan J,Fang K.Growth Curve Models and Statistical Diagnostics[M]. New York:Springer-Verlag,2002.

[4]Ratan D.Advances in Growth Curve Models:Topics From the Indian Statistical Institute[M].New York:Springer,2013.

[5]Weerahandi S,Berger VW.Exact Inference for Growth CurvesWith IntraclassCorrelation Structure[J].Biometrics,1999,（55）.

[6]Lin SH,Lee JC.Exact Tests in Simple Growth Curve Models and One-way ANOVAWith Equicorrelation Error Structure[J].Journal of Multivariate Analysis,2003,（84）.

[7]Krishnamoorthy K,Lu F,Mathew T.A Parametric Bootstrap Ap?proach for ANOVAWith Unequal Variances:Fixed and Random Mod?els[J].ComputationalStatistics&Data Analysis,2007,（51）.

[8]Krishnamoorthy K,Lu F.A Parametric Bootstrap Solution to the MANOVA Under Heteroscedasticity[J].J.Stat.Comput.Simul.2010,（80）.

[9]Xu LW,Yang FQ,Abula A,et al.A parametric Bootstrap Approach for Two-Way ANOVA In Presence of Possible InteractionsWith Un?equalVariances[J].JournalofMultivariate Analysis,2013.

[10]Rao CR.Least-Squares Theory Using an Estimated Dispersion Ma?trix and Its Application To Measurement of Signals[C].Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability,1967,(1).

（責(zé)任編輯/易永生）

0212

A

1002-6487（2016）19-0027-05

國家自然科學(xué)基金資助項(xiàng)目(11171002)；北京市屬高等學(xué)校高層次人才引進(jìn)與培養(yǎng)計(jì)劃項(xiàng)目(CIT&TCD201404002)；北京市自然科學(xué)基金資助項(xiàng)目(9144026)

徐禮文(1977—)，男，安徽滁州人，博士，教授，研究方向：復(fù)雜數(shù)據(jù)分析。