3_interpreting_literate_results_final.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "3_interpreting_literate_results_final.ipynb",
      "private_outputs": true,
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/dsilvestro/LiteRate/blob/master/3_interpreting_literate_results_final.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2Xr4GHRM4n5F"
      },
      "source": [
        "# Modeling the Dynamics of Cultural Diversification\n",
        "\n",
        "# 3. Interpreting LiteRate Results\n",
        "\n",
        "In the last tutorial, we introduced the LiteRate algorithm for discovering statistically-significant shifts in diversification rates over time. These rates can indicate the influence of historical events or evolutionary mechanisms on cultural lineages.  \n",
        "\n",
        "We also introduced our example dataset, Metal bands active from 1968 to 2000. The population is interesting in that bands can be understood as cultural lineages representing distinct musical, aesthetic and social ideas, as well as organizations competing with each other in a field for social capital. \n",
        "\n",
        "In this tutorial you will learn:\n",
        "\n",
        "* How to assess convergence of LiteRate MCMC chains\n",
        "* How to plot LiteRate results\n",
        "* How to interpret LiteRate results\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ASKSQJ0r5Jg8"
      },
      "source": [
        "# a. Setting up the Notebook"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "v-euiyJ_6dkd"
      },
      "source": [
        "We will once again clone the LiteRate repository "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "Z4vuculU6vH-",
        "cellView": "both"
      },
      "source": [
        "!git clone https://github.com/dsilvestro/LiteRate\n",
        "%cd 'LiteRate'"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KjUrZpRWLvF0"
      },
      "source": [
        "Now we'll assume that you ran LiteRate using the command:\n",
        "\n",
        "**`!python3 LiteRateForward.py -d ./example_data/metal_bands/metal_bands_1.tsv -s 10000 -p 10000000`**\n",
        "\n",
        "MCMCs achieve parameter estimates by random walking through the parameter space so they are not guaranteed to explore the entire space, even after 10,000,000 steps. **We recommend that you run LiteRate at least five times to ensure that your chains have mixed properly and give confidence in convergence.**\n",
        "\n",
        "---\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6qp2B07u-yc4"
      },
      "source": [
        "# b. MCMC Diagnostics\n",
        "\n",
        "## Setting up Tracer\n",
        "\n",
        "Before we interpret our results, we want to make sure that the MCMC samplers have explored the entire parameter space and reached a stationary distribution. The statistics we suggest looking at below are easy enough to script in Python or R, and there are packages to compute them as well. However, we use a standalone program called Tracer. Tracer can be downloaded [here](https://beast.community/tracer). Because it is java, the program should run without installation.\n",
        "\n",
        "Because we can't run Tracer in the notebook, the following diagnostic steps are just going to be shown with pictures. If you want, you can also [download](https://raw.githubusercontent.com/dsilvestro/LiteRate/master/example_data/metal_bands/single_run/metal_bands_1_mcmc.log) the log file ``metal_bands_1_mcmc.log`` from the Github on to your own computer and run it along with us.\n",
        "\n",
        "## Opening Tracer for the first time\n",
        "\n",
        "Click the \"+\" button to add files to Tracer. You can add multiple logs (e.g. for each chain) to Tracer to compare their results and give confidence in convergence. For now we have just one.\n",
        "\n",
        "There are four panels in Tracer.\n",
        "\n",
        "**In the top left** we have the list of logs, the number of states (i.e. steps) that the sampler explored and the amount of states you want to burn-in. Burn-in removes the first X% of steps in the random walk before the MCMC has (hopefully) reached its stationary distribution. Anywhere in the 10-20% range is fine. Tracer defaults to 10%.\n",
        "\n",
        "**In the bottom left**, you have the list of parameters in the log-file. LiteRate keeps track of the joint-posterior, the likelihood, the joint-prior and several meta parameters. The rates themselves are found in ``metal_bands_1_sp_rates.log`` and ``metal_bands_1_ex_rates.log`` and aren't readable by Tracer. For your interest the parameters listed here are described in the Deep Dive below. To assess convergence, we'll just focus on the posterior.\n",
        "\n",
        "<details>\n",
        "<summary><font size=\"4\" color=\"Indigo\"> Deep Dive: MCMC Log Parameters </font></summary>\n",
        "\n",
        "* **lambda_avg:** The average rate if you had just averaged over all birth rates (i.e. run a single birth-death process). This can be useful as a null hypothesis that there is no change in rates over time. \n",
        "* **mu_avg:** The average rate if you had just averaged over all death rates (i.e. run a single birth-death process). This can be useful as a null hypothesis that there is no change in rates over time.\n",
        "* **K_l:** The number of rate shifts in the birth rates.\n",
        "* **K_m:** The number of rate shifts in the death rates.\n",
        "* **root_age:** The first time unit of analysis.\n",
        "* **death_age:** The last time unit of analysis.\n",
        "* **gamma_rate_hp_BI:** This a hyperprior. Don't worry about this; you can read Silvestro et al. 2018 for details.\n",
        "* **gamma_rate_hp_D:** This a hyperprior. Don't worry about this; you can read Silvestro et al. 2018 for details.\n",
        "* **poisson_rate_hp:** This a hyperprior. Don't worry about this; you can read Silvestro et al. 2018 for details.\n",
        "* **corr_coeff:** For convenience, this is the correlation coefficient between the LiteRate rates and the empirical bin-wise rates. It is computed by concatenating the birth and death rates together before regressing the estimated rates on the empirical rates. \n",
        "* **r_squared:** The above, squared. \n",
        "* **gelman_squared:** This is the Bayesian $R^2$ described in [Gelman et al., 2018](http://www.stat.columbia.edu/~gelman/research/unpublished/bayes_R2.pdf). It should be very similar or identical to the above statistic.  \n",
        "</details>\n",
        "\n",
        "There are also columns for the **mean** and **effective sample size (ESS).** These are reproduced in the summary statistics panel described below. \n",
        "\n",
        "**In the top right,** we have the summary statistics panel. The **mean** averaged across all samples (after burnin) is your parameter estimate of interest. This panel also include 95\\% high-posterior density (HPD) interval. Any point within this interval has a higher density than any point outside of it. The last statistic of interest is the **effective sample size (ESS)**.\n",
        "\n",
        "**In the bottom right,** you have a plot window for which you can see histogram/density plots of your samples for each parameter, as well as **trace plots** to assess auto-correlation.\n",
        "\n",
        "<figure align=\"center\">\n",
        "<img src=\"https://drive.google.com/uc?id=1foc5pk63J60q3qEColw1_TJlOox8a90V\" alt=\"\" width=\"1000\" height=\"500\" border=\"0\">\n",
        "</figure>\n",
        "\n",
        "## Key Diagnostics\n",
        "\n",
        "### ESS\n",
        "In an MCMC chain, parameter estimates at each step are highly correlated because they are close to the previous step. In our example, we only sampled every 1000 steps with hopes that samples at this frequency are relatively uncorrelated. The ESS is a simple statistic that approximates the theoretical number of roughly uncorrelated samples we have in our chain, controlling for the actual correlation between our samples:\n",
        "$$N_{\\mathrm{eff}}= \\frac{N}{1 + 2 \\sum_{t = 1}^{\\infty} \\rho_t}$$\n",
        "where $t$ is the lag between samples and $\\rho_t$ is the autocorrelation between sample $t$ and $t+1$ [(STAN Manual)](https://mc-stan.org/docs/2_20/reference-manual/effective-sample-size-section.html). \n",
        "\n",
        "Intuitively, ESS just discounts samples that are highly correlated with the next sample. You can also see how the ESS is inversely related to the standard deviation of the mean. As a rule of thumb, you want an ESS of at least 200 for most of your parameter estimates. The best way to get a bigger ESS is to run the chain for a longer time.\n",
        "\n",
        "### Trace Plot\n",
        "\n",
        "The default window for each parameter is a histogram of samples. If the MCMC mixes well, this should look vaguely Gaussian for continuous parameters (although it could be skewed).\n",
        "\n",
        "The other nice diagnostic to check out is the \"Trace\" plot. To open up the trace plot click on the tab at the top of the plot window.\n",
        "<figure align=\"center\">\n",
        "<img src=\"https://drive.google.com/uc?id=1BmDMPOY43znUmLpe2YrQjrqT4vV4jA3h\" alt=\"\" width=\"1000\" height=\"500\" border=\"0\">\n",
        "</figure>\n",
        "\n",
        "The trace feature plots the estimates within the sample as a line chart. On the left, note how the burn-in is removed. The fact that the trace jumps up and down suggests that each individual sample is uncorrelated with the previous one. If we see trends in the plot over time, this indicates that the chain has not yet converged on its stationary distribution, and you should run it for longer. This chain looks good!   \n",
        "\n",
        "## Combining Chains\n",
        "\n",
        "You should check each of your chains (ESS and Trace plot) to confirm their convergence. If everything looks good, it makes sense to combine your log files (removing burn-in in each), giving you a bigger ESS and reducing variance in your estimates. You can do this using ``logCombiner.py``, a utility in the LiteRate github in a directory full of files you wish to combine. Alternatively, this functionality is built into ``plotRJForward.py`` demonstrated below.\n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "OOX-nNu77y8J"
      },
      "source": [
        "# c. Generating Results\n",
        "\n",
        "After you have confirmed the convergence of your chains, you can generate results using `plotRJforward.v3.py` in the LiteRate github. `plotRJforward.v3.py` collates the various log files, computes mean parameter estimates and 95\\% HPDs, and then generates a Rscript to plot the results. You can also open this script in R for further analysis of the rates, or you can run it to create a PDF. Let's look at the options for `plotRJforward.v3.py`:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_GO9ml8O-R3P"
      },
      "source": [
        "!python plotRJforward.v3.py -h"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "ntTQaGc1-n5g"
      },
      "source": [
        "- `combine` can be set to 0 or 1. When running MCMCs, it is common to practice to run several chains with different starting points to ensure that they have reached the same optima. If they do, you can combine the log files for greater statistical power.\n",
        "- `logT` plots rates on a log scale.\n",
        "-`burnin` It is common practice to drop the first 10-20% of samples from an MCMC chain before the chain has reached its stationary distribution.\n",
        "\n",
        "Now let's run it! Note that you only ran one chain, but these are combined results from 100 chains with different stochastic imputations."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "1V1ZC7BX7jCo"
      },
      "source": [
        "!python plotRJforward.v3.py ./example_data/metal_bands/combined_runs"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jSBI4LcH-d41"
      },
      "source": [
        "---\n",
        "\n",
        "# d. Visualizing Results\n",
        "\n",
        "This is what the resulting `metal_bands.100chains_plots.pdf` looks like:\n",
        "<html><img src='https://drive.google.com/uc?id=1OvKQj5OzVsj0HjlY-POcoORzA032oBKX', width=1000></html>\n",
        "\n",
        "The results are organized by colors and columns. With respect to color, blue indicates birth rates, red indicates death rates, and green indicates net (birth-death) rates/diversity. Let's now break it down by column.\n",
        "\n",
        "\\\\\n",
        "**Column 1** shows the distribution of the number of possible rate shifts (``K_j`` and ``K_m`` in the Tracer figure). This can be interesting just to see the number of shifts proposed and which one is most common. \n",
        "\n",
        "\\\\\n",
        "**Column 2** is the main result. This shows the average rates as a solid line and the 95% HPD in a single plot. Remember, the algorithm works by concatenating linear BD processes together, so individual samples will all have sharp rate shifts. But when you average over 10,000 samples, the rates appear rather smooth.\n",
        "\n",
        "\\\\\n",
        "**Column 3** shows the frequency with which rate shifts are observed in each time bin. These plots are especially interesting because LiteRate runs an MCMC simulation where rate shifts are sampled from their priors (no data). This simulation is used as a null hypothesis to compute Bayes factors for the statistical significance of rate shifts. The lower dotted line indicates the $2*log(Bayes\\ Factor)>2$. The upper dotted line indicates the $2*log(Bayes\\ Factor)>6$. We consider significant rate shifts to be those supported by $2\\log BF>2$ as ``Positive\" and $2\\log BF>6$ as ``Strong\", following  [Kass and  Raftery (1995)](https://doi.org/10.1080/01621459.1995.10476572). See [Silvestro 2019](https://doi.org/10.1017/pab.2019.23.) for more details.\n",
        "\n",
        "\\\\\n",
        "**Column 4** shows the net diversification rate and standing diversity. The net diversification rate is calculated by subtracting the death rate from the birth rate. You already know what the standing diversity is. \n",
        "\n",
        "**The second plot** integrates some of these plots, along with the empirical rates, in a pretty ``ggplot2`` plot. Remember, you can edit any of these plots on your own by modifying the ``.R`` file. An example is provided below.\n",
        "\n",
        "On the second page, you will find a **combined plot** that integrates empirical rates, estimated rates, and significant rate shifts. This is described below in the context of the Metal data. \n",
        "\n",
        "---"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2XE51MYw8loV"
      },
      "source": [
        "# e. Interpreting Results\n",
        "\n",
        "We now return to the Metal data from [Koch et al., in progress](https://osf.io/preprints/socarxiv/659bt) as an example of how to interpret these results. What can these plots tell us about how the Metal music genre changes over time? To clarifying what is happening, we plotted statistically significant rate shifts on top of these rates. We also plot the empirical rates for comparison.\n",
        "\n",
        "\\\\\n",
        "\n",
        "<html><img src='https://drive.google.com/uc?id=1w1hRJBYFENZSV6WIeiIszMaDH67wr6sX'</img>\n",
        "\n",
        "The LiteRate-estimated Metal band diversification rates largely support a multi-stage diversification trajectory also observed in the empirical rates. Broadly, we identified five phases: \n",
        "\n",
        "- Phase 1 (1968-1978): Experimentation by early/proto- Metal bands where birth rates start high, decline, and begin to rise again. This trend is observable in the empirical rates and, to a lesser extent, in the LiteRate-estimated rates. There are few bands active during this period, leading to high uncertainty around the LiteRate estimates as indicated by large shaded areas (95% highest posterior density intervals or HPDs). One possible interpretation is that rates start out high as the first-generation of bands emerges, but slow down as the sonic, aesthetic, and social parameters have not yet been concretized.  \n",
        "\n",
        "- Phase 2 (1978-1981): By 1978, the genre form solidifies (a key innovation) and there is sharp growth through 1981, coupled with increased turnover as bands rapidly explore the genre space (Phase 2).\n",
        "\n",
        "- Phase 3 (1981-1988): From 1981-1984, birth rates are stable and high, while death rates rise. In biology, the rapid churn of species early in a clade’s history has been interpreted as necessary for clades to build a \"standing diversity” of fit species that stabilize the clade’s existence. From 1984-1988, birth rates fall while death rates stabilize. The overall signature of converging birth and death rates in Phase 3, called *diversity dependent diversification*, is interpreted in biology as consistent with competition. We test this hypothesis in the next module (Interpreting Evolutionary Signatures). \n",
        "\n",
        "- Phase 4 (1988-1993): Birth and death rates stabilize, in between two periods of possible competition.\n",
        "- Phase 5 (1993-2000): This phase begins with another diversity dependence-like signature, before both birth and death rates slow, ossifying the stable of bands that make up the Metal genre. \n",
        "\n",
        "\n",
        "One of the most striking patterns in this data is that not only are birth rates greater than death rates *on average*, they are *always* greater than death rates. Moreover, while there is some range in death rates, they are not nearly as dynamic, suggesting that this is a birth driven process. \n",
        "\n",
        "A second puzzle is that the two periods of apparent competition (Phases 3, 5) are separated by a shelf with stable rates (Phase 4). The first possible interpretation of Phases 3-5 is that there is global competition operating from the genre's solidification in the early 1980s through 2000, and that Phase 4 is a momentary contextual deviation from this process. The second interpretation is that Phase 4 separates two independent periods of competition, potentially by two different groups of bands. \n",
        "\n",
        "It is interesting to note that Phase 4 also corresponds with Metal's peak in commercial appeal. This suggests that, at least for the majority of bands, the rise of Grunge music during this same time period did not cause a significant extinction in the Metal genre (although it may have for a subset of bands). However, the contraction in rates associated with the later Phase 5 coincides with what scholars call the Metal \"Dark Ages\" in the second half of the 1990s after the fall from pop music (Kahn-Harris 2006).\n",
        "\n",
        "---\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "018jPHTfMTM0"
      },
      "source": [
        "# Key Takeaways\n",
        "\n",
        "- **To assess convergence of LiteRate MCMC runs, make sure ESS is 200+ and traceplots of parameters look stable and stochastic over time.**\n",
        "\n",
        "\\\\\n",
        "\n",
        "- **Run multiple chains of LiteRate. If they all look good, combine them using ``logCombiner.py`` or the option in ``plotRJforward.py``.**\n",
        "\n",
        "\\\\\n",
        "\n",
        "- **``plotRJforward.py`` prints out average rates and timing of significant shifts in an ``.R`` files, and plots these results in a ``.pdf`` file.**\n",
        "\n",
        "\\\\\n",
        "\n",
        "- **The Metal example provides one interpretation of plotRJ results.**\n",
        "\n",
        "# Up next...\n",
        "\n",
        "One of our major claims throughout these modules is that diversification rates can give insight into the role of historical events and evolutionary mechanisms in cultural forms. One of the major claims in the Metal analysis is that the trajectory of this cultural form is driven by competition between cultural lineages. In the next tutorial, we explore theoretical evolutionary mechanisms at the macroevolutionary level as well as their corresponding rate signatures.\n",
        "\n",
        "\n",
        "---\n",
        "\n",
        "# References\n",
        "Kahn-Harris, Keith. *Extreme Metal: Music and Culture on the Edge.* Berg Publishers, Oxford 2006.\n",
        "\n",
        "Kass, Robert E., and Adrian E. Raftery. 'Bayes Factors’. *Journal of the American Statistical Association* 90, no. 430 (1995): 773–795.\n",
        "\n",
        "Silvestro, Daniele, Nicolas Salamin, Alexandre Antonelli, and Xavier Meyer. ‘Improved Estimation of Macroevolutionary Rates from Fossil Data Using a Bayesian Framework’. Paleobiology 45, no. 4 (2019): 546–70. https://doi.org/10.1017/pab.2019.23.\n",
        "\n",
        "Koch, Bernard, Daniele Silvestro, and Jacob G. Foster. n.d. “The Evolutionary Dynamics of Cultural Change (as Told Through the Birth and Brutal, Blackened Death of Metal Music).” SocArXiv. [osf.io/preprints/socarxiv/659bt](https://osf.io/preprints/socarxiv/659bt)."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "tb3Ec60Cwkv9"
      },
      "source": [
        "**License:** These tutorials are licensed under a [Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License](http://creativecommons.org/licenses/by-nc-sa/4.0/)."
      ]
    }
  ]
}