outbreaker2 / Git / [dfe06d] /docs/articles/customisation.html

Models:
MarcoTheBlack/
outbreaker2
Downloads: 1
[dfe06d]: / docs / articles / customisation.html
History
Download this file
1072 lines (1001 with data), 92.5 kB

<!DOCTYPE html>
<!-- Generated by pkgdown: do not edit by hand --><html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>Using custom priors, likelihood, or movements in outbreaker2 • outbreaker2</title>
<!-- jquery --><script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.4.1/jquery.min.js" integrity="sha256-CSXorXvZcTkaix6Yvo6HppcZGetbYMGWSFlBw8HfCJo=" crossorigin="anonymous"></script><!-- Bootstrap --><link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/css/bootstrap.min.css" integrity="sha256-bZLfwXAP04zRMK2BjiO8iu9pf4FbLqX6zitd+tIvLhE=" crossorigin="anonymous">
<script src="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/js/bootstrap.min.js" integrity="sha256-nuL8/2cJ5NDSSwnKD8VqreErSWHtnEP9E7AySL+1ev4=" crossorigin="anonymous"></script><!-- bootstrap-toc --><link rel="stylesheet" href="../bootstrap-toc.css">
<script src="../bootstrap-toc.js"></script><!-- Font Awesome icons --><link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/all.min.css" integrity="sha256-mmgLkCYLUQbXn0B1SRqzHar6dCnv9oZFPEC1g1cwlkk=" crossorigin="anonymous">
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/v4-shims.min.css" integrity="sha256-wZjR52fzng1pJHwx4aV2AO3yyTOXrcDW7jBpJtTwVxw=" crossorigin="anonymous">
<!-- clipboard.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script><!-- headroom.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/headroom.min.js" integrity="sha256-AsUX4SJE1+yuDu5+mAVzJbuYNPHj/WroHuZ8Ir/CkE0=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/jQuery.headroom.min.js" integrity="sha256-ZX/yNShbjqsohH1k95liqY9Gd8uOiE1S4vZc+9KQ1K4=" crossorigin="anonymous"></script><!-- pkgdown --><link href="../pkgdown.css" rel="stylesheet">
<script src="../pkgdown.js"></script><meta property="og:title" content="Using custom priors, likelihood, or movements in outbreaker2">
<meta property="og:description" content="outbreaker2">
<!-- mathjax --><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js" integrity="sha256-nvJJv9wWKEm88qvoQl9ekL2J+k/RWIsaSScxxlsrv8k=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/config/TeX-AMS-MML_HTMLorMML.js" integrity="sha256-84DKXVJXs0/F8OTMzX4UR909+jtl4G7SPypPavF+GfA=" crossorigin="anonymous"></script><!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]-->
</head>
<body data-spy="scroll" data-target="#toc">
    <div class="container template-article">
      <header><div class="navbar navbar-default navbar-fixed-top" role="navigation">
  <div class="container">
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar" aria-expanded="false">
        <span class="sr-only">Toggle navigation</span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
      <span class="navbar-brand">
        <a class="navbar-link" href="../index.html">outbreaker2</a>
        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">1.1.2</span>
      </span>
    </div>

    <div id="navbar" class="navbar-collapse collapse">
      <ul class="nav navbar-nav">
<li>
  <a href="../index.html">
    <span class="fas fa-home fa-lg"></span>
     
  </a>
</li>
<li>
  <a href="../reference/index.html">Reference</a>
</li>
<li class="dropdown">
  <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
    Articles
     
    <span class="caret"></span>
  </a>
  <ul class="dropdown-menu" role="menu">
<li>
      <a href="../articles/customisation.html">Using custom priors, likelihood, or movements in outbreaker2</a>
    </li>
    <li>
      <a href="../articles/introduction.html">Introduction to outbreaker2</a>
    </li>
    <li>
      <a href="../articles/overview.html">outbreaker2: package overview</a>
    </li>
    <li>
      <a href="../articles/Rcpp_API.html">outbreaker2: Rcpp API</a>
    </li>
  </ul>
</li>
<li>
  <a href="../news/index.html">Changelog</a>
</li>
      </ul>
<ul class="nav navbar-nav navbar-right"></ul>
</div>
<!--/.nav-collapse -->
  </div>
<!--/.container -->
</div>
<!--/.navbar -->

      

      </header><script src="customisation_files/header-attrs-2.6/header-attrs.js"></script><script src="customisation_files/accessible-code-block-0.0.1/empty-anchor.js"></script><div class="row">
  <div class="col-md-9 contents">
    <div class="page-header toc-ignore">
      <h1 data-toc-skip>Using custom priors, likelihood, or movements in outbreaker2</h1>
                        <h4 class="author">Thibaut Jombart</h4>
            
            <h4 class="date">2021-02-09</h4>
      
      
      <div class="hidden name"><code>customisation.Rmd</code></div>

    </div>

    
    
<p>In this vignette, we show how custom functions for priors, likelihood, or movement of parameters and augmented data can be used in <em>outbreaker2</em>. In all these functions, the process will be similar:</p>
<ol style="list-style-type: decimal">
<li>write your own function with the right arguments</li>
<li>pass this function as an argument to a <code>custom...</code> function</li>
<li>pass the result to <em>outbreaker2</em>
</li>
</ol>
<p>Note that 2-3 can be a single step if passing the function to the arguments of <em>outbreaker2</em> directly. Also note that <strong>all priors and likelihoods are expected on a log scale</strong>. Finally, also note that while the various <code>custom...</code> functions will try to some extent to check that the provided functions are valid, such tests are very difficult to implement. In short: you are using these custom features at your own risks - make sure these functions work before passing them to <em>outbreaker2</em>.</p>
<p><br></p>
<div id="customising-priors" class="section level1">
<h1 class="hasAnchor">
<a href="#customising-priors" class="anchor"></a>Customising priors</h1>
<p>Priors of <em>outbreaker2</em> must be a function of an <code>outbreaker_param</code> list (see <code><a href="../reference/create_param.html">?outbreaker_param</a></code>). Here, we decide to use a step function rather than the default Beta function as a prior for <em>pi</em>, the reporting probability, and a flat prior between 0 and 1 for the mutation rate (which is technically a probability in the basic genetic model used in <em>outbreaker2</em>).</p>
<p>We start by defining two functions: an auxiliary function <code>f</code> which returns values on the natural scale, and which we can use for plotting the prior distribution, and then a function <code>f_pi</code> which will be used for the customisation.</p>
<div class="sourceCode" id="cb1"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="va">f</span> <span class="op">&lt;-</span> <span class="kw">function</span><span class="op">(</span><span class="va">pi</span><span class="op">)</span> <span class="op">{</span>
    <span class="fu"><a href="https://rdrr.io/r/base/ifelse.html">ifelse</a></span><span class="op">(</span><span class="va">pi</span> <span class="op">&lt;</span> <span class="fl">0.8</span>, <span class="fl">0</span>, <span class="fl">5</span><span class="op">)</span>
<span class="op">}</span>

<span class="va">f_pi</span> <span class="op">&lt;-</span> <span class="kw">function</span><span class="op">(</span><span class="va">param</span><span class="op">)</span> <span class="op">{</span> 
    <span class="fu"><a href="https://rdrr.io/r/base/Log.html">log</a></span><span class="op">(</span><span class="fu">f</span><span class="op">(</span><span class="va">param</span><span class="op">$</span><span class="va">pi</span><span class="op">)</span><span class="op">)</span>
<span class="op">}</span>

<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">f</span>, type <span class="op">=</span> <span class="st">"s"</span>, col <span class="op">=</span> <span class="st">"blue"</span>, 
     xlab <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/expression.html">expression</a></span><span class="op">(</span><span class="va">pi</span><span class="op">)</span>, ylab <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/expression.html">expression</a></span><span class="op">(</span><span class="fu">p</span><span class="op">(</span><span class="va">pi</span><span class="op">)</span><span class="op">)</span>, 
     main <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/expression.html">expression</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/paste.html">paste</a></span><span class="op">(</span><span class="st">"New prior for "</span>, <span class="va">pi</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/f_pi-1.png" width="768"></p>
<p>While <code>f</code> is a useful function to visualise the prior, <code>f_pi</code> is the function which will be passed to <code>outbreaker</code>. To do so, we pass it to <code>custom_priors</code>:</p>
<div class="sourceCode" id="cb2"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="kw"><a href="https://rdrr.io/r/base/library.html">library</a></span><span class="op">(</span><span class="va">outbreaker2</span><span class="op">)</span>

<span class="va">f_mu</span> <span class="op">&lt;-</span> <span class="kw">function</span><span class="op">(</span><span class="va">param</span><span class="op">)</span> <span class="op">{</span>
  <span class="kw">if</span> <span class="op">(</span><span class="va">param</span><span class="op">$</span><span class="va">mu</span> <span class="op">&lt;</span> <span class="fl">0</span> <span class="op">||</span> <span class="va">param</span><span class="op">$</span><span class="va">mu</span> <span class="op">&gt;</span> <span class="fl">1</span><span class="op">)</span> <span class="op">{</span>
    <span class="kw"><a href="https://rdrr.io/r/base/function.html">return</a></span><span class="op">(</span><span class="op">-</span><span class="cn">Inf</span><span class="op">)</span>
  <span class="op">}</span> <span class="kw">else</span> <span class="op">{</span>
    <span class="kw"><a href="https://rdrr.io/r/base/function.html">return</a></span><span class="op">(</span><span class="fl">0.0</span><span class="op">)</span>
  <span class="op">}</span>
  
<span class="op">}</span>

<span class="va">priors</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/custom_priors.html">custom_priors</a></span><span class="op">(</span>pi <span class="op">=</span> <span class="va">f_pi</span>, mu <span class="op">=</span> <span class="va">f_mu</span><span class="op">)</span>
<span class="va">priors</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt;  ///// outbreaker custom priors ///</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; class: custom_priors list</span>
<span class="co">#&gt; number of items: 4 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; /// custom priors set to NULL (default used) //</span>
<span class="co">#&gt; $eps</span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $lambda</span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; /// custom priors //</span>
<span class="co">#&gt; $mu</span>
<span class="co">#&gt; function(param) {</span>
<span class="co">#&gt;   if (param$mu &lt; 0 || param$mu &gt; 1) {</span>
<span class="co">#&gt;     return(-Inf)</span>
<span class="co">#&gt;   } else {</span>
<span class="co">#&gt;     return(0.0)</span>
<span class="co">#&gt;   }</span>
<span class="co">#&gt;   </span>
<span class="co">#&gt; }</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $pi</span>
<span class="co">#&gt; function(param) { </span>
<span class="co">#&gt;     log(f(param$pi))</span>
<span class="co">#&gt; }</span></code></pre></div>
<p>Note that <code>custom_priors</code> does more than just adding the custom function to a list. For instance, the following customisations are all wrong, and rightfully rejected:</p>
<div class="sourceCode" id="cb3"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="co">## wrong: not a function</span>
<span class="co">## should be pi = function(x){0.0}</span>
<span class="fu"><a href="../reference/custom_priors.html">custom_priors</a></span><span class="op">(</span>pi <span class="op">=</span> <span class="fl">0.0</span><span class="op">)</span>
<span class="co">#&gt; Error in custom_priors(pi = 0): The following priors are not functions: pi</span>

<span class="co">## wrong: two arguments</span>
<span class="fu"><a href="../reference/custom_priors.html">custom_priors</a></span><span class="op">(</span>pi <span class="op">=</span> <span class="kw">function</span><span class="op">(</span><span class="va">x</span>, <span class="va">y</span><span class="op">)</span><span class="op">{</span><span class="fl">0.0</span><span class="op">}</span><span class="op">)</span>
<span class="co">#&gt; Error in custom_priors(pi = function(x, y) {: The following priors dont' have a single argument: pi</span></code></pre></div>
<p>We can now use the new priors to run <code>outbreaker</code> on the <code>fake_outbreak</code> data (see <a href="introduction.html"><em>introduction vignette</em></a>):</p>
<div class="sourceCode" id="cb4"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">dna</span> <span class="op">&lt;-</span> <span class="va">fake_outbreak</span><span class="op">$</span><span class="va">dna</span>
<span class="va">dates</span> <span class="op">&lt;-</span> <span class="va">fake_outbreak</span><span class="op">$</span><span class="va">sample</span>
<span class="va">w</span> <span class="op">&lt;-</span> <span class="va">fake_outbreak</span><span class="op">$</span><span class="va">w</span>
<span class="va">data</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/outbreaker_data.html">outbreaker_data</a></span><span class="op">(</span>dna <span class="op">=</span> <span class="va">dna</span>, dates <span class="op">=</span> <span class="va">dates</span>, w_dens <span class="op">=</span> <span class="va">w</span><span class="op">)</span>

<span class="co">## we set the seed to ensure results won't change</span>
<span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span><span class="op">(</span><span class="fl">1</span><span class="op">)</span>


<span class="va">res</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/outbreaker.html">outbreaker</a></span><span class="op">(</span>data <span class="op">=</span> <span class="va">data</span>, priors <span class="op">=</span> <span class="va">priors</span><span class="op">)</span></code></pre></div>
<p>We can check the results first by looking at the traces, and then by plotting the posterior distributions of <code>pi</code> and <code>mu</code>, respectively:</p>
<div class="sourceCode" id="cb5"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/traces_custom_priors-1.png" width="768"></p>
<div class="sourceCode" id="cb6"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res</span>, <span class="st">"pi"</span>, burnin <span class="op">=</span> <span class="fl">500</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/traces_custom_priors-2.png" width="768"></p>
<div class="sourceCode" id="cb7"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res</span>, <span class="st">"mu"</span>, burnin <span class="op">=</span> <span class="fl">500</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/traces_custom_priors-3.png" width="768"></p>
<div class="sourceCode" id="cb8"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res</span>, <span class="st">"pi"</span>, type <span class="op">=</span> <span class="st">"density"</span>, burnin <span class="op">=</span> <span class="fl">500</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/traces_custom_priors-4.png" width="768"></p>
<div class="sourceCode" id="cb9"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res</span>, <span class="st">"mu"</span>, type <span class="op">=</span> <span class="st">"hist"</span>, burnin <span class="op">=</span> <span class="fl">500</span><span class="op">)</span>
<span class="co">#&gt; `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.</span></code></pre></div>
<p><img src="figs-customisation/traces_custom_priors-5.png" width="768"></p>
<p>Note that we are using density and histograms here for illustrative purposes, but there is no reason to prefer one or the other for a specific parameter.</p>
<p>Interestingly, the trace of <code>pi</code> suggests that the MCMC oscillates between two different states, on either bound of the interval on which the prior is positive (it is <code>-Inf</code> outside (0.8; 1)). This may be a consequence of the step function, which causes sharp ‘cliffs’ in the posterior landscape. What shall one do to derive good samples from the posterior distribution in this kind of situation? There are several options, which in fact apply to typical cases of multi-modal posterior distributions:</p>
<ul>
<li><p>Avoid ‘cliffs’, i.e. sharp drops in the posterior landscape, typically created by using step-functions in likelihoods and in priors.</p></li>
<li><p>Use larger samples, i.e. run more MCMC iterations.</p></li>
<li><p>Use a different sampler, better than Metropolis-Hasting at deriving samples from multi-modal distributions.</p></li>
</ul>
<p>Because we know what the real transmission tree is for this dataset, we can assess how the new priors impacted the inference of the transmission tree.</p>
<div class="sourceCode" id="cb10"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span><span class="op">(</span><span class="va">res</span>, burnin <span class="op">=</span> <span class="fl">500</span><span class="op">)</span>
<span class="co">#&gt; $step</span>
<span class="co">#&gt;    first     last interval  n_steps </span>
<span class="co">#&gt;      550    10000       50      190 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $post</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;  -481.1  -467.5  -463.9  -464.7  -461.4  -456.7 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $like</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;  -482.7  -469.1  -465.5  -466.3  -463.0  -458.3 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $prior</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;   1.609   1.609   1.609   1.609   1.609   1.609 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $mu</span>
<span class="co">#&gt;      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. </span>
<span class="co">#&gt; 8.215e-05 1.254e-04 1.400e-04 1.419e-04 1.566e-04 2.159e-04 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $pi</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;  0.8031  0.9207  0.9556  0.9453  0.9813  0.9999 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $tree</span>
<span class="co">#&gt;    from to time     support generations</span>
<span class="co">#&gt; 1    NA  1   -1          NA          NA</span>
<span class="co">#&gt; 2     1  2    1 1.000000000           1</span>
<span class="co">#&gt; 3     2  3    3 1.000000000           1</span>
<span class="co">#&gt; 4    NA  4    3 0.005263158          NA</span>
<span class="co">#&gt; 5     3  5    4 0.994736842           1</span>
<span class="co">#&gt; 6     9  6    6 1.000000000           1</span>
<span class="co">#&gt; 7     4  7    5 0.994736842           1</span>
<span class="co">#&gt; 8     5  8    6 0.989473684           1</span>
<span class="co">#&gt; 9     4  9    5 0.968421053           1</span>
<span class="co">#&gt; 10    6 10    8 1.000000000           1</span>
<span class="co">#&gt; 11    7 11    7 0.689473684           1</span>
<span class="co">#&gt; 12    5 12    7 0.826315789           1</span>
<span class="co">#&gt; 13    9 13    7 1.000000000           1</span>
<span class="co">#&gt; 14    5 14    7 0.768421053           1</span>
<span class="co">#&gt; 15    5 15    7 0.752631579           1</span>
<span class="co">#&gt; 16    7 16    8 0.794736842           1</span>
<span class="co">#&gt; 17    7 17    7 0.626315789           1</span>
<span class="co">#&gt; 18    8 18    9 0.431578947           1</span>
<span class="co">#&gt; 19    9 19    8 1.000000000           1</span>
<span class="co">#&gt; 20   10 20   10 0.968421053           1</span>
<span class="co">#&gt; 21   11 21   10 0.973684211           1</span>
<span class="co">#&gt; 22   11 22   10 1.000000000           1</span>
<span class="co">#&gt; 23   13 23    9 1.000000000           1</span>
<span class="co">#&gt; 24   13 24    9 1.000000000           1</span>
<span class="co">#&gt; 25   13 25    9 1.000000000           1</span>
<span class="co">#&gt; 26   17 26    9 1.000000000           1</span>
<span class="co">#&gt; 27   17 27   10 1.000000000           1</span>
<span class="co">#&gt; 28   NA 28    9          NA          NA</span>
<span class="co">#&gt; 29   10 29   11 1.000000000           1</span>
<span class="co">#&gt; 30   13 30   10 1.000000000           1</span>
<span class="va">tree</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span><span class="op">(</span><span class="va">res</span>, burnin <span class="op">=</span> <span class="fl">500</span><span class="op">)</span><span class="op">$</span><span class="va">tree</span>

<span class="va">comparison</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html">data.frame</a></span><span class="op">(</span>case <span class="op">=</span> <span class="fl">1</span><span class="op">:</span><span class="fl">30</span>,
                         inferred <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/paste.html">paste</a></span><span class="op">(</span><span class="va">tree</span><span class="op">$</span><span class="va">from</span><span class="op">)</span>,
             true <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/paste.html">paste</a></span><span class="op">(</span><span class="va">fake_outbreak</span><span class="op">$</span><span class="va">ances</span><span class="op">)</span>,
             stringsAsFactors <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span>
             
<span class="va">comparison</span><span class="op">$</span><span class="va">correct</span> <span class="op">&lt;-</span> <span class="va">comparison</span><span class="op">$</span><span class="va">inferred</span> <span class="op">==</span> <span class="va">comparison</span><span class="op">$</span><span class="va">true</span>
<span class="va">comparison</span>
<span class="co">#&gt;    case inferred true correct</span>
<span class="co">#&gt; 1     1       NA   NA    TRUE</span>
<span class="co">#&gt; 2     2        1    1    TRUE</span>
<span class="co">#&gt; 3     3        2    2    TRUE</span>
<span class="co">#&gt; 4     4       NA   NA    TRUE</span>
<span class="co">#&gt; 5     5        3    3    TRUE</span>
<span class="co">#&gt; 6     6        9    4   FALSE</span>
<span class="co">#&gt; 7     7        4    4    TRUE</span>
<span class="co">#&gt; 8     8        5    5    TRUE</span>
<span class="co">#&gt; 9     9        4    6   FALSE</span>
<span class="co">#&gt; 10   10        6    6    TRUE</span>
<span class="co">#&gt; 11   11        7    7    TRUE</span>
<span class="co">#&gt; 12   12        5    8   FALSE</span>
<span class="co">#&gt; 13   13        9    9    TRUE</span>
<span class="co">#&gt; 14   14        5    5    TRUE</span>
<span class="co">#&gt; 15   15        5    5    TRUE</span>
<span class="co">#&gt; 16   16        7    7    TRUE</span>
<span class="co">#&gt; 17   17        7    7    TRUE</span>
<span class="co">#&gt; 18   18        8    8    TRUE</span>
<span class="co">#&gt; 19   19        9    9    TRUE</span>
<span class="co">#&gt; 20   20       10   10    TRUE</span>
<span class="co">#&gt; 21   21       11   11    TRUE</span>
<span class="co">#&gt; 22   22       11   11    TRUE</span>
<span class="co">#&gt; 23   23       13   13    TRUE</span>
<span class="co">#&gt; 24   24       13   13    TRUE</span>
<span class="co">#&gt; 25   25       13   13    TRUE</span>
<span class="co">#&gt; 26   26       17   17    TRUE</span>
<span class="co">#&gt; 27   27       17   17    TRUE</span>
<span class="co">#&gt; 28   28       NA   NA    TRUE</span>
<span class="co">#&gt; 29   29       10   10    TRUE</span>
<span class="co">#&gt; 30   30       13   13    TRUE</span>
<span class="fu"><a href="https://rdrr.io/r/base/mean.html">mean</a></span><span class="op">(</span><span class="va">comparison</span><span class="op">$</span><span class="va">correct</span><span class="op">)</span>
<span class="co">#&gt; [1] 0.9</span></code></pre></div>
<p><br></p>
</div>
<div id="customizing-likelihood" class="section level1">
<h1 class="hasAnchor">
<a href="#customizing-likelihood" class="anchor"></a>Customizing likelihood</h1>
<p>Likelihood functions customisation works identically to prior functions. The only difference is that custom functions will take two arguments (<code>data</code> and <code>param</code>) instead of one in the prior functions. The function used to specify custom likelihood is <code>custom_likelihoods</code>. Each custom function will correspond to a specific likelihood component:</p>
<div class="sourceCode" id="cb11"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="../reference/custom_likelihoods.html">custom_likelihoods</a></span><span class="op">(</span><span class="op">)</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt;  ///// outbreaker custom likelihoods ///</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; class: custom_likelihoods list</span>
<span class="co">#&gt; number of items: 5 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; /// custom likelihoods //</span>
<span class="co">#&gt; $genetic</span>
<span class="co">#&gt; $genetic[[1]]</span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $genetic[[2]]</span>
<span class="co">#&gt; [1] 0</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $reporting</span>
<span class="co">#&gt; $reporting[[1]]</span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $reporting[[2]]</span>
<span class="co">#&gt; [1] 0</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $timing_infections</span>
<span class="co">#&gt; $timing_infections[[1]]</span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $timing_infections[[2]]</span>
<span class="co">#&gt; [1] 0</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $timing_sampling</span>
<span class="co">#&gt; $timing_sampling[[1]]</span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $timing_sampling[[2]]</span>
<span class="co">#&gt; [1] 0</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $contact</span>
<span class="co">#&gt; $contact[[1]]</span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $contact[[2]]</span>
<span class="co">#&gt; [1] 0</span></code></pre></div>
<p>see <code><a href="../reference/custom_likelihoods.html">?custom_likelihoods</a></code> for details of these components, and see the section ‘Extending the model’ for new, other components. As for <code>custom_priors</code>, a few checks are performed by <code>custom_likelihoods</code>:</p>
<div class="sourceCode" id="cb12"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="co">## wrong: not a function</span>
<span class="fu"><a href="../reference/custom_likelihoods.html">custom_likelihoods</a></span><span class="op">(</span>genetic <span class="op">=</span> <span class="st">"fubar"</span><span class="op">)</span>
<span class="co">#&gt; Error in custom_likelihoods(genetic = "fubar"): The following likelihoods are not functions: genetic</span>

<span class="co">## wrong: only one argument</span>
<span class="fu"><a href="../reference/custom_likelihoods.html">custom_likelihoods</a></span><span class="op">(</span>genetic <span class="op">=</span> <span class="kw">function</span><span class="op">(</span><span class="va">x</span><span class="op">)</span><span class="op">{</span> <span class="fl">0.0</span> <span class="op">}</span><span class="op">)</span>
<span class="co">#&gt; Error in custom_likelihoods(genetic = function(x) {: The following likelihoods do not have arity two or three: genetic</span></code></pre></div>
<p>A trivial customisation is to disable some or all of the likelihood components of the model by returning a finite constant. Here, we apply this to two cases: first, we will disable all likelihood components as a sanity check, making sure that the transmission tree landscape is explored freely by the MCMC. Second, we will recreate the <a href="http://dx.doi.org/10.1093/aje/kwh255">Wallinga &amp; Teunis (1994)</a> model, by disabling specific components.</p>
<div id="a-null-model" class="section level2">
<h2 class="hasAnchor">
<a href="#a-null-model" class="anchor"></a>A null model</h2>
<div class="sourceCode" id="cb13"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">f_null</span> <span class="op">&lt;-</span> <span class="kw">function</span><span class="op">(</span><span class="va">data</span>, <span class="va">param</span><span class="op">)</span> <span class="op">{</span>
   <span class="kw"><a href="https://rdrr.io/r/base/function.html">return</a></span><span class="op">(</span><span class="fl">0.0</span><span class="op">)</span>
<span class="op">}</span>

<span class="va">null_model</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/custom_likelihoods.html">custom_likelihoods</a></span><span class="op">(</span>genetic <span class="op">=</span> <span class="va">f_null</span>,
                                 timing_sampling <span class="op">=</span> <span class="va">f_null</span>,
                                 timing_infections <span class="op">=</span> <span class="va">f_null</span>,
                                 reporting <span class="op">=</span> <span class="va">f_null</span>,
                                 contact <span class="op">=</span> <span class="va">f_null</span><span class="op">)</span>

<span class="va">null_model</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt;  ///// outbreaker custom likelihoods ///</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; class: custom_likelihoods list</span>
<span class="co">#&gt; number of items: 5 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; /// custom likelihoods //</span>
<span class="co">#&gt; $genetic</span>
<span class="co">#&gt; $genetic[[1]]</span>
<span class="co">#&gt; function(data, param) {</span>
<span class="co">#&gt;    return(0.0)</span>
<span class="co">#&gt; }</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $genetic[[2]]</span>
<span class="co">#&gt; [1] 2</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $reporting</span>
<span class="co">#&gt; $reporting[[1]]</span>
<span class="co">#&gt; function(data, param) {</span>
<span class="co">#&gt;    return(0.0)</span>
<span class="co">#&gt; }</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $reporting[[2]]</span>
<span class="co">#&gt; [1] 2</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $timing_infections</span>
<span class="co">#&gt; $timing_infections[[1]]</span>
<span class="co">#&gt; function(data, param) {</span>
<span class="co">#&gt;    return(0.0)</span>
<span class="co">#&gt; }</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $timing_infections[[2]]</span>
<span class="co">#&gt; [1] 2</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $timing_sampling</span>
<span class="co">#&gt; $timing_sampling[[1]]</span>
<span class="co">#&gt; function(data, param) {</span>
<span class="co">#&gt;    return(0.0)</span>
<span class="co">#&gt; }</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $timing_sampling[[2]]</span>
<span class="co">#&gt; [1] 2</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $contact</span>
<span class="co">#&gt; $contact[[1]]</span>
<span class="co">#&gt; function(data, param) {</span>
<span class="co">#&gt;    return(0.0)</span>
<span class="co">#&gt; }</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $contact[[2]]</span>
<span class="co">#&gt; [1] 2</span></code></pre></div>
<p>We also specify settings via the <code>config</code> argument to avoid detecting imported cases, reduce the number of iterations and sampling each of them:</p>
<div class="sourceCode" id="cb14"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">null_config</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></span><span class="op">(</span>find_import <span class="op">=</span> <span class="cn">FALSE</span>,
n_iter <span class="op">=</span> <span class="fl">500</span>,
sample_every <span class="op">=</span> <span class="fl">1</span><span class="op">)</span>

<span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span><span class="op">(</span><span class="fl">1</span><span class="op">)</span>

<span class="va">res_null</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/outbreaker.html">outbreaker</a></span><span class="op">(</span>data <span class="op">=</span> <span class="va">data</span>,
config <span class="op">=</span> <span class="va">null_config</span>,
likelihoods <span class="op">=</span> <span class="va">null_model</span><span class="op">)</span></code></pre></div>
<div class="sourceCode" id="cb15"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_null</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_null_model-1.png" width="768"></p>
<div class="sourceCode" id="cb16"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_null</span>, <span class="st">"pi"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_null_model-2.png" width="768"></p>
<div class="sourceCode" id="cb17"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_null</span>, <span class="st">"mu"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_null_model-3.png" width="768"></p>
<p>By typical MCMC standards, these traces look appaling, as they haven’t reach stationarity (i.e. same mean and variance over time), and are grossly autocorrelated in parts. Fair enough, as these are only the first 500 iterations of the MCMC, so that autocorrelation is expected. In fact, what we observe here literally is the random walk across the posterior landscape, which in this case is only impacted by the priors.</p>
<p>We can check that transmission trees are indeed freely explored:</p>
<div class="sourceCode" id="cb18"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_null</span>, type <span class="op">=</span> <span class="st">"alpha"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/null_trees-1.png" width="768"></p>
<p>Do <strong>not</strong> try to render the corresponding network using <code><a href="https://rdrr.io/r/graphics/plot.default.html">plot(..., type = "network")</a></code> as the force-direction algorithm will go insane. However, this network can be visualised using <em>igraph</em>, extracting the edges and nodes from the plot (without displaying it):</p>
<div class="sourceCode" id="cb19"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="co">## extract nodes and edges from the visNetwork object</span>
<span class="va">temp</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_null</span>, type <span class="op">=</span> <span class="st">"network"</span>, min_support <span class="op">=</span> <span class="fl">0</span><span class="op">)</span>
<span class="fu"><a href="https://rdrr.io/r/base/class.html">class</a></span><span class="op">(</span><span class="va">temp</span><span class="op">)</span>
<span class="co">#&gt; [1] "visNetwork" "htmlwidget"</span>
<span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span><span class="op">(</span><span class="va">temp</span><span class="op">$</span><span class="va">x</span><span class="op">$</span><span class="va">edges</span><span class="op">)</span>
<span class="co">#&gt;   from to value arrows   color</span>
<span class="co">#&gt; 1    1  1 0.004     to #CCDDFF</span>
<span class="co">#&gt; 2    1  2 0.026     to #CCDDFF</span>
<span class="co">#&gt; 3    1  3 0.036     to #CCDDFF</span>
<span class="co">#&gt; 4    1  4 0.040     to #CCDDFF</span>
<span class="co">#&gt; 5    1  5 0.020     to #CCDDFF</span>
<span class="co">#&gt; 6    1  6 0.042     to #CCDDFF</span>
<span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span><span class="op">(</span><span class="va">temp</span><span class="op">$</span><span class="va">x</span><span class="op">$</span><span class="va">nodes</span><span class="op">)</span>
<span class="co">#&gt;   id label value   color shape shaped</span>
<span class="co">#&gt; 1  1     1 0.916 #CCDDFF   dot   star</span>
<span class="co">#&gt; 2  2     2 0.940 #B2D9E3   dot   &lt;NA&gt;</span>
<span class="co">#&gt; 3  3     3 0.984 #98D6C7   dot   &lt;NA&gt;</span>
<span class="co">#&gt; 4  4     4 0.986 #7ED2AC   dot   &lt;NA&gt;</span>
<span class="co">#&gt; 5  5     5 0.900 #99CAA9   dot   &lt;NA&gt;</span>
<span class="co">#&gt; 6  6     6 1.066 #C2C0AD   dot   &lt;NA&gt;</span>

<span class="co">## make an igraph object</span>
<span class="kw"><a href="https://rdrr.io/r/base/library.html">library</a></span><span class="op">(</span><span class="va"><a href="https://igraph.org">igraph</a></span><span class="op">)</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Attaching package: 'igraph'</span>
<span class="co">#&gt; The following objects are masked from 'package:stats':</span>
<span class="co">#&gt; </span>
<span class="co">#&gt;     decompose, spectrum</span>
<span class="co">#&gt; The following object is masked from 'package:base':</span>
<span class="co">#&gt; </span>
<span class="co">#&gt;     union</span>

<span class="va">net_null</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/pkg/igraph/man/graph_from_data_frame.html">graph.data.frame</a></span><span class="op">(</span><span class="va">temp</span><span class="op">$</span><span class="va">x</span><span class="op">$</span><span class="va">edges</span>,
                             vertices <span class="op">=</span> <span class="va">temp</span><span class="op">$</span><span class="va">x</span><span class="op">$</span><span class="va">nodes</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">4</span><span class="op">]</span><span class="op">)</span>

<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">net_null</span>, layout <span class="op">=</span> <span class="va">layout.circle</span>,
     main <span class="op">=</span> <span class="st">"Null model, posterior trees"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/null_net-1.png" width="768"></p>
<p>We can derive similar diagnostics for the number of generations betweens cases (<code>kappa</code>), only constrained by default settings to be between 1 and 5, and for the infection dates (<em>t_inf</em>):</p>
<div class="sourceCode" id="cb20"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_null</span>, type <span class="op">=</span> <span class="st">"kappa"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_null_diag-1.png" width="768"></p>
<div class="sourceCode" id="cb21"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_null</span>, type <span class="op">=</span> <span class="st">"t_inf"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_null_diag-2.png" width="768"></p>
<p>Finally, we can verify that the distributions of <code>mu</code> and <code>pi</code> match their priors, respectively an exponential distribution with rate 1000 and a beta with parameters 10 and 1. Here, we get a qualitative assessment by comparing the observed distribution (histograms) to the densities of similar sized random samples from the priors:</p>
<div class="sourceCode" id="cb22"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/graphics/par.html">par</a></span><span class="op">(</span>xpd<span class="op">=</span><span class="cn">TRUE</span><span class="op">)</span>
<span class="fu"><a href="https://rdrr.io/r/graphics/hist.html">hist</a></span><span class="op">(</span><span class="va">res_null</span><span class="op">$</span><span class="va">mu</span>, prob <span class="op">=</span> <span class="cn">TRUE</span>, col <span class="op">=</span> <span class="st">"grey"</span>,
     border <span class="op">=</span> <span class="st">"white"</span>,
     main <span class="op">=</span> <span class="st">"Distribution of mu"</span><span class="op">)</span>

<span class="fu"><a href="https://rdrr.io/r/base/invisible.html">invisible</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/lapply.html">replicate</a></span><span class="op">(</span><span class="fl">30</span>,
     <span class="fu"><a href="https://rdrr.io/r/graphics/points.html">points</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/stats/density.html">density</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/stats/Exponential.html">rexp</a></span><span class="op">(</span><span class="fl">500</span>, <span class="fl">1000</span><span class="op">)</span><span class="op">)</span>, type <span class="op">=</span> <span class="st">"l"</span>, col <span class="op">=</span> <span class="st">"blue"</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_null_priors-1.png" width="768"></p>
<div class="sourceCode" id="cb23"><pre class="downlit sourceCode r">
<code class="sourceCode R">

<span class="fu"><a href="https://rdrr.io/r/graphics/hist.html">hist</a></span><span class="op">(</span><span class="va">res_null</span><span class="op">$</span><span class="va">pi</span>, prob <span class="op">=</span> <span class="cn">TRUE</span>, col <span class="op">=</span> <span class="st">"grey"</span>,
     border <span class="op">=</span> <span class="st">"white"</span>, main <span class="op">=</span> <span class="st">"Distribution of pi"</span><span class="op">)</span>

<span class="fu"><a href="https://rdrr.io/r/base/invisible.html">invisible</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/lapply.html">replicate</a></span><span class="op">(</span><span class="fl">30</span>,
     <span class="fu"><a href="https://rdrr.io/r/graphics/points.html">points</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/stats/density.html">density</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/stats/Beta.html">rbeta</a></span><span class="op">(</span><span class="fl">500</span>, <span class="fl">10</span>, <span class="fl">1</span><span class="op">)</span><span class="op">)</span>, type <span class="op">=</span> <span class="st">"l"</span>, col <span class="op">=</span> <span class="st">"blue"</span><span class="op">)</span><span class="op">)</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_null_priors-2.png" width="768"></p>
<p><br></p>
</div>
<div id="a-model-using-symptom-onset-only" class="section level2">
<h2 class="hasAnchor">
<a href="#a-model-using-symptom-onset-only" class="anchor"></a>A model using symptom onset only</h2>
<p>We can use data and likelihood customisation to change the default <em>outbreaker2</em> model into a <a href="http://dx.doi.org/10.1093/aje/kwh255">Wallinga &amp; Teunis (1994)</a> model. To do so, we need to:</p>
<ul>
<li><p>Remove the DNA sequences from the data; alternatively we could also specify a ‘null’ function (i.e. returning a finite constant, as above) for the genetic likelihood.</p></li>
<li><p>Disable all likelihood components other than <code>timing_infections</code> using <code>custom_likelihoods</code>. This means that the dates provided will be treated as dates of symptom onset, and the timing distribution <code>w</code> will be taken as the serial interval.</p></li>
<li>
<p>Disable the detection of imported cases, and forcing all <code>kappa</code> values to be</p>
<ol style="list-style-type: decimal">
<li>
</ol>
</li>
</ul>
<p>While these are fairly major changes, they are straightforward to implement. We first create the dataset and custom likelihood functions:</p>
<div class="sourceCode" id="cb24"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">onset_data</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/outbreaker_data.html">outbreaker_data</a></span><span class="op">(</span>dates <span class="op">=</span> <span class="va">fake_outbreak</span><span class="op">$</span><span class="va">onset</span>,
                          w_dens <span class="op">=</span> <span class="va">fake_outbreak</span><span class="op">$</span><span class="va">w</span><span class="op">)</span>

<span class="va">wt_model</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/custom_likelihoods.html">custom_likelihoods</a></span><span class="op">(</span>timing_sampling <span class="op">=</span> <span class="va">f_null</span>,
                               reporting <span class="op">=</span> <span class="va">f_null</span><span class="op">)</span></code></pre></div>
<p>To fix parameters or augmented data (here, fix all <code>kappa</code> values to 1), we set the initial states to the desired values and disable the corresponding moves:</p>
<div class="sourceCode" id="cb25"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">wt_config</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/create_config.html">create_config</a></span><span class="op">(</span>init_kappa <span class="op">=</span> <span class="fl">1</span>, move_kappa <span class="op">=</span> <span class="cn">FALSE</span>,
                           init_pi <span class="op">=</span> <span class="fl">1</span>, move_pi <span class="op">=</span> <span class="cn">FALSE</span>,
                           move_mu <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span></code></pre></div>
<p>We can now run the analyses for this new model:</p>
<div class="sourceCode" id="cb26"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span><span class="op">(</span><span class="fl">1</span><span class="op">)</span>
<span class="va">res_wt</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/outbreaker.html">outbreaker</a></span><span class="op">(</span>data <span class="op">=</span> <span class="va">onset_data</span>,
                     config <span class="op">=</span> <span class="va">wt_config</span>,
                     likelihoods <span class="op">=</span> <span class="va">wt_model</span><span class="op">)</span>
<span class="co">#&gt; Can't use seqTrack initialization with missing DNA sequences; using a star-like tree</span></code></pre></div>
<div class="sourceCode" id="cb27"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_wt</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_wt-1.png" width="768"></p>
<div class="sourceCode" id="cb28"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_wt</span>, burnin <span class="op">=</span> <span class="fl">500</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_wt-2.png" width="768"></p>
<div class="sourceCode" id="cb29"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_wt</span>, burnin <span class="op">=</span> <span class="fl">500</span>, type <span class="op">=</span> <span class="st">"alpha"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_wt-3.png" width="768"></p>
<div class="sourceCode" id="cb30"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span><span class="op">(</span><span class="va">res_wt</span><span class="op">)</span>
<span class="co">#&gt; $step</span>
<span class="co">#&gt;    first     last interval  n_steps </span>
<span class="co">#&gt;        1    10000       50      201 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $post</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt; -354.93  -35.56  -33.78  -35.23  -31.48  -26.79 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $like</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt; -357.23  -37.86  -36.08  -37.53  -33.78  -29.10 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $prior</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;   2.302   2.302   2.302   2.302   2.302   2.302 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $mu</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;   1e-04   1e-04   1e-04   1e-04   1e-04   1e-04 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $pi</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;       1       1       1       1       1       1 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $tree</span>
<span class="co">#&gt;    from to time    support generations</span>
<span class="co">#&gt; 1     4  1    1 0.05472637           1</span>
<span class="co">#&gt; 2    NA  2    1 0.05472637          NA</span>
<span class="co">#&gt; 3    25  3    1 0.06467662           1</span>
<span class="co">#&gt; 4    21  4    0 0.07462687           1</span>
<span class="co">#&gt; 5    26  5    0 0.08457711           1</span>
<span class="co">#&gt; 6    19  6   -1 0.06965174           1</span>
<span class="co">#&gt; 7    15  7    1 0.06467662           1</span>
<span class="co">#&gt; 8    11  8    0 0.06467662           1</span>
<span class="co">#&gt; 9    15  9   -1 0.04975124           1</span>
<span class="co">#&gt; 10   12 10   -1 0.05472637           1</span>
<span class="co">#&gt; 11   16 11    0 0.05970149           1</span>
<span class="co">#&gt; 12    2 12    0 0.07462687           1</span>
<span class="co">#&gt; 13   24 13    0 0.05970149           1</span>
<span class="co">#&gt; 14   10 14    1 0.06467662           1</span>
<span class="co">#&gt; 15    9 15    1 0.07462687           1</span>
<span class="co">#&gt; 16   22 16    0 0.05472637           1</span>
<span class="co">#&gt; 17    2 17    0 0.05472637           1</span>
<span class="co">#&gt; 18    1 18    1 0.05472637           1</span>
<span class="co">#&gt; 19   24 19    0 0.06965174           1</span>
<span class="co">#&gt; 20   28 20    0 0.05970149           1</span>
<span class="co">#&gt; 21   12 21   -1 0.05472637           1</span>
<span class="co">#&gt; 22   10 22   -1 0.05970149           1</span>
<span class="co">#&gt; 23    5 23    1 0.06965174           1</span>
<span class="co">#&gt; 24   11 24    0 0.05970149           1</span>
<span class="co">#&gt; 25   10 25    0 0.05472637           1</span>
<span class="co">#&gt; 26   22 26    2 0.07462687           1</span>
<span class="co">#&gt; 27    7 27    0 0.05970149           1</span>
<span class="co">#&gt; 28   17 28    2 0.05970149           1</span>
<span class="co">#&gt; 29    2 29    0 0.06965174           1</span>
<span class="co">#&gt; 30   22 30    2 0.07462687           1</span></code></pre></div>
<p>As before for the ‘null’ model, the transmission tree is very poorly resolved in this case. We use the same approach to visualise it: extract nodes and edges from the <code>visNetork</code> object, use this information to create an <code>igraph</code> object, and visualise the result using a circular layout:</p>
<div class="sourceCode" id="cb31"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="co">## extract nodes and edges from the visNetwork object</span>
<span class="va">temp</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_wt</span>, type <span class="op">=</span> <span class="st">"network"</span>, min_support <span class="op">=</span> <span class="fl">0.05</span><span class="op">)</span>
<span class="fu"><a href="https://rdrr.io/r/base/class.html">class</a></span><span class="op">(</span><span class="va">temp</span><span class="op">)</span>
<span class="co">#&gt; [1] "visNetwork" "htmlwidget"</span>
<span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span><span class="op">(</span><span class="va">temp</span><span class="op">$</span><span class="va">x</span><span class="op">$</span><span class="va">edges</span><span class="op">)</span>
<span class="co">#&gt;    from to      value arrows   color</span>
<span class="co">#&gt; 10    1 11 0.05472637     to #CCDDFF</span>
<span class="co">#&gt; 17    1 18 0.05472637     to #CCDDFF</span>
<span class="co">#&gt; 23    1 24 0.05472637     to #CCDDFF</span>
<span class="co">#&gt; 32    2  4 0.06965174     to #B2D9E3</span>
<span class="co">#&gt; 39    2 12 0.07462687     to #B2D9E3</span>
<span class="co">#&gt; 42    2 15 0.05970149     to #B2D9E3</span>
<span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span><span class="op">(</span><span class="va">temp</span><span class="op">$</span><span class="va">x</span><span class="op">$</span><span class="va">nodes</span><span class="op">)</span>
<span class="co">#&gt;   id label     value   color shape shaped</span>
<span class="co">#&gt; 1  1     1 0.9751244 #CCDDFF   dot   &lt;NA&gt;</span>
<span class="co">#&gt; 2  2     2 1.0696517 #B2D9E3   dot   star</span>
<span class="co">#&gt; 3  3     3 0.9402985 #98D6C7   dot   &lt;NA&gt;</span>
<span class="co">#&gt; 4  4     4 1.0845771 #7ED2AC   dot   &lt;NA&gt;</span>
<span class="co">#&gt; 5  5     5 0.9353234 #99CAA9   dot   &lt;NA&gt;</span>
<span class="co">#&gt; 6  6     6 0.9303483 #C2C0AD   dot   &lt;NA&gt;</span>

<span class="co">## make an igraph object</span>

<span class="va">net_wt</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/pkg/igraph/man/graph_from_data_frame.html">graph.data.frame</a></span><span class="op">(</span><span class="va">temp</span><span class="op">$</span><span class="va">x</span><span class="op">$</span><span class="va">edges</span>,
                             vertices <span class="op">=</span> <span class="va">temp</span><span class="op">$</span><span class="va">x</span><span class="op">$</span><span class="va">nodes</span><span class="op">[</span><span class="fl">1</span><span class="op">:</span><span class="fl">4</span><span class="op">]</span><span class="op">)</span>
                 
<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">net_wt</span>, layout <span class="op">=</span> <span class="va">layout.circle</span>,
     main <span class="op">=</span> <span class="st">"WT model, posterior trees"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/wt_net-1.png" width="768"></p>
<p><br></p>
</div>
</div>
<div id="customising-movements" class="section level1">
<h1 class="hasAnchor">
<a href="#customising-movements" class="anchor"></a>Customising movements</h1>
<p>Customising movements works in similar ways to priors and likelihoods. In practice, this type of customisation is more complex as it most likely will require evaluation of likelihoods and priors, and therefore require the user to know which functions to all, and how. These are documented in the <a href="Rcpp_API.html">API vignette</a>. In the following, we provide two examples:</p>
<ul>
<li><p>a (fake) Gibbs sampler for the movement of the mutation rate <code>mu</code></p></li>
<li><p>a new ‘naive’ movement of ancestries in which infectors are picked at random from all cases</p></li>
</ul>
<p>But before getting into these, we need to explicit how movements are happening in <em>outbreaker2</em>.</p>
<div id="movements-in-outbreaker2" class="section level2">
<h2 class="hasAnchor">
<a href="#movements-in-outbreaker2" class="anchor"></a>Movements in <em>outbreaker2</em>
</h2>
<p>At the core of the <code>outbreaker</code> function, movements are implemented as a list of functions, which are all evaluated in turn during every iteration of the MCMC. All movement functions must obey two rules:</p>
<ul>
<li><p>The first argument must be an <code>outbreaker_param</code> object (typically called <code>param</code> in the original code); see <code><a href="../reference/create_param.html">?create_param</a></code> for details.</p></li>
<li><p>All movement functions must return a valid, <code>outbreaker_param</code> object.</p></li>
</ul>
<p>However, a new difficulty compared to prior or likelihood customisation is that different movements may require different components of the model, and a different set of arguments after <code>param</code>. In fact, this can be seen by examining the arguments of all the default movement functions:</p>
<div class="sourceCode" id="cb32"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/base/lapply.html">lapply</a></span><span class="op">(</span><span class="fu"><a href="../reference/custom_moves.html">custom_moves</a></span><span class="op">(</span><span class="op">)</span>,  <span class="va">args</span><span class="op">)</span>
<span class="co">#&gt; $mu</span>
<span class="co">#&gt; function (param, data, config, custom_ll = NULL, custom_prior = NULL) </span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $pi</span>
<span class="co">#&gt; function (param, data, config, custom_ll = NULL, custom_prior = NULL) </span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $eps</span>
<span class="co">#&gt; function (param, data, config, custom_ll = NULL, custom_prior = NULL) </span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $lambda</span>
<span class="co">#&gt; function (param, data, config, custom_ll = NULL, custom_prior = NULL) </span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $alpha</span>
<span class="co">#&gt; function (param, data, list_custom_ll = NULL) </span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $swap_cases</span>
<span class="co">#&gt; function (param, data, list_custom_ll = NULL) </span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $t_inf</span>
<span class="co">#&gt; function (param, data, list_custom_ll = NULL) </span>
<span class="co">#&gt; NULL</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $kappa</span>
<span class="co">#&gt; function (param, data, config, list_custom_ll = NULL) </span>
<span class="co">#&gt; NULL</span></code></pre></div>
<p>To handle this difficulty, <em>outbreaker2</em> transforms every movement function before running the MCMC into a new function with a single parameter <code>param</code>, attaching all other required argument to the function’s environment. The function achieving this transformation is called <code>bind_moves</code>. This function ‘knows’ what these components are for known moves listed aboves. For new, unknown moves, it attaches the following components, provided they are used as arguments of the new function:</p>
<ul>
<li><p><code>data</code>: the processed data; see <code><a href="../reference/outbreaker_data.html">?outbreaker_data</a></code></p></li>
<li><p><code>config</code>: the configuration list; see <code>create_config</code></p></li>
<li><p><code>likelihoods</code>: a list of custom likelihood functions; see <code><a href="../reference/custom_likelihoods.html">?custom_likelihoods</a></code></p></li>
<li><p><code>priors</code>: a list of custom prior functions; see <code><a href="../reference/custom_priors.html">?custom_priors</a></code></p></li>
</ul>
<p>See examples in <code>?bind_moves</code> for details of how this process works.</p>
</div>
<div id="a-fake-gibbs-sampler-for-mu" class="section level2">
<h2 class="hasAnchor">
<a href="#a-fake-gibbs-sampler-for-mu" class="anchor"></a>A (fake) Gibbs sampler for <code>mu</code>
</h2>
<p>A Gibbs sampler supposes that the conditional distribution of a parameter is known and can directly be sampled from. Here, we use this principle to provide a toy example of custom movement for <code>mu</code>, assuming that this conditional distribution is always an Exponential distribution with a rate of 1000. This is easy to implement; to make sure that the function is actually used, we set a global variable changed when the function is called.</p>
<div class="sourceCode" id="cb33"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">move_mu</span> <span class="op">&lt;-</span> <span class="kw">function</span><span class="op">(</span><span class="va">param</span>, <span class="va">config</span><span class="op">)</span> <span class="op">{</span>

  <span class="va">NEW_MOVE_HAS_BEEN_USED</span> <span class="op">&lt;&lt;-</span> <span class="cn">TRUE</span>
  <span class="va">param</span><span class="op">$</span><span class="va">mu</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/stats/Exponential.html">rexp</a></span><span class="op">(</span><span class="fl">1</span>, <span class="fl">1000</span><span class="op">)</span>  
  <span class="kw"><a href="https://rdrr.io/r/base/function.html">return</a></span><span class="op">(</span><span class="va">param</span><span class="op">)</span>
  
<span class="op">}</span>

<span class="va">moves</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/custom_moves.html">custom_moves</a></span><span class="op">(</span>mu <span class="op">=</span> <span class="va">move_mu</span><span class="op">)</span>

<span class="va">quick_config</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></span><span class="op">(</span>n_iter <span class="op">=</span> <span class="fl">500</span>, sample_every <span class="op">=</span> <span class="fl">1</span>, find_import <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span></code></pre></div>
<p>Note that the new movement function <code>move_mu</code> has two arguments, and that we do not specify <code>config</code>. Internally, what happens is:</p>
<div class="sourceCode" id="cb34"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="co">## bind quick_config to function</span>
<span class="va">move_mu_intern</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/bind_to_function.html">bind_to_function</a></span><span class="op">(</span><span class="va">move_mu</span>, config <span class="op">=</span> <span class="va">quick_config</span><span class="op">)</span>

<span class="co">## new function has just one argument</span>
<span class="va">move_mu_intern</span>
<span class="co">#&gt; function (param) </span>
<span class="co">#&gt; {</span>
<span class="co">#&gt;     NEW_MOVE_HAS_BEEN_USED &lt;&lt;- TRUE</span>
<span class="co">#&gt;     param$mu &lt;- rexp(1, 1000)</span>
<span class="co">#&gt;     return(param)</span>
<span class="co">#&gt; }</span>
<span class="co">#&gt; &lt;environment: 0x0000000022775760&gt;</span>

<span class="co">## 'config' is in the function's environment</span>
<span class="fu"><a href="https://rdrr.io/r/base/names.html">names</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/environment.html">environment</a></span><span class="op">(</span><span class="va">move_mu_intern</span><span class="op">)</span><span class="op">)</span>
<span class="co">#&gt; [1] "config"</span>

<span class="co">## 'config' is actually 'quick_config'</span>
<span class="fu"><a href="https://rdrr.io/r/base/identical.html">identical</a></span><span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/environment.html">environment</a></span><span class="op">(</span><span class="va">move_mu_intern</span><span class="op">)</span><span class="op">$</span><span class="va">config</span>, <span class="va">quick_config</span><span class="op">)</span>
<span class="co">#&gt; [1] TRUE</span></code></pre></div>
<p>We perform a quick run using this new movement:</p>
<div class="sourceCode" id="cb35"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">NEW_MOVE_HAS_BEEN_USED</span> <span class="op">&lt;-</span> <span class="cn">FALSE</span>

<span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span><span class="op">(</span><span class="fl">1</span><span class="op">)</span>
<span class="va">res_move_mu</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/outbreaker.html">outbreaker</a></span><span class="op">(</span><span class="va">data</span>, <span class="va">quick_config</span>, moves <span class="op">=</span> <span class="va">moves</span><span class="op">)</span>
<span class="va">NEW_MOVE_HAS_BEEN_USED</span>
<span class="co">#&gt; [1] TRUE</span>

<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_move_mu</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/run_custom_move_mu-1.png" width="768"></p>
<div class="sourceCode" id="cb36"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_move_mu</span>, <span class="st">"pi"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/run_custom_move_mu-2.png" width="768"></p>
<div class="sourceCode" id="cb37"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_move_mu</span>, <span class="st">"mu"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/run_custom_move_mu-3.png" width="768"></p>
<p>This short, full trace, clearly hasn’t mixed well (as is to be expected). But while we see the effect of accept/reject sampling (Metropolis algorithm) for <code>pi</code> with a lot of autocorrelation, the trace of <code>mu</code> shows complete independence between successive values. While the Gibbs sampler used here is not correct, this result is: a Gibbs sampler will be more efficient than the classical Metropolis(-Hasting) algorithm for a given number a iterations.</p>
<p><br></p>
</div>
<div id="a-new-movement-of-ancestries" class="section level2">
<h2 class="hasAnchor">
<a href="#a-new-movement-of-ancestries" class="anchor"></a>A new movement of ancestries</h2>
<p>Moves of ancestries are done in two ways in outbreaker: by picking ancestors at random from any prior case, and by swapping cases from a transmission link. Here, we implement a new move, which will propose infectors which have been infected on the same day of the current infector. As before, we will use global variables to keep track of the resulting movements (see <code>N_ACCEPT</code> and <code>N_REJECT</code>).</p>
<div class="sourceCode" id="cb38"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">api</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/get_cpp_api.html">get_cpp_api</a></span><span class="op">(</span><span class="op">)</span>

<span class="va">new_move_ances</span> <span class="op">&lt;-</span> <span class="kw">function</span><span class="op">(</span><span class="va">param</span>, <span class="va">data</span>, <span class="va">custom_likelihoods</span> <span class="op">=</span> <span class="cn">NULL</span><span class="op">)</span> <span class="op">{</span>

<span class="kw">for</span> <span class="op">(</span><span class="va">i</span> <span class="kw">in</span> <span class="fl">1</span><span class="op">:</span><span class="va">data</span><span class="op">$</span><span class="va">N</span><span class="op">)</span> <span class="op">{</span>
    <span class="va">current_ances</span> <span class="op">&lt;-</span> <span class="va">param</span><span class="op">$</span><span class="va">alpha</span><span class="op">[</span><span class="va">i</span><span class="op">]</span>
    <span class="kw">if</span> <span class="op">(</span><span class="op">!</span><span class="fu"><a href="https://rdrr.io/r/base/NA.html">is.na</a></span><span class="op">(</span><span class="va">current_ances</span><span class="op">)</span><span class="op">)</span> <span class="op">{</span>
      <span class="co">## find cases infected on the same days</span>
      <span class="va">current_t_inf</span> <span class="op">&lt;-</span> <span class="va">param</span><span class="op">$</span><span class="va">t_inf</span><span class="op">[</span><span class="va">current_ances</span><span class="op">]</span>
      <span class="va">pool</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/which.html">which</a></span><span class="op">(</span><span class="va">param</span><span class="op">$</span><span class="va">t_inf</span> <span class="op">==</span> <span class="va">current_t_inf</span><span class="op">)</span>
      <span class="va">pool</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/sets.html">setdiff</a></span><span class="op">(</span><span class="va">pool</span>, <span class="va">i</span><span class="op">)</span>
      
      <span class="kw">if</span> <span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/base/length.html">length</a></span><span class="op">(</span><span class="va">pool</span><span class="op">)</span> <span class="op">&gt;</span> <span class="fl">0</span><span class="op">)</span> <span class="op">{</span>
        <span class="co">## propose new ancestor</span>
        <span class="va">current_ll</span> <span class="op">&lt;-</span> <span class="va">api</span><span class="op">$</span><span class="fu">cpp_ll_all</span><span class="op">(</span><span class="va">data</span>, <span class="va">param</span>, i <span class="op">=</span> <span class="va">i</span>, <span class="va">custom_likelihoods</span><span class="op">)</span>
        <span class="va">param</span><span class="op">$</span><span class="va">alpha</span><span class="op">[</span><span class="va">i</span><span class="op">]</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/sample.html">sample</a></span><span class="op">(</span><span class="va">pool</span>, <span class="fl">1</span><span class="op">)</span>
        <span class="va">new_ll</span> <span class="op">&lt;-</span> <span class="va">api</span><span class="op">$</span><span class="fu">cpp_ll_all</span><span class="op">(</span><span class="va">data</span>, <span class="va">param</span>, i <span class="op">=</span> <span class="va">i</span>, <span class="va">custom_likelihoods</span><span class="op">)</span>

        <span class="co">## likelihood ratio - no correction, move is symmetric</span>
        <span class="va">ratio</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Log.html">exp</a></span><span class="op">(</span><span class="va">new_ll</span> <span class="op">-</span> <span class="va">current_ll</span><span class="op">)</span>

        <span class="co">## accept / reject</span>
        <span class="kw">if</span> <span class="op">(</span><span class="fu"><a href="https://rdrr.io/r/stats/Uniform.html">runif</a></span><span class="op">(</span><span class="fl">1</span><span class="op">)</span> <span class="op">&lt;=</span> <span class="va">ratio</span><span class="op">)</span> <span class="op">{</span> <span class="co"># accept</span>
          <span class="va">N_ACCEPT</span> <span class="op">&lt;&lt;-</span> <span class="va">N_ACCEPT</span> <span class="op">+</span> <span class="fl">1</span>
        <span class="op">}</span> <span class="kw">else</span> <span class="op">{</span> <span class="co"># reject</span>
          <span class="va">N_REJECT</span> <span class="op">&lt;&lt;-</span> <span class="va">N_REJECT</span> <span class="op">+</span> <span class="fl">1</span>
          <span class="va">param</span><span class="op">$</span><span class="va">alpha</span><span class="op">[</span><span class="va">i</span><span class="op">]</span> <span class="op">&lt;-</span> <span class="va">current_ances</span>
        <span class="op">}</span>
      <span class="op">}</span>
    <span class="op">}</span>
  <span class="op">}</span>
  <span class="kw"><a href="https://rdrr.io/r/base/function.html">return</a></span><span class="op">(</span><span class="va">param</span><span class="op">)</span>
<span class="op">}</span>

<span class="va">moves</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/custom_moves.html">custom_moves</a></span><span class="op">(</span>new_move <span class="op">=</span> <span class="va">new_move_ances</span><span class="op">)</span></code></pre></div>
<p>We can now use this new move in our transmission tree reconstruction. We will use a shorter chain than the defaults as this new move is likely to be computer intensive.</p>
<div class="sourceCode" id="cb39"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="va">N_ACCEPT</span> <span class="op">&lt;-</span> <span class="fl">0</span>
<span class="va">N_REJECT</span> <span class="op">&lt;-</span> <span class="fl">0</span>

<span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span><span class="op">(</span><span class="fl">1</span><span class="op">)</span>
<span class="va">res_new_move</span> <span class="op">&lt;-</span> <span class="fu"><a href="../reference/outbreaker.html">outbreaker</a></span><span class="op">(</span><span class="va">data</span>, <span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></span><span class="op">(</span>move_kappa <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span>, moves <span class="op">=</span> <span class="va">moves</span><span class="op">)</span>

<span class="va">N_ACCEPT</span>
<span class="co">#&gt; [1] 150958</span>
<span class="va">N_REJECT</span>
<span class="co">#&gt; [1] 263986</span></code></pre></div>
<div class="sourceCode" id="cb40"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_new_move</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_new_move-1.png" width="768"></p>
<div class="sourceCode" id="cb41"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/graphics/plot.default.html">plot</a></span><span class="op">(</span><span class="va">res_new_move</span>, type <span class="op">=</span> <span class="st">"alpha"</span><span class="op">)</span></code></pre></div>
<p><img src="figs-customisation/res_new_move-2.png" width="768"></p>
<div class="sourceCode" id="cb42"><pre class="downlit sourceCode r">
<code class="sourceCode R"><span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span><span class="op">(</span><span class="va">res_new_move</span><span class="op">)</span>
<span class="co">#&gt; $step</span>
<span class="co">#&gt;    first     last interval  n_steps </span>
<span class="co">#&gt;        1    10000       50      201 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $post</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt; -1087.1  -447.3  -444.8  -448.5  -443.0  -436.8 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $like</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt; -1088.4  -449.5  -446.9  -450.6  -445.2  -439.0 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $prior</span>
<span class="co">#&gt;     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. </span>
<span class="co">#&gt; -0.01622  1.99783  2.14793  2.05799  2.23686  2.30230 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $mu</span>
<span class="co">#&gt;      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. </span>
<span class="co">#&gt; 8.573e-05 1.282e-04 1.421e-04 1.453e-04 1.619e-04 2.329e-04 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $pi</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;  0.7729  0.9667  0.9830  0.9737  0.9927  1.0000 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $tree</span>
<span class="co">#&gt;    from to time     support generations</span>
<span class="co">#&gt; 1    NA  1   -1 0.004975124          NA</span>
<span class="co">#&gt; 2     1  2    1 0.995024876           1</span>
<span class="co">#&gt; 3     2  3    3 0.995024876           1</span>
<span class="co">#&gt; 4    NA  4    3 0.014925373          NA</span>
<span class="co">#&gt; 5     3  5    4 0.985074627           1</span>
<span class="co">#&gt; 6     4  6    5 0.965174129           1</span>
<span class="co">#&gt; 7     4  7    5 0.985074627           1</span>
<span class="co">#&gt; 8     5  8    6 0.960199005           1</span>
<span class="co">#&gt; 9    13  9    7 0.995024876           1</span>
<span class="co">#&gt; 10    6 10    7 0.995024876           1</span>
<span class="co">#&gt; 11    7 11    7 0.681592040           1</span>
<span class="co">#&gt; 12    5 12    7 0.850746269           1</span>
<span class="co">#&gt; 13    6 13    6 0.995024876           1</span>
<span class="co">#&gt; 14    5 14    7 0.761194030           1</span>
<span class="co">#&gt; 15    5 15    7 0.800995025           1</span>
<span class="co">#&gt; 16    7 16    8 0.805970149           1</span>
<span class="co">#&gt; 17    7 17    7 0.626865672           1</span>
<span class="co">#&gt; 18    8 18    9 0.467661692           1</span>
<span class="co">#&gt; 19    9 19    9 1.000000000           1</span>
<span class="co">#&gt; 20   10 20   10 0.985074627           1</span>
<span class="co">#&gt; 21   11 21   10 0.970149254           1</span>
<span class="co">#&gt; 22   11 22   10 1.000000000           1</span>
<span class="co">#&gt; 23   13 23    9 1.000000000           1</span>
<span class="co">#&gt; 24   13 24    9 1.000000000           1</span>
<span class="co">#&gt; 25   13 25    9 1.000000000           1</span>
<span class="co">#&gt; 26   17 26    9 0.990049751           1</span>
<span class="co">#&gt; 27   17 27   10 1.000000000           1</span>
<span class="co">#&gt; 28   NA 28    9          NA          NA</span>
<span class="co">#&gt; 29   10 29   10 1.000000000           1</span>
<span class="co">#&gt; 30   13 30   10 1.000000000           1</span></code></pre></div>
<p>Results show a switch to a new mode at about 5000 iterations. Let us compare the consensus tree to the actual one (store in <code>fake_outbreak$ances</code>):</p>
<div class="sourceCode" id="cb43"><pre class="downlit sourceCode r">
<code class="sourceCode R">
<span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span><span class="op">(</span><span class="va">res_new_move</span>, burnin <span class="op">=</span> <span class="fl">5000</span><span class="op">)</span>
<span class="co">#&gt; $step</span>
<span class="co">#&gt;    first     last interval  n_steps </span>
<span class="co">#&gt;     5050    10000       50      100 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $post</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;  -452.5  -447.0  -444.5  -445.0  -443.0  -437.9 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $like</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;  -453.9  -449.2  -446.5  -447.1  -445.2  -440.2 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $prior</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;   1.212   2.020   2.170   2.086   2.238   2.302 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $mu</span>
<span class="co">#&gt;      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. </span>
<span class="co">#&gt; 0.0001006 0.0001281 0.0001416 0.0001450 0.0001620 0.0002308 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $pi</span>
<span class="co">#&gt;    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. </span>
<span class="co">#&gt;  0.8859  0.9691  0.9854  0.9765  0.9929  1.0000 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; $tree</span>
<span class="co">#&gt;    from to time support generations</span>
<span class="co">#&gt; 1    NA  1   -1      NA          NA</span>
<span class="co">#&gt; 2     1  2    1    1.00           1</span>
<span class="co">#&gt; 3     2  3    3    0.99           1</span>
<span class="co">#&gt; 4    NA  4    3    0.01          NA</span>
<span class="co">#&gt; 5     3  5    4    0.99           1</span>
<span class="co">#&gt; 6     4  6    5    0.96           1</span>
<span class="co">#&gt; 7     4  7    5    0.99           1</span>
<span class="co">#&gt; 8     5  8    6    0.95           1</span>
<span class="co">#&gt; 9    13  9    7    1.00           1</span>
<span class="co">#&gt; 10    6 10    7    1.00           1</span>
<span class="co">#&gt; 11    7 11    7    0.70           1</span>
<span class="co">#&gt; 12    5 12    7    0.88           1</span>
<span class="co">#&gt; 13    6 13    6    1.00           1</span>
<span class="co">#&gt; 14    5 14    7    0.79           1</span>
<span class="co">#&gt; 15    5 15    7    0.81           1</span>
<span class="co">#&gt; 16    7 16    8    0.85           1</span>
<span class="co">#&gt; 17    7 17    7    0.68           1</span>
<span class="co">#&gt; 18    8 18    9    0.47           1</span>
<span class="co">#&gt; 19    9 19    9    1.00           1</span>
<span class="co">#&gt; 20   10 20   10    0.97           1</span>
<span class="co">#&gt; 21   11 21   10    0.98           1</span>
<span class="co">#&gt; 22   11 22   10    1.00           1</span>
<span class="co">#&gt; 23   13 23    9    1.00           1</span>
<span class="co">#&gt; 24   13 24    9    1.00           1</span>
<span class="co">#&gt; 25   13 25    8    1.00           1</span>
<span class="co">#&gt; 26   17 26    9    1.00           1</span>
<span class="co">#&gt; 27   17 27   10    1.00           1</span>
<span class="co">#&gt; 28   NA 28    9      NA          NA</span>
<span class="co">#&gt; 29   10 29   10    1.00           1</span>
<span class="co">#&gt; 30   13 30   10    1.00           1</span>
<span class="va">tree2</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span><span class="op">(</span><span class="va">res_new_move</span>, burnin <span class="op">=</span> <span class="fl">5000</span><span class="op">)</span><span class="op">$</span><span class="va">tree</span>

<span class="va">comparison</span> <span class="op">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/data.frame.html">data.frame</a></span><span class="op">(</span>case <span class="op">=</span> <span class="fl">1</span><span class="op">:</span><span class="fl">30</span>,
                         inferred <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/paste.html">paste</a></span><span class="op">(</span><span class="va">tree2</span><span class="op">$</span><span class="va">from</span><span class="op">)</span>,
             true <span class="op">=</span> <span class="fu"><a href="https://rdrr.io/r/base/paste.html">paste</a></span><span class="op">(</span><span class="va">fake_outbreak</span><span class="op">$</span><span class="va">ances</span><span class="op">)</span>,
             stringsAsFactors <span class="op">=</span> <span class="cn">FALSE</span><span class="op">)</span>
             
<span class="va">comparison</span><span class="op">$</span><span class="va">correct</span> <span class="op">&lt;-</span> <span class="va">comparison</span><span class="op">$</span><span class="va">inferred</span> <span class="op">==</span> <span class="va">comparison</span><span class="op">$</span><span class="va">true</span>
<span class="va">comparison</span>
<span class="co">#&gt;    case inferred true correct</span>
<span class="co">#&gt; 1     1       NA   NA    TRUE</span>
<span class="co">#&gt; 2     2        1    1    TRUE</span>
<span class="co">#&gt; 3     3        2    2    TRUE</span>
<span class="co">#&gt; 4     4       NA   NA    TRUE</span>
<span class="co">#&gt; 5     5        3    3    TRUE</span>
<span class="co">#&gt; 6     6        4    4    TRUE</span>
<span class="co">#&gt; 7     7        4    4    TRUE</span>
<span class="co">#&gt; 8     8        5    5    TRUE</span>
<span class="co">#&gt; 9     9       13    6   FALSE</span>
<span class="co">#&gt; 10   10        6    6    TRUE</span>
<span class="co">#&gt; 11   11        7    7    TRUE</span>
<span class="co">#&gt; 12   12        5    8   FALSE</span>
<span class="co">#&gt; 13   13        6    9   FALSE</span>
<span class="co">#&gt; 14   14        5    5    TRUE</span>
<span class="co">#&gt; 15   15        5    5    TRUE</span>
<span class="co">#&gt; 16   16        7    7    TRUE</span>
<span class="co">#&gt; 17   17        7    7    TRUE</span>
<span class="co">#&gt; 18   18        8    8    TRUE</span>
<span class="co">#&gt; 19   19        9    9    TRUE</span>
<span class="co">#&gt; 20   20       10   10    TRUE</span>
<span class="co">#&gt; 21   21       11   11    TRUE</span>
<span class="co">#&gt; 22   22       11   11    TRUE</span>
<span class="co">#&gt; 23   23       13   13    TRUE</span>
<span class="co">#&gt; 24   24       13   13    TRUE</span>
<span class="co">#&gt; 25   25       13   13    TRUE</span>
<span class="co">#&gt; 26   26       17   17    TRUE</span>
<span class="co">#&gt; 27   27       17   17    TRUE</span>
<span class="co">#&gt; 28   28       NA   NA    TRUE</span>
<span class="co">#&gt; 29   29       10   10    TRUE</span>
<span class="co">#&gt; 30   30       13   13    TRUE</span>
<span class="fu"><a href="https://rdrr.io/r/base/mean.html">mean</a></span><span class="op">(</span><span class="va">comparison</span><span class="op">$</span><span class="va">correct</span><span class="op">)</span>
<span class="co">#&gt; [1] 0.9</span></code></pre></div>
</div>
</div>
  </div>

  <div class="col-md-3 hidden-xs hidden-sm" id="pkgdown-sidebar">

        <nav id="toc" data-toggle="toc"><h2 data-toc-skip>Contents</h2>
    </nav>
</div>

</div>



      <footer><div class="copyright">
  <p>Developed by Thibaut Jombart, Finlay Campbell, Rich Fitzjohn.</p>
</div>

<div class="pkgdown">
  <p>Site built with <a href="https://pkgdown.r-lib.org/">pkgdown</a> 1.6.1.</p>
</div>

      </footer>
</div>

  


  </body>
</html>