ltt.plot从进化树推断物种分化的历史：分类单元数与时间的关系

NM Hallinan原著

Nick Matzke增补

今天我们分析生物类群在时间上的多样性是如何变化的。我们将用到分类单元数-时间图（Lineages-through-time plot）,该图可以用来描述物种多样化的总体趋势。我们将真实的格局与零分布进行比较，从而推断其与期望的差别。同时，我们还将学习一下几种用来推断物种丰富度随时间变化的似然方法。当然，第一步是用R模拟生成进化树。
在R中进行进化树的模拟
在tree平衡的试验中，你已经做了一些这方面的分析。在R中进行进化树模拟有几种方法。我们下面就将看到几种。ape程序包提供了进化树模拟的两个函数，rtree和rcoal. 两个函数随机生成枝长相等的进化树，但是都不能满足我们的需求。

geiger程序包中的birthdeath.tree 用处更大些。 这一函数从两个分类单元开始，基于随机生灭过程（birth death BD）进行进化树的模拟。用户给出生成率，同时给出灭亡率， 以及该过程经历的时间或者分类单元数，就可完成进化树的模拟。如果让这一函数经历相同的时间，则它会给出不同的分类单元数，如果让它给出相同的分类单元数，则经历的时间会不同。下面首先要载入 geiger程序包。

1
2
library(geiger)
tree1  <- birthdeath.tree(0.3,0,10)

以上命令将建立一个λ=0.3, μ=0 的（物种生成和绝灭速率），经过十个时间单位的进化树。理论上，这将产生 2exp(0.310) 个分类单元，但是也可能多些或少些。我们看一下其结果:

图1代码

1
2
plot(tree1)
axisPhylo()

图 1
看起来怪怪的，像一只螃蟹？ 等一下，问题可以解决：

图2代码

1
2
3
tree1<-reorder(tree1)
plot(tree1)
axisPhylo()

图 2
你到进化树是否包括大约40个分类单元？ 是否具有10个时间单位？我们来建立另外一棵进化树，这棵进化树的灭绝速率不为0. 仍然设定λ=0.35, and μ=0.05，r仍然等于 0.3

图3代码

1
2
3
tree2 <- reorder(birthdeath.tree(0.35,0.05,10))
plot(tree2)
axisPhylo()

图3
请注意，该树既包括了现存的，也包括了已经灭绝的分类单元。上一个例子中，我们没有获得绝灭的分类单元，那是因为 μ 为 0。我们可以将已经灭绝的分类单元去除。

图4代码

1
2
3
tree2.pruned<-prune.extinct.taxa(tree2)
plot(tree2.pruned)
axisPhylo()

图 4

这就是上面的结果了。进化树中只包括了现存的分类单元。我们也可以生成一个分类单元数一定的进化树。

图5代码

1
2
3
tree3 <- reorder(birthdeath.tree(0.35,0.05,taxa.stop=40))
plot(tree3)
axisPhylo()

图 5
这就是想获得的40个分类单元。由于一些分支灭绝，可能获得的分类单元数多于40个，这一过程在生成40个分类单元的时候结束。这一过程需要的时间是随机分布的。可能你注意到了，最后一个时间段重复的很厉害。这是由于当前仅包括40个分类单元数引起的。这一模拟过程并不是完全可信的，因为之后你仍然能够获得40个分类单元。我们这里不考虑这一问题。

现在，我们将会用到我编写的函数，rdbtree.n. 这一函数的参数有 λ, μ, 时长,以及分类单元数，继而依据这些参数生成一棵随机进化树。与此同时，λ 和 μ随时间的变化可以发生变化。需要注意的是，该函数本函数将给出一系列分支长度，虽然你选择的速率可能产生你所设定的分类单元数可能不能用这些枝长来生成。因此，选择速率的时候，要尽量有所参考。从网站上下载 rbdtree.n3.R 文件，或者远程读取。
source(“http://ib.berkeley.edu/courses/ib200b/labs/lab12/rbdtree.n3.R“)
之后，运行函数 rbdtree.n：

图6代码

1
2
3
tree4<-rbdtree.n(40,10,0.35,0.05)
plot(tree4)
axisPhylo()

图6
很棒！ 现在的结果大体上能够符合我们检验所需要的零分布了。这是由于从分类单元数以及时间上，都设定了，而且，这两个因素对于分类单元数-时间图，都是有很大影响的。需要注意的是，可能存在着其他的假设更为合适的情况。

分类单元数-时间图

现在，我们就来做一些分类单元数时间图。该图为一个clade拥有的子代物种数量随时间变化的图。该过程的期望是，如果μ=0 而 λ为常数，则物种数随着时间指数性的增加。偏离于该期望时，就是说明多样化是随时间而变化的。我们将调用ape程序包的一些函数来绘制这些图形。

图7代码

1
ltt.plot(tree1)

图 7
按照经验，y轴要经过对数转换。我们应用R 绘图函数的log = “y”来实现相应的功能，该参数将改变所选轴数值的显示方式，但是并不会改变数值本身。

图8代码

1
ltt.plot(tree1,log="y")

图 8
对于对数转换的图形，期望值应该为一条直线，我们可以在图形上添加一条直线，用来比较其异同。

图9代码

1
lines(c(-10,0), c(2,40), lty=2)

图 9
结果看起来确实不错。你还可以继续添加物种数-时间曲线到现有的图形上

图10代码

1
ltt.lines(tree2,col=2)

图 10
不错，但是看起来有点儿奇怪，这是由于我们两棵进化树最后的分类单元数不同。我们也可以同时绘制多棵进化树的物种数-时间曲线。

图11代码

1
mltt.plot(tree1,tree2,tree3,tree4)

图 11
物种数-时间图的零分布
为了理解多样性系数会怎样影响物种数-时间图形的形状，我们需要了解这些图形在零假设的时候的形状。我们从一段给定时间的完全随机生灭过程的零假设开始。

1
trees.yule <- replicate(100, birthdeath.tree(0.3,0,10), simplify=FALSE)

以上函数将birthdeath.tree函数重复100次，创建了一个列表。我们用现有数据建立物种数-时间图

1
mltt.plot(trees.yule, legend=FALSE)

图12代码

图 12
问题来了，每棵进化树具有不同的分类单元数，你看到的大多数情况是，物种数比较多的图，其结尾分类单元数也比较多。如果你设定了物种数一定，那么将会有另外一个问题。让我们创建一些具有相同问题的进化树。

图13代码

1
2
trees.yule2<-replicate(100, rbdtree.n(40,10,0.3,0), simplify=FALSE)
mltt.plot(trees.yule2, legend=FALSE, log="y")

图 13
可以见到本图中的Yule进化树的分布，但是在与其他分布比较的时候会遇到困难，因为所有的线都是跳跃性的上升。因此，我编写了另一个函数，用来对任意参数值给出物种数-时间的零分布。

图14代码

1
ltt.null(40, 10, 0.3, 0)

图 14
该函数的参数与yule2进化树的例子是一致的。颜色的深浅代表了不同的p值。查看我们已经生成的yule2进化树的分布可以添加到这一图形中。

图15代码

1
x <- lapply(trees.yule2,ltt.lines)

图 15
我们将结果赋值给x, 因此，该函数就不会像屏幕中打印出一些无用的信息了。可以看出，大部分线都落在我们分布的中间部分。一部分跑到了95%置信区间以外。甚至有一条落在了99%置信区间，就像我们进行100次随机模拟所期望的。
现在看一下更高的灭绝速率下，这一分布会如何变化？ 我们将μ提高到0.1, λ提高到0.4, 因此净出生率仍然不发生变化，N==2exp(rt)仍然保持不变。在此之前，我们需要为零分布创建另外一个图，之后，用新参数模拟生成100棵进化树，最后，添加到我们的图上。

图16代码

1
2
3
ltt.null(40,10,0.3,0)
trees.bd<-replicate(100,rbdtree.n(40,10,0.4,0.1),simplify=FALSE)
x <- lapply(trees.bd,ltt.lines)

图 16
可以看出，这些分布并不相同。首先，灭绝速率较高时，起始部分延伸所需的时间更长。其次，图形变得更为凹陷，其中大部分分类单元低于Yule分布，甚至有些超出了99%的分布。这就意味着，lttplot的曲线越凹陷，则其μ值，也就是灭绝速率就越高。此时，我们倾向于认为，古老的分类单元更容易灭绝。但是大部分的分类单元仍然处于95%的置信区间内，因此，这一检验的效力是十分薄弱的。

进化树与零分布的比较

我们看一下一些实际的进化树，看其是否符合我们的零分布。laser程序包中已经收录了几棵进化树，我们以此为例。第一棵进化树包含了澳大利亚69种Agamid蜥蜴，相应的数据发表在Evolution上Rabosky (2006) Evolution 60:1152-1164。

图17代码

1
2
3
data(agamids)
agamid.tree<-read.tree(text=agamids)
plot(agamid.tree)

图 17
很好，看起来这棵进化树不错。我们把它和零分布做一下比较。我们从Yule过程开始，可以用极大似然法估计λ。

1
2
3
4
5
6
7
8
9
10
11
12
13
yule.agam<-yule(agamid.tree)
\> yule.agam
$lambda
[1] 6.807443

$se
[1] 0.8316616

$loglik
[1] 283.4636

attr(,"class")
[1] "yule"

同时，估计进化树的时间长度。即根的深度。

1
btimes.agam<-sort(branching.times(agamid.tree),decreasing=TRUE)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
\> btimes.agam
      70       71      131      132       72       73       74      133 
0.323855 0.322855 0.313397 0.307173 0.294147 0.270881 0.238577 0.238520 
     109      135       75      134       76       77       85      110 
0.222756 0.219480 0.217127 0.216939 0.213018 0.213018 0.213018 0.209613 
     112      113       86       93      116      117       94      111 
0.207659 0.202293 0.194939 0.183979 0.178225 0.172959 0.167084 0.162999 
      98      118      123      115      124      106      102      128 
0.162093 0.160592 0.157785 0.157288 0.154405 0.149574 0.144083 0.141278 
     125       95      114       87      136       99       91      119 
0.141127 0.140194 0.136160 0.133171 0.131567 0.124221 0.120235 0.118270 
     103      105      107       88      137       89      127       78 
0.112825 0.105969 0.104019 0.102561 0.087830 0.086799 0.085487 0.082827 
      92      126      120       97      104      129      121      100 
0.082356 0.074543 0.074123 0.073526 0.071486 0.069944 0.067407 0.063750 
     122       96      108      130      101       80       83       81 
0.062395 0.060306 0.051672 0.050408 0.049427 0.042118 0.035946 0.020464 
      79       90       82       84 
0.009766 0.008402 0.003883 0.001482 
\>

我们在零模型中可以运用这些值。

图18代码

1
ltt.null(69,btimes.agam[1],yule.agam$lambda,0)

图 18
之后，我们可以添加上进化树的ltt.

图19代码

1
ltt.lines(agamid.tree)

图 19

该进化树的ltt线超出了99%的置信区间。在一开始符合的很好，但是在中间逐渐超越了我们的零模型的预期，最后有回到了预期的置信区间。这与我们的预期中的零分布有较大的不同；当μ逐渐升高，零分布将变得越来越凹陷，而不是当前这样。从图中大约可以看出，在0.22时间左右，存在一个物种分化异常高的时期，同时也包括约0.15时期（遗憾的是，我并不知道时间单位）。我们可以尝试拟合一个BD模型，以查看其与灭绝速率大于0的结果是否一致。

1
birthdeath(agamid.tree)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
> birthdeath(agamid.tree)
Estimation of Speciation and Extinction Rates
            with Birth-Death Models
     Phylogenetic tree: agamid.tree 
        Number of tips: 69 
              Deviance: -566.9276 
        Log-likelihood: 283.4638 
   Parameter estimates:
      d / b = 0   StdErr = 0 
      b - d = 6.807467   StdErr = 0.5881344 
   (b: speciation rate, d: extinction rate)
   Profile likelihood 95% confidence intervals:
      d / b: [0, 0.1605283]
      b - d: [5.30515, 8.569748]

我们仍然得到灭绝速率为0，因此，不要保留那种想法。因为这种情况并不总是那么重要。在laser程序包中还有两个数据，并且只有物种分化时间，plethodon （highton and larson 1979 Syst. Zool. 28:579-599）以及 warblers （Lovette and Bermingham 1999 Proc. Roy. Soc. B. Lond. 266:1629-1636）。我们对这两套数据进行再分析。

图20代码

1
2
3
data(plethodon)
yule.pleth <- pureBirth(plethodon) ##This is like yule, but for vector of branching times
ltt.null(length(plethodon)+1, max(plethodon), yule.pleth$r1, 0)

图 20
由于只给出了物种数和分化时间，我们需要自己将信息添加到ltt图中。

图21代码

1
2
3
ntaxa.pleth <- 1:length(plethodon) + 1 ## Number of lineages corresponding to times
segments(-plethodon, ntaxa.pleth, -plethodon, ntaxa.pleth-1)
segments(-plethodon, ntaxa.pleth, -c(plethodon[-1],0), ntaxa.pleth)

图 21
这一结果好于Agamid进化树的结果。从结果看来，在15到30个时间单位上，很长时间都没有发生物种分化，但是这并不出乎我们的意料。但是，有可能某一个类群在30个时间单位以前已经有很多分类单元了。我们下面就来对此进行检验。

可以对分支时间数据拟合一个BD模型：

1
bd.pleth <- bd(plethodon)
1
2
3
4
5
6
7
\> bd.pleth 
$LH
[1] -45.28264
$r1
[1] 0.01098217
$a
[1] 0.8337962

aic [1] 94.56529 这一模型的拟合结果中，灭绝速率确实大于0，因此，我们设定如下的零分布。在此之前，我们必须先进行一个LRT(似然比检验)，来确保两者之间的差异是否显著。但是你已经知道了如何进行该检验，因此我们跳过这一步。该函数给出r和a=λ/μ,因此，我们需要首先计算λ和μ。 ## 图22代码 lam.pleth <- bd.plethaic [1] 94.56529 这一模型的拟合结果中，灭绝速率确实大于0，因此，我们设定如下的零分布。在此之前，我们必须先进行一个LRT(似然比检验)，来确保两者之间的差异是否显著。但是你已经知道了如何进行该检验，因此我们跳过这一步。该函数给出r和a=λ/μ,因此，我们需要首先计算λ和μ。

图22代码

1
2
3
4
lam.pleth <- bd.plethr1 / (1-bd.pletha)mu.pleth<−lam.pleth−bd.pletha)mu.pleth<−lam.pleth−bd.plethr1
ltt.null(length(plethodon)+1 ,max(plethodon), lam.pleth, mu.pleth)
segments(-plethodon, ntaxa.pleth, -plethodon, ntaxa.pleth-1)
segments(-plethodon, ntaxa.pleth, -c(plethodon[-1],0), ntaxa.pleth)

图 22
问题： 对warblers数据重复这一分析。你估计的lambda值是怎样的？该图与零分布有什么差异？（具体说明
对进化速率变化的进化树模拟

分类单元数-时间图与零分布不同究竟意味着什么？一般的假设认为这是由于多样化速率随时间变化。你可以用我提供的函数生成随机进化树，对这种零分布进行模拟。需要提供一个向量，该向量包含了每个时间段结束的时间。也可以设定出生和灭绝速率随时间的变化的向量，每个向量给出相应时间段的出生和灭绝速率。例如：

图23代码

1
2
tree5<-rbdtree.n(40, c(3,6,10), c(0.1,1,0.05), 0)
plot(tree5)

图 23
以上程序将生成一个40个分类单元的进化树，在时间段 0 - 3， 灭绝速率为0， μ= 0， λ is 0.1； 而在3 - 6时间段为1， 在第三个时间段，也就是6之后到结束的时候，λ 为 0.05. 可以看到，大部分分支出现在进化树的中间部分。 另外：

图24代码

1
2
tree6 <- mrbdtree.n(10, 40, c(3,6,10), c(0.1,1,0.05), c(0.05,0.05,0))
mltt.plot(tree6)

图 24
将生成10棵进化树，其中的μ随时间变化。应该能很清楚的看到lttplot期望的分布格局。

同时可以建立在某一个时刻具有一定物种数的进化树及类似的零分布，而不是对末端分类单元。 还记得plethodon图吗？

图25代码

1
2
3
ltt.null(length(plethodon)+1, max(plethodon), lam.pleth, mu.pleth)
segments(-plethodon, ntaxa.pleth, -plethodon, ntaxa.pleth-1)
segments(-plethodon, ntaxa.pleth, -c(plethodon[-1],0), ntaxa.pleth)

图 25
所有内容都在我们期望的范围内，但是有一段时间较为奇怪，该段时间没有新物种形成。我在怀疑，是否整条ltt图都超出了我们的零分布，虽然说每个时刻都处于我们的零分布的预期。让我们来看一下，当将20个时间单元设定到物种分布为11个分类单元时，零分布会变成什么样子。

图26代码

1
2
3
ltt.null(c(11, length(plethodon)+1), c(-20,0)+max(plethodon), lam.pleth, mu.pleth)
segments(-plethodon, ntaxa.pleth, -plethodon, ntaxa.pleth-1)
segments(-plethodon, ntaxa.pleth, -c(plethodon[-1],0), ntaxa.pleth)

图 26
在我看来，这些也是可以的。事实上，除了靠右面的第二部分被强制设定为11个种之外，我没有看出这些零分布的显著差别。

出生率的下降以及灭绝率的上升
好了，我们已经看到一些进化树的物种分化速率在初期似乎比后期要高。这一现象不能用常规的随机生灭过程解释。这一格局常常在大型进化树中遇到，通常被认为是适应辐射的标志。换句话说，这一格局可以解释为，演化早期拥有较高的物种分化速率， 而随着生态位被逐步填充， 物种分化速率降低。 但是， 也可以用灭绝速率升高来解释。

我们看一下，我们的猜想是否正确。我们将会建立一对主观模型（arbitrary model）. 例如，我们将以60个种结束，并且将进化树分成两个时间段；我们将第一阶段设定为一个时间单位，后一个阶段设定为9个时间单位。同时，我们认为多样化速率r在第一阶段（r1）在第一阶段是第二阶段（r2）的10倍。相应的物种数可以通过公式计算，N=2*exp(r1t1)*exp(r2t2),那么我们将其中的已知数填入，解方程，得到r2=log(30)/19=0.179, and r1=1.79.
下面，我们看一下随着 μ的上升或 λ的下降，我们的零分布将会怎样变化？假如，我们设定μ1=0.05, λ1=1.84.之后，我们的第一个模型λ下降，令μ保持不变，而令λ2=.229，对于模型二，我们令λ保持不变，而μ2=1.661。我们将这两种模型同时和Yule分布进行比较。Yule分布的λ=log(30)/10=0.34

图27代码

1
2
3
ltt.null(60, 10, 0.34, 0)
model.lambda <- mrbdtree.n(10, 60, c(1,10), c(1.84,0.229), 0.05)
lapply(model.lambda, ltt.lines) -> x

图 27

看起来非常不错，分类群在早期迅速增加，之后回到了预期的曲线中。或许两个时间段在长度上接近，我们会得到更为平滑的曲线。 我们现在看一下，如果灭绝速率升高，将会有什么发生：

图28代码

1
2
3
ltt.null(60, 10, 0.34, 0)
model.mu <- mrbdtree.n(10, 60, c(1,10), 1.84, c(0.05, 1.661))
lapply(model.mu, ltt.lines) -> x

图 28

Huh, 完全失效了。这表明，我们之间看到的格局（例如Agamid进化树），物种数在早期的增加，实际上是由于物种分化速率降低造成的。 我们无法提高灭绝速率，同时又生成一个凹陷的ltt图，虽然可以通过另一些参数的组合完成这一功能。
对这种格局的另外一种解释就是取样问题。我们认为随机取样能够降低最近的物种分化时间，但是对早期的分化则影响很小。这是由于古老的分支类群往往很大，从而在随机抽样的时候很难去除。在处理隐性物种（cryptic species）或者更高等级的分类单元时，这种情况变得尤为严重。我们下面就做一个80个物种的模拟数据，并且随机提取出20个种。

图29代码

1
2
3
ltt.null(60, 10, 0.34, 0)
make.samp <- function () prune.random.taxa(rbdtree.n(100, 10, 0.3912, 0), 40)
x <- replicate(10, ltt.lines(make.samp()))

图 29
我们看到了，上面的函数几乎不起什么作用。

用来探测分化速率变化的似然模型
laser程序包中有函数可以用来检验多样性速率的变化。第一个模型与Yule模型比较，其中的物种分化速率随时间是恒定的，而实际情况中，物种分化速率存在着发生者变化。（Rabosky 2006, Evolution 60:1152-1164）. 我们来检验一下模拟出的早期进化速率很高的那棵进化树。

首先，我们用单一速率检验该模型。所有这些模型都对分化时间进行拟合，而不是整个进化树，因此，我们需要创建一个物种分支时间的向量。

1
2
times.lambda <- branching.times(model.lambda[[1]])
pureBirth(times.lambda)
1
2
3
4
5
6
7
8
9
10
11
12
\> pureBirth(times.lambda)
$LH
      61 
36.49077

$aic
       61 
-70.98154

$r1
       61 
0.2117251

看起来非常合理。λ 比我们预想的要低一些，但这可能是由于曲线的凹凸度的后果。现在我们试验一下两个速率的模型。

1
2
yule2rate(times.lambda)
yule2rate(times.lambda)
1
2
        LH      st1       r1        r2       aic
6 45.92311 9.010908 1.882302 0.1887342 -85.84621

看一下这一模型的结果。LH是Log likehihood. stl是物种分户速率变化的转折点。但是这些结果与我们时间上所做的转折如何比较？r1 和 r2分别是每个时间段的物种分化速率。 AIC是Akaike’s information criterion, 这是一种比较似然率的另外一种方法；即使在模型为非巢式的时候，仍然很有效。 AIC定义为 2k-2log(likelihood)，需要选择AIC值最低的模型。那么按照这一标准，什么才是最好的模型？ 此外，也可以用似然比检验（Likelihood Ratio Test）来选择模型。

如果选择了两速率模型作为最好的，那么可以继续用 yule3rate, yule4rate, 或yule5rate. 同样可以用fitdAICrc来比较不同模型，虽然我没有看出在比较这些模型时的优势。
问题2。 Warbler进化树的最优模型是什么？ 该模型的参数是什么？

模拟物种分化速率连续变化的进化树

在Rabosky和Lovette 2008（Evolution 62 - 8: 1866 - 1875）中，作者提出了一个物种分化速率模型，该模型的多样性速率指数递减。他们提出了三种模型： λ(t)=λ0exp(-kt)， μ 对于 SPVAR是常数。λ 为常数， μ(t)=μ0(1-exp(-zt)) 对于EXVAR是常数，两者都是BOTHVAR的变量。我们估计成百上千的不同点上的速率来模拟这些数据。
对于SPVAR模型，期望的分类单元数为 N=2*exp(λ0[1-exp(-kt)]/k-μt)
因此，如果我们希望模拟的过程持续10个时间单位，产生60个物种，灭绝速率0.05，在开始的时候物种分化速率比后期高十倍，则满足 

$exp(-10k)=0.1$

$N=2λ0exp(-10k)/k-10μ)$

因此 k = 0.23, λ0= 0.997

如果我们将时间100等分，就将获得每个时刻的平均值 
λ0exp(-kt)[exp(-kΔt)-1]/kΔt=1.009exp(-0.23t)。 之后创建一个向量

1
2
lambdas <- 1.009*exp(-0.23*1:100/10)
lambdas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
  [1] 0.9860578 0.9636373 0.9417266 0.9203141 0.8993884 0.8789386 0.8589537
  [8] 0.8394232 0.8203368 0.8016844 0.7834561 0.7656422 0.7482334 0.7312205
 [15] 0.7145943 0.6983462 0.6824676 0.6669500 0.6517852 0.6369652 0.6224822
 [22] 0.6083285 0.5944966 0.5809792 0.5677692 0.5548596 0.5422434 0.5299142
 [29] 0.5178652 0.5060903 0.4945830 0.4833374 0.4723475 0.4616075 0.4511117
 [36] 0.4408546 0.4308306 0.4210346 0.4114613 0.4021057 0.3929628 0.3840278
 [43] 0.3752960 0.3667627 0.3584234 0.3502738 0.3423094 0.3345261 0.3269198
 [50] 0.3194865 0.3122222 0.3051230 0.2981853 0.2914053 0.2847794 0.2783043
 [57] 0.2719763 0.2657923 0.2597488 0.2538428 0.2480710 0.2424305 0.2369182
 [64] 0.2315313 0.2262668 0.2211221 0.2160943 0.2111809 0.2063792 0.2016866
 [71] 0.1971008 0.1926192 0.1882395 0.1839594 0.1797766 0.1756889 0.1716942
 [78] 0.1677903 0.1639752 0.1602468 0.1566032 0.1530424 0.1495626 0.1461619
 [85] 0.1428386 0.1395908 0.1364168 0.1333150 0.1302838 0.1273215 0.1244265
 [92] 0.1215973 0.1188325 0.1161306 0.1134900 0.1109096 0.1083877 0.1059233
 [99] 0.1035148 0.1011612

>

之后，可以建立一些进化树，并且绘制出来。

图30代码

1
2
3
trees.spvar <- mrbdtree.n(10, 60, 1:100/10, lambdas, 0.05)
ltt.null(60, 10, 0.34, 0)
lapply(trees.spvar, ltt.lines) -> x

图 30
这与我们想要解释的真实进化树极为相似，在早期都具有较高的分化速率。我们也可以建立一个相反的模型，在每个区间μ不断升高，而r保持不变，但是λ为常数。这可以通过 λ2=λ0 - μ1以及 μ2(t)=λ0 - λ1(t)实现。 我们可以建立这些进化树，并且与SPVAR进化树进行比较。

图31代码

1
2
3
mus <- 0.997-lambdas
trees.exvar <- mrbdtree.n(10, 60, 1:100/10, 0.947, mus)
lapply(trees.exvar,ltt.lines,col=3)->x

图 31
问题3 EXVAR模拟的进化树总体上应该是什么形状？ 它们看起来是否与我们说的早期物种分化速率较高的进化树相似？ 这些现象说明什么过程能够生成这些进化树?

用似然率连续变化模型的拟合
如果你注意到，上面提到的论文的作者，就是laser程序包的作者（现在是Huelsenbeck实验室的博士后， 很快就要成为Michigan大学的老师）， 那么在这个程序包中发现几个用来拟合物种分化速率连续变化的模型的函数就不奇怪了。该程序包中包含了拟合SPVAR, EXVAR和BOTHVAR的函数， 也就是我在前一段中提到的。同时，也有随着物种丰富度的升高，物种分化速率下降的模型。 我们仅简要介绍其中的内容，如果你有兴趣，可以查看其选项。
我们用SPVAR模型来拟合SPVAR模拟生成的进化树，同时看一下它是如何工作的。首先检验一下其物种分化速率为常量的零假设。

1
2
times.spvar <- branching.times(trees.spvar[[1]])
bd(times.spvar)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
$LH
     61 
35.7508
$r1
      61 
0.209041
$a
[1] 0
$aic
       61 
-67.50159

之后，检查SPVAR模型
fitSPVAR(times.spvar)
\> fitSPVAR(times.spvar)
$model
[1] "SPVAR"

$LH
[1] 47.27787

$aic
[1] -88.55573

$lam0
[1] 1.081366

$k
[1] 0.2593436

$mu0
[1] 0.001

你是如何估计模型参数的？ 你青睐那种模型？ 我们这里暂不涉及EXVAR和BOTHVAR，但是你应该了解到了如何使用它们。我们简要查看一下密度依赖多样化模型。

1
DDX(times.spvar)
1
2
3
4
5
6
7
8
9
10
11
12
13
\> DDX(times.spvar)
$LH
[1] 45.96522

$aic
[1] -87.93044

$r1
[1] 2.532487

$xparam
[1] 0.7286594

该模型拟合的效果好吗？ 如果想要理解这一模型的意思，则必须查看laser程序包的帮助文件。不仅如此，如果想要了解R任何函数的使用方法，都需要查看帮助文件。 我只介绍了其中很少一部分函数的用法，你可以参考更多的选项及功能。

References