Exercise 11: State dependent diversification analysis using diversitree

In the following we will use the MuSSE method - a multi-state character extension of the BiSSE model - to test various hypotheses about diversification in the fish family Terapontidae.

First, we can load the packages that we will be using today. These consist of ape, diversitree, phytools, and geiger.

library(ape)
library (diversitree)
library(phytools)
## Loading required package: maps
library(geiger)

Next, download the two files we will be using for today’s exercise:

  1. Terapontidae.tre (the tree file).
  2. SSE.csv (the data file).
phy<-read.tree("Terapontidae.tre")
phy
## 
## Phylogenetic tree with 38 tips and 37 internal nodes.
## 
## Tip labels:
##  Terapon_jarbua, Pelsartia_humeralis, Mesopristes_cancellatus, Mesopristes_argenteus, Rhyncopelates_oxyrhynchus, Terapon_puta, ...
## 
## Rooted; includes branch lengths.
plotTree(phy,ftype="i")

ltt(phy,show.tree=TRUE,lwd=3)

## Object of class "ltt" containing:
## 
## (1) A phylogenetic tree with 38 tips and 37 internal nodes.
## 
## (2) Vectors containing the number of lineages (ltt) and branching times (times) on the tree.
## 
## (3) A value for Pybus & Harvey's "gamma" statistic of 0.1285, p-value = 0.8977.

Now, let’s read in our data from file. Unfortunately, diversitree requires that states be numerically coded 1, 2, 3, etc.

xx <- read.csv("SSE.csv", row.names = 1)
xx
##                           Diet     Animal      Plant         PC1
## Amniataba_affinis            2  0.2818512 -0.2818512  0.29135223
## Amniataba_caudavittatus      1  2.5866893 -2.5866893 -0.63360293
## Amniataba_percoides          2 -0.2330491  0.2330491  0.19121023
## Bidyanus_bidyanus            3 -2.1972246  2.1972246  0.48236174
## Bidyanus_welchi              2  0.8472979 -0.8472979  0.28913974
## Hannia_greenwayi_1           2  1.1970524 -1.1970524 -0.18609133
## Helotes_octolineatus         3 -3.4094962  3.3776909  1.33077974
## Helotes_sexlineatus          3 -1.2656664  1.2656664  1.30024471
## Hephaestus_carbo             1  4.2545990 -4.2545990 -0.65802938
## Hephaestus_epirrhinos        1  1.3862944 -1.4051510 -1.39516149
## Hephaestus_fuliginosus       2 -0.7172447  0.7172447 -0.74690337
## Hephaestus_habbemai          1  2.5866893 -2.5866893 -0.64295076
## Hephaestus_jenkinsi          2 -0.1966311  0.1966311 -0.56697125
## Hephaestus_transmontanus     1  5.5174529 -5.5174529 -0.30074728
## Hephaestus_tulliensis        2 -1.1914483  1.1914483  0.01499809
## Leiopotherapon_aheneus       2 -0.7583712  0.7583712  0.02378382
## Leiopotherapon_plumbeus      1  2.5866893 -2.5866893 -0.80360103
## Leiopotherapon_unicolor      1  2.1972246 -2.1751973 -1.06302329
## Mesopristes_argenteus        1  2.9444390 -2.8438517 -0.14004817
## Mesopristes_cancellatus      2  1.0986123 -1.0986123  0.02272593
## Pelates_quadrilineatus       1  4.8202816 -4.8202816  0.27274089
## Pelates_sexlineatus          1  3.9441335 -3.9441335  0.49450293
## Pelsartia_humeralis          1  3.1780538 -3.1780538 -0.45071433
## Pingalla_gilberti            3 -1.6434218  1.6434218  1.48379891
## Pingalla_lorentzi            2 -0.6722197  0.6766921  0.99779512
## Pingalla_midgleyi            3 -3.2044128  3.2044128  1.35889133
## Rhyncopelates_oxyrhynchus    1  6.9067548 -6.9067548 -0.38469018
## Scortum_barcoo               3 -2.9873640  2.9873640  0.50486583
## Scortum_hillii               3 -3.8918203  3.8918203  0.52885797
## Scortum_ogilbyi              3 -2.4423470  2.4423470  0.56911574
## Scortum_parviceps            3 -3.4760987  3.4760987  0.48617250
## Syncomistes_butleri          3 -6.2126061  6.2126061  1.33305640
## Syncomistes_rastellus        3 -2.4838238  2.4838238  1.03094100
## Syncomistes_trigonicus       3 -3.5845472  3.5845472  1.36530955
## Terapon_jarbua               1  4.5951198 -4.5951198 -1.61420290
## Terapon_puta                 1  3.2875724 -3.2875724 -0.19368748
## Terapon_theraps              1  6.9067548 -6.9067548 -0.61994670
## Variichthys_lacustris        2 -0.1201443  0.1201443 -0.49414189
##                                   PC2
## Amniataba_affinis         -0.21526981
## Amniataba_caudavittatus    0.15468468
## Amniataba_percoides        0.07985197
## Bidyanus_bidyanus         -0.78346644
## Bidyanus_welchi           -0.77832681
## Hannia_greenwayi_1        -0.29919632
## Helotes_octolineatus       0.98925031
## Helotes_sexlineatus        0.85397920
## Hephaestus_carbo           0.22459829
## Hephaestus_epirrhinos     -0.26280957
## Hephaestus_fuliginosus     0.07468081
## Hephaestus_habbemai        0.39513150
## Hephaestus_jenkinsi        0.13175860
## Hephaestus_transmontanus  -0.03448754
## Hephaestus_tulliensis      0.03155424
## Leiopotherapon_aheneus    -0.02685418
## Leiopotherapon_plumbeus   -0.34694022
## Leiopotherapon_unicolor   -0.13183385
## Mesopristes_argenteus     -0.22825584
## Mesopristes_cancellatus   -0.34447877
## Pelates_quadrilineatus    -0.24877258
## Pelates_sexlineatus       -0.59433403
## Pelsartia_humeralis        0.07651360
## Pingalla_gilberti          0.56462925
## Pingalla_lorentzi          0.60814139
## Pingalla_midgleyi          0.09966254
## Rhyncopelates_oxyrhynchus -0.25877534
## Scortum_barcoo            -0.01681409
## Scortum_hillii            -0.11574660
## Scortum_ogilbyi           -0.07753850
## Scortum_parviceps         -0.20918208
## Syncomistes_butleri        0.40985322
## Syncomistes_rastellus     -1.04410803
## Syncomistes_trigonicus    -0.51744194
## Terapon_jarbua             1.41833043
## Terapon_puta              -0.26355821
## Terapon_theraps           -0.62775691
## Variichthys_lacustris      0.40517581
Diet <- setNames(xx$Diet,rownames(xx))
Diet
##         Amniataba_affinis   Amniataba_caudavittatus 
##                         2                         1 
##       Amniataba_percoides         Bidyanus_bidyanus 
##                         2                         3 
##           Bidyanus_welchi        Hannia_greenwayi_1 
##                         2                         2 
##      Helotes_octolineatus       Helotes_sexlineatus 
##                         3                         3 
##          Hephaestus_carbo     Hephaestus_epirrhinos 
##                         1                         1 
##    Hephaestus_fuliginosus       Hephaestus_habbemai 
##                         2                         1 
##       Hephaestus_jenkinsi  Hephaestus_transmontanus 
##                         2                         1 
##     Hephaestus_tulliensis    Leiopotherapon_aheneus 
##                         2                         2 
##   Leiopotherapon_plumbeus   Leiopotherapon_unicolor 
##                         1                         1 
##     Mesopristes_argenteus   Mesopristes_cancellatus 
##                         1                         2 
##    Pelates_quadrilineatus       Pelates_sexlineatus 
##                         1                         1 
##       Pelsartia_humeralis         Pingalla_gilberti 
##                         1                         3 
##         Pingalla_lorentzi         Pingalla_midgleyi 
##                         2                         3 
## Rhyncopelates_oxyrhynchus            Scortum_barcoo 
##                         1                         3 
##            Scortum_hillii           Scortum_ogilbyi 
##                         3                         3 
##         Scortum_parviceps       Syncomistes_butleri 
##                         3                         3 
##     Syncomistes_rastellus    Syncomistes_trigonicus 
##                         3                         3 
##            Terapon_jarbua              Terapon_puta 
##                         1                         1 
##           Terapon_theraps     Variichthys_lacustris 
##                         1                         2

Just as with the birth-death model we fit before, we can specify sampling fractions that are less than 100%. Here, though, we need to specify the sampling fraction for each character state:

samplingf<-c(0.58, 0.78, 0.85)

The MuSSE method also comes with a helpfing function called starting.point.musse which returns “reasonable” starting values for the parameters in the analysis.

p <- starting.point.musse(phy, k=3)
p
##    lambda1    lambda2    lambda3        mu1        mu2        mu3 
## 0.08725939 0.08725939 0.08725939 0.01013885 0.01013885 0.01013885 
##        q12        q13        q21        q23        q31        q32 
## 0.01542411 0.01542411 0.01542411 0.01542411 0.01542411 0.01542411

If we check, we will see that the birth & death starting values are the MLEs by Nee et al.’s method:

obj<-birthdeath(phy)
bd(obj)
##          b          d 
## 0.08726392 0.01014791

Now we’re ready to starting fitting our models. Remember, in diversitree we have to make a likelihood function, and then optimize it.

To start, we can write down our MuSSE likelihood function as follows:

lik.musse<-make.musse(phy,Diet,3,sampling.f=samplingf)
lik.musse
## MuSSE likelihood function:
##   * Parameter vector takes 12 elements:
##      - lambda1, lambda2, lambda3, mu1, mu2, mu3, q12, q13, q21,
##        q23, q31, q32
##   * Function takes arguments (with defaults)
##      - pars: Parameter vector
##      - condition.surv [TRUE]: Condition likelihood on survial?
##      - root [ROOT.OBS]: Type of root treatment
##      - root.p [NULL]: Vector of root state probabilities
##      - intermediates [FALSE]: Also return intermediate values?
##   * Phylogeny with 38 tips and 37 nodes
##      - Taxa: Terapon_jarbua, Pelsartia_humeralis, ...
##   * Reference:
##      - FitzJohn (submitted)
## R definition:
## function (pars, condition.surv = TRUE, root = ROOT.OBS, root.p = NULL, 
##     intermediates = FALSE)

This is our general likelihood function. We can constrain it in different ways to fit the models that we are interested in.

To start, we will fit a model - the “null” model in a sense - in which all birth & death rates are equal between states, the character evolution is ordered (that is 1<->2<->3), and there is a single character transition rate:

lik.null<-constrain(lik.musse,lambda2 ~ lambda1, lambda3 ~ lambda1, 
    mu2 ~ mu1, mu3 ~ mu1, q13 ~ 0, q21 ~ q12, q23 ~ q12, q31 ~ 0, 
    q32 ~ q12)
lik.null
## MuSSE likelihood function:
##   * Parameter vector takes 3 elements:
##      - lambda1, mu1, q12
##   * Function constrained (original took 12 elements):
##      - lambda2 ~ lambda1
##      - lambda3 ~ lambda1
##      - mu2 ~ mu1
##      - mu3 ~ mu1
##      - q13 ~ 0
##      - q21 ~ q12
##      - q23 ~ q12
##      - q31 ~ 0
##      - q32 ~ q12
##   * Function takes arguments (with defaults)
##      - pars: Parameter vector
##      - ...: Additional arguments to underlying function
##      - pars.only [FALSE]: Return full parameter vector?
##   * Phylogeny with 38 tips and 37 nodes
##      - Taxa: Terapon_jarbua, Pelsartia_humeralis, ...
##   * Reference:
##      - FitzJohn (submitted)
## R definition:
## function (pars, ..., pars.only = FALSE)

To fit, we need to supply starting conditions, x.init. We can do this most easily using:

p[argnames(lik.null)]
##    lambda1        mu1        q12 
## 0.08725939 0.01013885 0.01542411

from our call to starting.point.musse earlier.

## fit the model
fit.null<-find.mle(lik.null,x.init=p[argnames(lik.null)])
fit.null
## $par
##    lambda1        mu1        q12 
## 0.12196704 0.04332135 0.04160960 
## 
## $lnLik
## [1] -159.3551
## 
## $counts
## [1] 167
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
## NULL
## 
## $method
## [1] "subplex"
## 
## $par.full
##    lambda1    lambda2    lambda3        mu1        mu2        mu3 
## 0.12196704 0.12196704 0.12196704 0.04332135 0.04332135 0.04332135 
##        q12        q13        q21        q23        q31        q32 
## 0.04160960 0.00000000 0.04160960 0.04160960 0.00000000 0.04160960 
## 
## $func.class
## [1] "constrained" "musse"       "dtlik"       "function"   
## 
## attr(,"func")
## MuSSE likelihood function:
##   * Parameter vector takes 3 elements:
##      - lambda1, mu1, q12
##   * Function constrained (original took 12 elements):
##      - lambda2 ~ lambda1
##      - lambda3 ~ lambda1
##      - mu2 ~ mu1
##      - mu3 ~ mu1
##      - q13 ~ 0
##      - q21 ~ q12
##      - q23 ~ q12
##      - q31 ~ 0
##      - q32 ~ q12
##   * Function takes arguments (with defaults)
##      - pars: Parameter vector
##      - ...: Additional arguments to underlying function
##      - pars.only [FALSE]: Return full parameter vector?
##   * Phylogeny with 38 tips and 37 nodes
##      - Taxa: Terapon_jarbua, Pelsartia_humeralis, ...
##   * Reference:
##      - FitzJohn (submitted)
## R definition:
## function (pars, ..., pars.only = FALSE)  
## attr(,"class")
## [1] "fit.mle.musse" "fit.mle"

Next, let’s fit the most complex model in which all rates of speciation & extinction depend on the character state for our multi-state character:

fit.full<-find.mle(lik.musse,x.init=p[argnames(lik.musse)])
fit.full
## $par
##      lambda1      lambda2      lambda3          mu1          mu2 
## 3.333387e-02 1.565995e-01 2.277772e-01 1.344360e-02 7.284552e-02 
##          mu3          q12          q13          q21          q23 
## 1.254970e-01 5.912277e-07 6.072168e-03 1.054161e-01 1.122270e-08 
##          q31          q32 
## 4.499693e-07 6.969455e-02 
## 
## $lnLik
## [1] -151.4568
## 
## $counts
## [1] 2954
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
## NULL
## 
## $method
## [1] "subplex"
## 
## $func.class
## [1] "musse"    "dtlik"    "function"
## 
## attr(,"func")
## MuSSE likelihood function:
##   * Parameter vector takes 12 elements:
##      - lambda1, lambda2, lambda3, mu1, mu2, mu3, q12, q13, q21,
##        q23, q31, q32
##   * Function takes arguments (with defaults)
##      - pars: Parameter vector
##      - condition.surv [TRUE]: Condition likelihood on survial?
##      - root [ROOT.OBS]: Type of root treatment
##      - root.p [NULL]: Vector of root state probabilities
##      - intermediates [FALSE]: Also return intermediate values?
##   * Phylogeny with 38 tips and 37 nodes
##      - Taxa: Terapon_jarbua, Pelsartia_humeralis, ...
##   * Reference:
##      - FitzJohn (submitted)
## R definition:
## function (pars, condition.surv = TRUE, root = ROOT.OBS, root.p = NULL, 
##     intermediates = FALSE)  
## attr(,"class")
## [1] "fit.mle.musse" "fit.mle"

There are also lots of other intermediate models between a completely flexible model, and our simplest model. In an empirical study we should identify these models a priori rather than dredging our data through each model to identify one with good fit. In this case, we have decided to fit models in which only the speciation rate (λ) varies between states, only the extinction rate (μ) varies, and finally, one in which neither λ nor μ vary, but the transition rates differ between types of transitions (e.g., ordered, unordered, etc.).

First, the speciation rate, λ:

## create likelihood function
lik.lambda <- constrain(lik.musse, mu2 ~ mu1, mu3 ~ mu1, q13 ~ 0, q21 ~ q12, 
    q23 ~ q12, q31 ~ 0, q32 ~ q12)
## fit it
fit.lambda <- find.mle(lik.lambda, p[argnames(lik.lambda)])
fit.lambda
## $par
##      lambda1      lambda2      lambda3          mu1          q12 
## 7.640919e-02 1.021870e-01 1.728302e-01 2.224129e-07 3.773011e-02 
## 
## $lnLik
## [1] -157.5085
## 
## $counts
## [1] 285
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
## NULL
## 
## $method
## [1] "subplex"
## 
## $par.full
##      lambda1      lambda2      lambda3          mu1          mu2 
## 7.640919e-02 1.021870e-01 1.728302e-01 2.224129e-07 2.224129e-07 
##          mu3          q12          q13          q21          q23 
## 2.224129e-07 3.773011e-02 0.000000e+00 3.773011e-02 3.773011e-02 
##          q31          q32 
## 0.000000e+00 3.773011e-02 
## 
## $func.class
## [1] "constrained" "musse"       "dtlik"       "function"   
## 
## attr(,"func")
## MuSSE likelihood function:
##   * Parameter vector takes 5 elements:
##      - lambda1, lambda2, lambda3, mu1, q12
##   * Function constrained (original took 12 elements):
##      - mu2 ~ mu1
##      - mu3 ~ mu1
##      - q13 ~ 0
##      - q21 ~ q12
##      - q23 ~ q12
##      - q31 ~ 0
##      - q32 ~ q12
##   * Function takes arguments (with defaults)
##      - pars: Parameter vector
##      - ...: Additional arguments to underlying function
##      - pars.only [FALSE]: Return full parameter vector?
##   * Phylogeny with 38 tips and 37 nodes
##      - Taxa: Terapon_jarbua, Pelsartia_humeralis, ...
##   * Reference:
##      - FitzJohn (submitted)
## R definition:
## function (pars, ..., pars.only = FALSE)  
## attr(,"class")
## [1] "fit.mle.musse" "fit.mle"

Next, the extinction rate, μ:

lik.mu <- constrain(lik.musse, lambda2 ~ lambda1, lambda3 ~ lambda1, 
    q13 ~ 0, q21 ~ q12, q23 ~ q12, q31 ~ 0, q32 ~ q12)
fit.mu <- find.mle(lik.mu, p[argnames(lik.mu)])
fit.mu
## $par
##      lambda1          mu1          mu2          mu3          q12 
## 1.219280e-01 2.641176e-02 1.484539e-01 2.356453e-06 4.967633e-02 
## 
## $lnLik
## [1] -158.1715
## 
## $counts
## [1] 530
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
## NULL
## 
## $method
## [1] "subplex"
## 
## $par.full
##      lambda1      lambda2      lambda3          mu1          mu2 
## 1.219280e-01 1.219280e-01 1.219280e-01 2.641176e-02 1.484539e-01 
##          mu3          q12          q13          q21          q23 
## 2.356453e-06 4.967633e-02 0.000000e+00 4.967633e-02 4.967633e-02 
##          q31          q32 
## 0.000000e+00 4.967633e-02 
## 
## $func.class
## [1] "constrained" "musse"       "dtlik"       "function"   
## 
## attr(,"func")
## MuSSE likelihood function:
##   * Parameter vector takes 5 elements:
##      - lambda1, mu1, mu2, mu3, q12
##   * Function constrained (original took 12 elements):
##      - lambda2 ~ lambda1
##      - lambda3 ~ lambda1
##      - q13 ~ 0
##      - q21 ~ q12
##      - q23 ~ q12
##      - q31 ~ 0
##      - q32 ~ q12
##   * Function takes arguments (with defaults)
##      - pars: Parameter vector
##      - ...: Additional arguments to underlying function
##      - pars.only [FALSE]: Return full parameter vector?
##   * Phylogeny with 38 tips and 37 nodes
##      - Taxa: Terapon_jarbua, Pelsartia_humeralis, ...
##   * Reference:
##      - FitzJohn (submitted)
## R definition:
## function (pars, ..., pars.only = FALSE)  
## attr(,"class")
## [1] "fit.mle.musse" "fit.mle"

Next, we can fit a model in which λ & μ vary, but in which character transitions are ordered:

lik.lambda.mu <- constrain(lik.musse, q13 ~ 0, q21 ~ q12, q23 ~ q12, 
    q31 ~ 0, q32 ~ q12)
fit.lambda.mu <- find.mle(lik.lambda.mu, p[argnames(lik.lambda.mu)])
fit.lambda.mu
## $par
##      lambda1      lambda2      lambda3          mu1          mu2 
## 4.919620e-02 2.753739e-01 3.595936e-01 1.138420e-07 2.137111e-01 
##          mu3          q12 
## 3.217308e-01 4.928915e-02 
## 
## $lnLik
## [1] -155.2664
## 
## $counts
## [1] 961
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
## NULL
## 
## $method
## [1] "subplex"
## 
## $par.full
##      lambda1      lambda2      lambda3          mu1          mu2 
## 4.919620e-02 2.753739e-01 3.595936e-01 1.138420e-07 2.137111e-01 
##          mu3          q12          q13          q21          q23 
## 3.217308e-01 4.928915e-02 0.000000e+00 4.928915e-02 4.928915e-02 
##          q31          q32 
## 0.000000e+00 4.928915e-02 
## 
## $func.class
## [1] "constrained" "musse"       "dtlik"       "function"   
## 
## attr(,"func")
## MuSSE likelihood function:
##   * Parameter vector takes 7 elements:
##      - lambda1, lambda2, lambda3, mu1, mu2, mu3, q12
##   * Function constrained (original took 12 elements):
##      - q13 ~ 0
##      - q21 ~ q12
##      - q23 ~ q12
##      - q31 ~ 0
##      - q32 ~ q12
##   * Function takes arguments (with defaults)
##      - pars: Parameter vector
##      - ...: Additional arguments to underlying function
##      - pars.only [FALSE]: Return full parameter vector?
##   * Phylogeny with 38 tips and 37 nodes
##      - Taxa: Terapon_jarbua, Pelsartia_humeralis, ...
##   * Reference:
##      - FitzJohn (submitted)
## R definition:
## function (pars, ..., pars.only = FALSE)  
## attr(,"class")
## [1] "fit.mle.musse" "fit.mle"

Finally, we can fit a model in which λ & μ are constant, but in which the transition process is full flexible (i.e., unordered):

lik.unordered <- constrain(lik.musse, lambda2 ~ lambda1, lambda3 ~ lambda1, 
    mu2 ~ mu1, mu3 ~ mu1)
fit.unordered <- find.mle(lik.unordered, p[argnames(lik.unordered)])
fit.unordered
## $par
##      lambda1          mu1          q12          q13          q21 
## 1.182333e-01 3.697392e-02 1.781070e-06 3.555106e-03 1.014136e-01 
##          q23          q31          q32 
## 1.164671e-06 9.953547e-08 7.091252e-02 
## 
## $lnLik
## [1] -156.6166
## 
## $counts
## [1] 651
## 
## $convergence
## [1] 0
## 
## $message
## NULL
## 
## $hessian
## NULL
## 
## $method
## [1] "subplex"
## 
## $par.full
##      lambda1      lambda2      lambda3          mu1          mu2 
## 1.182333e-01 1.182333e-01 1.182333e-01 3.697392e-02 3.697392e-02 
##          mu3          q12          q13          q21          q23 
## 3.697392e-02 1.781070e-06 3.555106e-03 1.014136e-01 1.164671e-06 
##          q31          q32 
## 9.953547e-08 7.091252e-02 
## 
## $func.class
## [1] "constrained" "musse"       "dtlik"       "function"   
## 
## attr(,"func")
## MuSSE likelihood function:
##   * Parameter vector takes 8 elements:
##      - lambda1, mu1, q12, q13, q21, q23, q31, q32
##   * Function constrained (original took 12 elements):
##      - lambda2 ~ lambda1
##      - lambda3 ~ lambda1
##      - mu2 ~ mu1
##      - mu3 ~ mu1
##   * Function takes arguments (with defaults)
##      - pars: Parameter vector
##      - ...: Additional arguments to underlying function
##      - pars.only [FALSE]: Return full parameter vector?
##   * Phylogeny with 38 tips and 37 nodes
##      - Taxa: Terapon_jarbua, Pelsartia_humeralis, ...
##   * Reference:
##      - FitzJohn (submitted)
## R definition:
## function (pars, ..., pars.only = FALSE)  
## attr(,"class")
## [1] "fit.mle.musse" "fit.mle"

Note that there are other models we might have fit, depending on our a priori biological hypotheses. The fact that there are so many possible models makes identifying those that are sensible, from the perspective of our data and question, very important.

It is straightforward to compare our fitted models. We can do that firstly using the generic function anova as follows:

AnovaResults <- anova(fit.null, 
    all.different=fit.full, 
    free.lambda=fit.lambda, 
    free.mu=fit.mu, 
    free.lambda.mu=fit.lambda.mu, 
    free.q=fit.unordered)
AnovaResults
##                Df   lnLik    AIC   ChiSq Pr(>|Chi|)  
## minimal         3 -159.35 324.71                     
## all.different  12 -151.46 326.91 15.7966    0.07125 .
## free.lambda     5 -157.51 325.02  3.6931    0.15778  
## free.mu         5 -158.17 326.34  2.3671    0.30618  
## free.lambda.mu  7 -155.27 324.53  8.1773    0.08529 .
## free.q          8 -156.62 329.23  5.4771    0.36047  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This approach uses a Χ2-test with the likelihood-ratio compared to our “null” model; however we can also use information theory and AIC to assess the relative fit of our different models, penalizing for their parameterization. We can compute the relative weight of evidence in favor of each of our different hypotheses using AIC weights:

aicw(setNames(AnovaResults$AIC,rownames(AnovaResults)))
##                     fit     delta          w
## minimal        324.7102 0.1773376 0.26116649
## all.different  326.9136 2.3807677 0.08678580
## free.lambda    325.0171 0.4842196 0.22401591
## free.mu        326.3431 1.8101975 0.11543737
## free.lambda.mu 324.5329 0.0000000 0.28538150
## free.q         329.2331 4.7002695 0.02721293

From this we see that the preferred model is one in which λ & μ differ, and the character transitions are ordered (although just barely, though let’s ignore this for the time being):

Just as we learned for the birth-death model, the flexible diversitree package also allows us to set up & run a Bayesian MCMC with our model. Here, we will focus on doing this with our best-fitting model (although in general we would have probably chosen to either use Bayesian or Likelihood for all our analyses a priori).

Bayesian Analysis

This exercise illustrates how one conducts a Bayesian analysis in Diversitree. We won’t go over it in class but this is a good resource if you decide to try this kind of analysis yourself.

Because this is a Bayesian analysis, we start by setting a prior distribution for the variables in our model:

prior <- make.prior.exponential(1/2)
prior
## function (pars) 
## sum(dexp(pars, r, log = TRUE))
## <environment: 0x7f88b0c1d318>

The prior distribution is just a probability density function that will be used internally by mcmc.

The other thing we need to do, because of the type of MCMC used by diversitree, is to set a control parameter w. This parameter affects how many function evaluations are required between updates to the MCMC chain; however we need to do some preliminary sampling from the posterior distribution to set it:

prelim <- mcmc(lik.lambda.mu, fit.lambda.mu$par, nsteps=100, prior=prior, 
    w=1, print.every=0)
head(prelim)
##   i     lambda1   lambda2   lambda3          mu1       mu2       mu3
## 1 1 0.092789083 0.2090567 0.4418774 0.0039896788 0.1511720 0.4236700
## 2 2 0.099020496 0.2319716 0.4497358 0.0501746848 0.2247217 0.4551220
## 3 3 0.122357587 0.1740923 0.5025504 0.0194961663 0.1838763 0.4568868
## 4 4 0.071654259 0.2818159 0.4743468 0.0103490478 0.2372909 0.4148908
## 5 5 0.016866308 0.2640114 0.4230620 0.0004929538 0.1705742 0.4339026
## 6 6 0.009789571 0.2438943 0.4464046 0.0031362104 0.1268606 0.4386646
##          q12         p
## 1 0.05042854 -162.1751
## 2 0.06755596 -163.5020
## 3 0.05082352 -163.2312
## 4 0.05174632 -161.6295
## 5 0.04159170 -162.8161
## 6 0.04305355 -164.3945
## the following is how Fitzjohn recommends setting w
w <- diff(sapply(prelim[2:(ncol(prelim)-1)], quantile, c(0.05, 0.95)))
w
##       lambda1   lambda2  lambda3        mu1       mu2       mu3      q12
## 95% 0.1350635 0.3562925 0.460306 0.09688166 0.4778757 0.5525677 1.195777

Now we are ready to run our MCMC!

mcmc.lambda.mu <- mcmc(lik.lambda.mu, p[colnames(w)], nsteps=1000, 
    prior=prior, w=w, print.every=100)
## 100: {0.1030, 0.4697, 0.8083, 0.0072, 0.6692, 0.6317, 0.0659} -> -168.58348
## 200: {0.1014, 0.1858, 0.3184, 0.0205, 0.2350, 0.2197, 0.0820} -> -162.50068
## 300: {0.1855, 0.3457, 0.5126, 0.1264, 0.2190, 0.4471, 0.0474} -> -168.25737
## 400: {0.0557, 0.1663, 0.7498, 0.0012, 0.0488, 0.6797, 0.2092} -> -168.27216
## 500: {0.0520, 0.4634, 0.9053, 0.0065, 0.3980, 0.9118, 0.1490} -> -168.13784
## 600: {0.0510, 0.2222, 0.8157, 0.0613, 0.2532, 0.6821, 0.1695} -> -168.22404
## 700: {0.1202, 0.3162, 0.3324, 0.0747, 0.5001, 0.2615, 0.0299} -> -168.47286
## 800: {0.1324, 0.1890, 0.6572, 0.0346, 0.2645, 0.6475, 0.1263} -> -165.69212
## 900: {0.0539, 0.1777, 0.5375, 0.0178, 0.1544, 0.4851, 0.0706} -> -162.30841
## 1000: {0.1267, 0.5071, 0.3322, 0.0009, 0.6285, 0.2588, 0.1314} -> -165.43472
head(mcmc.lambda.mu)
##   i    lambda1   lambda2   lambda3        mu1         mu2       mu3
## 1 1 0.07776414 0.2137392 0.3029475 0.05238125 0.001800773 0.2081768
## 2 2 0.03575911 0.1285946 0.3064494 0.09618955 0.110016845 0.3031415
## 3 3 0.09820649 0.1685717 0.3642669 0.09897791 0.022640871 0.3457640
## 4 4 0.09283533 0.2160474 0.1367783 0.01305970 0.102383855 0.2782124
## 5 5 0.13558836 0.2080410 0.1707534 0.03632909 0.015177519 0.3478259
## 6 6 0.15310330 0.0773500 0.2499495 0.01429214 0.012949337 0.2522562
##         q12         p
## 1 0.2508332 -171.3811
## 2 0.2049265 -170.6124
## 3 0.2800079 -168.8064
## 4 0.3805234 -170.6065
## 5 0.5073986 -170.9392
## 6 0.3354549 -168.7375

We can start with a likelihood-profile, showing how the likelihood changed through the MCMC run:

plot(mcmc.lambda.mu$i, mcmc.lambda.mu$p, type="l", xlab="generation",
    ylab="log(L)")

Let’s cut out the first 20% of our samples as burn-in. Since we started in a reasonable area of parameter space, this might be overly conservative, but this will often not be the case:

mcmc.lambda.mu. <- mcmc.lambda.mu[601:3000,]

We can compute the mean of our posterior sample, and we can also compute things like the density or the 95% HPD interval, just as we would with any standard Bayesian MCMC.

colMeans(mcmc.lambda.mu)[2:ncol(mcmc.lambda.mu)]
##       lambda1       lambda2       lambda3           mu1           mu2 
##    0.10079659    0.34803264    0.55570550    0.04494997    0.38926039 
##           mu3           q12             p 
##    0.51157601    0.09442285 -166.41022325
colors<-setNames(c("yellow","red","green"),1:3)
profiles.plot(mcmc.lambda.mu[,grep("lambda",colnames(mcmc.lambda.mu))], 
    col.line=colors, las=1, legend.pos="topright")

profiles.plot(mcmc.lambda.mu[,grep("mu",colnames(mcmc.lambda.mu))], 
    col.line=colors, las=1, legend.pos="topright")

We can do something similar with net diversification, that is, the difference between λ & μ:

net.div<-mcmc.lambda.mu[,grep("lambda",colnames(mcmc.lambda.mu))]-
    mcmc.lambda.mu[,grep("mu",colnames(mcmc.lambda.mu))]
colnames(net.div)<-paste("lambda-mu(",1:3,")",sep="")
profiles.plot(net.div, 
    xlab="Net diversification rate", ylab="Probability density",
    legend.pos="topleft",col.line=setNames(colors,colnames(net.div)),
    lty=1)

Modified by Liam J. Revell from a script written by Ricardo Betancur. Last updated 20 December 2017.