馬尾松主要蛀干害蟲林間誘捕序列的似乎不相關分析

陳繪畫 張汝忠

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

徐志宏

(浙江農林大學)

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

1 材料與方法

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

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

1.2 似乎不相關線性回歸分析

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

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

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

令 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個林間種群的監(jiān)測序列{Ni0、Ni1、Ni2、…、Niq，i=1、2、…、n}，則有:

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

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

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

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

亦即:

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

也就是E(e)=0。

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

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

2 結果與分析

2.1 模型自變量的確定

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

2.2 似乎不相關線性回歸模型的建立

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

分別將使用似乎不相關線性回歸方法和最小二乘法擬合的參數(shù)結果列于表2。由表2的結果看出，短角幽天牛林間種群數(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)成正比。最小二乘法的回歸結果與似乎不相關的一致，差別僅體現(xiàn)在回歸系數(shù)的不同。

2.3 似乎不相關線性回歸模型的檢驗

2.3.1 是否為對角陣的檢驗

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

表1 4個分模型間殘差的相關系數(shù)

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

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

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

表2 似乎不相關和最小二乘法回歸結果

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

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

2.3.4 殘差分析

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

分模型殘差的獨立性檢驗:殘差的獨立性檢驗通常采用Durbin－Watson方法，其計算公式為:

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

2.4 似乎不相關法與最小二乘法估計結果的差異

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

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

3 結論與討論

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

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

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

［1］張建軍，張潤志，陳京元.松材線蟲媒介昆蟲種類及其擴散能力［J］.浙江林學院學報，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］池田俊彌.マツのマダラカシキリの誘引物質とその利用［J］.森林防疫，1986，6(11):95 －99.

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

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

［8］崔相富，陳繪畫，趙錦年，等.括蒼山北麓松林間鞘翅目主要昆蟲種群動態(tài)研究［J］.林業(yè)科學研究，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ù)學雜志，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)非線性建模理論與應用［J］.生物數(shù)學學報，2001，16(4):439 －447.

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

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

［17］傅軍，丁晶，鄧育仁.嘉陵江流域形態(tài)及流量過程分維研究［J］.成都科技大學學報，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］梁洪川.似乎不相關回歸模型及其在醫(yī)學中的應用［D］.太原:山西醫(yī)科大學，2006.

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

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

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