模型選擇與多模型推斷

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

斯幸峰

三墩職業技術學校空想資本主義學院理論想像學研究所, 杭州 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語言中有現成的命令(stat包中的AIC命令，及stats包中的extractAIC命令)。如果自己動手算，也可以：假設條件是模型的誤差服從獨立正態分布，n為觀察數, RSS為殘差平方和，則

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

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

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

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

$\hat{c}$ 為方差膨脹係數(VIF)或者過度離散係數(overdispersion coefficient)。如果 $\hat{c}$ 大於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<2，第一個模型無法「代表」其他模型，於是所有模型都得參與進來。假設 Y^ 值為預測值(鳥類種數或墓葬出現概率)，則平均預測值為：

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

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

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

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

其中 $\hat{\bar{θ}}$ 是模型的平均參數估計，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. 所以，在樣本量較大的前提下，最後參數的置信區間為

實戰演練

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

install.packages("glmulti")

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

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

導入 glmulti包

library(glmulti)

###Loading required package: rJava

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

tilbird <- read.table("http://sixf.org/files/code/2014/tilbird.txt", h = T)  #找到 'tilbird.txt'文件並打開str(tilbird)  ##檢查`til.bird`的數據結構

###'data.frame': 40 obs. of  9 variables:### $ birdspp  : int  43 34 35 32 31 27 30 33 24 24 ...### $ area     : num  1289.2 143.2 109 55.1 46.4 ...### $ isolation: num  897 1415 965 954 730 ...### $ plants   : int  36 50 88 86 65 68 45 49 45 31 ...### $ habitats : int  3 6 3 3 3 3 3 7 4 4 ...### $ Pe       : num  105965 17465 12022 7570 10444 ...### $ PAR      : num  82.2 122 110.3 137.4 225.2 ...### $ SI       : num  832 412 325 288 433 ...### $ elev     : num  298 251 227 198 174 ...

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

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

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

cor.sig = function(test) {    res.cor = cor(test)    res.sig = res.cor    res.sig[abs(res.sig) > 0] = NA    nx = dim(test)[2]    for (i in 1:nx) {        for (j in 1:nx) {            res.cor1 = as.numeric(cor.test(test[, i], test[, j])$est)            res.sig1 = as.numeric(cor.test(test[, i], test[, j])$p.value)            if (res.sig1 <= 0.001) {                sig.mark = "***"            }            if (res.sig1 <= 0.01 & res.sig1 > 0.001) {                sig.mark = "** "            }            if (res.sig1 <= 0.05 & res.sig1 > 0.01) {                sig.mark = "*  "            }            if (res.sig1 > 0.05) {                sig.mark = "   "            }            if (res.cor1 > 0) {                res.sig[i, j] = paste(" ", as.character(round(res.cor1, 3)),                  sig.mark, sep = "")            } else {                res.sig[i, j] = paste(as.character(round(res.cor1, 3)), sig.mark,                  sep = "")            }        }    }    as.data.frame(res.sig)}

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

cor.sig(tilbird[, 2:9])  #第一列不算，那是鳥類物種數，即Y值。

###               area isolation    plants  habitats        Pe       PAR###area           1*** -0.115    -0.139    -0.064     0.996*** -0.429** ###isolation -0.115         1*** -0.101      -0.1    -0.117     0.299   ###plants    -0.139    -0.101         1***  -0.16    -0.138    -0.048   ###habitats  -0.064      -0.1     -0.16         1*** -0.057    -0.035   ###Pe         0.996*** -0.117    -0.138    -0.057         1*** -0.481** ###PAR       -0.429**   0.299    -0.048    -0.035    -0.481**       1***###SI         0.857*** -0.045    -0.167    -0.034     0.898*** -0.619***###elev       0.726*** -0.127    -0.039    -0.032     0.775*** -0.803***###                 SI      elev###area       0.857***  0.726***###isolation -0.045    -0.127   ###plants    -0.167    -0.039   ###habitats  -0.034    -0.032   ###Pe         0.898***  0.775***###PAR       -0.619*** -0.803***###SI             1***  0.888***###elev       0.888***      1***

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

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

global.model <- lm(birdspp ~ area + isolation + plants + habitats, data = tilbird)

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

bird.model <- glmulti(global.model, level = 1, crit = "aicc")  #選用AICc進行評判模型

###Initialization...###TASK: Exhaustive screening of candidate set.###Fitting...###Completed.

summary(bird.model)

###$name###[1] "glmulti.analysis"######$method###[1] "h"######$fitting###[1] "lm"######$crit###[1] "aicc"######$level###[1] 1######$marginality###[1] FALSE######$confsetsize###[1] 100######$bestic###[1] 223.7######$icvalues### [1] 223.7 223.8 225.0 225.7 226.2 226.7 228.6 228.7 243.6 244.1 244.5###[12] 244.5 246.0 246.7 246.8 247.0######$bestmodel###[1] "birdspp ~ 1 + area + habitats"######$modelweights### [1] 2.871e-01 2.708e-01 1.455e-01 1.049e-01 8.123e-02 6.348e-02 2.432e-02### [8] 2.259e-02 1.332e-05 1.036e-05 8.538e-06 8.461e-06 3.984e-06 2.794e-06###[15] 2.652e-06 2.496e-06######$includeobjects###[1] TRUE

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

lm9 <- lm(birdspp ~ area + habitats, data = tilbird)summary(lm9)

######Call:###lm(formula = birdspp ~ area + habitats, data = tilbird)######Residuals:###   Min     1Q Median     3Q    Max ###-6.606 -2.107 -0.263  1.911  8.705 ######Coefficients:###            Estimate Std. Error t value Pr(>|t|)    ###(Intercept) 20.69295    2.08432    9.93  5.6e-12 ***###area         0.01564    0.00289    5.41  3.9e-06 ***###habitats     1.29893    0.55652    2.33    0.025 *  ###---###Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1######Residual standard error: 3.67 on 37 degrees of freedom###Multiple R-squared:  0.474, Adjusted R-squared:  0.445 ###F-statistic: 16.6 on 2 and 37 DF,  p-value: 7e-06

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

summary(bird.model)$icvalue

### [1] 223.7 223.8 225.0 225.7 226.2 226.7 228.6 228.7 243.6 244.1 244.5###[12] 244.5 246.0 246.7 246.8 247.0

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

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

16個可能模型一起平均，

library(MuMIn)lm.ave <- model.avg(lm1, lm2, lm3, lm4, lm5, lm6, lm7, lm8, lm9, lm10, lm11,    lm12, lm13, lm14, lm15, lm16)summary(lm.ave)

######Call:###model.avg.default(object = lm1, lm2, lm3, lm4, lm5, lm6, lm7, ###    lm8, lm9, lm10, lm11, lm12, lm13, lm14, lm15, lm16)######Component models:###       df logLik  AICc Delta Weight###12      4 -107.3 223.7  0.00   0.29###123     5 -106.0 223.8  0.12   0.27###124     5 -106.6 225.0  1.36   0.15###1234    6 -105.6 225.7  2.01   0.10###13      4 -108.5 226.2  2.53   0.08###1       3 -110.0 226.7  3.02   0.06###134     5 -108.4 228.6  4.94   0.02###14      4 -109.8 228.7  5.08   0.02###3       3 -118.5 243.6 19.96   0.00###23      4 -117.5 244.1 20.46   0.00###(Null)  2 -120.1 244.5 20.85   0.00###2       3 -118.9 244.5 20.86   0.00###34      4 -118.4 246.0 22.37   0.00###234     5 -117.5 246.7 23.08   0.00###4       3 -120.1 246.8 23.18   0.00###24      4 -118.9 247.0 23.31   0.00######Term codes:###     area  habitats isolation    plants ###        1         2         3         4 ######Model-averaged coefficients: ###             Estimate Std. Error Adjusted SE z value Pr(>|z|)    ###(Intercept) 22.023011   3.272613    3.327045    6.62   <2e-16 ***###area         0.015423   0.002945    0.003044    5.07   <2e-16 ***###habitats     1.287322   0.560392    0.579478    2.22    0.026 *  ###isolation   -0.001104   0.000728    0.000753    1.47    0.143    ###plants       0.019405   0.021519    0.022247    0.87    0.383    ###---###Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1######Full model-averaged coefficients (with shrinkage): ### (Intercept)      area  habitats isolation    plants###   22.023011  0.015422  1.040595 -0.000531  0.005770######Relative variable importance:###     area  habitats isolation    plants ###     1.00      0.81      0.48      0.30

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

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

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

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

dredge(global.model)

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

pred.mat <- matrix(NA, ncol = 16, nrow = 40, dimnames = list(paste("isl", 1:40,

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

bird.pred <- pred.mat %*% summary(lm.ave)$msTable$weight[c(4,14,3,2,7,16,10,13,1,8,5,6,9,15,12,11)]round(t(bird.pred),2) #把矩陣換方向，給頁面省點空間，跟分析無關

###      isl1  isl2  isl3  isl4  isl5  isl6  isl7  isl8  isl9 isl10###[1,] 44.76 30.01 26.82 25.98 25.85 24.97 24.87 29.67 26.19 25.46###     isl11 isl12 isl13 isl14 isl15 isl16 isl17 isl18 isl19 isl20###[1,] 25.11 25.06 24.81 23.79 25.03 24.41    26 25.03 24.92  26.6###     isl21 isl22 isl23 isl24 isl25 isl26 isl27 isl28 isl29 isl30###[1,] 25.63 25.12 24.32  24.1 25.02 24.62 24.15 26.52 27.03 23.69###     isl31 isl32 isl33 isl34 isl35 isl36 isl37 isl38 isl39 isl40###[1,] 25.59  24.5 24.27 26.11 25.51  26.5 24.88 24.95 28.28 24.91

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

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

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

tiltomb <- read.table("http://sixf.org/files/code/2014/tiltomb.txt", h = T)  #讀取墓的虛擬數據 'tiltomb.txt'cor.sig(tiltomb[, -1])

###              area    plants  habitats        SI      elev    convex###area          1*** -0.139    -0.064     0.857***  0.726***  0.041   ###plants   -0.139         1***  -0.16    -0.167    -0.039    -0.107   ###habitats -0.064     -0.16         1*** -0.034    -0.032    -0.187   ###SI        0.857*** -0.167    -0.034         1***  0.888***  0.237   ###elev      0.726*** -0.039    -0.032     0.888***      1***  0.307   ###convex    0.041    -0.107    -0.187     0.237     0.307         1***###slope     0.248     0.193     0.115     0.247     0.322*    0.264   ###aspect   -0.114    -0.081     0.141     0.069     0.278     0.075   ###Al        0.088    -0.066     0.193     0.223     0.326*     0.06   ###Si        0.055    -0.099     0.101      0.14     0.243     0.081   ###sand     -0.207     0.197     0.111     -0.22    -0.191    -0.311   ###pH       -0.194    -0.132    -0.204     -0.17    -0.228     0.018   ###             slope    aspect        Al        Si      sand        pH###area      0.248    -0.114     0.088     0.055    -0.207    -0.194   ###plants    0.193    -0.081    -0.066    -0.099     0.197    -0.132   ###habitats  0.115     0.141     0.193     0.101     0.111    -0.204   ###SI        0.247     0.069     0.223      0.14     -0.22     -0.17   ###elev      0.322*    0.278     0.326*    0.243    -0.191    -0.228   ###convex    0.264     0.075      0.06     0.081    -0.311     0.018   ###slope         1*** -0.093     0.412**   0.334*    0.111     -0.53***###aspect   -0.093         1***  0.075     0.101    -0.086     0.198   ###Al        0.412**   0.075         1***  0.887***  0.615*** -0.756***###Si        0.334*    0.101     0.887***      1***  0.504*** -0.646***###sand      0.111    -0.086     0.615***  0.504***      1*** -0.598***###pH        -0.53***  0.198    -0.756*** -0.646*** -0.598***      1***

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

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

cor.sig(tiltomb[, c("plants", "habitats", "SI", "convex", "aspect", "Al", "sand")])

###            plants  habitats        SI    convex    aspect        Al###plants        1***  -0.16    -0.167    -0.107    -0.081    -0.066   ###habitats  -0.16         1*** -0.034    -0.187     0.141     0.193   ###SI       -0.167    -0.034         1***  0.237     0.069     0.223   ###convex   -0.107    -0.187     0.237         1***  0.075      0.06   ###aspect   -0.081     0.141     0.069     0.075         1***  0.075   ###Al       -0.066     0.193     0.223      0.06     0.075         1***###sand      0.197     0.111     -0.22    -0.311    -0.086     0.615***###              sand###plants    0.197   ###habitats  0.111   ###SI        -0.22   ###convex   -0.311   ###aspect   -0.086   ###Al        0.615***###sand          1***

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

global.model.tomb <- glm(tomb ~ plants + habitats + SI + convex + aspect + Al +    sand, family = binomial("logit"), data = tiltomb)tomb.model <- glmulti(global.model.tomb, level = 1, crit = "aicc")

###Initialization...###TASK: Exhaustive screening of candidate set.###Fitting...######After 50 models:###Best model: tomb~1+SI###Crit= 57.9820910321992###Mean crit= 64.0858355584437

######After 100 models:###Best model: tomb~1+SI###Crit= 57.9820910321992###Mean crit= 64.9421343165768

######After 150 models:###Best model: tomb~1+SI###Crit= 57.9820910321992###Mean crit= 64.5619346833708

###Completed.

summary(tomb.model)

###$name###[1] "glmulti.analysis"######$method###[1] "h"######$fitting###[1] "glm"######$crit###[1] "aicc"######$level###[1] 1######$marginality###[1] FALSE######$confsetsize###[1] 100######$bestic###[1] 57.98######$icvalues###  [1] 57.98 58.33 59.61 60.17 60.22 60.30 60.39 60.46 60.54 60.75 60.82### [12] 60.94 62.15 62.17 62.18 62.24 62.39 62.44 62.67 62.70 62.77 62.77### [23] 62.79 62.87 62.89 62.91 62.95 62.98 63.01 63.09 63.19 63.27 63.32### [34] 63.51 63.53 63.60 63.66 64.05 64.43 64.52 64.53 64.60 64.62 64.77### [45] 64.80 64.87 64.88 64.91 64.92 64.92 64.95 64.96 65.08 65.12 65.14### [56] 65.40 65.40 65.45 65.52 65.54 65.58 65.65 65.66 65.67 65.69 65.69### [67] 65.72 65.81 65.82 65.88 66.02 66.08 66.14 66.14 66.22 66.32 66.47### [78] 66.56 66.64 66.72 66.82 66.85 66.91 66.92 66.93 67.22 67.35 67.37### [89] 67.38 67.64 67.65 67.66 67.67 67.80 67.83 67.83 67.86 67.87 68.02###[100] 68.17######$bestmodel###[1] "tomb ~ 1 + SI"######$modelweights###  [1] 0.1201201 0.1011728 0.0531133 0.0401592 0.0392729 0.0377621 0.0359927###  [8] 0.0348480 0.0333513 0.0300502 0.0290503 0.0273379 0.0149678 0.0148350### [15] 0.0147370 0.0143245 0.0132729 0.0129413 0.0115051 0.0113448 0.0109770### [22] 0.0109687 0.0108420 0.0104403 0.0103423 0.0102346 0.0100197 0.0098576### [29] 0.0097113 0.0093423 0.0088919 0.0085364 0.0083084 0.0075580 0.0074968### [36] 0.0072460 0.0070090 0.0057865 0.0047700 0.0045613 0.0045489 0.0043911### [43] 0.0043397 0.0040282 0.0039725 0.0038450 0.0038208 0.0037599 0.0037378### [50] 0.0037365 0.0036946 0.0036696 0.0034468 0.0033923 0.0033552 0.0029419### [57] 0.0029405 0.0028648 0.0027779 0.0027422 0.0026895 0.0025989 0.0025888### [64] 0.0025730 0.0025506 0.0025475 0.0025079 0.0023998 0.0023883 0.0023173### [71] 0.0021582 0.0020990 0.0020310 0.0020308 0.0019495 0.0018537 0.0017263### [78] 0.0016468 0.0015814 0.0015211 0.0014462 0.0014228 0.0013853 0.0013789### [85] 0.0013704 0.0011849 0.0011078 0.0011006 0.0010928 0.0009612 0.0009542### [92] 0.0009487 0.0009451 0.0008861 0.0008726 0.0008723 0.0008583 0.0008546### [99] 0.0007955 0.0007370######$includeobjects###[1] TRUE

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

tomb7=tiltomb[, c("plants", "habitats", "SI", "convex", "aspect", "Al", "sand","tomb")]npar=7modPar=c("plants", "habitats", "SI", "convex", "aspect", "Al", "sand","tomb")unit=c(1,0)parEst=rep(unit,each=2^(npar-1))for (i in 2:npar){  unit=c(i,0)  parEst.tmp=rep(rep(unit,each=2^(npar-i)),2^(i-1))  parEst=cbind(parEst,parEst.tmp)}parMat=cbind(parEst[,1:npar],1)dimnames(parMat)=list(1:(2^npar),modPar)allModel=list()for (i in 1:(dim(parMat)[1]-1)) {    tomb7.tmp=tomb7[,parMat[i,]!=0]    allModel[[i]]=glm(tomb~.,family = binomial("logit"),data=tomb7.tmp)}modelC=glm(tomb~1,family = binomial("logit"),data=tomb7)lm.ave <- model.avg(allModel,modelC)summary(lm.ave)

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

### Full model-averaged coefficients (with shrinkage): ### (Intercept)          SI          Al        sand      aspect    habitats      convex      plants### -2.96791424  0.01986124 -0.09423029 -0.08413135  0.00217080 -0.04294219  0.00446755  0.00043138######Relative variable importance:###      SI   aspect       Al     sand habitats   convex   plants ###   0.98     0.27     0.27     0.26     0.24     0.23     0.22

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

dredge(global.model.tomb)

結果

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

討論

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

PS: 以下是娛樂時間。

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

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

cor.test(tilbird[, 1], tiltomb[, 1])

###### Pearson's product-moment correlation######data:  tilbird[, 1] and tiltomb[, 1]###t = 3.256, df = 38, p-value = 0.002378###alternative hypothesis: true correlation is not equal to 0###95 percent confidence interval:### 0.1821 0.6797###sample estimates:###   cor ###0.4671

結果是顯著正相關(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.

往期精彩文章：

第一個被認為「科學家」的人:泰勒斯

數學思維比數學運算更重要

二十世紀的十大科學騙局

瞎扯現代數學的基礎

x背後的軼聞趣事

主宰這個世界的10大算法

16個讓你燒腦讓你暈的悖論

機器學習中距離和相似性度量方法

傳說中的快排是怎樣的

玻璃秘史：一個人 改變了全世界

程序人生的四個象限和兩條主線

比特幣的原理及運作機制

概率論公式，你值得擁有

分類算法之樸素貝葉斯算法

採樣定理：有限個點構建出整個函數

長按二維碼識別

歡迎大家把獨特見解分享出來，投稿郵箱：

[e-mail]