周益林, 聶 曉, 劉 偉, 范潔茹, 梁革梅, 陸宴輝
(中國農業科學院植物保護研究所,植物病蟲害生物學國家重點實驗室, 北京 100193)
數理統計學是隨著人類社會發展和管理的需要而從應用統計學發展起來的一個分支學科,它是一門通過運用多種模型與技術分析、社會調查與統計分析來收集數據與處理數據的方法,對科技前沿、國民經濟中復雜或者重大的問題以及社會和政府中的大量問題進行定量分析的學科。此學科從卡爾·皮爾遜(Karl Pearson)[1]于20世紀初創立以來,作為一個實用性極強的學科已廣泛應用于社會各個領域,其在植物保護學科領域中的大量應用始于20世紀60年代,特別是70年代計算機及其相關工程技術的快速發展,大大促進了數理統計在植保學科特別是昆蟲生態學、植物病害流行學及病蟲害監測預警等領域中的廣泛應用。郭予元院士作為數理統計學在植保學科的應用研究先驅者之一,從20世紀50-80年代就開始在這方面做了大量工作,為數理統計在植物保護學科中應用研究以及教育普及工作作出了重大貢獻。
20世紀60年代郭予元院士在寧夏永寧農校工作期間,計算器/計算機還不普及,大量的數理統計的計算工作還需要依賴手工計算,其中相關分析是在病蟲害調查數據分析中經常使用的重要方法之一。簡單的相關系數比較容易計算,但多重相關分析由于公式復雜,其手工計算過程比較麻煩,為了簡化這一過程,郭予元院士根據分解線性回歸方程的原理,創造性地提出了一種緊湊計算法。這種方法不論有多少要分析的因素,都無需查找特殊的公式就可簡捷地把各種相關系數算出來[2]。在此基礎上郭予元院士于1985年又進一步提出了緊湊消去法解多元回歸聯立方程的計算方法,這兩種計算相關系數和回歸系數的緊湊消去法的特別之處在于把幾階矩陣的中間計算過程串聯在一起,減少了手工抄寫,大大提高了計算速度,節省了計算時間,而且涉及的因子越多,其優點更加突出,同時計算過程中便于檢驗,能避免出錯[3]。郭予元院士(1981年)通過對科技期刊上已發表的求致死中量的若干實例數據研究分析發現,在統計分析殺蟲劑的致死劑量中采用多項式配線精確度一般優于常用的機值分析法,特別是在數據點偏離劑量-死亡線的情況下,機值分析法誤差大且結果很不可靠。但是采用多項式配線的最大問題是運算過程太麻煩,遠遠比不上機值分析法簡便快捷,盡管英國著名統計學家羅納德·費希爾(Ronald Fisher)[4]在1938年提供的正交多項式表可加快系數變換后的配線速度,但是這種方法解決不了系數復原的簡化問題,進行高次方程的系數復原非常費事。郭予元院士根據正交多項式編制了系數復原表,使整個配線過程全部表格化,成功解決了這一問題[5-6]。同時他還通過進一步數據分析提出了以上兩種方法適用的數據特征和各自的局限性。
在20世紀80年代初,在比較高端、有統計功能的計算器如Casio fx-180p、Casio fx-3600p等還不普及的情況下,郭予元院士(1983年)就開發出基于單變量函數式普通電子計算器直線回歸分析程序,由于這類早期型號的計算器(如Casio fx-120)只帶有一個存儲器,還不能實現真正意義的統計學分析如直線回歸分析,他巧妙地利用這單一存儲器創造性地通過對該類型號計算器編程實現快速連續運算,約1 min可直接求出直線回歸的截距、斜率和回歸誤差[7]。這一絕妙的思路和方法成為他以后利用Casio fx-180p、Casio fx-3600p等型號計算器進行統計學分析方法編程的開端,由此可以看出郭予元院士深厚的數理統計功底和功夫。隨著Casio fx-180p系列型號計算器的推出和普及,郭予元院士基于此系列型號計算器開發出幾乎涵蓋植保學科領域常用的數理統計學方法或模型的編碼程序,如分組資料的平均數和方差計算、方差同質性分析、自相關和偏相關系數計算、機值法求LD50及卡方檢驗、Leslie矩陣最大特征值和特征向量、邏輯斯蒂曲線擬合及K值求法、3~5次方程求根、病蟲害各種空間分布型擬合(正態分布、二項分布、泊松分布、負二項分布、奈曼A分布、泊松-二項分布等)、諧波分析、生態位和病蟲危害指數等參數計算。從郭予元院士這些編碼程序可以看出,他把此型號系列計算器的2程序存儲區共38步可用的功能發揮到極致,只可惜這些程序編碼當時都沒有在公開期刊或雜志發表,只保留在郭先生的工作筆記和課程教案中,或者參加過他的數理統計培訓班或上過他的昆蟲生態學課程的學生或學員的筆記中有部分保存,這也是非常遺憾的事情,本文在此展示兩個例子(正態分布和諧波分析),以體會下這些編碼程序的精妙思路及它的美學之處。
例子1的正態分布由法國數學家亞伯拉罕·棣莫弗(Abraham De Moivre)于1733年首次提出,后由德國數學家Gauss率先將其應用于天文學研究,故正態分布又叫高斯分布,該分布是許多統計學方法的基礎,應用范圍廣泛。其公式如下:
上式中π為常數,μ為數學期望,σ2為方差,隨機變量x服從參數為μ、σ2的正態分布或高斯(Gauss)分布。

例子2的諧波分析在18世紀和19世紀已經奠定了良好的基礎,傅里葉等提出的諧波分析方法仍被廣泛應用。盡管此分析方法起源于電磁學領域,但也可用于對病蟲害時間序列的分析和模擬預測,滿足一定條件(Dirichlet條件)的、以T為周期的時間周期函數f(x),在連續點處,可用下述三角函數的線性組合(傅里葉級數)來表示:

上式稱為f(x)的傅里葉級數,a0為直流分量,r為諧波次數,ar和br為諧波幅度,L為最大諧波次數,x為等距值。
Casio fx-180p求a0,ar,br: 1 Kin 5 MODE 0 INV PCL P1Kout 6 M+ MR sin*ENT Kin+3 Kin 4=Kin+1 MR cos*Kout 4=Kin+2 INV RTN INV P2Kout 5÷÷Kout 2=INV HLT Kout 1=MODE·INV KAC INV Min 2πm/n Kin 6 INV M-n Kin 5 P1y1RUNy2RUN····AC INV P2→arRUN→brKout 3÷N→a0。
郭予元院士早期在數理統計學理論和方法方面的貢獻中其“正交多項式配線求殺蟲劑的致死中量”和“多元回歸分析的因子相關選擇法”兩項成果分別獲得了1983年寧夏技術改進獎二等獎和寧夏優秀科技成果獎三等獎。
20世紀50年代中期郭予元院士在寧夏回族自治區試驗農場做稻瘟病(病原菌Pyriculariaoryzae)發生規律、預測預報和防治研究的時候,就開始在試驗研究中應用相關的數理統計學方法。他從調查稻瘟病發生規律入手,在田間調查和病菌孢子捕捉技術收集數據基礎上,采用多元回歸分析方法對多年稻瘟病病情和氣象資料進行分析,通過對分析結果的實地驗證,證實關鍵時期雨量和雨日的乘積與稻瘟病發生程度極顯著相關,由此制定了預測稻瘟病病情和確定最佳防治時間的中長期預測模型,此模型可準確測報稻瘟病的發生,有效地指導對該病的防治。當時,采用數理統計方法定量化研究病蟲害規律并建立病蟲害預測模型的研究在國內是非常少見的,只可惜限于當時客觀條件,這些有重要意義的結果大多都沒能在相關學術期刊上正式發表[8]。

1982年底郭予元院士調入中國農業科學院植物保護研究所,此后在植物保護研究所工作的35年間,在他的引領下他和團隊的老師及學生使用數理統計的方法解決了不少植保學科數據分析中的難題、障礙和科學問題,取得許多豐碩的研究成果,為數理統計在植保學科中的應用作出了重大貢獻。此文鑒于篇幅所限,在此只介紹郭予元先生在害蟲生命表和種群動態及多病蟲危害產量損失和防治指標兩個方面的研究工作。
20世紀80年代開始郭予元院士和團隊的研究人員就開展了麥田黏蟲Mythimnaseparata生命表的研究,在對黏蟲發育階段致死因子K值和種群趨勢指數分析的基礎上,首次采用了歐式距離相似程度的方法分析其關鍵因子,篩選出了影響黏蟲種群的關鍵因子為蛹的被寄生和1~3齡幼蟲的被捕食,這2個因子是造成一代黏蟲年度間變動的主要原因。此研究采用的歐式距離分析法比前人常用的圖解法、b值分析法、r2值分析方法等所得結果更精確,還可避免采用逐步回歸分析存在的缺陷,即不能排除選中高負值r的因子可能性(這種情況是不符合生物學邏輯,即由此因子造成的死亡率越高,最終殘蟲率也越高,是不合理的)[11]。1989年郭予元院士和團隊的研究人員通過田間試驗首次建立了麥長管蚜Sitobionavenae(目前已更名為荻草谷網蚜Sitobionmiscanthi)田間自然種群的生命表,此研究在對麥長管蚜各代種群趨勢指數分析的基礎上,采用聚類分析中的距離系數方法對影響麥長管蚜種群的影響因子K值進行了分析,明確了影響種群波動的主要因素是天氣條件和天敵,低齡若蚜主要受風雨影響,高齡若蚜和成蚜則受天敵和風雨共同作用[12]。隨后還建立了棉鈴蟲Helicoverpaarmigera自然種群生命表,通過關鍵因子分析明確了影響全年棉鈴蟲種群數量變動的關鍵因子是2齡期被寄生,年度間關鍵因子是5齡期被捕食,明確了化學防治對種群I值的影響,給出了相應的化學防治策略[13]。1994年通過對田間麥無網長管蚜Metopolophiumdirhodum系統調查,分析了該蚜蟲的時空分布格局和動態變化,通過對生態位的定量分析,明確了麥無網長管蚜、禾谷縊管蚜Rhopalosiphumpadi和麥長管蚜的種間競爭關系,在用相關分析和通徑分析明確麥無網長管蚜種群消長關鍵因子的基礎上,通過多元回歸分析法建立了基于關鍵因子的預測模型[14]。
20世紀80年代開始郭予元院士以及團隊的科研人員還開展了棉鈴蟲對棉花的危害和防治指標的研究,組建了不同年份、不同肥料水平及不同代棉鈴蟲為害棉花的產量損失率模型,在這項研究過程中采用協方差分析成功地解決了模型之間的斜率和截距差異性比較,明確了不同模型的適用性范圍和穩定性差異,最終分別建立了適用于三種情況,即高產二代或中產二代或高產、中產、低產三代及低產二代的棉鈴蟲損失率估計模型[15]。1988年通過田間自然蟲源輔助人工控制培養研究了4種麥蚜混合種群對小麥產量的影響,此研究采用偏相關分析方法,成功地解決了穗數對小麥產量的影響及不同蚜種對小麥產量的影響問題,從而明確了4種蚜蟲混合種群中麥長管蚜和禾谷縊管蚜對產量有明顯影響,并通過多元回歸分析方法獲得了復合蚜量的動態防治指標[16]。郭予元院士等于1994年通過田間開放式小區試驗,巧妙地采用了回歸最優設計的“206”田間試驗設計,對河南省二代棉鈴蟲和葉螨(朱砂葉螨Tetranychuscinnabarinus和截型葉螨T.truncatus的復合種群)為害棉花的復合防治指標進行了研究,在數據分析中用偏相關分析發現黃萎病對葉螨有一定干擾作用,因此通過協方差方法排除了黃萎病的干擾,并用主成分分析明確了棉鈴蟲和葉螨的一次和二次型數值,尤其是二者的交互作用,建立了棉鈴蟲和葉螨復合為害的二次曲線方程,并組建了一系列復合防治指標[17]。另外在棉鈴蟲預測預報研究方面,在華北棉區用掃網法和目測查蟲法調查麥田一代棉鈴蟲幼蟲發生,用比較法和生命表分析法預測二代棉鈴蟲在棉田的發生程度和蟲口密度,經1991年-1994年應用檢驗,其預測結果均與實際發生情況吻合[18]。
郭予元院士值得一書的還有他在植保數理統計學的教學和普及方面的重要貢獻。多年來他通過在全國各地舉辦植保數理統計培訓班和研究生授課,大大提升了我國從事植保學科的一代又一代科技人員在試驗設計、科學地統計和分析試驗數據的水平,解決了一些科研和生產及相關的數理統計中的重大問題和復雜問題,為科學地預測和防治病蟲害以及提高我國植保的科學水平都起到重要的作用。
20世紀80年代初郭先生在實際研究和工作中發現,由于我國教育體系受文化大革命的影響,植保領域的研究人員特別是從事病蟲害預測預報等方面的年輕工作者普遍存在數理統計知識薄弱的問題,嚴重制約植保學科領域科學研究的發展和研究水平的提升。他把自己多年來堅持自學的數理統計知識和經驗總結整理出適用于植保科學試驗研究的一套統計分析方法—《植保數理統計學方法》,并自己刻板油印,先后在我國7個省20多個農業科研和推廣單位舉辦多次數理統計講習班,不但使當時許多植保科研和推廣人員的研究論文數據分析的科學性得到明顯提高,更重要的是培養了一大批懂數理統計的植保領域的專家和學者,這些人員至今還活躍在各自的研究領域并發揮著重要的作用。
郭予元院士油印本的《植保數理統計學方法》內容廣泛,不僅涵蓋了當時植保學科領域所有能用到的基本統計學方法,如回歸分析、相關分析和通徑分析、曲線方程擬合、病情發生趨勢預測方法、列聯表分析法、聚類分析和判別分析、t檢驗和方差分析、正交試驗、卡方分析、協方差分析、空間分布型和抽樣技術、病蟲危害特性和防治指標等,還包括了當時一些剛剛出現的數學或統計新興領域的方法如模糊數學、灰色理論等,還有一些不太常見的統計學方法或當時比較新的田間試驗設計方法如題總相關分析、嶺回歸方法、灰色突變長期預測模型、回歸正交旋轉設計等也在培訓班中多有涉及。一些學員或學生在學習后就很快在自己研究中應用,如有學員1989年在小麥白粉病(病原菌:Blumeriagraminisf.sp.tritici)田間藥劑防治試驗中就采用了回歸正交旋轉設計,這種試驗設計的優點是試驗處理少,獲得的信息多,精確度高,而且計算簡便,從藥劑試驗結果看,它不但能獲得每種藥劑的防效大小和藥劑之間交互作用的有無,更重要的是可以獲得每種藥劑用量與病情之間的關系模型,從而在防治時可根據實際情況,計算出最適用藥量[19]。
郭予元院士授課通俗實用,所用的實例都是植保領域田間試驗或病蟲害發生流行的相關數據實例,每次上課都根據學員的實際情況舉例,而且在課堂上實時采用Casio fx-3600p、Casio fx-180p計算器編程,進行實際操作計算,使學員能學到數理統計的真本領,尤其重要的是他在授課時不但教各種統計學方法的使用,還特別強調使用時需要注意的問題,為此他還專門整理油印了相關的材料如《病蟲測報中的幾個問題》《棉蟲試驗的幾種統計分析方法》等供學員學習。比如在《病蟲測報中的幾個問題》中就強調了預報因子在選擇時不能只憑生物學邏輯或簡單相關系數挑選因子,一定還要注意因子之間的自相關性和主次關系;在對有重復使用的數據進行回歸分析時,不但要檢驗回歸的顯著性,還要檢驗重復之間失擬的顯著性;在《棉蟲試驗的幾種統計分析方法》中強調了數據轉換的問題以及數據轉換與數據分布型的關系,例如在對棉蚜Aphisgossypii田間種群量數據進行方差分析時,應先把數據轉換成對數值,這是因為棉蚜在田間是聚集分布的。統計分析數據轉換和數據分布型在當前依然是植保數理統計中的薄弱環節或被忽視的問題。
郭予元院士不但通過舉辦培訓班傳授數理統計知識,國內植保同行或學生也經常通過信件向他請教相關問題,郭先生都是有問必答,甚至親自幫助不少老師處理和分析他們的試驗數據。如1986年前后新疆植保站的榮麗君老師就數次給郭先生寫信請教多年度害蟲為害產量損失數據的分析和模型的組建問題,郭先生都一一回信做了解答,而且還親自對他的試驗數據進行了統計分析,通過采用協方差的分析方法,解決了黃地老虎Agrotissegetum、冬麥地老虎Caradrinaauguroides、麥長管蚜、棉鈴蟲和牧草盲蝽Lyguspratensis等5種害蟲不同年份模型差異性的檢驗問題,并建立了適用不同情況或條件下相關蟲害為害產量損失的損失估計模型,相關的文章也在學術期刊上及時發表[20-22]。
郭予元院士10多年來一直擔任中國農業科學院研究生院《昆蟲生態學》課程的主講老師,他的課程授課內容豐富、涵蓋范圍廣,包括昆蟲種群的時間動態、空間模式、種間關系以及群落生態系統等諸多方面,其課程內容的最大特點是每個講授內容或章節都有實際的例子,而且如果涉及數理統計內容,都有課堂實時分析計算,使同學們不但學到了昆蟲生態學的理論知識,而且還學到數理統計的實戰方法,并在自己的科研中很快就能得到應用。特別是在講授種群的空間圖式(正態分布、泊松分布、均勻分布、二項分布和負二項分布、奈曼A型分布、泊松-二項分布等P-E分布型)、種群增長模型(邏輯斯蒂曲線、指數增長曲線、Gompertz曲線、冪函數曲線、差數指數曲線、差數冪函數曲線、雙曲線、麥肯齊模型等)、種間關系(競爭模型、捕食模型等)以及這些模型的擬合和檢驗時,郭先生都有實際的植保例子供學習和分析,最有特色的是在課堂上用他開發的基于計算器編程的程序進行實時計算,這種教學方法和模式,大大提升了同學們昆蟲生態學知識和相應的數理統計學知識水平。
郭予元院士在植物保護研究領域特別是發展棉鈴蟲綜合防治策略和技術體系、有效解決生產抗性棉鈴蟲的種群治理問題和促進農作物病蟲害綜合防治理論的創新等方面做出重要貢獻,在他的帶動和引領下,年輕的一代植保科研工作者,尤其是研究團隊繼承和發揚他們這一輩老科學家思想和精神,青出于藍而勝于藍,在棉花害蟲研究領域做出了世界水平研究成果,在這些科學研究中也大量地使用了數理數據統計的方法,其研究結果先后在《科學》《自然》等國際著名期刊發表。例如棉鈴蟲研究團隊通過對我國北方6省過去10年100多個樣點棉鈴蟲在多種作物上卵(多代)和幼蟲發生數據的回歸分析發現,我國華北地區大規模種植Bt棉,不僅降低了棉花上棉鈴蟲的數量,而且減輕了周邊其他非轉基因農作物(玉米、大豆和花生)上的棉鈴蟲為害程度。這項研究結果是世界首個針對轉基因作物抗蟲性進行的大規模長期跟蹤研究,該論文曾作為2008年《科學》雜志的封面文章[23];進一步的研究工作又發現雖然大范圍種植Bt棉大大降低了棉鈴蟲發生和危害,但另一類害蟲盲蝽種群在棉花中快速上升,通過對長期多點收集的數據進行模型擬合,結果顯示,Bt棉種植比例與該害蟲的種群密度呈線性正相關,與盲蝽殺蟲劑的使用次數呈非線性指數增長關系,通過對種群密度與多因子包括殺蟲劑的使用次數、主要氣象因子溫度、降雨等的回歸和相關分析表明,該害蟲的種群密度增長主要是種植Bt棉帶來的棉鈴蟲殺蟲劑打藥次數減少造成的。另外研究團隊對我國華北大范圍種植Bt棉后該害蟲在棉花、蘋果、葡萄、桃、梨、棗等作物上的發生危害嚴重度多年多點發生數據的模型模擬發現,其危害嚴重度與Bt棉的種植比例也呈指數增長關系,該研究論文2010年再次登上《科學》雜志[24]。2012年研究團隊又闡明了Bt棉種植后由于打藥次數減少導致的天敵昆蟲種群和棉蚜種群等的變化趨勢,論文發表于《自然》雜志[25],同年該研究入選“中國科學十大進展”。這篇論文在數據分析中也采用了數理統計的多種方法如多因子方差分析、逐步回歸、相關分析、直線或曲線模型擬合等方法。
郭予元院士不僅在植保科研上做出了巨大的貢獻,而且在數理統計學的理論和方法研究、數理統計在植保研究領域的應用以及數理統計學知識的教學普及和推廣方面做出了重要的貢獻。近年來隨著數理統計學本身的發展,如廣義線性模型、加性模型、混合效應模型、隨機系數模型等發展,同時一些重要數理支持軟件的發展和改進如SAS、SPSS等,特別是開源軟件R軟件的快速發展,隨著大數據、機器學習、AI人工智能技術時代的來臨,相信郭予元院士奠基和引領的植保數理統計學必將在植保領域發揮越來越重要的作用,并產出更多更大的研究成果。