陳繪畫 張汝忠

(浙江省仙居縣林業(yè)局，仙居，317300)

徐志宏

(浙江農(nóng)林大學(xué))

松墨天牛(Monochamus alternatus Hope)、短角幽天牛(Spondylis buprestoides(L.))等是我國馬尾松(Pinus massoniana)林的主要鉆蛀害蟲，同時(shí)還能攜帶松材線蟲［1］，尤其是松墨天牛為我國松材線蟲病最重要的傳播媒介。由于蛀干類害蟲成蟲體壁及鞘翅堅(jiān)厚、耐藥力較強(qiáng)、羽化期不整齊，在日本人發(fā)現(xiàn)從松樹伐倒木上收集的揮發(fā)性物質(zhì)對松墨天牛成蟲具有強(qiáng)烈的引誘作用后［2－5］，用引誘劑誘殺蛀干類害蟲的方法便得到較快地推廣［6－8］。蛀干類害蟲引誘劑在誘殺林間松墨天牛、短角幽天牛、立毛角脛象(Shirahoshizo erectus Chen)和松瘤象(Hyposipalus gigas Linnaeus)等松樹主要蛀干害蟲成蟲的同時(shí)，連續(xù)使用還會顯著降低林間蛀干類害蟲的蟲口密度［9］，并且可以連續(xù)監(jiān)測這些蛀干害蟲的發(fā)生時(shí)間、發(fā)生量和發(fā)生規(guī)律，為監(jiān)測、預(yù)報(bào)、控制松樹主要蛀干害蟲提供直接、客觀的依據(jù)，節(jié)省人力、物力。筆者利用連續(xù)的林間松墨天牛、短角幽天牛、立毛角脛象和松瘤象蛀干害蟲誘捕資料，用似乎不相關(guān)線性回歸(SUR)方法分析林間各害蟲種群的變化及其相互關(guān)系，為及時(shí)、全面、客觀地了解其林間動態(tài)和監(jiān)測預(yù)報(bào)提供技術(shù)支撐。

1 材料與方法

1.1 馬尾松蛀干類主要害蟲種群數(shù)量動態(tài)觀測

蛀干類害蟲引誘劑(原稱M－99)，規(guī)格為300 mL/瓶，以及配套的小型折疊式誘捕器由浙江省林業(yè)有害生物防治檢疫局、中國林業(yè)科學(xué)研究院亞熱帶林業(yè)研究所研制，寧波中化化學(xué)品有限公司生產(chǎn)。試驗(yàn)于2006—2010年每年4月至11月進(jìn)行。試驗(yàn)地設(shè)在浙江省仙居縣大北地溪林場馬尾松人工純林內(nèi)，面積 40.2 hm2，海拔66 m，林齡59 a，平均樹高約31 m，平均胸徑約22 cm。將誘捕器架在距地面1.5 m、兩相鄰松樹間的鉛絲上，共設(shè)誘捕器22只。每3—4周更換1次引誘劑，每周收蟲1次。在室內(nèi)進(jìn)行分類鑒定，列出種類名稱及數(shù)量，共得觀測數(shù)據(jù)131個(gè)。

1.2 似乎不相關(guān)線性回歸分析

似乎不相關(guān)回歸也稱半相依回歸、似不相關(guān)回歸，簡記為SUR［10］。SUR是一種特殊的聯(lián)立方程模型，模型里只有方程間的誤差相關(guān)而沒有聯(lián)立性，同時(shí)各方程影響因變量的自變量數(shù)量可以不同，由Zellner［11－12］于 1962 年首先提出。

設(shè){N0、N1、N2、…、Nq}為某林間種群監(jiān)測序列，描述單個(gè)監(jiān)測序列動態(tài)最簡單的模型為離散差分方程，其3階滯后的一般形式為:

式中:Nt為t時(shí)誘集到的某種群密度;t是監(jiān)測次數(shù)(t=0、1、2、…、q);f是描述某次誘集到的某種群密度影響下次誘集密度的函數(shù);Et是均值為0的隨機(jī)誤差。

令 Xt=logNt－1，Yt=logNt－2，Zt=logNt－3，rt=log(Nt/Nt－1)，則(1)式化為線性模型為:

式中:a、b、c、d、e、f、g、h、j、k 是待定系數(shù)。

對于n個(gè)林間種群的監(jiān)測序列{Ni0、Ni1、Ni2、…、Niq，i=1、2、…、n}，則有:

式中:i=1、2、…、n。

將(4)式化為一般形式，則變?yōu)閚個(gè)似乎不相關(guān)方程:

式中:yi為(q×1)向量，對應(yīng)于yi的q個(gè)觀測值;xi為(q×ki)設(shè)計(jì)矩陣;βi為(ki×1)未知參數(shù)向量;(ei=e1i，e2i，…，eqi)為(q ×1)隨機(jī)誤差向量。

將(5)式寫成矩陣形式:

亦即:

假定因變量y的q個(gè)觀測間相互獨(dú)立，各方程的誤差項(xiàng) ei間不獨(dú)立，則有:E(ei)=0，E(ei，ei)=σiiIq，E(ei，ej)= σijIq，i、j=1、2、…、n。

也就是E(e)=0。

式中:V為模型誤差項(xiàng)的方差——協(xié)方差矩陣;∑是同一觀測點(diǎn)誤差項(xiàng)的方差——協(xié)方差矩陣;?是叉乘(克羅內(nèi)克乘積)符號。從形式上看，模型的n個(gè)方程為各自獨(dú)立的自變量所解釋，各因變量間似乎是不相關(guān)的。不過n個(gè)方程的同一觀測的模型誤差項(xiàng)ei的各分量相關(guān)，也就是∑為非對角陣，說明這n個(gè)方程實(shí)際上是通過同一觀測的模型誤差項(xiàng)ei的各分量相互聯(lián)系的。

似乎不相關(guān)回歸模型參數(shù)估計(jì)的方法及步驟詳見文獻(xiàn)［15］，［16］106－111。文中的計(jì)算及圖形繪制均在MATLAB上完成。

2 結(jié)果與分析

2.1 模型自變量的確定

對2006—2010年松墨天牛、短角幽天牛、立毛角脛象和松瘤象4種蛀干害蟲誘捕數(shù)序列的年際間連接做近似處理，也就是連接各蛀干害蟲的羽化期，將非羽化期的空間軌線概化為一點(diǎn)［17］。使用偏自相關(guān)函數(shù)法可以解決模型的滯后代數(shù)難以確定的問題［18］，由 Bartlett的 2 /法則確定短角幽天牛、松墨天牛、松瘤象和立毛角脛象4種馬尾松主要蛀干害蟲林間監(jiān)測序列模型自變量的最大滯后代數(shù)為3代。根據(jù)響應(yīng)面建模原理，對于生態(tài)時(shí)間序列在確定自變量滯后代數(shù)后，通常取(2)式那樣的滯后2次型［14］，這樣以4種蛀干害蟲中任一種群為因變量、該種群自身及其余3種害蟲種群作自變量，4種蛀干害蟲任一種群模型有自變量36個(gè)，為防止自變量間的作用重疊，采用逐步回歸法選取自變量，短角幽天牛最終入選的自變量為自身滯后3步序列的觀測值(x13)、松墨天牛滯后1步序列觀測值的平方(x24);松墨天牛為自身滯后1步序列(x21)、自身滯后3步序列(x23)、自身滯后1步序列的平方(x24)、立毛角脛象滯后1步序列的平方(x44);松瘤象為自身滯后1步序列(x31)、自身滯后2步序列的平方(x35);立毛角脛象為短角幽天牛滯后3步序列的觀測值(x13)、自身滯后1步序列(x41)、自身滯后1步序列的平方(x44)。為使自變量序列的長度與因變量序列長度相等，將某種群全部自變量的平均值作為該種群自變量序列的第1個(gè)變量值、第2個(gè)變量值或第3個(gè)變量值。

2.2 似乎不相關(guān)線性回歸模型的建立

表1列出了分別用最小二乘法對短角幽天牛、松墨天牛、松瘤象和立毛角脛象4種林間種群模型進(jìn)行獨(dú)立估計(jì)時(shí)4個(gè)模型間殘差的相關(guān)系數(shù)矩陣。根據(jù)表1數(shù)據(jù)可知，其中3個(gè)分模型的殘差在置信水平為0.01下相關(guān)性顯著，表明短角幽天牛、松墨天牛、松瘤象和立毛角脛象4種害蟲林間種群存在著一定的相關(guān)性，即:cov(ei，ei)≠0，因而對這4個(gè)方程的參數(shù)進(jìn)行似乎不相關(guān)聯(lián)立估計(jì)能夠提高參數(shù)估計(jì)的有效性和一致性。

分別將使用似乎不相關(guān)線性回歸方法和最小二乘法擬合的參數(shù)結(jié)果列于表2。由表2的結(jié)果看出，短角幽天牛林間種群數(shù)量的某次觀測值隨自身前3次的觀測值(x13)或松墨天牛前1次觀測值的平方(x24)的增加而降低。松墨天牛林間種群數(shù)量與自身前1次的觀測值(x21)成反比，與自身前3次的觀測值(x23)、自身前1次觀測值的平方(x24)或立毛角脛象前1次觀測值的平方(x44)成正比。松瘤象林間種群數(shù)量則隨自身前1次觀測值(x31)的增加而減少、隨自身前2次觀測值的平方(x35)的增加而增加;立毛角脛象林間種群數(shù)量與短角幽天牛前3次的觀測值(x13)、自身前1次的觀測值(x41)成反比，與自身前1次觀測值的平方(x44)成正比。最小二乘法的回歸結(jié)果與似乎不相關(guān)的一致，差別僅體現(xiàn)在回歸系數(shù)的不同。

2.3 似乎不相關(guān)線性回歸模型的檢驗(yàn)

2.3.1 是否為對角陣的檢驗(yàn)

為進(jìn)一步確認(rèn)表1的模型殘差間相關(guān)系數(shù)所揭示的可以用似乎不相關(guān)方法對4種蛀干害蟲林間種群模型進(jìn)行聯(lián)立估計(jì)，用模型間的方差—協(xié)方差矩陣∑是否為對角陣來判斷。若∑為對角陣，模型中參數(shù)的最小二乘估計(jì)量與廣義最小二乘估計(jì)量相等，而最小二乘估計(jì)量就是模型的最優(yōu)無偏估計(jì)，∑是否為對角陣的具體判定方法見文獻(xiàn)［16］110，［19］。經(jīng)計(jì)算，拉格朗日乘子統(tǒng)計(jì)量λ=76.97＞=9.49，因此拒絕∑為對角陣的假設(shè)。

表1 4個(gè)分模型間殘差的相關(guān)系數(shù)

2.3.2 模型擬合優(yōu)度的檢驗(yàn)

模型對數(shù)據(jù)的整體擬合評價(jià)用確定系數(shù)R2來判定，其計(jì)算公式為［19］:

經(jīng)計(jì)算，R2=0.295 2，因此整體模型的相關(guān)系數(shù)在置信水平為0.01下顯著。同時(shí)對4種蛀干害蟲的分模型進(jìn)行F檢驗(yàn)，結(jié)果也表明極顯著(表2)。

表2 似乎不相關(guān)和最小二乘法回歸結(jié)果

2.3.3 模型回歸系數(shù)的檢驗(yàn)

建立的似乎不相關(guān)模型通過了確定系數(shù)顯著性檢驗(yàn)并不代表每個(gè)回歸方程的自變量Xi都對相應(yīng)的因變量Yi有顯著影響，也許模型中某個(gè)方程的自變量Xi對相應(yīng)的因變量Yi影響不顯著，因而需作模型回歸系數(shù)的顯著性t檢驗(yàn)［20］40，t檢驗(yàn)的結(jié)果見表2。從表2可知，用似乎不相關(guān)聯(lián)立估計(jì)法得到的各方程回歸系數(shù)均呈極顯著差異，說明各回歸方程均拒絕回歸系數(shù)bi=0的假設(shè)，4種蛀干害蟲林間種群模型的回歸系數(shù)均通過t檢驗(yàn)。

2.3.4 殘差分析

殘差的正態(tài)性檢驗(yàn):對整體模型的殘差進(jìn)行四元正態(tài)分布以及4個(gè)分模型的殘差進(jìn)行兩兩配對的二元正態(tài)分布檢驗(yàn)［21］表明，各模型的殘差均符合四元正態(tài)分布或二元正態(tài)分布，檢驗(yàn)方法及評判標(biāo)準(zhǔn)見文獻(xiàn)［22］406－408。

分模型殘差的獨(dú)立性檢驗(yàn):殘差的獨(dú)立性檢驗(yàn)通常采用Durbin－Watson方法，其計(jì)算公式為:

計(jì)算結(jié)果為 D短角幽天牛=2.37、D松墨天牛=1.93、D松瘤象=1.92、D立毛角脛象=1.94，各分模型的 D 值約等于2，故各分模型的殘差不存在一階自相關(guān)，也不存在多階自相關(guān)，評判標(biāo)準(zhǔn)見文獻(xiàn)［20］46－47，［22］364。對各分模型的殘差進(jìn)行自相關(guān)計(jì)算，其滯后1階自相關(guān)系數(shù)和2階自相關(guān)系數(shù)均顯示相關(guān)不顯著，這與Durbin－Watson方法的檢驗(yàn)結(jié)果相符。

2.4 似乎不相關(guān)法與最小二乘法估計(jì)結(jié)果的差異

從表2看出，似乎不相關(guān)聯(lián)立估計(jì)的各自變量回歸系數(shù)與最小二乘參數(shù)估計(jì)值間均存在偏差，似乎不相關(guān)回歸聯(lián)立估計(jì)參數(shù)的標(biāo)準(zhǔn)誤均小于最小二乘估計(jì)，因此似乎不相關(guān)聯(lián)立估計(jì)的精度高于最小二乘估計(jì)?？疾?個(gè)分模型間因變量的相關(guān)系數(shù)(表3)，顯示各方程間不獨(dú)立，未將方程間的相關(guān)信息納入考慮范圍正是最小二乘法估計(jì)精度低的原因，同時(shí)似乎不相關(guān)模型的殘差極差也小于最小二乘法，由此說明用似乎不相關(guān)法對4種蛀干害蟲種群的監(jiān)測序列進(jìn)行建模，所建立的模型在盡可能地反映4種害蟲種群監(jiān)測序列各自的動態(tài)變化及相互間影響的同時(shí)，也取得了較為滿意的擬合精度。

表3 4個(gè)分模型因變量間的相關(guān)系數(shù)

3 結(jié)論與討論

似乎不相關(guān)是一種特殊的聯(lián)立方程模型，它以模型方程間的誤差相關(guān)性為前提，在模型的擬合過程中根據(jù)方程間的誤差調(diào)整各方程的回歸系數(shù)，使模型的整體誤差趨于最小。通過在林間懸掛誘捕器連續(xù)監(jiān)測馬尾松主要蛀干害蟲種群動態(tài)，由逐步回歸法分別篩選出影響各蛀干害蟲種群動態(tài)變化的因子，然后用似乎不相關(guān)法進(jìn)行聯(lián)合估計(jì)模型的各回歸系數(shù)，所建立的模型既反映各蛀干害蟲種群自身林間的動態(tài)變化，又根據(jù)模型各方程間的誤差調(diào)整模型參數(shù)，從而間接地反映4種害蟲種群間的相互影響與相互聯(lián)系，并使模型的整體殘差達(dá)到最小，擬合結(jié)果令人較為滿意。

由于似乎不相關(guān)法聯(lián)合估計(jì)模型回歸參數(shù)的過程中考慮了方程間的相關(guān)信息，故比傳統(tǒng)的最小二乘法估計(jì)模型參數(shù)更合理。不過當(dāng)各方程的因變量間相關(guān)性不大(即∑為對角陣)時(shí)，利用似乎不相關(guān)法的回歸估計(jì)結(jié)果與最小二乘法的估計(jì)結(jié)果相同;在樣本數(shù)量較少的情況下，有時(shí)還可能出現(xiàn)最小二乘的估計(jì)結(jié)果優(yōu)于似乎不相關(guān)線性回歸的估計(jì)結(jié)果，尤其是當(dāng)方程間的殘差相關(guān)性較弱時(shí)更可能出現(xiàn)上述情況，因此需要在使用似乎不相關(guān)線性回歸估計(jì)前進(jìn)行Σ是否為對角陣的檢驗(yàn)，經(jīng)檢驗(yàn)后再決定采用哪種估計(jì)方法［16］124。

似乎不相關(guān)法允許模型內(nèi)各方程的自變量數(shù)量不同，因而它比多元多重(也稱多對多)線性回歸方法具有更大的靈活性;另外，似乎不相關(guān)法運(yùn)用方程間殘差的相關(guān)信息矩陣來提高參數(shù)的估計(jì)效率，構(gòu)成了似乎不相關(guān)線性回歸法的兩個(gè)主要特點(diǎn)。在昆蟲種群生態(tài)中非獨(dú)立的數(shù)據(jù)是非常普遍的，忽略觀測數(shù)據(jù)間的相關(guān)性會使統(tǒng)計(jì)結(jié)果偏離真實(shí)情況，這為似乎不相關(guān)線性回歸方法的應(yīng)用提供了廣闊的舞臺。

［1］張建軍，張潤志，陳京元.松材線蟲媒介昆蟲種類及其擴(kuò)散能力［J］.浙江林學(xué)院學(xué)報(bào)，2007，24(3):350 －356.

［2］Ikeda T，Oda K，Yamane A，et al.Volatiles from pine logs as the attractant for the Japanese pine sawyer，Monochamus alternatus Hope(Coleoptera:Cerambycidae)［J］.Jap For Soc，1980，62(4):150－152.

［3］Ikeda T，Enda N，Yamane A，et al.Attractants for the Japanese pine sawyer，Monochamus alternatus Hope(Coleoptera:Cerambycidae)［J］.Appl Ent Zool，1980，15(3):358 －361.

［4］Ikeda T，Yamane A，Enda N，et al.Attractiveness of chemicaltreated pine trees for Monochamus alternatus Hope(Coleoptera:Cerambycidae)［J］.J Jap For Soc，1981，63(6):201 －207.

［5］池田俊彌.マツのマダラカシキリの誘引物質(zhì)とその利用［J］.森林防疫，1986，6(11):95 －99.

［6］王玉嬿，舒超然，孫永春.松褐天牛引誘試驗(yàn)初報(bào)［J］.林業(yè)科學(xué)，1991，27(2):186 －189.

［7］趙錦年，林長春，姜禮元，等.M99－1引誘劑誘捕松墨天牛等松甲蟲的研究［J］.林業(yè)科學(xué)研究，2001，14(5):523 －529.

［8］崔相富，陳繪畫，趙錦年，等.括蒼山北麓松林間鞘翅目主要昆蟲種群動態(tài)研究［J］.林業(yè)科學(xué)研究，2008，21(3):340 －345.

［9］Yoichi K.The pine wood nematode and the Japanese pine sawyer［M］.Tokyo:Thomas Co Ltd，1995:208－209.

［10］張莉莉，史建紅.半相依線性回歸模型的影響分析［J］.數(shù)學(xué)雜志，2010，30(1):137 －144.

［11］Zellner A.An efficient method of estimating seemingly unrelated regressions and tests for aggregation bias［J］.J Am Statist Assoc，1962，57:348 －368.

［12］Zellner A.Estimators for seemingly unrelated regression equations:some exact finite sample results［J］.J Am Statist Assoc，1963，58:977 －992.

［13］Turchin P，Taylor A D.Complex dynamics in ecological time series［J］.Ecology，1992，73(1):289 －305.

［14］趙中華，沈佐銳.昆蟲種群動態(tài)非線性建模理論與應(yīng)用［J］.生物數(shù)學(xué)學(xué)報(bào)，2001，16(4):439 －447.

［15］梁洪川，韓宏，郎素平，等.似乎不相關(guān)回歸模型及其在老年認(rèn)知功能研究中的應(yīng)用［J］.中國衛(wèi)生統(tǒng)計(jì)，2005，22(6):362－364，372.

［16］唐守正，李勇.生物數(shù)學(xué)模型的統(tǒng)計(jì)學(xué)基礎(chǔ)［M］.北京:科學(xué)出版社，2002.

［17］傅軍，丁晶，鄧育仁.嘉陵江流域形態(tài)及流量過程分維研究［J］.成都科技大學(xué)學(xué)報(bào)，1995(1):74 －80.

［18］Berryman A，Turchin P.Identifying the density-dependent structure underlying ecological time series［J］.Oikos，2001，92(2):265－270.

［19］梁洪川.似乎不相關(guān)回歸模型及其在醫(yī)學(xué)中的應(yīng)用［D］.太原:山西醫(yī)科大學(xué)，2006.

［20］向東進(jìn)，李宏偉，劉小雅.實(shí)用多元統(tǒng)計(jì)分析［M］.武漢:中國地質(zhì)大學(xué)出版社，2005.

［21］Dennis B，Kemp W P，Taper M L.Joint density dependence［J］.Ecology，1998，79(2):426 －441.

［22］蘇金明，阮沈勇.MATLAB 6.1實(shí)用指南:下冊［M］.北京:電子工業(yè)出版社，2002.