revbayes_mcmc_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### `RevBayes` with jupyter\n",
    "This notebook demonstrates how to run a simple `RevBayes` analysis using jupyter. Clicking a cell will allow you to modify its contents. Note that some cells contain `RevBayes` code and others contain `Markdown` code. Pressing '`Shift+Enter`' will execute or render the code within a cell. The prompt on the left hand side reading e.g. '`In [1]:`' indicates the sequence of executed cells (where '`[2]`' is executed *after* '`[1]`').\n",
    "\n",
    "First, we'll create filepath and filename variables for IO."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# IO\n",
    "dat_fp = \"example/data/\"\n",
    "dat_fn = dat_fp + \"primates_cytb.nex\"\n",
    "out_fp = \"example/output/\"\n",
    "out_fn = out_fp + \"primates\"\n",
    "print(\"path to data: \" + dat_fn)\n",
    "print(\"path to output: \" + out_fn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, create helper variables to configure the MCMC analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# MCMC settings\n",
    "mvi = 1\n",
    "mni = 1\n",
    "n_gen = 1e3\n",
    "sample_freq = floor(n_gen/1e2)\n",
    "print(\"n_gen: \" + n_gen)\n",
    "print(\"sample_freq: \" + sample_freq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the NEXUS file stored in `dat_fn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read alignment\n",
    "dat = readDiscreteCharacterData(dat_fn)\n",
    "dat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the alignment's dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data dimensions\n",
    "taxa = dat.taxa()\n",
    "n_sites = dat.nchar()\n",
    "n_taxa = taxa.size()\n",
    "print(\"taxon names:\")\n",
    "print(taxa)\n",
    "print(\"n_sites: \" + n_sites)\n",
    "print(\"n_taxa: \" + n_taxa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's create a simple pure birth process with unit height and birth rate, $\\lambda \\sim  \\text{Exp}(1)$. The three moves---`mvNNI`, `mvFNPR`, and `mvNodeTimeSlideUniform`---will be used to instruct MCMC how to explore, and thus integrate over, tree space."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create unit Yule tree\n",
    "birth ~ dnExp(1)\n",
    "phy ~ dnBDP(lambda=birth, mu=0., rootAge=1, taxa=taxa)\n",
    "\n",
    "# tree moves\n",
    "mv[mvi++] = mvNNI(phy, weight=n_taxa)\n",
    "mv[mvi++] = mvFNPR(phy, weight=n_taxa/2)\n",
    "mv[mvi++] = mvNodeTimeSlideUniform(phy, weight=n_taxa)\n",
    "\n",
    "# print the tree's value\n",
    "phy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll assume characters evolve according to a generalized time-reversible substitution process, with a rate matrix, `Q`, determined by the base frequencies, `pi`, and the exchangeability rates, `er`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# base frequencies\n",
    "pi ~ dnDirichlet(rep(1,4))\n",
    "mv[mvi++] = mvSimplexElementScale(pi, alpha=10, weight=4)\n",
    "\n",
    "# exchangeability rates\n",
    "er ~ dnDirichlet(rep(1,6))\n",
    "mv[mvi++] = mvSimplexElementScale(er, alpha=10, weight=6)\n",
    "\n",
    "# GTR rate matrix\n",
    "Q := fnGTR(exchangeRates=er,\n",
    "           baseFrequencies=pi)\n",
    "    \n",
    "# print the matrix's value\n",
    "Q"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we'll create our phylogenetic substitution model, `seq`, conditioned on the sequence data, `dat`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# build phylogenetic CTMC\n",
    "seq ~ dnPhyloCTMC(tree=phy,\n",
    "                  Q=Q,\n",
    "                  branchRates=1.,\n",
    "                  nSites=n_sites,\n",
    "                  type=\"DNA\")\n",
    "\n",
    "# treat the simulated data as 'observed'\n",
    "seq.clamp(dat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample posterior tree and parameter estimates from MCMC every `sample_freq` iterations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create monitors\n",
    "mn[mni++] = mnScreen(pi, printgen=sample_freq)\n",
    "mn[mni++] = mnModel(filename=out_fn+\".model.log\",\n",
    "                    printgen=sample_freq)\n",
    "mn[mni++] = mnFile(phy,\n",
    "                   filename=out_fn+\".tre\",\n",
    "                   printgen=sample_freq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build the `Model` and `Mcmc` objects from the model graph."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# create MCMC\n",
    "mdl = model(phy)\n",
    "ch = mcmc(mdl,mv,mn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the MCMC chain, `ch`, for `n_gen` iterations. Samples from `mnScreen` will be pushed to the jupyter console, but also saved in the `output` directory relative to this notebook's directory.\n",
    "\n",
    "Highlight this final cell and press '`Shift+Enter`' to start the MCMC analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# run MCMC\n",
    "ch.run(n_gen)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "RevBayes",
   "language": "bash",
   "name": "revbayes_kernel"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "r"
   },
   "file_extension": ".Rev",
   "help_links": [
    {
     "text": "RevBayes",
     "url": "https://revbayes.org"
    },
    {
     "text": "RevBayes Kernel",
     "url": "https://github.com/sdwfrost/revbayes_kernel"
    },
    {
     "text": "MetaKernel Magics",
     "url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
    }
   ],
   "mimetype": "text/x-rsrc",
   "name": "RevBayes",
   "pygments_lexer": "R"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}