Skip to content

Instantly share code, notes, and snippets.

@sixf

sixf/TIL-model-selection

Last active April 5, 2024 16:12
Show Gist options
  • Save sixf/9504429 to your computer and use it in GitHub Desktop.
Save sixf/9504429 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
TIL-model-selection
---
title: 千岛湖鸟类与墓群的决定因素分析
tags: [多模型推断,模型选择,千岛湖,逐步回归,盗墓]
categories: [统计,娱乐]
layout: post
comments: yes
---
<script src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML" type="text/javascript"></script>
#### －基于模型选择和多模型推断的方法
[**斯幸峰**](http://sixf.org/cn)
三墩职业技术学校空想资本主义学院理论空间学研究所, 杭州 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米水位时)，乃名副其实的“千岛湖”，是一个得天独厚的路桥岛屿天然实验场所。[本研究团队](http://mypage.zju.edu.cn/personnelCard/pingding)自2002年开始千岛湖地区的鸟类调查，到目前已经逐渐拓展到蜘蛛、蜥蜴、青蛙、蛇、猴子、昆虫、兽类、蝴蝶以及植物等各项业务，欢迎广大生态爱好者和有志之士前来参观与洽谈。撰写本文的起因是早先跟本团队中的“蜘蛛侠”吴博士尝试探讨鸟类多样性与风水的关系，加上近日刚好看了一些有关模型选择和多模型推断(model selection and multimodel inference)的文献，采用“先进”的AIC(Akaike information criterion)技术，探讨该学术问题的可能性。
本文主要探讨的问题包括两部分：1) AIC是啥？莫非是美国国际大学(American International College)吗？ 2) 模型选择的操作步骤；3) 千岛湖岛屿上鸟类和墓葬分布的机理。
## 材料与方法
### 研究地点与岛屿参数
按照面积和隔离度，进行分层随机(stratified random sampling)在千岛湖选取40个岛屿。自2002年起开始实地考察并详细并测量了跟鸟类多样性相关的各种岛屿参数：面积、隔离度、 植被物种数、生境数目、周长、周长面积比、形状指数、海拔，并于昨晚想像了各种与及盗墓可能性相关的岛屿参数：凹凸度、坡度、朝向、铝和硅的含量，沙土指数和pH值。
其中铝和硅的含量是白膏泥的主要组成元素。由于白膏泥防水性能好，是墓葬出没的指标。沙土指数反映了建墓的可能性，即如果沙土含量过多，土质不夯实，容易测漏。pH值，跟墓葬中的有机体“发酵”程度相关。形状指数、凹凸度、坡度和朝向是判断风水优劣的关键，因为圆山、朝南、土层厚及石头少的生境是墓葬出现的高发区。
### AIC
AIC(Akaika Information Criterion)即赤池信息量准则，是评估统计模型的复杂度和衡量统计模型拟合优良性的一种标准。最早由日本统计学家赤池弘次创立和发展，由此得名。
AIC在一般情况下，可以表示为:
$$AIC = 2k-2ln(L)$$
其中: k是参数的数量, L是似然函数(likelihood function)。这是公式，知道就可以，R语言中有现成的命令(`stat`包中的`AIC`命令，及`stats`包中的`extractAIC`命令)。如果自己动手算，也可以：假设条件是模型的误差服从独立正态分布，n为观察数, RSS为残差平方和，则
$$AIC = 2k+n ln(\frac{RSS}{n})$$
增加了自由参数提高了拟合的优良性，即AIC鼓励数据的优良性但是尽量避免出现过度拟合(overfitting)的情况，所以优先考虑的模型是AIC值最小的那一只。
其中在样本小的情况下 $\left( \frac{n}{k} <40 \right)$ ，AIC 转变成AICc(corrected AIC)，即：
$$AICc = AIC + 2k \frac {k+1}{n-k-1}$$
当n增加时，AICc收敛成AIC。所以AICc可以应用于任何样本大小的情况下(注: 这部分内容主要抄自[维基百科](http://zh.wikipedia.org/wiki/赤池信息量准则)，不过维基百科中的文献引用有个小错误，即参考书是 Burham & Anderson(2002)，而不是2004)
如果数据有过度离散(overdispersion)的影响，则需要考虑Q版的AIC，即
$$QAIC = \frac {-2ln(L)}{\hat{c}} + 2k$$
$\hat{c}$ 为方差膨胀系数(VIF)或者过度离散系数(overdispersion coefficient)。如果 $\hat{c}$ 大于1，则需要采用QAIC。当然，Q版的，也有QAICc，道理同上。一般可以参数进入模型前，只要保证参数的独立性，则可以避免过度离散的情况。
### 计算模型权重
得到各个模型的AIC值后，按照AIC从小到大排列，然后每个模型的AIC值与最小的AIC值相减，得到ΔAIC。
通过得到的ΔAIC，计算各个模型的模型权重，即Akaika weight($w_i$)。其中第 $i$ 个模型的模型权重为：
$$w_i = \frac {e^{(- \frac{1}{2} \Delta AIC_i)}}{\sum_{r=1}^R e^{(- \frac {1}{2} \Delta AIC_r)}}$$
公式不复杂，而且R中有现成的命令计算 $w_i$ 。 $w_i$ 在0至1之间，并且所有模型权重之和为1。模型权重越大，表示该模型是真实模型的可能性就越大。比如第二个模型的 $w_2$ 为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，第一个模型无法“代表”其他模型，于是所有模型都得参与进来。假设 $\hat {Y_i}$ 值为预测值(鸟类种数或墓葬出现概率)，则平均预测值 $\hat {\bar {Y}}$ 为：
$$\hat {\bar {Y}} = \sum_{r=1}^R (w_i \hat {Y_i})$$
啥意思？假设有十个可能模型，则有十个模型的权重，以及可以计算出十个预测值。如今，平均预测值就是预测值分别乘以权重后的和，比如
$$w_1 \hat{Y_1} + w_2 \hat{Y_2} + ... + w_{10} \hat{Y_{10}}$$
既然预测值 $\hat {Y}$ 需要模型平均，参数估计值也得平均，道理跟估计预测值相似。假设参数 $i$ 的参数估计为 $\theta_i$ ，本来当ΔAIC>2时只要直接采用最小AIC模型的 $\theta_i$ 值即可，现在则需要把含有参数 $i$ 的所有模型列出来，进行相似的模型平均：
$$\hat {\bar {\theta}} = \sum_{r=1}^R (w_i \hat {\theta_i})$$
同理，计算参数估计的方差时，也得进行模型平均，得到非条件方差估计(unconditional variance estimate)，详见(Burnham & Anderson, 2002, p162):
$$\hat {var}\left(\hat {\bar {\theta}}\right) = \left[\sum_{i=1}^R w_i \sqrt{\hat {var}(\hat {\theta_i}|g_i)+(\hat {\theta_i}-\hat {\bar {\theta}})^2} \right]^2$$
Anderson大神似乎对这个公式也不是很满意，建议更新为Anderson(2008)第111页的公式，其实结果相差不多：
$$var\left(\hat {\bar {\theta}}\right) = \sum_{i=1}^R w_i\left[var(\hat {\theta_i}|g_i)+(\hat {\theta_i}-\hat {\bar {\theta}})^2 \right]$$
其中 $\hat {\bar {\theta}}$ 是模型的平均参数估计， $w_i$ 是模型权重，以及 $g_i$ 表示第 $i$ 个模型。简言之，非条件方差估计就是包括两部分：一部分是本身的取样方差 $var(\hat {\theta_i}|g_i)$ ，另外一部分是由于模型选择不确定导致的方差 $(\hat {\theta_i}-\hat {\bar {\theta}})^2$ 。所以，把后者考虑进去以后，最后的方差估计不会由于模型的不确定性而降低准确性。我怕表达不准，列出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. 所以，最后参数的置信区间为
$$\hat {\bar {\theta}} \pm 1.96 \times \sqrt{var\left({\hat {\bar {\theta}}}\right)}$$
### 实战演练
演练开始之前，请确保已经安装下列软件包：`glmulti`, `MuMIn`, `bbmle`。网速给力的情况下，最简单的方法是直接在R语言操作界面输入
`install.packages("glmulti")`
否则，从R的镜像网站下载压缩包后再本地安装。
#### 演练一：千岛湖鸟类多样性的决定因素
导入 `glmulti`包
```r
library(glmulti)
```
```
## Loading required package: rJava
```
导入千岛湖鸟类和岛屿数据(注：这个数据是真实的，只是我把数据随机调换顺序了)
```r
tilbird <- read.table("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`，代码为：
```r
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)
}
```
所有岛屿参数进行相关分析，
```r
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)，即包括所有参数：
```r
global.model <- lm(birdspp ~ area + isolation + plants + habitats, data = tilbird)
```
然后利用`glmulti`包中的函数`glmulti`对所有可能模型中来选择最优模型。此处由于是4个参数，则共有2^4=16个可能模型(**此处不考虑交互效应**)。
```r
bird.model <- glmulti(global.model, level = 1, crit = "aicc") #选用AICc进行评判模型
```
```
## Initialization...
## TASK: Exhaustive screening of candidate set.
## Fitting...
## Completed.
```
```r
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
```
结果出来了，最优模型只包括面积和生境的参数，看看：
```r
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结果：
```r
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为225.6429-225.4588=0.1841。坑爹啊！如果此时ΔAICc>2，则模型选择到此结束，即最优模型为第一个模型。可是，现实比较残忍，继续模型平均，列出所有可能模型：
```r
lm1 <- lm(birdspp ~ area + isolation + plants + habitats, data = tilbird)
lm2 <- lm(birdspp ~ isolation + plants + habitats, data = tilbird)
lm3 <- lm(birdspp ~ area + plants + habitats, data = tilbird)
lm4 <- lm(birdspp ~ area + isolation + habitats, data = tilbird)
lm5 <- lm(birdspp ~ area + isolation + plants, data = tilbird)
lm6 <- lm(birdspp ~ plants + habitats, data = tilbird)
lm7 <- lm(birdspp ~ isolation + habitats, data = tilbird)
lm8 <- lm(birdspp ~ isolation + plants, data = tilbird)
lm9 <- lm(birdspp ~ area + habitats, data = tilbird)
lm10 <- lm(birdspp ~ area + plants, data = tilbird)
lm11 <- lm(birdspp ~ area + isolation, data = tilbird)
lm12 <- lm(birdspp ~ area, data = tilbird)
lm13 <- lm(birdspp ~ isolation, data = tilbird)
lm14 <- lm(birdspp ~ plants, data = tilbird)
lm15 <- lm(birdspp ~ habitats, data = tilbird)
lm16 <- lm(birdspp ~ 1, data = tilbird)
```
看着比较壮观，但是碰到十个参数，共 $2^10=1024$ 个可能模型的时候就比较麻烦了。没事，可以再编个程序循环一下就行，此处暂时不提。
16个可能模型一起平均，
```r
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.22**，即为真实模型的可能性为22%(其实是非常低的，一般达到0.6-0.7就很不错了，当然，这里使用的数据是被我随机化过的，所以结果没有实际参考价值)
其中的第4部分，'Full model-averaged coefficients'，即使平均参数估计， $\hat {\bar {\theta}}$ 。
第5部分，'Relative variable importance'，即是各个参数的重要值。最大为1，可见该例子中，面积是最重要的，次之是生境。至于隔离度和植物数量，则在模型中贡献不大。
此时如果打算计算岛1预测的鸟类物种数，则可以如下进行模型平均：
```r
pred.mat <- matrix(NA, ncol = 16, nrow = 40, dimnames = list(paste("isl", 1:40,
sep = ""), paste("lm", 1:16, sep = ""))) #建立一个空矩阵，放40个岛的各16各模型预测值，如下所示
pred.mat[, 1] <- predict(lm1)
pred.mat[, 2] <- predict(lm2)
pred.mat[, 3] <- predict(lm3)
pred.mat[, 4] <- predict(lm4)
pred.mat[, 5] <- predict(lm5)
pred.mat[, 6] <- predict(lm6)
pred.mat[, 7] <- predict(lm7)
pred.mat[, 8] <- predict(lm8)
pred.mat[, 9] <- predict(lm9)
pred.mat[, 10] <- predict(lm10)
pred.mat[, 11] <- predict(lm11)
pred.mat[, 12] <- predict(lm12)
pred.mat[, 13] <- predict(lm13)
pred.mat[, 14] <- predict(lm14)
pred.mat[, 15] <- predict(lm15)
pred.mat[, 16] <- predict(lm16)
# 输出40个岛屿的平均预测值，即上述的 hat-bar(Y)
bird.pred <- pred.mat %*% summary(lm.ave)$summary$Weight
t(bird.pred) #把矩阵换方向，给页面省点空间，跟分析无关
```
```
## isl1 isl2 isl3 isl4 isl5 isl6 isl7 isl8 isl9 isl10 isl11
## [1,] 37.64 29.45 26.74 26.21 26.15 24.7 24.62 30.58 26.76 25.43 25.56
## isl12 isl13 isl14 isl15 isl16 isl17 isl18 isl19 isl20 isl21 isl22
## [1,] 25.51 24.19 24.07 25.48 24.31 26.46 25.51 24.43 26.67 26.58 25.67
## isl23 isl24 isl25 isl26 isl27 isl28 isl29 isl30 isl31 isl32 isl33
## [1,] 24.18 23.78 25.49 24.75 23.87 26.52 27.5 23.02 26.37 24.45 24.08
## isl34 isl35 isl36 isl37 isl38 isl39 isl40
## [1,] 26.61 26.38 26.48 25.21 25.34 28.94 25.29
```
还有一点是非条件方差估计，这个，有点麻烦，等以后再说。计算方法其实跟上述的 $\hat {\bar {Y}}$ 类似。
#### 实战演练2: 千岛湖墓群的决定因素
这个分析就跟上述方法相似了，按部就班：
```r
tiltomb <- read.table("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去除，估计过了上千年，墓葬中的有机质早化成泥土了。
再看看选取参数后的结果，
```r
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***
```
后续步骤跟演练1类似，不同的是，此处的应变量为二元结构，即presence-absence数据，得用广义线性模型中的逻辑斯帝回归(logistic regreθ)。其他注解省略，直接上程序，
```r
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
```
![plot of chunk unnamed-chunk-15](figure/unnamed-chunk-151.png)
```
##
## After 100 models:
## Best model: tomb~1+SI
## Crit= 57.9820910321992
## Mean crit= 64.9421343165768
```
![plot of chunk unnamed-chunk-15](figure/unnamed-chunk-152.png)
```
##
## After 150 models:
## Best model: tomb~1+SI
## Crit= 57.9820910321992
## Mean crit= 64.5619346833708
```
![plot of chunk unnamed-chunk-15](figure/unnamed-chunk-153.png)
```
## Completed.
```
```r
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$ 个可能模型，那个循环程序还没写好，所以就此为止。
## 结果
千岛湖鸟类多样性主要取决于岛屿面积和生境多样性，而墓葬可能性取决于岛屿的生境多样性。
## 讨论
听说统计上有一个更牛的利器是[随机森林模型](http://blog.sciencenet.cn/blog-661364-615921.html)(random forest model)，可以无视参数是否独立，直接进入模型而且可以精确预测。哪天有兴趣琢磨琢磨。
PS: 以下是娱乐时间。
生境丰富的山头是墓葬的首选，所以，各位看官以后到千岛湖旅游，不要去什么猴岛蛇岛，选择生境种类最多的岛，才是王道！怎么算生境种类，建议看看[群落生态学](http://en.wikipedia.org/wiki/Community_(ecology))。
最后检验一下鸟类多样性跟墓葬出现的相关性分析：
```r
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)。墓葬的出现，表示该岛风水还不错，所以最后证实本文的最初假设，即跟蜘蛛侠讨论时所做的预测：鸟类多样性与风水有显著的相关性。至于机理的讨论，还得依赖于路径分析(pathway analyses)，即是不是中间由于生境多样性搭了桥的原因？对于有关该科学问题的进一步分析，这不是本篇论文能够解决的。请听下回分解。
## 致谢
谢谢看官的一路捧场，浏览完这块又长又臭的博文。谢谢实验室提供的平台和提供的支助，给于了我想像的空间，以及岛屿的数据。有关墓葬的生境数据，来自[古田山大样地](http://blog.sciencenet.cn/blog-267448-463699.html)，我想像着搬到千岛湖了，在此致谢。分析方法部分参考于<a href="https://gist.github.com/sixf/9488518">此</a>。本文的源代码及数据可以点击[此处](https://github.com/sixf/TIL-model-selection/archive/master.zip)下载。看官就是reviewer(评审员)，若有任何reviews，请在留言处指出，谢谢！
## 参考文献
1. Anderson, David R. (2008) *Model based inference in the life sciences: a primer on evidence*. New York: Springer.
1. Burnham, Kenneth P., and David R. Anderson. (2002) *Model selection and multimodel inference: a practical information-theoretic approach*. Springer.
2. Symonds, Matthew RE, and Adnan Mouθalli. (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
3. Whittingham, Mark J., et al. (2006) Why do we still use stepwise modelling in ecology and behaviour?. *Journal of animal ecology*, **75**: 1182-1189.
Raw
tilbird
"birdspp" "area" "isolation" "plants" "habitats" "Pe" "PAR" "SI" "elev"
"1" 43 1289.23 897.41 198 7 105964.7686 82.19229 832.4723 298.5
"2" 34 143.19 1415.09 99 6 17464.9831 121.9707 411.7433 251.1
"3" 35 109.03 964.97 86 6 12022.2003 110.26507 324.7969 226.9
"4" 32 55.08 953.95 59 5 7569.6441 137.42999 287.7235 198
"5" 31 46.37 729.8 51 5 10444.3768 225.23996 432.6807 174
"6" 27 35.64 2110.41 49 5 7318.3728 205.3416 345.8136 175.2
"7" 30 32.29 1936.95 57 5 5891.5018 182.45592 292.4811 196.7
"8" 33 5.69 21.85 69 3 1489.0858 261.7022 176.1085 109
"9" 24 3.42 583 74 4 1164.0401 340.36259 177.5585 131.1
"10" 24 2.9 1785.3 85 3 1241.7205 428.17949 205.6961 126.3
"11" 30 2.83 1238.14 86 4 900.8172 318.30996 151.0635 125.5
"12" 29 2.29 973.85 65 4 956.9893 417.89926 178.395 130.5
"13" 26 2.23 3261.96 53 3 1341.0224 601.35534 253.3318 113.5
"14" 23 2 1042.38 45 3 794.5238 397.26191 158.4913 126
"15" 24 1.93 888.05 50 4 658.0149 340.94035 133.62 125.3
"16" 26 1.74 2293.25 100 3 604.3852 347.34779 129.2567 126.2
"17" 28 1.54 711.04 88 3 706.2847 458.62641 160.5491 125.5
"18" 24 1.52 849.88 40 3 604.3504 397.599 138.2827 123.5
"19" 23 1.52 2849.99 53 3 715.7862 470.91194 163.7761 126.6
"20" 23 1.4 1760.34 49 3 722.5082 516.07727 172.2591 115.9
"21" 27 1.26 54.86 65 3 427.3114 339.1361 107.3828 138.5
"22" 26 1.2 657.72 56 3 505.729 421.44082 130.2297 125.1
"23" 25 1.2 2128.52 68 3 551.9353 459.94607 142.1382 116.7
"24" 22 1.17 2453.37 69 3 608.5181 520.10096 158.6975 116.6
"25" 23 1.15 847.12 33 3 503.0566 437.44052 132.3302 121
"26" 18 1.03 1458.81 36 3 488.1191 473.90207 135.6801 115.6
"27" 20 1.01 2437.85 29 3 667.929 661.31584 187.4852 116.7
"28" 27 1.01 2103.85 36 3 492.2143 487.34092 138.1583 124
"29" 29 0.97 938.85 70 3 486.7272 501.78064 139.4073 120.1
"30" 21 0.96 3133.96 50 3 595.3203 620.1253 171.4 115
"31" 23 0.91 1339.71 50 4 412.1105 452.86873 121.8633 122.2
"32" 25 0.86 2321.51 56 3 428.6627 498.445 130.3991 119.5
"33" 22 0.83 2298.5 50 3 432.4952 521.07852 133.9143 123
"34" 20 0.83 1098.58 45 4 397.6094 479.04748 123.1127 120.7
"35" 23 0.8 102.6 68 3 487.9542 609.9427 153.8995 114
"36" 29 0.8 2097.52 80 2 450.65 563.3125 142.1382 117.8
"37" 23 0.73 1320.4 31 3 526.9358 721.82989 173.973 106.7
"38" 18 0.67 1139.87 39 3 419.6842 626.39428 144.6331 109.8
"39" 25 0.59 640.53 42 3 313.8599 531.96593 115.2726 113
"40" 26 0.59 1018.42 55 3 407.6879 690.99651 149.725 110
Raw
tiltomb
"tomb" "area" "plants" "habitats" "SI" "elev" "convex" "slope" "aspect" "Al" "Si" "sand" "pH"
"1" 1 1289.23 198 7 832.4723 298.5 -1.74375 46.78939 114.01775 3.334376 0.2147654 3.694624 4.713993
"2" 1 143.19 99 6 411.7433 251.1 -2.215 45.69906 263.08339 4.202102 0.2451958 4.373228 4.77312
"3" 1 109.03 86 6 324.7969 226.9 3.84075 44.46079 256.71699 3.833674 0.226052 4.289754 4.935478
"4" 1 55.08 59 5 287.7235 198 4.326 39.24653 237.81106 3.079796 0.2017985 3.95893 5.087159
"5" 1 46.37 51 5 432.6807 174 1.779 37.27584 110.11573 2.887753 0.1963663 3.644821 5.167524
"6" 1 35.64 49 5 345.8136 175.2 4.287 35.67033 185.92226 3.678946 0.2422317 4.624101 4.815154
"7" 1 32.29 57 5 292.4811 196.7 1.848 33.06503 201.07674 4.654444 0.2853413 5.090561 4.678613
"8" 0 5.69 69 3 176.1085 109 -2.68175 37.39011 184.01256 4.710466 0.2841518 5.814682 4.57159
"9" 1 3.42 74 4 177.5585 131.1 -1.253 45.07064 150.5588 4.292108 0.23286 5.726786 4.682153
"10" 0 2.9 85 3 205.6961 126.3 -6.42925 45.64218 131.79889 3.798971 0.2010399 5.481689 4.765959
"11" 1 2.83 86 4 151.0635 125.5 -1.68575 42.62368 122.48938 3.232953 0.1929789 5.132699 4.599959
"12" 1 2.29 65 4 178.395 130.5 2.06225 39.33636 99.11854 3.050695 0.1883907 4.651636 4.825771
"13" 1 2.23 53 3 253.3318 113.5 1.78125 41.37146 179.76741 2.946751 0.2007715 4.551206 5.091734
"14" 1 2 45 3 158.4913 126 1.32875 33.53033 157.16585 2.796017 0.2025215 4.615628 5.075911
"15" 0 1.93 50 4 133.62 125.3 0.6375 23.12768 158.78088 2.787408 0.2189087 5.072189 4.969378
"16" 0 1.74 100 3 129.2567 126.2 -2.58625 40.5828 144.53658 2.777316 0.2293555 5.304275 4.946582
"17" 1 1.54 88 3 160.5491 125.5 -2.41025 29.9744 179.66997 2.490265 0.2131747 4.70214 5.118917
"18" 1 1.52 40 3 138.2827 123.5 2.4145 40.31063 203.94778 2.083769 0.1587637 3.804712 5.184095
"19" 0 1.52 53 3 163.7761 126.6 2.19825 33.43398 204.3187 1.762225 0.1288696 3.330006 5.196133
"20" 0 1.4 49 3 172.2591 115.9 0.03225 30.76756 168.02526 1.684568 0.1291059 3.678364 5.111234
"21" 0 1.26 65 3 107.3828 138.5 -1.1075 51.28462 106.23779 3.82252 0.2846606 3.936857 4.666953
"22" 0 1.2 56 3 130.2297 125.1 4.1321 48.21877 179.71606 4.180629 0.2936661 4.701744 4.696333
"23" 0 1.2 68 3 142.1382 116.7 2.10755263 42.36636 175.59017 3.401895 0.2254496 4.537955 5.03401
"24" 0 1.17 69 3 158.6975 116.6 -2.86494737 32.92537 175.87571 2.763922 0.1969102 4.040445 5.261254
"25" 1 1.15 33 3 132.3302 121 -5.73885526 22.06466 153.19076 2.575541 0.1751494 3.476163 5.404008
"26" 0 1.03 36 3 135.6801 115.6 -4.0025921 33.3992 168.35525 3.226501 0.2148751 4.377678 4.900335
"27" 0 1.01 29 3 187.4852 116.7 -1.96175 37.61729 136.47969 4.043538 0.2474189 5.213403 4.649036
"28" 0 1.01 36 3 138.1583 124 -3.22177632 37.75938 127.03075 4.253066 0.2570321 5.68848 4.848047
"29" 1 0.97 70 3 139.4073 120.1 -3.94242105 39.29461 141.44578 3.933928 0.238971 5.440202 4.870482
"30" 0 0.96 50 3 171.4 115 -4.28860526 39.66588 125.71218 3.392143 0.2099845 5.362482 4.874699
"31" 0 0.91 50 4 121.8633 122.2 -4.57039474 44.01055 126.50366 2.914111 0.181405 5.688733 4.821245
"32" 0 0.86 56 3 130.3991 119.5 -5.86013158 31.17616 125.20316 2.707486 0.1625752 5.352435 4.889009
"33" 0 0.83 50 3 133.9143 123 -8.67830263 29.8614 193.29856 2.628504 0.1731227 4.979624 5.016178
"34" 1 0.83 45 4 123.1127 120.7 -4.86335526 48.07563 215.95427 2.65283 0.1862995 4.49594 4.992888
"35" 1 0.8 68 3 153.8995 114 -7.35480263 27.46834 205.26209 2.778418 0.1907243 4.622172 4.996512
"36" 1 0.8 80 2 142.1382 117.8 -11.94153947 36.0866 192.72386 2.96296 0.2320561 4.717983 4.892908
"37" 0 0.73 31 3 173.973 106.7 -7.30802632 29.04003 191.99504 2.768074 0.213337 4.550363 5.071504
"38" 0 0.67 39 3 144.6331 109.8 -4.15302632 31.72337 159.91387 2.290645 0.1665244 3.531768 5.079191
"39" 0 0.59 42 3 115.2726 113 -0.42901667 38.08879 110.85501 2.017465 0.1412382 3.137476 5.000523
"40" 1 0.59 55 3 149.725 110 2.73575 40.82871 109.24371 2.047162 0.1664801 3.387708 5.047433
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment