基于多重比較檢驗的函數(shù)型數(shù)據(jù)方差分析

白永昕 ，田茂再 ,2a，2b

（1.蘭州財經(jīng)大學(xué) 統(tǒng)計學(xué)院，蘭州 730020；2.中國人民大學(xué)a.應(yīng)用統(tǒng)計科學(xué)研究中心；b.統(tǒng)計學(xué)院，北京 100972）

0 引言

函數(shù)型方差分析是函數(shù)型數(shù)據(jù)分析的對象之一。方差分析的目的在于研究諸多因素在不同狀態(tài)下對因變量的取值是否有顯著影響。因此，對函數(shù)型的研究對象在某方面的差異進(jìn)行統(tǒng)計分析時往往需要用到函數(shù)型方差分析。函數(shù)型方差分析模型擬合和參數(shù)估計的研究相對較多[1-3]。在方差分析中，若控制變量確實對觀測變量的顯著影響存在，那么本文需要進(jìn)一步研究控制變量的不同水平對觀測變量的影響程度。單純的通過兩獨立樣本的t檢驗會增大犯棄真的概率，對此可通過多重比較檢驗解決[4-7]。

甘肅是一個地理分布極其特殊的省份。李棟梁[8]根據(jù)氣候區(qū)劃分指標(biāo)將甘肅省分為8個氣候區(qū)。從東南到西北包括了北亞熱帶濕潤區(qū)到高寒區(qū)、半干旱區(qū)和干旱區(qū)的各種氣候類型。年平均氣溫在0～14℃之間,但是各個氣候帶的氣溫差異很大。為此，甘肅省不同氣候帶對氣溫的影響很有研究的必要。本文對函數(shù)型方差分析模型的參數(shù)估計進(jìn)行介紹，并將多元中方差分析的多重比較檢驗應(yīng)用到函數(shù)型方差分析中，分析了甘肅省不同氣候帶對氣溫的影響。

1 模型及方法

1.1 模型

一般地，假定控制變量A有m個水平，每個水平均有n個樣本。那么在水平Aj下的第i次連續(xù)的觀測結(jié)果為：

其中，μ(t)為觀測變量總的期望均值函數(shù)；aj(t)'j=1'…'m是控制變量水平 j對觀測變量產(chǎn)生的附加影響,稱為水平 j對觀測變量產(chǎn)生的效應(yīng),且是殘差函數(shù)，假定每個εij(t)都是均值為零的獨立的高斯過程。

定義一個(mn)×(m+1)的設(shè)計陣Z,每一行代表單獨的一個觀測。符號(ij)表示第 j個水平下第i個觀測相應(yīng)的行,此行的第一列是1,第 j+1列是1,剩余全為0。z(ij)k是設(shè)計陣Z中對應(yīng)第(ij)行第k列的值。同時，定義一個相應(yīng)的函數(shù)型向量 β1=μ'β2=a1,以此類推 βm+1=am,因此得到函數(shù)型向量β=(μ'…'am),則模型（1）可以表達(dá)為：

假設(shè)響應(yīng)變量是離散數(shù)據(jù)平滑的結(jié)果,因此在實際應(yīng)用中,一般直接用原始觀測數(shù)據(jù)Y來估計函數(shù)型參數(shù)。為了簡化,其模型（2）的向量表達(dá)形式為：

其中，Y是包含mn個氣溫的函數(shù)型向量,β(t)是m+1維的函數(shù)型參數(shù)向量,ε(t)是m+1維的函數(shù)型殘差向量。

1.2 參數(shù)估計

1.2.1 參數(shù)B的估計

下面開始討論模型中未知函數(shù)型參數(shù)β(t)的估計。對上述函數(shù)型線性模型中的β(t)用基函數(shù)展開,則有：

其中,B是5×Kβ的系數(shù)矩陣,θ(t)是長度為Kβ的基函數(shù)向量,Kβ是基函數(shù)展開的個數(shù)。對于B的懲罰函數(shù),引入一個線性算子L來定義：

關(guān)于微分算子的計算將在下文給出。

如果模型是標(biāo)準(zhǔn)的一般線性模型，β(t)則可以通過標(biāo)準(zhǔn)的最小二乘估計得到。將最小二乘方法拓展到函數(shù)型數(shù)據(jù)，只需以適當(dāng)?shù)姆绞街匦露x殘差平方和。此時Y-ZBθ是一個函數(shù)型的向量，由方程（3）、方程（4）和方程（5）得到懲罰最小二乘為：

運(yùn)用Kronecker積最小化上式可以得到B的估計值為：

1.2.2 微分算子的估計

定義線性微分算子為：

與之對應(yīng)的微分方程為：

其中，C=(C0'…'CP-1)表示的是微分算子系數(shù)；θ=(1'sin(wt)'…)T是傅里葉基函數(shù)。

令W=(θ'…'Dp-1θ)，則方程（6）可以表達(dá)：

求解以上方程得到向量C的估計值C?為：

由上式則可得到：

進(jìn)一步，R的估計值R?為：

2 方差分析

為了研究水平Aj對觀測變量產(chǎn)生的效應(yīng)，上文中估計了不同水平下的觀測變量曲線。接下來需要進(jìn)一步分析控制變量在固定的時間t上是否對觀測變量函數(shù)產(chǎn)生了顯著影響。為此開始進(jìn)行函數(shù)型方差分析。與普通的方差相同，函數(shù)型方差分析需要考慮殘差平方和函數(shù)LMSE和F比值函數(shù)FRATIO。不同于標(biāo)準(zhǔn)多元情況的是在特定時間這些量是相互獨立的。在上述模型中,殘差平方和函數(shù)為：

同時，令：

類似與單因素方差分析，可以進(jìn)一步計算殘差的均方函數(shù)為：

MSE=SSE/df(error)

其中,df(error)是誤差的自由度。

同樣,回歸的均方函數(shù)是SSY(t)與SSE(t)的差值除以回歸的自由度。因此：

其中,df(regression)表示兩個模型中誤差自由度的差。最后,可以得到F檢驗函數(shù)為：

3 假設(shè)檢驗

方差分析只能判斷控制變量是否對觀測變量產(chǎn)生了顯著影響。如果控制變量確實對觀測變量產(chǎn)生了顯著影響，那么感興趣的是：控制變量的不同水平對觀測變量的影響程度如何；其中哪個水平的作用明顯異于其他水平；哪個水平的作用不顯著等。盡管這些問題可以通過兩獨立樣本的t檢驗解決，但是進(jìn)行多次比較會增大犯第一類錯誤的概率。對此，本文通過多元方差分析中的多重比較檢驗來解決。

Tukey的HSD檢驗是較為常用的多重比較檢驗。它將控制變量k個水平下觀測變量的總體均值做兩兩對比檢驗。原假設(shè)為第i個總體與第 j個總體的均值函數(shù)無顯著差異。即：

H0:μ1=…=μm

H1:μ1≠…≠μm

檢驗統(tǒng)計量為q統(tǒng)計量。第i個總體與第 j個總體均值對比檢驗的檢驗統(tǒng)計量定義為：

其中，LMSE為殘差平方和函數(shù)；ni和nj分別為第i和第 j水平下的樣本量。不同于LSD檢驗的t統(tǒng)計量，q統(tǒng)計量近似服從有n-k個自由度的學(xué)生化極差q分布。假設(shè)檢驗的顯著性通過q統(tǒng)計量的p值判斷。

4 實證分析

4.1 數(shù)據(jù)來源

數(shù)據(jù)來自中國氣象局提供的天水（站號57006）、宕昌（站號56095）、蘭州（站號52889）、靖遠(yuǎn)（站號52895）、平?jīng)觯ㄕ咎?3915）、正寧（站號53935）、張掖（站號52625）、金昌（站號52675）、武威（站號52679）、敦煌（站號52418）、碌曲（站號56071）、卓尼（站號56082）和舟曲（站號56094）共13個地面氣候觀測站常規(guī)4次觀測得到的2015年逐日日最高溫度、日最低溫度。日平均氣溫由每日四個時點氣溫值求均值得到。每個站所有缺失數(shù)據(jù)均少于7天,所缺失的數(shù)據(jù)由相鄰前后5天氣溫的平均值代替。

4.2 模型

根據(jù)甘肅省特殊的地理特征,本文把甘肅省劃分為隴南暖溫帶濕潤區(qū)、隴中溫帶半干旱區(qū)、河西暖溫帶干旱區(qū)和甘南高寒濕潤區(qū)四個不同的氣候帶。同時,在這些氣候帶中選出13個可靠性較高的氣象觀測站,它們分別是隴南暖溫帶濕潤區(qū)的天水和宕昌站、隴中溫帶半干旱區(qū)的蘭州、靖遠(yuǎn)、平?jīng)龊驼龑幷?、河西暖溫帶干旱區(qū)的張掖、金昌、武威和敦煌站以及甘南高寒濕潤區(qū)的碌曲、卓尼和舟曲站。

每個氣候帶有不同數(shù)量的觀測站，假設(shè)第g個氣候帶中第m個氣溫函數(shù)的模型為：

Ymg(t)=μ(t)+ag(t)+εmg

其中,Ymg表示第g個氣候帶中第m個氣溫觀測站的溫度函數(shù)。均值函數(shù)μ(t)表示甘肅省平均氣溫，即所有氣溫觀測站的平均溫度。ag(t)表示第g個氣候帶對氣溫的影響函數(shù)，并且滿足殘差函數(shù)εmg表示第g個氣候帶中第m個氣溫觀測站不可解釋的變化。

4.3 方差分析

從上述模型出發(fā)，使用R軟件進(jìn)行數(shù)據(jù)處理，得到了甘肅省4個氣候帶的氣溫曲線圖（見圖1）以及不同氣候帶中總體均值的95%置信帶（見圖2）。

圖1 實線表示不同氣候帶的氣溫曲線，虛線表示甘肅省平均氣溫

圖2 不同氣候帶中總體均值的95%置信帶

從圖1中可以看到：隴南地區(qū)的全年氣溫均高于甘肅省平均氣溫。春季和秋季的氣溫與均值的差異較小，在夏季的七月份則差異較為明顯。冬季氣溫普遍高出全省均值4℃左右，冬季的差異最明顯。隴中地區(qū)全年的溫度最接近甘肅省平均氣溫。只有冬季氣溫略高于平均氣溫，但差異較小。河西走廊地區(qū)相比于全省平均氣溫則變化較為多樣性。春季和秋季氣溫略低于平均氣溫，夏季則高于平均氣溫，但是與平均氣溫的差異都不太明顯。到了冬季，氣溫則明顯下降，比平均氣溫要低3～4℃。甘南地區(qū)的氣溫變化情況與河西走廊地區(qū)恰恰相反。甘南地區(qū)在春、夏、秋三個季節(jié)的氣溫都低于全省平均氣溫。其中，春季和秋季與均值的差異較小，在夏季達(dá)到最大，比平均氣溫要高出3～4℃。冬季氣溫則要高于平均氣溫，并且差異相對較小。

方差分析的F比值函數(shù)如圖3所示。

圖3 函數(shù)型方差分析的F比值函數(shù)

從圖3方差分析的F比值函數(shù)可以看出，F(xiàn)比值均高于5%顯著水平下的2.92，表明4個氣候帶對氣溫的影響有顯著的差異。

以上分析結(jié)果可以表明4個氣候帶對氣溫的影響有顯著的差異，但是并不能表明每兩個氣候帶之間的差異都顯著，也不能表明哪些氣候帶的影響異于其他氣候帶的影響，哪些的影響不顯著等。因此，本文用單因素方差分析中的HSD多重比較檢驗進(jìn)行兩兩氣候帶之間的比較。同時，將整個函數(shù)型觀測變量的值域分成春季（3月至5月）、夏季（6月至8月）、秋季（9月至11月）和冬季（12月至次年2月）四個區(qū)間。多重比較檢驗的結(jié)果見表1所示。在顯著性數(shù)值上用“*”表示均值之間的對比在0.05的置信水平上有顯著差異。

表1 HSD多重比較表

多重比較檢驗的結(jié)果分析如下：在春季，隴南和隴中（p=0.3578）、隴南和河西走廊（p=0.1018）以及隴中和河西走廊（p=0.9147）的對比都沒有顯著差異。隴南和甘南（p=0）、隴中和甘南（p=0.0069）以及河西走廊（p=0.0480）的對比出現(xiàn)了顯著差異。在夏季，除了隴中和河西走廊（p=0.9986）的對比無顯著差異。其余氣候帶的兩兩對比都存在顯著差異。同樣在秋季，顯著差異出現(xiàn)在隴南和河西走廊（p=0.3895）、隴南和甘南（p=0.0051）以及河西走廊和甘南（p=0.0002）的對比中。冬季同夏季相似，除了隴中和甘南（p=0.8677）的對比無顯著差異，其余氣候帶的兩兩對比均存在顯著差異。

最后，感興趣的是在整個區(qū)間上氣候帶對氣溫影響的成對比較（見圖4）。從圖中可以看到：隴南和隴中對氣溫的影響只有在夏季（p=0.0161）和冬季（p=0）有顯著差異；隴南和河西走廊對比，除了春季的其他季節(jié)均有顯著差異；隴南和甘南則在四個季節(jié)均出現(xiàn)了顯著性的差異；而隴中與河西走廊對比，只有冬季有顯著差異；相比之下，隴中和甘南的顯著性差異出現(xiàn)在春季（p=0.069）和夏季（p=0）；類似于隴南和甘南，河西走廊和甘南在四個季節(jié)也都出現(xiàn)了顯著差異。

圖4 不同氣候帶顯著性水平的兩兩對比

4.4 分析結(jié)果

分析結(jié)果可以概括如下：

從季節(jié)來看,秋季氣溫受氣候帶的影響較小,一些地區(qū)秋季氣溫幾乎等于平均氣溫；春季氣溫受氣候帶的影響次之；夏季和冬季氣溫受氣候帶影響最大。在夏季，由于受到氣候帶的影響，甘南的氣溫要比全省平均氣溫高5℃左右，隴南的氣溫也要高于全省平均氣溫2～4℃。在冬季，除了甘南的其他地區(qū)氣溫均受到氣候帶的影響。隴南冬季的氣溫比平均氣溫高5℃左右，而河西走廊冬季氣溫則相反，比平均氣溫低5℃左右?？傮w來看，氣溫在秋季和春季的波動較小，在夏季和冬季的波動相對較大。

從空間變異來看,隴中地區(qū)全年氣溫受氣候帶的影響最小。除了冬季氣溫稍有波動，全年氣溫曲線幾乎與全省平均氣溫曲線重合。河西走廊地區(qū)次之，除了冬季氣溫有較大波動，全年氣溫與平均氣溫相似。隴南地區(qū)和甘南地區(qū)最大。由于受到氣候帶的影響,隴南地區(qū)和甘南全年氣溫較平均氣溫有較大的波動，而且全年氣溫差異性較大。總體來看,甘肅省西北部受氣候帶影響較大,東南部則較小。

5 結(jié)論

本文詳細(xì)介紹了函數(shù)型方差分析的參數(shù)估計和多重比較檢驗方法，并通過對實例進(jìn)行分析，給出了函數(shù)型方差分析的方法和結(jié)果。將多元方差分析中的多重比較檢驗應(yīng)用到函數(shù)型方差分析中，得到不同氣候帶兩兩之間顯著性差異的情況，從而確定甘肅省不同氣候帶對氣溫的具體影響，為甘肅省氣溫的研究提供參考。

[1]Wahba G,Klein B.Smoothing Spline ANOVA for Exponential Families,With Application to the Wisconsin Epidemiological Study of Diabetic Retinopathy[J].Annals of Statistics,1995,23(6).

[2]Stone C J,Hansen M H,Kooperberg C,et al.Polynomial Splines and their Tensor Products in Extended Linear Modeling[J].Annals of Statistics,1997,25(4).

[3]Huang J Z.Projection Estimation in Multiple Regression With Application to Functional Anova Models[J].Annals of Statistics,1998,26(1).

[4]Tallarida R J,Murray R B.Newman-Keuls Test[M].Manual of Pharmacologic Calculations,1987.

[5]Keselman H J,Burt H,Cribbie R A.Multiple Comparison Procedures[M].New York:Wiley,1987.

[6]Tuncer D.The Influence of Blood and/or Hemostatic Agent Contamination on Micro-TBS to Dentin[J].Journal of Adhesion Science&Technology,2015,29(13).

[7]Wang T L,Tseng Y K.Do Thinking Styles Matter for Science Achievement and Attitudes Toward Science Class in Male and Female Elementary School Students in Taiwan?[J].International Journal of Science&Mathematics Education,2015,(13).

[8]李棟梁.中國夏季月平均氣溫異常研究[J].高原氣象,1995,14(2).