<!DOCTYPE html>
<html lang="" xml:lang="">
<head>

  <meta charset="utf-8" />
  <meta http-equiv="X-UA-Compatible" content="IE=edge" />
  <title>8 Estimating birth-death models | CRABS: Congruent Rate Analyses in Birth-death Scenarios</title>
  <meta name="description" content="8 Estimating birth-death models | CRABS: Congruent Rate Analyses in Birth-death Scenarios" />
  <meta name="generator" content="bookdown 0.36 and GitBook 2.6.7" />

  <meta property="og:title" content="8 Estimating birth-death models | CRABS: Congruent Rate Analyses in Birth-death Scenarios" />
  <meta property="og:type" content="book" />
  
  
  

  <meta name="twitter:card" content="summary" />
  <meta name="twitter:title" content="8 Estimating birth-death models | CRABS: Congruent Rate Analyses in Birth-death Scenarios" />
  
  
  

<meta name="author" content="Sebastian Höhna" />
<meta name="author" content="Bjørn T. Kopperud" />
<meta name="author" content="Andrew F. Magee" />
<meta name="author" content="Jérémy Andréoletti" />


<meta name="date" content="2023-10-23" />

  <meta name="viewport" content="width=device-width, initial-scale=1" />
  <meta name="apple-mobile-web-app-capable" content="yes" />
  <meta name="apple-mobile-web-app-status-bar-style" content="black" />
  
  
<link rel="prev" href="exploring-the-posterior.html"/>
<link rel="next" href="references.html"/>
<script src="lib/jquery-3.6.0/jquery-3.6.0.min.js"></script>
<script src="https://cdn.jsdelivr.net/npm/fuse.js@6.4.6/dist/fuse.min.js"></script>
<link href="lib/gitbook-2.6.7/css/style.css" rel="stylesheet" />
<link href="lib/gitbook-2.6.7/css/plugin-table.css" rel="stylesheet" />
<link href="lib/gitbook-2.6.7/css/plugin-bookdown.css" rel="stylesheet" />
<link href="lib/gitbook-2.6.7/css/plugin-highlight.css" rel="stylesheet" />
<link href="lib/gitbook-2.6.7/css/plugin-search.css" rel="stylesheet" />
<link href="lib/gitbook-2.6.7/css/plugin-fontsettings.css" rel="stylesheet" />
<link href="lib/gitbook-2.6.7/css/plugin-clipboard.css" rel="stylesheet" />








<link href="lib/anchor-sections-1.1.0/anchor-sections.css" rel="stylesheet" />
<link href="lib/anchor-sections-1.1.0/anchor-sections-hash.css" rel="stylesheet" />
<script src="lib/anchor-sections-1.1.0/anchor-sections.js"></script>


<style type="text/css">
pre > code.sourceCode { white-space: pre; position: relative; }
pre > code.sourceCode > span { display: inline-block; line-height: 1.25; }
pre > code.sourceCode > span:empty { height: 1.2em; }
.sourceCode { overflow: visible; }
code.sourceCode > span { color: inherit; text-decoration: inherit; }
pre.sourceCode { margin: 0; }
@media screen {
div.sourceCode { overflow: auto; }
}
@media print {
pre > code.sourceCode { white-space: pre-wrap; }
pre > code.sourceCode > span { text-indent: -5em; padding-left: 5em; }
}
pre.numberSource code
  { counter-reset: source-line 0; }
pre.numberSource code > span
  { position: relative; left: -4em; counter-increment: source-line; }
pre.numberSource code > span > a:first-child::before
  { content: counter(source-line);
    position: relative; left: -1em; text-align: right; vertical-align: baseline;
    border: none; display: inline-block;
    -webkit-touch-callout: none; -webkit-user-select: none;
    -khtml-user-select: none; -moz-user-select: none;
    -ms-user-select: none; user-select: none;
    padding: 0 4px; width: 4em;
    color: #aaaaaa;
  }
pre.numberSource { margin-left: 3em; border-left: 1px solid #aaaaaa;  padding-left: 4px; }
div.sourceCode
  {   }
@media screen {
pre > code.sourceCode > span > a:first-child::before { text-decoration: underline; }
}
code span.al { color: #ff0000; font-weight: bold; } /* Alert */
code span.an { color: #60a0b0; font-weight: bold; font-style: italic; } /* Annotation */
code span.at { color: #7d9029; } /* Attribute */
code span.bn { color: #40a070; } /* BaseN */
code span.bu { } /* BuiltIn */
code span.cf { color: #007020; font-weight: bold; } /* ControlFlow */
code span.ch { color: #4070a0; } /* Char */
code span.cn { color: #880000; } /* Constant */
code span.co { color: #60a0b0; font-style: italic; } /* Comment */
code span.cv { color: #60a0b0; font-weight: bold; font-style: italic; } /* CommentVar */
code span.do { color: #ba2121; font-style: italic; } /* Documentation */
code span.dt { color: #902000; } /* DataType */
code span.dv { color: #40a070; } /* DecVal */
code span.er { color: #ff0000; font-weight: bold; } /* Error */
code span.ex { } /* Extension */
code span.fl { color: #40a070; } /* Float */
code span.fu { color: #06287e; } /* Function */
code span.im { } /* Import */
code span.in { color: #60a0b0; font-weight: bold; font-style: italic; } /* Information */
code span.kw { color: #007020; font-weight: bold; } /* Keyword */
code span.op { color: #666666; } /* Operator */
code span.ot { color: #007020; } /* Other */
code span.pp { color: #bc7a00; } /* Preprocessor */
code span.sc { color: #4070a0; } /* SpecialChar */
code span.ss { color: #bb6688; } /* SpecialString */
code span.st { color: #4070a0; } /* String */
code span.va { color: #19177c; } /* Variable */
code span.vs { color: #4070a0; } /* VerbatimString */
code span.wa { color: #60a0b0; font-weight: bold; font-style: italic; } /* Warning */
</style>

<style type="text/css">
  
  div.hanging-indent{margin-left: 1.5em; text-indent: -1.5em;}
</style>
<style type="text/css">
/* Used with Pandoc 2.11+ new --citeproc when CSL is used */
div.csl-bib-body { }
div.csl-entry {
  clear: both;
}
.hanging div.csl-entry {
  margin-left:2em;
  text-indent:-2em;
}
div.csl-left-margin {
  min-width:2em;
  float:left;
}
div.csl-right-inline {
  margin-left:2em;
  padding-left:1em;
}
div.csl-indent {
  margin-left: 2em;
}
</style>

</head>

<body>



  <div class="book without-animation with-summary font-size-2 font-family-1" data-basepath=".">

    <div class="book-summary">
      <nav role="navigation">

<ul class="summary">
<li class="chapter" data-level="1" data-path="overview.html"><a href="overview.html"><i class="fa fa-check"></i><b>1</b> Overview</a></li>
<li class="chapter" data-level="2" data-path="installation.html"><a href="installation.html"><i class="fa fa-check"></i><b>2</b> Installation</a></li>
<li class="chapter" data-level="3" data-path="features.html"><a href="features.html"><i class="fa fa-check"></i><b>3</b> Features</a>
<ul>
<li class="chapter" data-level="3.1" data-path="features.html"><a href="features.html#overview-of-functions-in-crabs"><i class="fa fa-check"></i><b>3.1</b> Overview of functions in <code>CRABS</code></a></li>
<li class="chapter" data-level="3.2" data-path="features.html"><a href="features.html#samplebasic"><i class="fa fa-check"></i><b>3.2</b> <code>sample.basic.models</code></a></li>
<li class="chapter" data-level="3.3" data-path="features.html"><a href="features.html#samplejoint"><i class="fa fa-check"></i><b>3.3</b> <code>sample.basic.models.joint</code></a></li>
</ul></li>
<li class="chapter" data-level="4" data-path="creating-the-congruence-class.html"><a href="creating-the-congruence-class.html"><i class="fa fa-check"></i><b>4</b> Creating the congruence class</a>
<ul>
<li class="chapter" data-level="4.1" data-path="creating-the-congruence-class.html"><a href="creating-the-congruence-class.html#rate-transformations"><i class="fa fa-check"></i><b>4.1</b> Rate transformations</a></li>
<li class="chapter" data-level="4.2" data-path="creating-the-congruence-class.html"><a href="creating-the-congruence-class.html#slopes-and-delta-t"><i class="fa fa-check"></i><b>4.2</b> Slopes and <span class="math inline">\(\Delta t\)</span></a></li>
<li class="chapter" data-level="4.3" data-path="creating-the-congruence-class.html"><a href="creating-the-congruence-class.html#the-model-object"><i class="fa fa-check"></i><b>4.3</b> The model object</a></li>
</ul></li>
<li class="chapter" data-level="5" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html"><i class="fa fa-check"></i><b>5</b> Testing alternative hypotheses</a>
<ul>
<li class="chapter" data-level="5.1" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#alt-ext"><i class="fa fa-check"></i><b>5.1</b> Testing alternative extinction rates</a>
<ul>
<li class="chapter" data-level="5.1.1" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#constant-extinction-rates"><i class="fa fa-check"></i><b>5.1.1</b> Constant extinction rates</a></li>
<li class="chapter" data-level="5.1.2" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#linearly-increasing-extinction-rates"><i class="fa fa-check"></i><b>5.1.2</b> Linearly increasing extinction rates</a></li>
<li class="chapter" data-level="5.1.3" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#exponential-increasing-extinction-rates"><i class="fa fa-check"></i><b>5.1.3</b> Exponential increasing extinction rates</a></li>
<li class="chapter" data-level="5.1.4" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#two-epoch-extinction-rates"><i class="fa fa-check"></i><b>5.1.4</b> Two-epoch extinction rates</a></li>
<li class="chapter" data-level="5.1.5" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#reversed-trend-extinction-rates"><i class="fa fa-check"></i><b>5.1.5</b> Reversed trend extinction rates</a></li>
</ul></li>
<li class="chapter" data-level="5.2" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#testing-alternative-speciation-rates"><i class="fa fa-check"></i><b>5.2</b> Testing alternative speciation rates</a>
<ul>
<li class="chapter" data-level="5.2.1" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#constant-speciation-rates"><i class="fa fa-check"></i><b>5.2.1</b> Constant speciation rates</a></li>
<li class="chapter" data-level="5.2.2" data-path="testing-alternative-hypotheses.html"><a href="testing-alternative-hypotheses.html#exponential-increasingdecreasing-speciation-rates"><i class="fa fa-check"></i><b>5.2.2</b> Exponential increasing/decreasing speciation rates</a></li>
</ul></li>
</ul></li>
<li class="chapter" data-level="6" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html"><i class="fa fa-check"></i><b>6</b> Exploring the congruence class</a>
<ul>
<li class="chapter" data-level="6.1" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#assumptions-on-rates"><i class="fa fa-check"></i><b>6.1</b> Assumptions on rates</a></li>
<li class="chapter" data-level="6.2" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#independent-or-joint-rate-sampling"><i class="fa fa-check"></i><b>6.2</b> Independent or joint rate sampling</a></li>
<li class="chapter" data-level="6.3" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#sample-alternative-extinction-rates"><i class="fa fa-check"></i><b>6.3</b> Sample alternative extinction rates</a>
<ul>
<li class="chapter" data-level="6.3.1" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#temporally-autocorrelated"><i class="fa fa-check"></i><b>6.3.1</b> Temporally autocorrelated</a></li>
<li class="chapter" data-level="6.3.2" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#exponentially-decreasing"><i class="fa fa-check"></i><b>6.3.2</b> Exponentially decreasing</a></li>
<li class="chapter" data-level="6.3.3" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#linearly-increasing"><i class="fa fa-check"></i><b>6.3.3</b> Linearly increasing</a></li>
<li class="chapter" data-level="6.3.4" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#episodic"><i class="fa fa-check"></i><b>6.3.4</b> Episodic</a></li>
<li class="chapter" data-level="6.3.5" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#uniform-independent"><i class="fa fa-check"></i><b>6.3.5</b> Uniform independent</a></li>
<li class="chapter" data-level="6.3.6" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#brownian-motion-with-a-trend"><i class="fa fa-check"></i><b>6.3.6</b> Brownian motion with a trend</a></li>
</ul></li>
<li class="chapter" data-level="6.4" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#sample-alternative-speciation-rates"><i class="fa fa-check"></i><b>6.4</b> Sample alternative speciation rates</a>
<ul>
<li class="chapter" data-level="6.4.1" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#temporally-autocorrelated-1"><i class="fa fa-check"></i><b>6.4.1</b> Temporally autocorrelated</a></li>
</ul></li>
<li class="chapter" data-level="6.5" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#jointly-sample-alternative-speciation-and-extinction-rates"><i class="fa fa-check"></i><b>6.5</b> Jointly sample alternative speciation and extinction rates</a>
<ul>
<li class="chapter" data-level="6.5.1" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#constant-joint-trajectory"><i class="fa fa-check"></i><b>6.5.1</b> Constant joint trajectory</a></li>
<li class="chapter" data-level="6.5.2" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#hsmrf-joint-trajectory"><i class="fa fa-check"></i><b>6.5.2</b> HSMRF joint trajectory</a></li>
<li class="chapter" data-level="6.5.3" data-path="exploring-the-congruence-class.html"><a href="exploring-the-congruence-class.html#hsmrf-joint-trajectory-with-global-regularity-filter"><i class="fa fa-check"></i><b>6.5.3</b> HSMRF joint trajectory with global regularity filter</a></li>
</ul></li>
</ul></li>
<li class="chapter" data-level="7" data-path="exploring-the-posterior.html"><a href="exploring-the-posterior.html"><i class="fa fa-check"></i><b>7</b> Exploring the posterior</a>
<ul>
<li class="chapter" data-level="7.1" data-path="exploring-the-posterior.html"><a href="exploring-the-posterior.html#reading-the-posterior"><i class="fa fa-check"></i><b>7.1</b> Reading the posterior</a></li>
<li class="chapter" data-level="7.2" data-path="exploring-the-posterior.html"><a href="exploring-the-posterior.html#sampling-congruent-models-in-the-posterior"><i class="fa fa-check"></i><b>7.2</b> Sampling congruent models in the posterior</a></li>
</ul></li>
<li class="chapter" data-level="8" data-path="estimating-birth-death-models.html"><a href="estimating-birth-death-models.html"><i class="fa fa-check"></i><b>8</b> Estimating birth-death models</a>
<ul>
<li class="chapter" data-level="8.1" data-path="estimating-birth-death-models.html"><a href="estimating-birth-death-models.html#reading-the-phylogenetic-tree"><i class="fa fa-check"></i><b>8.1</b> Reading the phylogenetic tree</a></li>
<li class="chapter" data-level="8.2" data-path="estimating-birth-death-models.html"><a href="estimating-birth-death-models.html#tess"><i class="fa fa-check"></i><b>8.2</b> TESS Analysis</a></li>
<li class="chapter" data-level="8.3" data-path="estimating-birth-death-models.html"><a href="estimating-birth-death-models.html#treepar"><i class="fa fa-check"></i><b>8.3</b> TreePar</a></li>
<li class="chapter" data-level="8.4" data-path="estimating-birth-death-models.html"><a href="estimating-birth-death-models.html#revbayes"><i class="fa fa-check"></i><b>8.4</b> The HSMRF model in RevBayes</a></li>
</ul></li>
<li class="chapter" data-level="9" data-path="references.html"><a href="references.html"><i class="fa fa-check"></i><b>9</b> References</a></li>
</ul>

      </nav>
    </div>

    <div class="book-body">
      <div class="body-inner">
        <div class="book-header" role="navigation">
          <h1>
            <i class="fa fa-circle-o-notch fa-spin"></i><a href="./">CRABS: Congruent Rate Analyses in Birth-death Scenarios</a>
          </h1>
        </div>

        <div class="page-wrapper" tabindex="-1" role="main">
          <div class="page-inner">

            <section class="normal" id="section-">
<div id="estimating-birth-death-models" class="section level1 hasAnchor" number="8">
<h1><span class="header-section-number">8</span> Estimating birth-death models<a href="estimating-birth-death-models.html#estimating-birth-death-models" class="anchor-section" aria-label="Anchor link to header"></a></h1>
<p>One of the main features and aims of <code>CRABS</code> is to explore the congruence class for the diversification rates that you estimated using another software or R package.
Here we provide you some examples using standard software for estimating diversification rates through time, but if your preferred software is not included then you can either specify the diversification rate functions manually, or try to follow our examples in loading the output.</p>
<div id="reading-the-phylogenetic-tree" class="section level2 hasAnchor" number="8.1">
<h2><span class="header-section-number">8.1</span> Reading the phylogenetic tree<a href="estimating-birth-death-models.html#reading-the-phylogenetic-tree" class="anchor-section" aria-label="Anchor link to header"></a></h2>
<p>First, we assume that you have a phylogenetic tree which you can load into R.
There are several tutorials on how to do so (e.g., <a href="https://cran.r-project.org/web/packages/TESS/vignettes/Bayesian_Diversification_Rate_Analysis.pdf">TESS</a> or <a href="https://cran.r-project.org/web/packages/ape/ape.pdf">APE</a>, and we simply refer to these tutorials here).
As an example, we used the <em>Primates</em> phylogeny from <span class="citation">Springer et al. (2012)</span>.
You can download the data from our <a href="https://revbayes.github.io/tutorials/divrate/ebd.html">RevBayes tutorial</a>.
First we need to load the APE package <span class="citation">(Paradis &amp; Schliep, 2019)</span>.</p>
<div class="sourceCode" id="cb75"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb75-1"><a href="estimating-birth-death-models.html#cb75-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(ape)</span></code></pre></div>
<p>Now we can read our phylogenetic tree</p>
<div class="sourceCode" id="cb76"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb76-1"><a href="estimating-birth-death-models.html#cb76-1" aria-hidden="true" tabindex="-1"></a>tree <span class="ot">&lt;-</span> <span class="fu">read.tree</span>(<span class="at">file=</span><span class="st">&quot;data/primates.tre&quot;</span>)</span></code></pre></div>
<p>Next, read the software examples to fit birth-death models to the tree, and estimate diversification rates.</p>
</div>
<div id="tess" class="section level2 hasAnchor" number="8.2">
<h2><span class="header-section-number">8.2</span> TESS Analysis<a href="estimating-birth-death-models.html#tess" class="anchor-section" aria-label="Anchor link to header"></a></h2>
<p>First, load the TESS package.</p>
<div class="sourceCode" id="cb77"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb77-1"><a href="estimating-birth-death-models.html#cb77-1" aria-hidden="true" tabindex="-1"></a><span class="fu">library</span>(TESS)</span></code></pre></div>
<p>For some of our plotting, we need the speciation times.
So we extract these from the tree.</p>
<div class="sourceCode" id="cb78"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb78-1"><a href="estimating-birth-death-models.html#cb78-1" aria-hidden="true" tabindex="-1"></a>btimes <span class="ot">&lt;-</span> <span class="fu">sort</span>( <span class="fu">as.numeric</span>( <span class="fu">branching.times</span>( tree ) ) )</span>
<span id="cb78-2"><a href="estimating-birth-death-models.html#cb78-2" aria-hidden="true" tabindex="-1"></a>max_t <span class="ot">&lt;-</span> <span class="fu">max</span>(btimes)</span></code></pre></div>
<div class="sourceCode" id="cb79"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb79-1"><a href="estimating-birth-death-models.html#cb79-1" aria-hidden="true" tabindex="-1"></a><span class="fu">tess.analysis</span>(tree,</span>
<span id="cb79-2"><a href="estimating-birth-death-models.html#cb79-2" aria-hidden="true" tabindex="-1"></a>              <span class="at">empiricalHyperPriors =</span> <span class="cn">TRUE</span>,</span>
<span id="cb79-3"><a href="estimating-birth-death-models.html#cb79-3" aria-hidden="true" tabindex="-1"></a>              <span class="at">samplingProbability =</span> <span class="fl">1.0</span>,</span>
<span id="cb79-4"><a href="estimating-birth-death-models.html#cb79-4" aria-hidden="true" tabindex="-1"></a>              <span class="at">numExpectedRateChanges =</span> <span class="dv">2</span>,</span>
<span id="cb79-5"><a href="estimating-birth-death-models.html#cb79-5" aria-hidden="true" tabindex="-1"></a>              <span class="at">MAX_ITERATIONS =</span> <span class="dv">100000</span>,</span>
<span id="cb79-6"><a href="estimating-birth-death-models.html#cb79-6" aria-hidden="true" tabindex="-1"></a>              <span class="at">MAX_TIME =</span> <span class="dv">2</span><span class="sc">*</span><span class="dv">60</span><span class="sc">*</span><span class="dv">60</span>,</span>
<span id="cb79-7"><a href="estimating-birth-death-models.html#cb79-7" aria-hidden="true" tabindex="-1"></a>              <span class="at">MIN_ESS =</span> <span class="dv">625</span>,</span>
<span id="cb79-8"><a href="estimating-birth-death-models.html#cb79-8" aria-hidden="true" tabindex="-1"></a>              <span class="at">dir =</span> <span class="st">&quot;primates_comet&quot;</span>)</span></code></pre></div>
<pre><code>## Estimating empirical hyper-parameters.
## 
## Burning-in the chain ...
## 0--------25--------50--------75--------100
## ==========================================
## Finished burnin period!
## 
## Running the chain ...
## 0--------25--------50--------75--------100
## ==========================================
## Finished MCMC!
## 
## Parameter | delta | Acceptance Probability
## ==========================================
## diversification      | 0.152     | 0.487
## turnover     | 0.031     | 0.406
## 
## Performing CoMET analysis.
## 
## Burning-in the chain ...
## 0--------25--------50--------75--------100
## ==========================================
## 
## Running the chain ... 
## 0--------25--------50--------75--------100
## ==========================================</code></pre>
<div class="sourceCode" id="cb81"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb81-1"><a href="estimating-birth-death-models.html#cb81-1" aria-hidden="true" tabindex="-1"></a>tess.output <span class="ot">&lt;-</span> <span class="fu">tess.process.output</span>(<span class="st">&quot;primates_comet&quot;</span>,</span>
<span id="cb81-2"><a href="estimating-birth-death-models.html#cb81-2" aria-hidden="true" tabindex="-1"></a>                              <span class="at">numExpectedRateChanges =</span> <span class="dv">2</span>)</span></code></pre></div>
<div class="sourceCode" id="cb82"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb82-1"><a href="estimating-birth-death-models.html#cb82-1" aria-hidden="true" tabindex="-1"></a>est_speciation <span class="ot">&lt;-</span> <span class="fu">apply</span>(tess.output[[<span class="st">&quot;speciation rates&quot;</span>]], <span class="dv">2</span>, median)</span>
<span id="cb82-2"><a href="estimating-birth-death-models.html#cb82-2" aria-hidden="true" tabindex="-1"></a>est_extinction <span class="ot">&lt;-</span> <span class="fu">apply</span>(tess.output[[<span class="st">&quot;extinction rates&quot;</span>]], <span class="dv">2</span>, median)</span></code></pre></div>
<p>Let us make a data frame representation of the rates. Since TESS is parameterized in forward time (i.e. with <span class="math inline">\(t = 0\)</span> being the root of the phylogeny, increasing toward the present), we want to reverse the rates to be compatible with <code>CRABS</code>.</p>
<div class="sourceCode" id="cb83"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb83-1"><a href="estimating-birth-death-models.html#cb83-1" aria-hidden="true" tabindex="-1"></a>times <span class="ot">&lt;-</span> <span class="fu">seq</span>(<span class="dv">0</span>, max_t, <span class="at">length.out =</span> <span class="dv">100</span>)</span>
<span id="cb83-2"><a href="estimating-birth-death-models.html#cb83-2" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb83-3"><a href="estimating-birth-death-models.html#cb83-3" aria-hidden="true" tabindex="-1"></a>primates_ebd_tess <span class="ot">&lt;-</span> <span class="fu">tibble</span>(</span>
<span id="cb83-4"><a href="estimating-birth-death-models.html#cb83-4" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;time&quot;</span> <span class="ot">=</span> times,</span>
<span id="cb83-5"><a href="estimating-birth-death-models.html#cb83-5" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;lambda&quot;</span> <span class="ot">=</span> <span class="fu">rev</span>(est_speciation),</span>
<span id="cb83-6"><a href="estimating-birth-death-models.html#cb83-6" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;mu&quot;</span> <span class="ot">=</span> <span class="fu">rev</span>(est_extinction),</span>
<span id="cb83-7"><a href="estimating-birth-death-models.html#cb83-7" aria-hidden="true" tabindex="-1"></a>)</span>
<span id="cb83-8"><a href="estimating-birth-death-models.html#cb83-8" aria-hidden="true" tabindex="-1"></a><span class="fu">head</span>(primates_ebd_tess)</span></code></pre></div>
<pre><code>## # A tibble: 6 × 3
##    time lambda    mu
##   &lt;dbl&gt;  &lt;dbl&gt; &lt;dbl&gt;
## 1 0      0.194 0.337
## 2 0.721  0.592 0.337
## 3 1.44   0.593 0.355
## 4 2.16   0.592 0.421
## 5 2.88   0.586 0.477
## 6 3.60   0.569 0.491</code></pre>
<p>Now let’s have a look at the estimated rate functions</p>
<div class="sourceCode" id="cb85"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb85-1"><a href="estimating-birth-death-models.html#cb85-1" aria-hidden="true" tabindex="-1"></a>p <span class="ot">&lt;-</span> primates_ebd_tess <span class="sc">%&gt;%</span></span>
<span id="cb85-2"><a href="estimating-birth-death-models.html#cb85-2" aria-hidden="true" tabindex="-1"></a>  tidyr<span class="sc">::</span><span class="fu">gather</span>(<span class="st">&quot;rate&quot;</span>, <span class="st">&quot;value&quot;</span>, <span class="sc">-</span>time) <span class="sc">%&gt;%</span></span>
<span id="cb85-3"><a href="estimating-birth-death-models.html#cb85-3" aria-hidden="true" tabindex="-1"></a>  <span class="fu">ggplot</span>(<span class="fu">aes</span>(<span class="at">x =</span> time, <span class="at">y =</span> value, <span class="at">color =</span> rate)) <span class="sc">+</span></span>
<span id="cb85-4"><a href="estimating-birth-death-models.html#cb85-4" aria-hidden="true" tabindex="-1"></a>  <span class="fu">geom_line</span>() <span class="sc">+</span></span>
<span id="cb85-5"><a href="estimating-birth-death-models.html#cb85-5" aria-hidden="true" tabindex="-1"></a>  <span class="fu">theme_classic</span>() <span class="sc">+</span></span>
<span id="cb85-6"><a href="estimating-birth-death-models.html#cb85-6" aria-hidden="true" tabindex="-1"></a>  <span class="fu">scale_x_reverse</span>() <span class="sc">+</span></span>
<span id="cb85-7"><a href="estimating-birth-death-models.html#cb85-7" aria-hidden="true" tabindex="-1"></a>  <span class="fu">theme</span>(<span class="at">legend.position =</span> <span class="fu">c</span>(<span class="fl">0.2</span>, <span class="fl">0.2</span>)) <span class="sc">+</span></span>
<span id="cb85-8"><a href="estimating-birth-death-models.html#cb85-8" aria-hidden="true" tabindex="-1"></a>  <span class="fu">labs</span>(<span class="at">y =</span> <span class="st">&quot;rate&quot;</span>, <span class="at">x =</span> <span class="st">&quot;time before present&quot;</span>)</span>
<span id="cb85-9"><a href="estimating-birth-death-models.html#cb85-9" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>(p)</span></code></pre></div>
<div class="figure" style="text-align: center"><span style="display:block;" id="fig:plot-tess-rates"></span>
<img src="figures/plot-tess-rates-1.png" alt="The estimated posterior median results from an episodic birth-death model fitted in RevBayes." width="480" />
<p class="caption">
Figure 8.1: The estimated posterior median results from an episodic birth-death model fitted in RevBayes.
</p>
</div>
<p>The data frame <code>primates_ebd_tess</code> is what we get when we call the <code>data("primates_ebd_tess")</code> function in <code>CRABS</code>. We use the RevBayes estimated rates in our vignette examples, but it is easy to substitute those for the TESS analysis results to construct the congruence class.</p>
</div>
<div id="treepar" class="section level2 hasAnchor" number="8.3">
<h2><span class="header-section-number">8.3</span> TreePar<a href="estimating-birth-death-models.html#treepar" class="anchor-section" aria-label="Anchor link to header"></a></h2>
<p>Let us walk here through a very simple TreePar analysis <span class="citation">(Stadler, 2011)</span>.
For a real publication worthy analysis, you may want to do a more thorough analysis by allowing more possible shifts and a finer grid for the search space.</p>
<p>First, load the TreePar package.</p>
<div class="sourceCode" id="cb86"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb86-1"><a href="estimating-birth-death-models.html#cb86-1" aria-hidden="true" tabindex="-1"></a>  <span class="fu">library</span>(TreePar)</span></code></pre></div>
<p>TreePar needs the ages of the phylogeny to perform it’s likelihood computation.
We use the TreePar builtin functionality for this.</p>
<div class="sourceCode" id="cb87"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb87-1"><a href="estimating-birth-death-models.html#cb87-1" aria-hidden="true" tabindex="-1"></a>x     <span class="ot">&lt;-</span> <span class="fu">sort</span>(<span class="fu">getx</span>(tree),<span class="at">decreasing=</span><span class="cn">TRUE</span>)</span>
<span id="cb87-2"><a href="estimating-birth-death-models.html#cb87-2" aria-hidden="true" tabindex="-1"></a>max_t <span class="ot">&lt;-</span> <span class="fu">max</span>(x)</span></code></pre></div>
<p>For our toy analysis, we setup some simpler than normal search values.
We specify that there are at most three rate shifts, and that there are ten possible shift points.
For your own analysis, you should most likely increase these values.
We also specify an interval offset of 0.05, which means that in the most recent 5% of the tree age there was no rate shift event.
This is simply for improving the parameter estimation.</p>
<div class="sourceCode" id="cb88"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb88-1"><a href="estimating-birth-death-models.html#cb88-1" aria-hidden="true" tabindex="-1"></a>MAX.NUM.SHIFTS  <span class="ot">=</span> <span class="dv">4</span></span>
<span id="cb88-2"><a href="estimating-birth-death-models.html#cb88-2" aria-hidden="true" tabindex="-1"></a>NUM.INTERVALS   <span class="ot">=</span> <span class="dv">20</span></span>
<span id="cb88-3"><a href="estimating-birth-death-models.html#cb88-3" aria-hidden="true" tabindex="-1"></a>INTERVAL.OFFSET <span class="ot">=</span> <span class="fl">0.05</span></span></code></pre></div>
<p>Now we run the full TreePar maximum likelihood search.</p>
<div class="sourceCode" id="cb89"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb89-1"><a href="estimating-birth-death-models.html#cb89-1" aria-hidden="true" tabindex="-1"></a>res   <span class="ot">&lt;-</span> <span class="fu">bd.shifts.optim</span>(x,</span>
<span id="cb89-2"><a href="estimating-birth-death-models.html#cb89-2" aria-hidden="true" tabindex="-1"></a>                         <span class="at">sampling=</span><span class="fu">rep</span>(<span class="dv">1</span>,MAX.NUM.SHIFTS<span class="sc">+</span><span class="dv">1</span>),</span>
<span id="cb89-3"><a href="estimating-birth-death-models.html#cb89-3" aria-hidden="true" tabindex="-1"></a>                         <span class="at">grid=</span>max_t<span class="sc">/</span>NUM.INTERVALS,</span>
<span id="cb89-4"><a href="estimating-birth-death-models.html#cb89-4" aria-hidden="true" tabindex="-1"></a>                         <span class="at">start=</span>INTERVAL.OFFSET<span class="sc">*</span>max_t,</span>
<span id="cb89-5"><a href="estimating-birth-death-models.html#cb89-5" aria-hidden="true" tabindex="-1"></a>                         <span class="at">end=</span>max_t,</span>
<span id="cb89-6"><a href="estimating-birth-death-models.html#cb89-6" aria-hidden="true" tabindex="-1"></a>                         <span class="at">survival=</span><span class="dv">1</span>,</span>
<span id="cb89-7"><a href="estimating-birth-death-models.html#cb89-7" aria-hidden="true" tabindex="-1"></a>                         <span class="at">posdiv=</span><span class="cn">TRUE</span>)[[<span class="dv">2</span>]]</span></code></pre></div>
<p>Now we have the maximum likelihood estimates for the different number of shifts.
We want to test if <span class="math inline">\(i\)</span> shifts explain the tree significantly better than <span class="math inline">\(i-1\)</span> shifts.</p>
<div class="sourceCode" id="cb90"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb90-1"><a href="estimating-birth-death-models.html#cb90-1" aria-hidden="true" tabindex="-1"></a>best <span class="ot">&lt;-</span> <span class="dv">0</span></span>
<span id="cb90-2"><a href="estimating-birth-death-models.html#cb90-2" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> (i <span class="cf">in</span> <span class="dv">1</span><span class="sc">:</span>(<span class="fu">length</span>(res)<span class="sc">-</span><span class="dv">1</span>)) {</span>
<span id="cb90-3"><a href="estimating-birth-death-models.html#cb90-3" aria-hidden="true" tabindex="-1"></a>    test <span class="ot">&lt;-</span> <span class="fu">pchisq</span>(<span class="dv">2</span><span class="sc">*</span>(res[[i]][<span class="dv">1</span>]<span class="sc">-</span>res[[i<span class="sc">+</span><span class="dv">1</span>]][<span class="dv">1</span>]),<span class="dv">3</span>)</span>
<span id="cb90-4"><a href="estimating-birth-death-models.html#cb90-4" aria-hidden="true" tabindex="-1"></a>    <span class="co">#if test&gt;0.95 then i shifts is significantly</span></span>
<span id="cb90-5"><a href="estimating-birth-death-models.html#cb90-5" aria-hidden="true" tabindex="-1"></a>    <span class="co"># better than i-1 shifts at a 5% error</span></span>
<span id="cb90-6"><a href="estimating-birth-death-models.html#cb90-6" aria-hidden="true" tabindex="-1"></a>    <span class="cf">if</span> ( test <span class="sc">&gt;</span> <span class="fl">0.95</span> ) {</span>
<span id="cb90-7"><a href="estimating-birth-death-models.html#cb90-7" aria-hidden="true" tabindex="-1"></a>        best <span class="ot">&lt;-</span> i</span>
<span id="cb90-8"><a href="estimating-birth-death-models.html#cb90-8" aria-hidden="true" tabindex="-1"></a>    } <span class="cf">else</span> {</span>
<span id="cb90-9"><a href="estimating-birth-death-models.html#cb90-9" aria-hidden="true" tabindex="-1"></a>        <span class="cf">break</span></span>
<span id="cb90-10"><a href="estimating-birth-death-models.html#cb90-10" aria-hidden="true" tabindex="-1"></a>    }</span>
<span id="cb90-11"><a href="estimating-birth-death-models.html#cb90-11" aria-hidden="true" tabindex="-1"></a>}</span></code></pre></div>
<div class="sourceCode" id="cb91"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb91-1"><a href="estimating-birth-death-models.html#cb91-1" aria-hidden="true" tabindex="-1"></a>best_result <span class="ot">&lt;-</span> res[[best<span class="sc">+</span><span class="dv">1</span>]]</span>
<span id="cb91-2"><a href="estimating-birth-death-models.html#cb91-2" aria-hidden="true" tabindex="-1"></a>speciation <span class="ot">&lt;-</span> <span class="fu">c</span>()</span>
<span id="cb91-3"><a href="estimating-birth-death-models.html#cb91-3" aria-hidden="true" tabindex="-1"></a>extinction <span class="ot">&lt;-</span> <span class="fu">c</span>()</span>
<span id="cb91-4"><a href="estimating-birth-death-models.html#cb91-4" aria-hidden="true" tabindex="-1"></a><span class="cf">for</span> (i <span class="cf">in</span> <span class="dv">1</span><span class="sc">:</span>(best<span class="sc">+</span><span class="dv">1</span>)) {</span>
<span id="cb91-5"><a href="estimating-birth-death-models.html#cb91-5" aria-hidden="true" tabindex="-1"></a>    speciation[i] <span class="ot">&lt;-</span> best_result[<span class="dv">1</span><span class="sc">+</span>i<span class="sc">+</span>(best<span class="sc">+</span><span class="dv">1</span>)] <span class="sc">/</span></span>
<span id="cb91-6"><a href="estimating-birth-death-models.html#cb91-6" aria-hidden="true" tabindex="-1"></a>                     (<span class="dv">1</span><span class="sc">-</span>best_result[<span class="dv">1</span><span class="sc">+</span>i])</span>
<span id="cb91-7"><a href="estimating-birth-death-models.html#cb91-7" aria-hidden="true" tabindex="-1"></a>    extinction[i] <span class="ot">&lt;-</span> best_result[<span class="dv">1</span><span class="sc">+</span>i<span class="sc">+</span>(best<span class="sc">+</span><span class="dv">1</span>)] <span class="sc">/</span></span>
<span id="cb91-8"><a href="estimating-birth-death-models.html#cb91-8" aria-hidden="true" tabindex="-1"></a>                     (<span class="dv">1</span><span class="sc">/</span>best_result[<span class="dv">1</span><span class="sc">+</span>i]<span class="sc">-</span><span class="dv">1</span>)</span>
<span id="cb91-9"><a href="estimating-birth-death-models.html#cb91-9" aria-hidden="true" tabindex="-1"></a>}</span>
<span id="cb91-10"><a href="estimating-birth-death-models.html#cb91-10" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb91-11"><a href="estimating-birth-death-models.html#cb91-11" aria-hidden="true" tabindex="-1"></a><span class="cf">if</span> ( best <span class="sc">&gt;</span> <span class="dv">0</span> ) {</span>
<span id="cb91-12"><a href="estimating-birth-death-models.html#cb91-12" aria-hidden="true" tabindex="-1"></a>    shift_times <span class="ot">&lt;-</span> best_result[ (<span class="fu">length</span>(best_result)<span class="sc">-</span>best<span class="sc">+</span><span class="dv">1</span>) <span class="sc">:</span></span>
<span id="cb91-13"><a href="estimating-birth-death-models.html#cb91-13" aria-hidden="true" tabindex="-1"></a>                                <span class="fu">length</span>(best_result)]</span>
<span id="cb91-14"><a href="estimating-birth-death-models.html#cb91-14" aria-hidden="true" tabindex="-1"></a>    tmp_lambda <span class="ot">&lt;-</span> <span class="cf">function</span>(x) {</span>
<span id="cb91-15"><a href="estimating-birth-death-models.html#cb91-15" aria-hidden="true" tabindex="-1"></a>       index <span class="ot">&lt;-</span> <span class="fu">findInterval</span>(x, shift_times);</span>
<span id="cb91-16"><a href="estimating-birth-death-models.html#cb91-16" aria-hidden="true" tabindex="-1"></a>       <span class="fu">return</span> ( speciation[index<span class="sc">+</span><span class="dv">1</span>] )</span>
<span id="cb91-17"><a href="estimating-birth-death-models.html#cb91-17" aria-hidden="true" tabindex="-1"></a>    }</span>
<span id="cb91-18"><a href="estimating-birth-death-models.html#cb91-18" aria-hidden="true" tabindex="-1"></a>    tmp_mu     <span class="ot">&lt;-</span> <span class="cf">function</span>(x) {</span>
<span id="cb91-19"><a href="estimating-birth-death-models.html#cb91-19" aria-hidden="true" tabindex="-1"></a>       index <span class="ot">&lt;-</span> <span class="fu">findInterval</span>(x, shift_times);</span>
<span id="cb91-20"><a href="estimating-birth-death-models.html#cb91-20" aria-hidden="true" tabindex="-1"></a>       <span class="fu">return</span> ( extinction[index<span class="sc">+</span><span class="dv">1</span>] )</span>
<span id="cb91-21"><a href="estimating-birth-death-models.html#cb91-21" aria-hidden="true" tabindex="-1"></a>    }</span>
<span id="cb91-22"><a href="estimating-birth-death-models.html#cb91-22" aria-hidden="true" tabindex="-1"></a>    lambda <span class="ot">&lt;-</span> <span class="fu">Vectorize</span>(tmp_lambda)</span>
<span id="cb91-23"><a href="estimating-birth-death-models.html#cb91-23" aria-hidden="true" tabindex="-1"></a>    mu     <span class="ot">&lt;-</span> <span class="fu">Vectorize</span>(tmp_mu)</span>
<span id="cb91-24"><a href="estimating-birth-death-models.html#cb91-24" aria-hidden="true" tabindex="-1"></a>} <span class="cf">else</span> {</span>
<span id="cb91-25"><a href="estimating-birth-death-models.html#cb91-25" aria-hidden="true" tabindex="-1"></a>    lambda <span class="ot">&lt;-</span> <span class="cf">function</span>(x) { <span class="fu">return</span> ( <span class="fu">rep</span>(speciation[<span class="dv">1</span>],<span class="fu">length</span>(x)) ) }</span>
<span id="cb91-26"><a href="estimating-birth-death-models.html#cb91-26" aria-hidden="true" tabindex="-1"></a>    mu     <span class="ot">&lt;-</span> <span class="cf">function</span>(x) { <span class="fu">return</span> ( <span class="fu">rep</span>(extinction[<span class="dv">1</span>],<span class="fu">length</span>(x)) ) }</span>
<span id="cb91-27"><a href="estimating-birth-death-models.html#cb91-27" aria-hidden="true" tabindex="-1"></a>}</span></code></pre></div>
<p>Let us make a data frame representation of the rate functions</p>
<div class="sourceCode" id="cb92"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb92-1"><a href="estimating-birth-death-models.html#cb92-1" aria-hidden="true" tabindex="-1"></a>times <span class="ot">&lt;-</span> <span class="fu">seq</span>(<span class="dv">0</span>, max_t, <span class="at">length.out =</span> <span class="dv">100</span>)</span>
<span id="cb92-2"><a href="estimating-birth-death-models.html#cb92-2" aria-hidden="true" tabindex="-1"></a></span>
<span id="cb92-3"><a href="estimating-birth-death-models.html#cb92-3" aria-hidden="true" tabindex="-1"></a>primates_ebd_treepar <span class="ot">&lt;-</span> <span class="fu">tibble</span>(</span>
<span id="cb92-4"><a href="estimating-birth-death-models.html#cb92-4" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;time&quot;</span> <span class="ot">=</span> times,</span>
<span id="cb92-5"><a href="estimating-birth-death-models.html#cb92-5" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;lambda&quot;</span> <span class="ot">=</span> <span class="fu">lambda</span>(times),</span>
<span id="cb92-6"><a href="estimating-birth-death-models.html#cb92-6" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;mu&quot;</span> <span class="ot">=</span> <span class="fu">mu</span>(times),</span>
<span id="cb92-7"><a href="estimating-birth-death-models.html#cb92-7" aria-hidden="true" tabindex="-1"></a>)</span>
<span id="cb92-8"><a href="estimating-birth-death-models.html#cb92-8" aria-hidden="true" tabindex="-1"></a><span class="fu">head</span>(primates_ebd_treepar)</span></code></pre></div>
<pre><code>## # A tibble: 6 × 3
##    time lambda         mu
##   &lt;dbl&gt;  &lt;dbl&gt;      &lt;dbl&gt;
## 1 0      0.392 0.00000181
## 2 0.721  0.392 0.00000181
## 3 1.44   0.392 0.00000181
## 4 2.16   0.392 0.00000181
## 5 2.88   0.392 0.00000181
## 6 3.60   0.243 0.216</code></pre>
<p>Now let’s have a look at the estimated rate functions</p>
<p>The data frame <code>primates_ebd_treepar</code> is what we get when we call the <code>data("primates_ebd_treepar")</code> function in <code>ACDC</code>. We use the RevBayes estimated rates in our vignette examples, but it is easy to substitute those for the two-epoch TreePar results to construct the congruence class.</p>
</div>
<div id="revbayes" class="section level2 hasAnchor" number="8.4">
<h2><span class="header-section-number">8.4</span> The HSMRF model in RevBayes<a href="estimating-birth-death-models.html#revbayes" class="anchor-section" aria-label="Anchor link to header"></a></h2>
<p>RevBayes <span class="citation">(Höhna et al., 2016)</span> provides several methods to estimate diversification rates through time.
Here we used the <a href="RevBayes%20tutorial">https://revbayes.github.io/tutorials/divrate/ebd.html</a> using the episodic-birth-death model with a Horseshoe-Markov-Random-Field (HSMRF) prior distribution to estimate diversification rates <span class="citation">(Magee et al., 2020)</span>.
We assume here that you followed the RevBayes tutorial and have the output stored as a log-file in the <em>output</em> folder.
If you have not run a RevBayes analysis and wish to explore functionality of the package, we have provided a thinned output file, accessible with <code>data("primates_ebd_log")</code>.</p>
<p>For some of our plotting, we need the speciation times.
So we extract these from the tree.</p>
<div class="sourceCode" id="cb94"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb94-1"><a href="estimating-birth-death-models.html#cb94-1" aria-hidden="true" tabindex="-1"></a>btimes <span class="ot">&lt;-</span> <span class="fu">sort</span>( <span class="fu">as.numeric</span>( <span class="fu">branching.times</span>( tree ) ) )</span>
<span id="cb94-2"><a href="estimating-birth-death-models.html#cb94-2" aria-hidden="true" tabindex="-1"></a>max_t <span class="ot">&lt;-</span> <span class="fu">max</span>(btimes)</span></code></pre></div>
<p>Then we read in the RevBayes log-file as a simple tab-delimited file.</p>
<div class="sourceCode" id="cb95"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb95-1"><a href="estimating-birth-death-models.html#cb95-1" aria-hidden="true" tabindex="-1"></a><span class="co">#samples &lt;- read.table(file=&quot;output/HSMRFBDP_primates.log&quot;,</span></span>
<span id="cb95-2"><a href="estimating-birth-death-models.html#cb95-2" aria-hidden="true" tabindex="-1"></a><span class="co">#                      stringsAsFactors=FALSE,</span></span>
<span id="cb95-3"><a href="estimating-birth-death-models.html#cb95-3" aria-hidden="true" tabindex="-1"></a><span class="co">#                      header=TRUE)</span></span>
<span id="cb95-4"><a href="estimating-birth-death-models.html#cb95-4" aria-hidden="true" tabindex="-1"></a><span class="fu">data</span>(<span class="st">&quot;primates_ebd_log&quot;</span>)</span>
<span id="cb95-5"><a href="estimating-birth-death-models.html#cb95-5" aria-hidden="true" tabindex="-1"></a>samples <span class="ot">&lt;-</span> primates_ebd_log</span></code></pre></div>
<p>Then, from this log-file we want to automatically extract the speciation and extinction rates.
Be careful that you used the same names as we did.
Otherwise you need to modify this code-snippet accordingly.</p>
<div class="sourceCode" id="cb96"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb96-1"><a href="estimating-birth-death-models.html#cb96-1" aria-hidden="true" tabindex="-1"></a>par_speciation <span class="ot">&lt;-</span> samples[, <span class="fu">startsWith</span>(<span class="fu">names</span>(samples), <span class="st">&quot;speciation_rate.&quot;</span>)]</span>
<span id="cb96-2"><a href="estimating-birth-death-models.html#cb96-2" aria-hidden="true" tabindex="-1"></a>par_extinction <span class="ot">&lt;-</span> samples[, <span class="fu">startsWith</span>(<span class="fu">names</span>(samples), <span class="st">&quot;extinction_rate.&quot;</span>)]</span></code></pre></div>
<p>In principle we can use now the full posterior estimates.
However, in this example we are only interested in the posterior median.
We strongly encourage to use the medians and not the means for HSMRF models, because the posterior median is a more robust estimate for rates, see <span class="citation">Magee et al. (2020)</span>.</p>
<div class="sourceCode" id="cb97"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb97-1"><a href="estimating-birth-death-models.html#cb97-1" aria-hidden="true" tabindex="-1"></a>est_speciation <span class="ot">&lt;-</span> <span class="fu">apply</span>(par_speciation, <span class="dv">2</span>, median)</span>
<span id="cb97-2"><a href="estimating-birth-death-models.html#cb97-2" aria-hidden="true" tabindex="-1"></a>est_extinction <span class="ot">&lt;-</span> <span class="fu">apply</span>(par_extinction, <span class="dv">2</span>, median)</span>
<span id="cb97-3"><a href="estimating-birth-death-models.html#cb97-3" aria-hidden="true" tabindex="-1"></a>times <span class="ot">&lt;-</span> <span class="fu">seq</span>(<span class="dv">0</span>, max_t, <span class="at">length.out =</span> <span class="fu">length</span>(est_speciation))</span></code></pre></div>
<p>We can also reorganize the rates in a data frame</p>
<div class="sourceCode" id="cb98"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb98-1"><a href="estimating-birth-death-models.html#cb98-1" aria-hidden="true" tabindex="-1"></a>primates_ebd_revbayes <span class="ot">&lt;-</span> <span class="fu">tibble</span>(</span>
<span id="cb98-2"><a href="estimating-birth-death-models.html#cb98-2" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;lambda&quot;</span> <span class="ot">=</span> est_speciation,</span>
<span id="cb98-3"><a href="estimating-birth-death-models.html#cb98-3" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;mu&quot;</span> <span class="ot">=</span> est_extinction,</span>
<span id="cb98-4"><a href="estimating-birth-death-models.html#cb98-4" aria-hidden="true" tabindex="-1"></a>  <span class="st">&quot;time&quot;</span> <span class="ot">=</span> times,</span>
<span id="cb98-5"><a href="estimating-birth-death-models.html#cb98-5" aria-hidden="true" tabindex="-1"></a>)</span></code></pre></div>
<p>Now let’s have a look at the estimated rate functions</p>
<div class="sourceCode" id="cb99"><pre class="sourceCode r"><code class="sourceCode r"><span id="cb99-1"><a href="estimating-birth-death-models.html#cb99-1" aria-hidden="true" tabindex="-1"></a>p <span class="ot">&lt;-</span> primates_ebd_revbayes <span class="sc">%&gt;%</span></span>
<span id="cb99-2"><a href="estimating-birth-death-models.html#cb99-2" aria-hidden="true" tabindex="-1"></a>  tidyr<span class="sc">::</span><span class="fu">gather</span>(<span class="st">&quot;rate&quot;</span>, <span class="st">&quot;value&quot;</span>, <span class="sc">-</span>time) <span class="sc">%&gt;%</span></span>
<span id="cb99-3"><a href="estimating-birth-death-models.html#cb99-3" aria-hidden="true" tabindex="-1"></a>  <span class="fu">ggplot</span>(<span class="fu">aes</span>(<span class="at">x =</span> time, <span class="at">y =</span> value, <span class="at">color =</span> rate)) <span class="sc">+</span></span>
<span id="cb99-4"><a href="estimating-birth-death-models.html#cb99-4" aria-hidden="true" tabindex="-1"></a>  <span class="fu">geom_line</span>() <span class="sc">+</span></span>
<span id="cb99-5"><a href="estimating-birth-death-models.html#cb99-5" aria-hidden="true" tabindex="-1"></a>  <span class="fu">theme_classic</span>() <span class="sc">+</span></span>
<span id="cb99-6"><a href="estimating-birth-death-models.html#cb99-6" aria-hidden="true" tabindex="-1"></a>  <span class="fu">scale_x_reverse</span>() <span class="sc">+</span></span>
<span id="cb99-7"><a href="estimating-birth-death-models.html#cb99-7" aria-hidden="true" tabindex="-1"></a>  <span class="fu">theme</span>(<span class="at">legend.position =</span> <span class="fu">c</span>(<span class="fl">0.2</span>, <span class="fl">0.6</span>)) <span class="sc">+</span></span>
<span id="cb99-8"><a href="estimating-birth-death-models.html#cb99-8" aria-hidden="true" tabindex="-1"></a>  <span class="fu">labs</span>(<span class="at">y =</span> <span class="st">&quot;rate&quot;</span>, <span class="at">x =</span> <span class="st">&quot;time before present&quot;</span>)</span>
<span id="cb99-9"><a href="estimating-birth-death-models.html#cb99-9" aria-hidden="true" tabindex="-1"></a><span class="fu">plot</span>(p)</span></code></pre></div>
<div class="figure" style="text-align: center"><span style="display:block;" id="fig:plot-revbayes-rates"></span>
<img src="figures/plot-revbayes-rates-1.png" alt="The estimated posterior median results from an episodic birth-death model fitted in RevBayes." width="480" />
<p class="caption">
Figure 8.2: The estimated posterior median results from an episodic birth-death model fitted in RevBayes.
</p>
</div>
<p>The data frame <code>primates_ebd_revbayes</code> is essentially equal to the data frame we get when we call the <code>data("primates_ebd")</code> function in <code>CRABS</code>. We use these estimated rates in our vignette examples, but we note that it is possible to use any kind of rate functions to construct the congruence class.</p>
</div>
</div>
            </section>

          </div>
        </div>
      </div>
<a href="exploring-the-posterior.html" class="navigation navigation-prev " aria-label="Previous page"><i class="fa fa-angle-left"></i></a>
<a href="references.html" class="navigation navigation-next " aria-label="Next page"><i class="fa fa-angle-right"></i></a>
    </div>
  </div>
<script src="lib/gitbook-2.6.7/js/app.min.js"></script>
<script src="lib/gitbook-2.6.7/js/clipboard.min.js"></script>
<script src="lib/gitbook-2.6.7/js/plugin-search.js"></script>
<script src="lib/gitbook-2.6.7/js/plugin-sharing.js"></script>
<script src="lib/gitbook-2.6.7/js/plugin-fontsettings.js"></script>
<script src="lib/gitbook-2.6.7/js/plugin-bookdown.js"></script>
<script src="lib/gitbook-2.6.7/js/jquery.highlight.js"></script>
<script src="lib/gitbook-2.6.7/js/plugin-clipboard.js"></script>
<script>
gitbook.require(["gitbook"], function(gitbook) {
gitbook.start({
"sharing": {
"github": false,
"facebook": true,
"twitter": true,
"linkedin": false,
"weibo": false,
"instapaper": false,
"vk": false,
"whatsapp": false,
"all": ["facebook", "twitter", "linkedin", "weibo", "instapaper"]
},
"fontsettings": {
"theme": "white",
"family": "sans",
"size": 2
},
"edit": {
"link": null,
"text": null
},
"history": {
"link": null,
"text": null
},
"view": {
"link": null,
"text": null
},
"download": null,
"search": false,
"toc": {
"collapse": "subsection"
}
});
});
</script>

<!-- dynamically load mathjax for compatibility with self-contained -->
<script>
  (function () {
    var script = document.createElement("script");
    script.type = "text/javascript";
    var src = "true";
    if (src === "" || src === "true") src = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.9/latest.js?config=TeX-MML-AM_CHTML";
    if (location.protocol !== "file:")
      if (/^https?:/.test(src))
        src = src.replace(/^https?:/, '');
    script.src = src;
    document.getElementsByTagName("head")[0].appendChild(script);
  })();
</script>
</body>

</html>
