徐禮文,瞿開毅
(1.中國人民大學(xué)統(tǒng)計學(xué)院,北京100872;2.北方工業(yè)大學(xué)理學(xué)院,北京100144)
帶有等相關(guān)誤差結(jié)構(gòu)生長曲線模型的參數(shù)boo tstrap檢驗
徐禮文1,2,瞿開毅2
(1.中國人民大學(xué)統(tǒng)計學(xué)院,北京100872;2.北方工業(yè)大學(xué)理學(xué)院,北京100144)
文章研究了具有等相關(guān)誤差結(jié)構(gòu)的生長曲線模型回歸系數(shù)的檢驗問題,構(gòu)造了參數(shù)bootstrap(PB)檢驗統(tǒng)計量和相應(yīng)的PB檢驗,并與已有的廣義p值(GP)檢驗進行了比較。模擬研究表明,PB方法和GP方法在單處理組情形下的表現(xiàn)趨于一致,均能很好的控制第一類錯誤率;在多處理組情形下,GP方法在一些情形下不能很好地控制犯第一類錯誤的概率,而PB方法則在很好地保證檢驗名義水平的前提下,同時也具有良好的勢表現(xiàn)。
生長曲線模型;重復(fù)觀測;bootstrap重抽樣;廣義p值
生長曲線模型在生物學(xué)、醫(yī)藥學(xué)、社會經(jīng)濟學(xué)、心理學(xué)等領(lǐng)域具有非常廣泛的應(yīng)用,許多學(xué)者對該模型進行了大量研究[1-4]。由于生長曲線模型用于建模重復(fù)觀測的數(shù)據(jù),且模型協(xié)方差矩陣因包含方差分量通常未知,傳統(tǒng)檢驗方法大多基于近似理論,難以得到精確檢驗。Weerahandi和Berger[5]提出運用廣義p值(GP)方法構(gòu)造了具有誤差獨立結(jié)構(gòu)的簡單生長曲線模型回歸系數(shù)的精確檢驗。Lin和Lee[6]進一步研究了等相關(guān)誤差結(jié)構(gòu)的生長曲線模型回歸系數(shù)的廣義p值檢驗。其中文獻[5,6]所說精確性其實是指廣義p值的計算公式是有精確表達式的。但廣義p值和經(jīng)典p值的定義卻有很大差異,有時無法保證經(jīng)典p值檢驗的優(yōu)良性。
大量的研究表明,bootstrap方法在檢驗問題中具有維持名義水平的優(yōu)良性質(zhì),在許多情形下要優(yōu)于GP方法[7-9]。而鮮有學(xué)者研究生長曲線模型回歸系數(shù)檢驗問題中兩種方法的優(yōu)良性比較。因此本文針對回歸系數(shù)檢驗問題構(gòu)造了參數(shù)bootstrap(PB)檢驗,并與GP方法進行比較研究。
單處理情形生長曲線模型的一般形式為:
其中,Yit表示因變量在第i個體第t個時刻的觀測,Xt是協(xié)變量向量,β為相應(yīng)的回歸系數(shù),αi和εit分別表示個體隨機效應(yīng)和隨機誤差項。其矩陣表達式為:
該協(xié)方差矩陣具有等相關(guān)結(jié)構(gòu)。在該協(xié)方差陣結(jié)構(gòu)下,可以證明回歸系數(shù)的廣義最小二乘估計和普通最小二乘估計是相同的[10],也是此模型下最大似然估計。模型(2)的最小二乘估計為
1.1GP檢驗
(bk,ssw,ssb)表示的一組樣本觀測,TGP的樣本觀測值T0=bk-βk,與除檢驗參數(shù)外的其他討厭參數(shù)無關(guān)。因此,廣義p值可表示為:
1.3PB-S檢驗
上文討論了只包含單個處理組時生長曲線模型回歸系數(shù)的檢驗問題,下面考慮包含多個處理組時生長曲線模型回歸系數(shù)的相等性檢驗,即H0:β1=…=βm.考慮的模型一般形式為:
αij和εijt分別為個體隨機效應(yīng)和隨即干擾項。其矩陣形式為:
2.1廣義F(GF)檢驗
對于模型(8),稱協(xié)方差矩陣Σi不相等的情形為異方差。模型的殘差平方和為其中
且有:
eij~N(0,IT).令為原假設(shè)下的標(biāo)準(zhǔn)化殘差平方和,為備擇假設(shè)下的標(biāo)準(zhǔn)化殘差平方和。對于檢驗H0:β1=…=βm的廣義檢驗變量為:
其中,
2.2PB-F檢驗
于是,PB-F檢驗統(tǒng)計量為:
為了比較參數(shù)bootstrap方法與廣義p值方法在上述不同情形下的優(yōu)良性,我們采用Monte Carlo方法進行模擬研究。具體地,分別模擬兩種檢驗方法對于檢驗問題的第一類錯誤概率和勢函數(shù)。為了便于模擬,對模型作如下假定
3.1單處理組回歸系數(shù)的模擬檢驗
考慮β2及線性組合β2+β3的單邊檢驗。設(shè)定回歸系數(shù)的真值為β=(10,2,1)′,模擬模型在不同樣本量n和不同組合情形下的第一類錯誤率和勢函數(shù)的表現(xiàn)?,F(xiàn)給出單個β2的PB-S檢驗的p值算法如下:
內(nèi)循環(huán):l=1 to L;
(4)利用(6)式計算TBS.如果TBS>T0,記countl=1;否則記為0;
結(jié)束內(nèi)循環(huán);
檢驗的勢可通過類似的算法得到??紤]到該模擬涉及到內(nèi)外兩層循環(huán),本文設(shè)定內(nèi)循環(huán)和外循環(huán)的次數(shù)為(L,M)=(3000,2000)。計算過程是在SAS軟件環(huán)境下進行的。第一類錯誤概率和檢驗勢的模擬結(jié)果分別見表1和表2。
表1 單處理組情形下回歸系數(shù)的第一類錯誤概率的估計值
從表1的第一類錯誤率估計中可以看出,兩種方法在各種情形下都能很好的控制第一類錯誤概率;而從表2可以看出,兩種方法的檢驗勢表現(xiàn)基本一致。值得注意的是,在不同的組合情形下,兩種方法的檢驗勢差別較大。
表2 單處理組情形下回歸系數(shù)的檢驗勢的估計值
3.2多處理組的回歸系數(shù)的相等性檢驗?zāi)M比較
對于含3個處理組的生長曲線模型,我們模擬比較GF和PB-FS三種方法在不同樣本組合和不同協(xié)方差矩陣情形下的第一類錯誤率和檢驗勢表現(xiàn)。PB-FS方法的p值算法如下:
循環(huán):for i=1 to L;
(3)根據(jù)(13)式計算TPB-F,若TPB-F>T0,記counti=1;否則記為0;
(4)即為p值的估計值;
與單處理組類似,假定回歸系數(shù)向量的真值為β1=β2=β3=(10,2,1)′,設(shè)定內(nèi)循環(huán)和外循環(huán)的次數(shù)為(L,M)=(3000,2000).并且在檢驗勢的模擬中,為了便于比較,保持回歸系數(shù)向量βi的第一和第三分量不變,僅改變第二分量。檢驗的模擬結(jié)果如表3和表4。
從表3可看出,當(dāng)樣本量非單調(diào)、三個處理組隨機誤差項的方差相差較大時,GF方法的第一類錯誤率顯著大于0.05,尤其當(dāng)樣本量=(8,4,6)時,GF的第一類錯誤率達到0.08,表現(xiàn)出一定的隨意性;而PB方法在各種情形下均表現(xiàn)穩(wěn)定。
表3 多處理組回歸系數(shù)相等性檢驗的第一類錯誤率估計
表4 多處理組回歸系數(shù)相等性檢驗勢估計
從表4可看出GF方法在樣本量非單調(diào)、三個處理組隨機誤差項的方差相差較大時,GF的檢驗勢顯著大于PF-F方法,其原因是此情形下具有較高的第一類錯誤率。在其他情形下,兩種方法的檢驗勢趨于一致。綜上分析,PB-F方法表現(xiàn)較為穩(wěn)定,而GF方法在一些情形下表現(xiàn)過于自由。
當(dāng)考慮多個處理組的生長曲線模型回歸系數(shù)相等性比較時,廣義p值檢驗無法控制第一類錯誤的問題沒有受到足夠重視。本文從單處理組和多處理組兩方面分別比較了參數(shù)bootstrap方法(PB)和廣義p值方法(GP)在生長曲線模型回歸系數(shù)檢驗的功效表現(xiàn)。模擬結(jié)果表明,兩種方法在單處理組情形下的表現(xiàn)是較為一致的,均能很好的控制第一類錯誤率;而在多處理組情形下,當(dāng)樣本量非單調(diào)、處理組隨機誤差項的方差相差較大時,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é)基金資助項目(11171002);北京市屬高等學(xué)校高層次人才引進與培養(yǎng)計劃項目(CIT&TCD201404002);北京市自然科學(xué)基金資助項目(9144026)
徐禮文(1977—),男,安徽滁州人,博士,教授,研究方向:復(fù)雜數(shù)據(jù)分析。