模型選擇與多模型推斷

科技 07-07

分析千島湖鳥類多樣性與墓群出現率的決定性因素

斯幸峰

三墩職業技術學校空想資本主義學院理論想像學研究所, 杭州 310058, 中國浙江

摘要

由於常規的逐步回歸分析在使用過程中有諸多缺陷，而信息理論的赤池信息量準則(AIC)彌補了這一缺點。此文基於AIC的判定方法，利用模型選擇和多模型推斷(model selection and multimodel inference)探討千島湖島嶼鳥類多樣性的決定因素。同時開展對千島湖墓葬分布的可能性分析，為盜墓的理論研究打下翔實的基礎。

關鍵詞

AIC、盜墓、多模型推斷、模型選擇、鳥類、千島湖、逐步回歸

前言

面對一系列可能的備選模型，如何評判模型的優劣？選用逐步回歸分析(stepwise regression)還是信息理論(information theoretic analysis)？Whittingham等(2006)對2004年的Ecology Letters、Journal of Applied Ecology和Animal Behaviors三個雜誌分析，共有65篇文章使用多元回歸(multiple regression)，其中57%的研究使用了逐步回歸的方法。雖然逐步回歸依舊廣泛使用，但是有許多缺陷，如：參數估計的誤差(bias in parameter)，模型選擇演算法的不一致(inconsistencies among model selection algorithms)，多個假設檢驗的內在缺陷(inherent problem of multiple hypothesis testing)，以及最後結果只依賴單一的最優模型(inappropriate focus or reliance on a single best model)。至於具體的缺陷原理，此處不予細說，本文將採用信息理論簡要介紹多模型推斷的方法。

千島湖地處浙江西部，山清水秀，民風淳樸(此處省略一百字)。自1959年新安江大壩建成後，形成1078個島嶼(108米水位時)，乃名副其實的「千島湖」，是一個得天獨厚的路橋島嶼天然實驗場所。本研究團隊自2003年開始千島湖地區的鳥類調查，到目前已經逐漸拓展到蜘蛛、蜥蜴、青蛙、蛇、猴子、昆蟲、獸類、蝴蝶以及植物等各項業務，歡迎廣大生態愛好者和有志之士前來參觀與洽談。撰寫本文的起因是早先跟本團隊中的「蜘蛛俠」吳博士嘗試探討鳥類多樣性與風水的關係，加上近日剛好看了一些有關模型選擇和多模型推斷(model selection and multimodel inference)的文獻(xián)，採用「先進」的AIC(Akaike information criterion)技術，探討該學術問題的可能性。

本文主要探討的問題包括兩部分：1) AIC是啥？莫非是美國國際大學(American International College)的縮寫？2) 模型選擇的操作步驟；3) 千島湖島嶼上鳥類和墓葬分布的機理。

材料與方法

研究地點與島嶼參數

按照面積和隔離度，利用分層隨機抽樣法(stratified random sampling)在千島湖選取40個島嶼。自2002年開始實地考察並詳細並測量了跟鳥類多樣性相關的各種島嶼參數：面積、隔離度、植被物種數、生境種類、周長、周長面積比、形狀指數、海拔，並於昨晚想像了各種與盜墓可能相關的島嶼參數：凹凸度、坡度、朝向、鋁和硅的含量，沙土指數和pH值。

其中鋁和硅的含量是白膏泥的主要組成元素。由於白膏泥防水性能好，是墓葬出沒的指標。沙土指數反映了建墓的可能性，即如果沙土含量過多，土質不夯實，容易測漏。pH值，跟墓葬中的有機體「發酵」程度相關。形狀指數、凹凸度、坡度和朝向是判斷風水優劣的關鍵，因為圓山、朝南、土層厚及石頭少的生境是墓葬出現的高發區。

AIC

AIC(Akaika Information Criterion)即赤池信息量準則，是評估統計模型的複雜度和衡量統計模型擬合優良性的一種標準。最早由日本統計學家赤池弘次創立和發展，由此得名。

AIC在一般情況下，可以表示為

其中: k是參數的數量, L是似然函數(likelihood function)。這是公式，知道就可以，R語言中有現成的命令(包中的命令，及包中的命令)。如果自己動手算，也可以：假設條件是模型的誤差服從獨立正態分布，n為觀察數, RSS為殘差平方和，則

增加了自由參數提高了擬合的優良性，即AIC鼓勵數據的優良性但是盡量避免出現過度擬合(overfitting)的情況，所以優先考慮的模型是AIC值最小的那一隻。

其中在小樣本的情況下(n/k < 40)，AIC 轉變成AICc (corrected AIC)，即：

當n增加時，AICc收斂成AIC。所以AICc可以應用於任何樣本大小的情況下(注: 這部分內容主要抄自維基百科，不過維基百科的該頁中文文獻引用有個小錯誤，即參考書是 Burham & Anderson(2002)，而不是2004)

如果數據有過度離散(overdispersion)的影響，則需要考慮Q版的AIC，即

$hat$ 為方差膨脹係數(VIF)或者過度離散係數(overdispersion coefficient)。如果 $hat$ 大於1，則需要採用QAIC。當然，Q版的，也有QAICc，道理同上。一般在參數進入模型前，只要保證參數的獨立性，則可以避免過度離散的情況。

計算模型權重

得到各個模型的AIC值後，按照AIC從小到大排列，然後每個模型的AIC值與最小的AIC值相減，得到ΔAIC。

通過得到的ΔAIC，計算各個模型的模型權重，即Akaika weight(wi)。其中第i個模型的模型權重為：

公式不複雜，而且R中有現成的命令計算wi。wi在0至1之間，並且所有模型權重之和為1。模型權重越大，表示該模型是真實模型的可能性就越大。比如第二個模型的w2為0.31，則表示這個模型為真實模型(best possible model)的可能性為31%。

通過模型權重還可以計算各個參數的重要值(importance)。方法很簡單，比如參數1，則挑出含參數1的所有模型，然後把這些模型的權重相加，即是該參數的權重。各個參數的權重值一比，就知道哪個參數最重要了。

模型選擇的不確定性和多模型推斷

其實現實一般不會這麼完美的，上述所有結論都建立在ΔAIC>2的基礎上，即第二個模型的AIC值比最小模型的AIC值差值大於2。如果小於2，則說明第一個模型跟第二個模型(或者連續前四五個模型)為真實模型的可能性差不多，無法決定優劣。咋么辦？終極武器：模型平均(model averaging)。

曾經ΔAIC>2是條金科玉律(Burnham & Anderson, 2002)，但是Anderson大神在2008版的書中似乎把ΔAIC>2給降級了(Andersion, 2008)，建議不要輕信這條規律，而是建議把所有模型統統進行模型平均，也就是不要隨便剔除一些看似不可能模型，哪怕這些模型的權重都小得接近於零。如果ΔAIC>2，通過最優模型，代入實際島嶼參數測量值，就可以計算出預測的鳥類種數或者存在墓葬的可能性。現在由於ΔAIC

啥意思？假設有九個可能模型，則有九個模型的權重，以及可以計算出九個預測值。如今，平均預測值就是預測值分別乘以權重後的和，比如

既然預測值Y^需要模型平均，參數估計值也得平均，道理跟估計預測值相似。假設參數i的參數估計為θi，本來當ΔAIC>2時只要直接採用最小AIC模型的θi值即可，現在則需要把含有參數i的所有模型列出來，進行模型平均：

同理，計算參數估計的方差時，也得進行模型平均，得到非條件方差估計(unconditional variance estimate)，詳見(Burnham & Anderson, 2002, p.162):

Anderson大神似乎對這個公式也不是很滿意，建議更新為Anderson (2008)第111頁的公式，其實計算結果相差不多：

其中 $hat{ar}$ 是模型的平均參數估計，wi是模型權重，以及gi表示第i個模型。簡言之，非條件方差估計就是包括兩部分：根號內的前部分是本身的取樣方差，另外一部分是由於模型選擇不確定導致的方差。所以，把後者考慮進去以後，最後的方差估計不會由於模型的不確定性而降低準確性。我怕表達有所不準，列出Anderson(2008)第111頁的原文: an estimator of the variance of parameter estimater esimates that incorporates both sampling variance, given a model, and a variance component for model selection uncertainty. 所以，在樣本量較大的前提下，最後參數的置信區間為

實戰演練

演練開始之前，請確保已經安裝下列軟體包：,,。網速給力的情況下，最簡單的方法是直接在R語言操作界面中輸入

否則，得從R的鏡像網站下載壓縮包後再本地安裝。

演練一：千島湖鳥類多樣性的決定因素

導入包

導入千島湖鳥類和島嶼數據(註：這個數據是真實的，只是我把數據的順序隨機調換了)

數據中第一列為鳥類物種數，其餘八列為島嶼參數，分別為：面積、隔離度、植物物種數、生境類別數、島嶼周長、周長面積比(越大表示邊緣越多)、形狀指數(完全的圓形，則形狀指數為1)和海拔。

模型開始之前得進行島嶼參數的獨立性檢驗。其中方法可以使用相關分析(correlation test)，方差膨脹係數(VIF)和主成份分析(PCA)，這裡採用常用的相關分析。

相關分析的R語言命令是，這是兩兩檢驗。是多個參數一起檢驗，可以多個參數一起檢驗的時候，結果不給出p值，於是我寫了一個小函數，就是多個參數檢驗的時候也同時給出p值。命令名稱為，代碼為：

所有島嶼參數進行相關分析，

結果表明，面積跟周長、周長面積比、形狀指數和海拔呈顯著相關。考慮到這些因素的生物學意義，很明顯，除去其他顯著相關的參數而保留面積是合理的，因為在島嶼生物地理學框架下，面積是極為重要的參數，且這裡的其他參數都可能由於面積而產生。比如海拔，由於是島嶼，在坡度相似的情況下，面積越大，海拔越高。所以，最後進入模型的是四個參數：面積、隔離度、植物數和生境數。

權且採用最常見的線性模型(linear model)，創建總模型(global model)，即包括所有參數：

然後利用包中的函數對所有可能模型中來選擇最優模型。此處由於是4個參數，則共有2^4=16個可能模型(此處不考慮交互效應)。

結果出來了，最優模型只包括面積和生境的參數，看看：

但再看看剛才的模型的AICc結果：

發現第二個模型的ΔAICc為223.8-223.7=0.1。坑爹啊！如果此時ΔAICc>2，則模型選擇到此結束，即最優模型為第一個模型。可是，現實比較殘忍，繼續模型平均，列出所有可能模型：

看著比較壯觀，但是碰到十個參數，共 2^10=1024 個可能模型的時候就比較麻煩了。沒事，可以再編個程序循環一下就行，後文會再次提及此問題。

16個可能模型一起平均，

結果中的第一部分，』Component models』，即列出了所有模型的自由度(df)，對數似然函數(logLik)，AICc值，ΔAICc值和模型權重。比如最優模型的模型權重為0.29，即為真實模型的可能性為29%(其實是非常低的，一般達到0.6-0.7就很不錯了，當然，這裡使用的數據是被我隨機化過的，所以結果沒有實際參考價值)

其中的第4部分，』Full model-averaged coefficients』，即是平均參數估計， $hat{ar}$ 。

第5部分，』Relative variable importance』，即是各個參數的重要值。最大為1，可見該例子中，面積是最重要的，次之是生境。至於隔離度和植物數量，則在模型中貢獻不大。

包裡面其實有現成的命令得到上述的部分結果，不信，試試輸入以下一句命令：

此時如果打算計算各島的預鳥類物種數，則可以如下進行模型平均：

pred.mat

sep=""),paste("lm",1:16,sep="")))#建立一個空矩陣，放40個島的各16各模型預測值，如下所示

還有一點是非條件方差估計，這個，有點麻煩，等以後再說。計算方法其實跟上述的 $hat{ar}$ 類似。

實戰演練二: 千島湖墓群的決定因素

這個分析就跟上述方法相似了，按部就班：

結果發現面積、形狀指數和海拔顯著相關。考慮島實際因素，島嶼面積或者說千島湖以前的山頭大小估計不會是墓葬考慮的因素，而這個山頭圓不圓，這關乎風水的事，應該是主要因素，所以剔除面積和海拔。再看發現形狀指數跟沙土也有正相關，可是考慮沙土多少是決定建不建墓的關鍵因素，予以保留，何況不是非常強烈的正相關(coef. = 0.373)。再看發現鋁、硅和坡度有相關，可以確信鋁和硅，其中之一是冗餘的，因為白膏泥富含鋁和硅。白膏泥相對鋁含量較多，此處選擇去除硅，以及另外的坡度。pH跟沙土相關，看來得把pH去除，估計過了上千年，墓葬中的有機質早化成泥土了。

再看看選取參數後的結果，

後續步驟跟演練一類似，不同的是，此處的應變數為二元結構，即presence-absence數據，得用廣義線性模型中的邏輯斯帝回歸(logistic regression)。其他註解省略，直接上程序，

結果一看，最優模型只包括形狀指數，看來理論想像的數據也不錯嘛，雖然煩人的ΔAICc依舊小於2，因此還得繼續模型平均了。因為 2^7=128 個可能模型，手動輸入運算則是比較折騰了，所以得寫個循環程序讓電腦來運算。

其中128個模型平均後的部分結果為：

再次放個大招吧，一次性列出模型平均後的部分結果：

結果

千島湖鳥類多樣性主要取決於島嶼面積和生境多樣性，而墓葬可能性取決於島嶼的形狀指數。

討論

聽說統計上有一個更牛的利器是隨機森林模型(random forest model)，可以無視參數是否獨立，直接進入模型而且可以精確預測。哪天有興趣琢磨琢磨。

PS: 以下是娛樂時間。

圓山頭是墓葬的首選，所以，各位看官以後到千島湖旅遊，不要去什麼猴島蛇島，選擇山頭比較圓的島，才是王道！

最後檢驗一下鳥類多樣性跟墓葬出現的相關性分析：

結果是顯著正相關(t = 3.2562, df = 38, p-value = 0.002378)。墓葬的出現，表示該島風水還不錯，所以最後證實本文的最初假設，即跟蜘蛛俠討論時所做的預測：鳥類多樣性與風水有顯著的相關性。至於機理等科學問題的討論，不是本篇論文能夠解決的。請聽下回分解。

致謝

謝謝看官的一路捧場，瀏覽完這塊又長又臭的博文。謝謝實驗室提供的平台和提供的支助，給於了我想像的空間，以及島嶼的數據。有關墓葬的生境數據，來自古田山大樣地，我想像著搬到千島湖了，在此一併致謝。分析方法部分參考於此。本文的源代碼及數據可以點擊此處(已更新)或此處(已過期)下載。看官就是reviewer(評審員)，若有任何reviews，請盡請留言，謝謝！

參考文獻

Anderson, David R. (2008)Model based inference in the life sciences: a primer on evidence. New York: Springer.

Burnham, Kenneth P. and David R. Anderson. (2002)Model selection and multimodel inference: a practical information-theoretic approach. Springer.

Symonds, Matthew RE and Adnan Moussalli. (2011) A brief guide to model selection, multimodel inference and model averaging in behavioural ecology using Akaike』s information criterion.Behavioral Ecology and Sociobiology,65: 13-21. APA

Whittingham, Mark J.et al.(2006) Why do we still use stepwise modelling in ecology and behaviour?.Journal of Animal Ecology,75: 1182-1189.

本文引用格式：Si X., Pimm S.L., Russell G.J. & P. Ding. (2014)Turnover of breeding bird communities on islands in an inundated lake.Journal of Biogeography, 41, 2283–2292.

來源：演算法與數學之美

史上最系統的量化投資學習！

量化投資背景闡述

常見量化投資模型構建

策略特點與交易行為分析

atlab工具箱的使用及Bar內交易和日內交易策略介紹

CTA之趨勢策略

高頻交易、支持向量機、小波分析……

1.第一部分：（7月8日14:00-18:00）

2.第二部分：（7月9日9：:00-18:00）

3.第三部分：（7月15日14：:00-7月16日18:00）

上課地點：深圳福田