ipynb/balance_trees.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# __Balance Trees__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "We will introduce the concept of balances here.  Balances were originally introduced in the geological sciences with the introduction of the isometric log ratio (ilr) transform.  Here we will refer to the visualizations representing balances as Balance Trees.\n",
    "\n",
    "Balances are an alternative way to interpret relative abundance information, without worrying about dependency properties between features inherient in multinomial distributions.  Often times, we cannot count all of the individuals within a community, and we are stuck analyzing representative samples.  If we have a community of 2 species, we can visualize the ratio or *balance* between these two species.  Below is a more intuitive explanation of a balance.\n",
    "\n",
    "<img src=\"images/Slide1.jpg\">\n",
    "\n",
    "Here we have two species, Red and Green that have equal counts.  The balance is taken as the log ratio between them.  In our convention, we will always take the ratio of the left species over the right species.  So in this case the balance is calculated as:\n",
    "\n",
    "$$\\log \\frac{Red}{Green} = \\log \\frac{5}{5} = \\log 1 = 0$$\n",
    "\n",
    "If the number of red individuals outnumber the number of green individuals 5 to 1, you'll end up with a positive balance of:\n",
    "\n",
    "$$\\log \\frac{Red}{Green} = \\log \\frac{5}{1} = \\log 5$$\n",
    "\n",
    "<img src=\"images/Slide2.jpg\">\n",
    "\n",
    "Likewise, if the number of green individuals outnumber the number of red individuals 5 to 1, you'll end up with a negative balance of\n",
    "\n",
    "$$\\log \\frac{Red}{Green} = \\log \\frac{1}{5} = - \\log 5$$\n",
    "\n",
    "<img src=\"images/Slide3.jpg\">\n",
    "\n",
    "Now, this scheme of balances can be generalized to any arbituary bifurcating trees.  But rather than comparing species abundances, we would be comparing mean log abundances of entire subtrees.  We'll refer to the concept of visualizing balances on bifurcating trees as __Balance Trees__. This concept has come under several names, such as sequential binary partitions, or CoDa dendrograms in the geological sciences.\n",
    "\n",
    "Below is how we would make this calculation on a tree with 4 leaves.\n",
    "\n",
    "<img src=\"images/Slide4.jpg\">\n",
    "\n",
    "The full equation on how to estimate balances is as follows\n",
    "$$\n",
    "b_i = \\sqrt{\\frac{l_i r_i}{l_i + r_i}} \n",
    " \\log \\big( \\frac{(\\prod_{x_n \\in L}{x_n})^{1/l_i}}\n",
    " {(\\prod_{x_m \\in R}{x_m})^{1/r_i}} \\big)\n",
    "$$\n",
    "\n",
    "where $b_i$ is the $i$th balance, $l_i$ is the number of leaves in the left subtree, $r_i$ is the number of leaves in the right subtree  $x_n$ are the abundances of the species in the left subtree, and $x_m$ are the abundances of species in the right subtree.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from skbio.stats.composition import ilr\n",
    "from skbio import TreeNode\n",
    "from gneiss.balances import balanceplot, balance_basis\n",
    "from gneiss.layouts import barchart_layout\n",
    "import scipy\n",
    "from functools import partial\n",
    "from ete3 import Tree, TreeStyle, faces, AttrFace, BarChartFace, CircleFace\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To illustrate how to use these modules, consider the following tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                    /-a\n",
      "          /e-------|\n",
      "         |          \\-b\n",
      "-r-------|\n",
      "         |          /-c\n",
      "          \\f-------|\n",
      "                    \\-d\n"
     ]
    }
   ],
   "source": [
    "tree = TreeNode.read([u\"((a,b)e, (c,d)f)r;\"])\n",
    "print(tree.ascii_art())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And consider the following abundances for the given sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "table = pd.Series([10, 20, 10, 10],\n",
    "                  index=['a','b','c','d'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given this data, we can visualize the ratios between the abundances according to this tree.\n",
    "Here, we will be visualizing the log ratios between b/a, d/c and f/e\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "basis, nodes = balance_basis(tree)\n",
    "ilr_coords = ilr(table, basis=basis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can visualize these log ratios as balances on the tree.  Since the c, d have the same abundances, their ratio is 1, and their log ratio is zero. So that balance on the c,d subtree is zero.\n",
    "\n",
    "On the other hand, the abundance of b is double the abundance of a.  So the ratio of b/a is 2, making the log ratio log(2). \n",
    "\n",
    "The logratios of subtrees can be visualized as balances on the tree as shown below.  In this scenario, the left subtree is on the bottom, and the right subtree is on the top."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAADECAYAAAA71CydAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAN1wAADdcBQiibeAAAFKBJREFUeJzt3X1sU9fdB/DvdRJiCmRNSkZSg4likzZkJI7iGlooiYY6\nQR9VE6PrSqEU1gqtgyFKlf6xaqJR24mKMNjatBNrljVIQXW8ojXt2m5dcSClYSKMZE0nh5DZeWlg\nD2KVkqXYvr7n+SNP3Jq8jMR2jnPz/Uj5w74vv2Mhfzn33ONzFSGEABGRRAbZDSAiYhARkXQMIiKS\njkFERNIxiIhIOgYREUnHICIi6RhERCQdg4iIpGMQEZF0DCKKmWAwCL/fP2vqUuwky24AzXwejwdN\nTU1ISkqCoigIBAJYsWIFVq1apcu6NLby8nIcPHhwSscyiChq58+fx44dO2AwfNXBdrlc8Pv9SE1N\n1V1diqSqKvbu3Yu6ujoGEckjhEAwGAx/+UOhEFRVhaIouqxLkQYGBvC9730PPp9vyudQuAwIRau7\nuxuNjY3hcZrk5GQ4HA4sX75cl3VpbA888AAaGhqmdCyDiIhiIpog4qUZRc3r9cLtdiMYDAIADAYD\n7HY7ioqKdFmXJnagphdHXZ9j10PZePqxJXj11VdRVVWFzz//HHPmzIHNZsOBAwdQXFwcPoY9Iora\n8ePH8eCDDyIlJQUAoGka6uvrsXHjRsyZMydh6lZUVMStLTSs4/JtcLXYEAhqSE0x4MTh5ej7x1sI\nhUL41re+hVAohJdffhkff/wxenp6kJSUBIA9IooBIUTEnStFUWAwGBDv/+MmU7eiogInT55EWVlZ\nXNtEAqFQCMDwDQNFAZ544onw1kAggDlz5sDlcuHSpUvIy8sDwCCiGHA4HKitrYWmaeE7VkVFRXG/\nhT7ZumVlZXjuuefi2iYCSmr7UFXfhz2bTVh/TwZ6e3vx7LPP4uTJk+jp6cEtt9wCABgcHAwfwyCi\nqFmtVlitVgghRvVS9FiXJrZvmwn7tpkADM8xcjgcKC4uxu9//3sUFBRgaGgImZmZEcfwX45iZuTS\nyO/3o7m5Wfd16b+7cOEC+vv78atf/Qp33XVXuDd0IwYRxZwQAoFAYNbUpfHl5+cjLS0Nx48fxxdf\nfIG2tjbs2bNn1H4MIoo5o9GItWvXzpq6NL558+bh5z//OX75y18iPT0d9957L5YsWTJqP96+p7ho\naWlBSUlJXGt0dXXh1KlT4XlEiqKgpKQkYn7KiIqKCgghOFgtSSAQwMWLF2GxWGA0Gkdt52A1Re3c\nuXMYGhqKeM/j8cQ9iM6ePYstW7aE5xEJIeB0OlFQUBDX+Us0eXPmzEFBQcG42xlEFLUFCxagp6cn\nYkZzb29v3OsKIUb9wHU65i9R7DGIKGpWqxVz586F2WwOv6dpWtzrrly5EseOHUMoFIKiKFAUBcXF\nxVwCZAZiEFHUkpKSIkIIGA6neLNYLLBYLADG7h3RzMG7ZqQLN4ZQeXm5pJbQVDCISFdUVcXu3btR\nXV0tuyk0Cbw0I12JxWqBNP3YIyJdSU9Px7e//W3ZzaBJYhARkXS8NKMp8/l8EEIgJydHdlPG5PP5\n8KqzH/VuA5alL8Hd1m7ZTaJxMIgk0cNqgV6vFwASMojsdjue2HMY7g4LVGMRui9bkJn2H9nNonEw\niCTQ02qBI2sBJSYBLTQ8sVIRw68pMTGIJNHDaoFOpxP9/f2w2WyymzJK62uv4YF19yB7cQo+7jYi\nL/1TWL95TXazaBwMIpoyk8mE9vZ2XLhwQXZTwjRNQ+eJE7jS0oKCO+7Azs33oXaNAxUV7yFhO27E\nIKKpU1UV6enpMX98zxtv5Ez52JDfj+u3rMK12z5DRsYgVFWNWbsofhhEFBWbzRbzsa7GxujP8Y9/\n/hM2mzWBx6/o6ziPiIikYxARkXQMItKl/Ecekd0EmgQGERFJxyAiIul0G0TBYFDKM66CwSD8fv+0\n1yWayXR3+97j8aCpqQlJSUlQFAWBQAArVqzAqlWrdFmXKFGUl5fj4MGDUzpWd0F0/vx57NixI+I5\n6C6XC36/P66LqsuqSySbqqrYu3cv6urqphxEurs0E0KEH7gHAKFQCKqqxn1hdVl1iWQbWRVz9erV\nY25/6aWXcPz48QnPobse0Zo1a+B0OsPjNMnJybjrrrvi/sC9seo6HA4+6I90b2RVzMOHD4+5/c9/\n/jNsNhs2b9487jl0F0RmsxmPPvrorKlLpAe6CyKv1wu32x2+TDIYDLDb7TH/YWai1CVKNAdqenHU\n9Tl2PZSNpx9bAmD48u2pp57CW2+9hUAggMcffxzPP/98eOhCd0F05swZPPLII+FLIk3TUF9fj/z8\n/LheJn3yyScRz2Gfrrp08xobG3WxMmaistvt2PLkK3C12BAIani2yofl1nkAgN/+9rd4/PHH8Ytf\n/ALvvvsuXnrpJZjNZuzcuROADoMIGH7y6AhFUableehCiIg7ZtNVl27O/v37E3w1Sb0QCIVCAIZ7\nOiP3anbt2oUjR44AADZt2oQvvvgC1dXV+g0ih8OB2tpaaJoW7vYVFRXF/Ra6rLp082b6ipgzRUlt\nH6rq+7Bnswnr78lAJYZv3nzdypUrUVlZGX6tuyCyWq2wWq3h//2+3kvRY12iRLNvmwn7tpkm3Gfh\nwoURU1t0+20ZuTTy+/1obm7WfV2imaStrQ35+fnh17rrEd1ICCHlN2ey6k631tbWmE/a9HqXxuQ8\nra2tKCwsjMm5KDqdnZ24evUqUlNT8bvf/Q7V1dWoqakJb9d9EBmNRqxdu3bW1J1OpaWlcRkAfuwx\nb4zOVKiLRzbpQVtbGxYvXhxeQ/xnP/sZvv/974e36z6IAKClpQUlJSVxrdHV1YVTp06F5xEpioKS\nkhIUFxfHta5s/KLTf/Phhx8CAAYHB+Hz+ZCTk4N58+ZF7KO7IDp37hyGhoYi3vN4PHEPorNnz0bM\nIxJCwOl0oqCggPOIiADMnz8fBQUFY27TXRAtWLAAPT09ETOae3t7415XCDFqrITziIhuju6CyGq1\nYu7cuTCbzeH3NE2Le92VK1fi2LFjCIVCUBQFiqKguLiY84iIboLugigpKSkihIDhcIo3i8UCi8UC\nYOzeERGNT7fziGS6MYTKy8sltYRoZmAQxZGqqti9ezeqq6tlN4Uooenu0iyRjKxc5/P5ZDeFKKGx\nRxRHIyvXEdHEGEREJJ1uLs18Ph+EEMjJyZHdlDH5fD686uxHvduAZelLcLe1W3aTiBJGMgBdrFrn\n9XoBICGDyG6344k9h+HusEA1FqH7sgWZaf+R3SyihJEMQDezfxN7BT4BLTQ8sVIRw6+JaFgyoI+V\n65xOJ/r7+2Gz2WQ3ZZTW117DA+vuQfbiFHzcbURe+qewfvOa7GYRJQzdjBGZTCa0t7fjwoULspsS\npmkaOk+cwJWWFhTccQd2br4PtWscqKh4DwnbcSOSQDdBpKoq0tPTE+rxPYGBAZg2boRHVXHmrxvQ\n3Z2Fv/wFcLtLAQATDc3t3z9NjSRKALoJIgCw2WwJuT7OiY8+QtfCLNnNIEpYnEdERNIxiIhIOgbR\nNNjY0CC7CUQJjUFERNIxiIhIOgYRxUwwGITf7581dSl2dHX7nuTweDxoampCUlISFEVBIBDAihUr\nsGrVKl3WpbGVl5fj4MGDUzqWQURRO3/+PHbs2AGD4asOtsvlgt/vj+vDA2TVpUiqqmLv3r2oq6tj\nEJE8QggEg8Hwlz8UCkFV1bg/QEBWXYoUi5VIGUQUtTVr1sDpdIbHaZKTk+FwOOL+YElZdSnSyEqk\nhw8fnvI5GEQUNbPZjEcffXTW1KXYYxBR1LxeL9xuN4LBIIDhJ9za7fa4/wBZVl2a2IGaXhx1fY5d\nD2Xj6ceWAABef/11HD16FF1dXTCbzdiyZQuefvrp8DEMIoraJ598gi1btiAlJQXA8PIn9fX1yM/P\nj+tl0mTr6mEl0kRmt9ux5clX4GqxIRDU8GyVD8ut89B6+ihefPFF/PCHP0RlZSXa2tpQU1PDIKLY\nEkJE3LlSFAUGgyHuq2VOpm5FRQVOnjyZkKsz6ItAKBQCMHzDIOD/Es8//zyeeuopvPDCCwCAtWvX\nYufOnRFHMYgoag6HA7W1tdA0LXzHqqioKO630Cdbt6ysTBerkSa6kto+VNX3Yc9mE5Zm/C+GhoZw\n//33R+xzY4+VQURRs1qtsFqt4TXDv95L0WNdmti+bSbs22YCADQ39wAYvrM2Ef7LUcyMXBr5/X40\nNzfrvi79d3l5eQCA999/f8L9dNUjam1tTdjJbFeuGLFo0SLZzZgWQggEAoFZU5fGl5GRgYcffhgv\nv/wysrOzsWnTJng8Hrzyyiv49a9/Hd5PN0FUWlqa0I8TeuaZL1FWlvP/rxohhMD+/WUSWxQ/RqMR\na9eunTV1aWJVVVXYvXs3tm7diq1bt0JRFDz88MMR++gmiADwjkgCaWlpQUlJSVxrdHV14dSpU+F5\nRIqioKSkBMXFxXGtS5OTkZGBuro61NTUoK+vD7fffjuMRmPEProKIpLj3LlzGBoainjP4/HEPYjO\nnj0bMY9ICAGn04mCggL+zCMBpaamIjc3d8xtDCKK2oIFC9DT0xMxo7m3tzfudYUQo8YEp2P+EsUe\ng4iiZrVaMXfuXJjN5vB7mqbFve7KlStx7NgxhEIhKIoCRVFQXFzMJUBmIAYRRS0pKSkihIDhcIo3\ni8UCi8UCYOzeEc0cnEdEunBjCJWXl0tqCU0Fg4h0RVVV7N69G9XV1bKbQpPASzPSlVisFkjTjz0i\n0pWR1QJpZmEQEZF0vDSjKfP5fBBCICcnR3ZTxuTz+fCqsx/1bgOWpS/B3dZu2U2icTCIJGlsbJzx\nKwZ6vV4ASMggstvteGLPYbg7LFCNRei+bEFm2n9kN4vGwSCSYP/+/Qn9A93JSOzPIaCFhidWKmL4\nNSUmBpEkelgp0Ol0or+/HzabTXZTRml97TU8sO4eZC9OwcfdRuSlfwrrN6/JbhaNg0FEU2YymdDe\n3o4LFy7IbkqYpmnoPHECV1paUHDHHdi5+T7UrnGgouI9JGzHjRhENHWqqiI9PT2hHt8TGBiAaeNG\neFQVaRkZUFVVdpPoJjCIKCo2my0h14E68dFHyLDZEnj8ir6O84iISDoGERFJxyAiXdrY0CC7CTQJ\nDCIiko5BRETSMYiISDoGERHFRDSrYjKIiCgqsVgVk0FERFEZWRVz9erVN33MT3/6U/zxj38Mv2YQ\nEVFUprIqZkNDA9rb28OvGUREJB2DiIhi6kBNL3L/56849EYPgOE1qyorK2G1WrFgwQKsW7cOV65c\niTiGP3qlWUMPq2ImMrvdji1PvgJXiw2BoIZnq3xYbp2HtqbfoKKiAtu3b8cPfvADdHZ24plnnok4\nlkFEs4KeVsVMbAKhUAjA8AMvNU3Fiy++iH379uGFF14AAJSWluLIkSMRRzGIaNbQw6qYM0FJbR+q\n6vuwZ7MJlsxrGBgYwHe+850Jj2EQEVFM7dtmwr5tJgBAc3MzACA7O3vCYzhYTURxYzabAQBnz56d\ncD/2iCgqra2tUBRFdjPG1NraisLCQtnNmNVuv/12rF69GgcPHsSiRYuwbNkyvPHGG+jo6IjYjz0i\nmrLS0lIUFRWFB4ET7a+wsDAhl7GdbY4ePRoeJ8rNzUVXVxdMJlPEPorgbQQiijNN09DZ2YmsrCyk\npaWN2s4gIiLpeGlGRNIxiIhIOgYREUnHICIi6RhERCQdg4iIpGMQEZF0DCIiko5BRETSMYiISDpD\nY2Mj3G637HZETS+fg2g2Mrjdbl18gfXyOYhmI16aEZF0DCIiko5BRETSJXu9Xni93hn/vCe3242c\nnBzZzSCiKUgGoJvnPenhMxDNRsk5OTlYunTpjH/mU0VFBYOIaIbiGBERSccgIiLpGEREJF1yWVlZ\nxNiK3+9Hampq1CceHBxEcnIyjEZj1Oe6GTd+DiKaOQylpaUoKyvD+fPnsXr1aqSlpSEvLw8ffPDB\nuAeZTCbceuut4b+77747vO3KlSvYsGEDli5divnz52Pbtm3T8Tkw8jmIaOZJBoDr169j69at2L59\nO9xuN9599108+OCD+Oyzz7BkyZKIAy5fvoxgMIgzZ86E3/t6r+fJJ5+E1WrFO++8g08//RQbNmxA\nU1MT1qxZM00fiYhmHCGE+NOf/iQyMzNFKBQSI/Ly8sShQ4fEjU6fPi3uvffeUe8LIURHR4cAILq7\nu8Pv/e1vfxOXLl0ac38iIiGEMACAx+OBzWaDwfDV2HVeXh46OjpGBVdnZyeGhoawfv16FBQUYPv2\n7ejr6wMAdHR0IDc3Fx988AHWrVuHH/3oR/jyyy+Rm5s7PalKRDOSAQC6urpw6623RmxYunQpBgcH\nRx3g9XoxMDCAn/zkJ6isrITP58O6detw/fp1XLx4Ef39/aivr8ePf/xjJCcno6ysDO+///70fBoi\nmpGSAWDevHn417/+FbFhcHAQd95556gDnnvuuYhZ2MuWLcOyZcvQ3NyMUCiEYDCI2tpaLFq0CJs2\nbUJvby/efPNNrF+/Pr6fhIhmLAMAWCwWXLp0KWLDxYsXUVhYOOqADz/8EN3d3eHXWVlZSE1NRVZW\nFsxmM77xjW9g4cKF4e3Z2dnw+/3xaj8R6YABAL773e/i6tWrcLlcEEKgrq4Of//73+FwOAAAb775\nJtrb2wEAp0+fxsaNG6GqKvx+Pw4dOoT8/HzceeeduO+++6BpGl5//XVomobOzk68/fbbuP/+++V9\nQiJKfCOj1n/4wx/E3LlzRXZ2tpg/f75wuVzhEe3FixeLI0eOCCGEuHbtmigrKxOLFy8WmZmZIisr\nS5w8eTK8b0NDg0hPTxfZ2dnCYDCIXbt2CU3Tpm30nYhmHkWIr6Yj+/1+eDwe5OfnIyUlZcIAu3r1\nKv79738jNzcXSUlJEdtGzmMymXDbbbfFJ0GJSDcigoiISAb+6JWIpGMQEZF0DCIiko5BRETSMYiI\nSDoGERFJxyAiIukYREQkHYOIiKRjEBGRdAwiIpKOQURE0jGIiEg6BhERSccgIiLpGEREJB2DiIik\nYxARkXQMIiKSjkFERNIxiIhIOgYREUnHICIi6RhERCQdg4iIpGMQEZF0DCIiko5BRETS/R/IDtMu\nDN60JgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr, ts = balanceplot(ilr_coords, tree, mode='r', layout=barchart_layout) \n",
    "ts.branch_vertical_margin = 20  # adjust the spacing between leaves.\n",
    "tr.render(file_name='%%inline', tree_style=ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cool thing about these sorts of plots is statistical tests can be run directly on the balances.  \n",
    "\n",
    "For example, below, we have a table of 6 samples $s1$, $s2$, $s3$, $s4$, $s5$, $s6$ and 7 features $a$, $b$, $d$, $f$, $h$, $j$, $l$, which are related with a given tree.  Suppose that samples $s1$, $s2$, $s3$ are within a control group, and samples $s4$, $s5$, $s6$ are within a treatment group.  We would like to see if there is a difference between the sample and treatment groups. \n",
    "We can do this by running ANOVAs on each balance, and visualize the p-values on each internal node.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>a</th>\n",
       "      <th>b</th>\n",
       "      <th>d</th>\n",
       "      <th>f</th>\n",
       "      <th>h</th>\n",
       "      <th>j</th>\n",
       "      <th>l</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>s1</th>\n",
       "      <td>10</td>\n",
       "      <td>11</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s2</th>\n",
       "      <td>10</td>\n",
       "      <td>12</td>\n",
       "      <td>11</td>\n",
       "      <td>11</td>\n",
       "      <td>11</td>\n",
       "      <td>11</td>\n",
       "      <td>13</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s3</th>\n",
       "      <td>11</td>\n",
       "      <td>11</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s4</th>\n",
       "      <td>20</td>\n",
       "      <td>21</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s5</th>\n",
       "      <td>20</td>\n",
       "      <td>21</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "      <td>10</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>s6</th>\n",
       "      <td>21</td>\n",
       "      <td>21</td>\n",
       "      <td>10</td>\n",
       "      <td>9</td>\n",
       "      <td>10</td>\n",
       "      <td>11</td>\n",
       "      <td>12</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     a   b   d   f   h   j   l\n",
       "s1  10  11  10  10  10  10  10\n",
       "s2  10  12  11  11  11  11  13\n",
       "s3  11  11  10  10  10  10  10\n",
       "s4  20  21  10  10  10  10  10\n",
       "s5  20  21  10  10  10  10  10\n",
       "s6  21  21  10   9  10  11  12"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "table = pd.DataFrame([\n",
    "            [10, 11, 10, 10, 10, 10, 10],\n",
    "            [10, 12, 11, 11, 11, 11, 13],\n",
    "            [11, 11, 10, 10, 10, 10, 10],\n",
    "            [20, 21, 10, 10, 10, 10, 10],\n",
    "            [20, 21, 10, 10, 10, 10, 10],\n",
    "            [21, 21, 10, 9,  10, 11, 12]],\n",
    "            index=['s1', 's2', 's3', 's4', 's5', 's6'],\n",
    "            columns=['a', 'b', 'd', 'f', 'h', 'j', 'l'])\n",
    "cats = pd.Series(['Control', 'Control', 'Control', \n",
    "                  'Treatment', 'Treatment', 'Treatment'],\n",
    "                index=['s1', 's2', 's3', 's4', 's5', 's6'])\n",
    "tree = TreeNode.read([u\"((((a,b)5, (f, d)6)3, (h, j)4)2, l)1;\"])\n",
    "table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The tree is structure as shown below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                        /-a\n",
      "                              /5-------|\n",
      "                             |          \\-b\n",
      "                    /3-------|\n",
      "                   |         |          /-f\n",
      "                   |          \\6-------|\n",
      "          /2-------|                    \\-d\n",
      "         |         |\n",
      "         |         |          /-h\n",
      "-1-------|          \\4-------|\n",
      "         |                    \\-j\n",
      "         |\n",
      "          \\-l\n"
     ]
    }
   ],
   "source": [
    "print(tree.ascii_art())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As before, we can transform these abundances into log ratios as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "basis, nodes = balance_basis(tree)\n",
    "balances = ilr(table.values, basis=basis)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can calculate ANOVAs on each balance, by splitting up the balances by the control and treatment groups."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pvalues = []\n",
    "for b in range(balances.shape[1]):\n",
    "    groups = [balances[(cats == x).values, b] for x in set(cats)]\n",
    "    f, p = scipy.stats.f_oneway(*groups)\n",
    "    pvalues.append(-np.log(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now we can plot the balances on the tree.  As expected, the only difference is in the balance with a,b vs f,d.\n",
    "The effect of that balance diminishes on higher levels of the tree."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPsAAACkCAYAAABRqliEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAN1wAADdcBQiibeAAAEfxJREFUeJzt3X9Q1XW+x/EnHFDkpydFERIJxUDqrlop5bS59nNqW0er\n2w+bakaud11ZK29qYzlFv1YGK1qtHBp07d7VNU2bWNtp1ZRp07QhjXQ1RUIFyzJFDx4Ez4/7x6cv\ncACRgHO+53s+78fMd45yfszn6Lz4fL+fz/f9+YR5vV4vQoiQF252A4QQgSFhF8JC5s6d2+33StiF\nsACXy0VeXh4lJSXd/oyIXmyPEMJPHA4HU6dO5ciRI93+DOnZhbAAu93OpEmTevQZEnYhNCFhF8KC\nFq2oIf2uXby68hgAb731FtnZ2djtdgYPHsztt9/O7t27fd4TJvPsQnf5+flmN+EXOfj9ANaVj6bp\ngoe+keFseH0UtfvX43a7ueqqq3C73SxZsoTPPvuMY8eOYbPZABmgE5rLz89n69atTJw40eym/AJe\n3G43EAZAWBjk5uY2P9vU1ESfPn1Yt24dhw8fZuTIkYCEXWhm7ty5FBYW+vxs4sSJPP/88+Y0qJuu\nebeWN9fWMvvBFO644TJqamp45pln2Lp1K8eOHSM6OhqA+vr65vfINbvQQm/MUweTOY+kcLh0HI8/\nlILL5WLcuHGcPHmS999/n3PnznU4RSc9u9BCb8xTB6s9e/bw3Xff8emnnzJ8+HAAnE5nu9dJzy60\n0Bvz1MEqKyuL+Ph4Vq9eTV1dHRUVFcyePbvd6yTsQlhcTEwMr7zyCm+88QZ2u50bb7yRoUOHtnud\nhF1oJePMGaiubp6n3n6ofSisaNasWdTW1rJ3715OnDhBQUEBXq+XsWPHNr9G5tlFj1llnjqlpoYb\nams5647moRPT+DZsALZwLw+Mr+D/3vqj2c3zO+nZRY8Y89Rerzf4D4Bz5+D0GdweDwBhXgA9+jsZ\njRc9ZqV56i8ffZScF1/k8W023lxby0j7XkYMOmV2s7qso/sEukp6dqEFj8vFlrw8tpeWQmpq8zx1\nzvBjZjetS6SeXYguanI4yJg6lbMWnWe/1H0CBQUFpKam8uCDD170M6RnF1qIsttJtfA8+6XuE9i0\naRPl5eWdfoaEXbT44QeoroamJrNbIvxATuN198EHsGsXfPst9OkD/fpBXR0kJsLVV8NvfwvJyWa3\nUrSxaEUNxeuOM+s/h/A/j6p7BRwOB08++STr16+nqamJ6dOn8+KLLxIWpqrjJOy6+vprKChQoU5I\ngAkTVNjDw8HrhbNnoaoKFiyAW2+FadPMbnGvmFJa2u5nZWVllrlX4Nprr2XazKXN9ezPvHmEUSNi\nAFi+fDnTp0/ntddeY+PGjc3X8TNmzAAk7Hr64ANYvx7GjoXoaOjbVx19+oDNpgqkBwyAlBRwOuGL\nL2DnTigqUr8MQshzzz3XPA9vHe3r2UHdRVdUVATAPffcQ11dHSUlJRJ2ba1dC3//O1xzjQp6v37q\niIpSgY+MVIEHde3udMJ116lr+ZkzmRsf3+153mBllXsEWmtbz74YiIjwjfP48eNZvHhx898l7Dqp\nqoING2DcuJaQR0f7ht7o4cPD4cIFdcdZZCQer5cdGzfi+uorCLGwW9GcR1KY80hKp68ZOHBg8/U6\nSNj1UlAA2dkqzMape1SUOtqGvk8f8HjUtXtYGI0OB8m33srk776DQ4cgI8PsbyMuoaKigqysrOa/\nh9YFmLi47dvVY2ysOlWPjFSBNh6Nwwh+fLwavOvfH+Li6JeYyBXZ2ZwaOBDWrDH3u4gOVVZWcvLk\nSRwOB0uWLKGkpIS8vLzm56Vn18WePSrINhtERKjHtn+22dTpu3Hq5/Go541fBH37cj42FvbtM/e7\niA5VVFRw+eWX43K5AFi4cCH33Xdf8/MSdl0cOgR2e0uYjUfjADXl5vGAywWNjepnTU3qZ2FhLb8c\nwsJYUrSb1zddIMM+lOtHHDXvewkANm/eDKgFJo8cOUJaWhoxMTE+r5GwBwF/z/GGeTxM+eQTfho6\nFGJiWq7P216ntx6gs9lU+Bsbob5e3Wjz008kp6TwWXklf63YzbeRIzj6/XAS48/5tf2i62JjY8nO\nzu7wOblmN1kg6sGhVcW219vSg3s84Harw+VSvXhjIzQ0qCm3c+fU4/nz6jm3W70Hfp7n1ase3Oqk\nZw8CAakHP3u2ZdCt7fy6McduDNgZN9Z4PCrkDQ2qdz9zhur33iPt6lRyr7+OU/9sslw9uNVJPbu4\ntIwMdSpu9OIul5pHv3ChpUc/f963Vzd6dqcTr9NJ1YYN/LBzJ3g8/PcfR1uqHtzqeqOeXcKui9Gj\nVaCNsLcNetuwtw58QwPuujoGpKdjHzBAFciIgDLq2SdMmHDR1+zfv5+nn36ae++9lx07drR7XsKu\nixtuUKfm9fUdB72hoX3QjT87nUR4vSQkJRF1+jR0skCC8I9L1bMfOHCA0aNHU1dXx/jx49vdOgty\nza6X+fNVFdvYsS0FLcZgnTHf3vrnbrf6pWD0+ocOcbJ/f4amp5v3HUSHVq9ezaRJk1i2bNlFXyNh\n10laGtx7r7o//le/8p1Xj4hQQTfCbozUG3PuJ06A2011aipjTP0SAnzr2aPqP2Tp0qVERkaSk5ND\nYmIipR2U8krYdTN1qhpxX7NG3ScfFeXbq4eFtUzPGT37N9+okfslS5gSYiWuYJ1176Hjevai2WPI\nzMwkKiqKvLw8oqKiOnyvhF1Hd98N6emqei0uTh0JCeqXgMHhgNOn1XHnnXD//ea1149CYX/2tPRM\nUlJSiI2NZcqUKRd9l4RdV9nZ8Je/QGmpWpzi669Vrx4drabohgxRo+533QVJSWa3tteE6v7sy4su\n/R4Ju+7uvlsdACdPqtH35GR1ah9CXC4XTzzxBKtWrQqJxTe6Us/eVmj9j4qeGTjQ7Bb4TSjvz95V\noTfaIkQHQnl/9q6SsAuhCTmNF1qxNzbC8eMs+thD8brjIVOP/957713yNRL23lJdDYcPw/79aqAr\nNlYVn2RlwciRZrfOr6yy7npMfT3znE52PPY4b//4MEexa1WPL2HvDW+/DXv3qrvN7HZVRupwwJYt\n8I9/wGWXwVNPqccQY7V11z0uF566OtxuD9j0qseXsPeEwwF/+AMMGqTmrdsuxezxqFtNa2th9mxY\nsIC5K1aExNRPa1aao94xaxYTFi5kzj/dsj+76CKnE2bMUPeYX3mlugOtf3/Vew8cCIMHqxtTkpMh\nMxPPpEl8ed997CguNrvlWjL2Zy9fvRqSkrTcn13C3l1/+pO65TQ2Vt11FhPTctup3a4Cb4Q9KYnG\nhATsU6eyKC7O7JZrydifPaWTevBg1pV69rYWLFjARx991Px3CXt3bN8OP/2kbiM1lnVqvcNKbKzq\n5RMTVdgHDaJfcjJXjBnD6UGD4J13zP4G2gn1/dk7Ulpayr5Wy35L2Lvjk09UmI3NFlofrdZYp18/\n1dvHxjafATguuwzKysz+BkJDEvbuMNZgN8pCjcfW5Z9utxqcczrVQF14eMsvA2BJ0W7S79rF9kND\nTfoSwsoWragh/a5dvLpSjTl4vV4WL17MiBEjiIuL4+abb+bEiRM+7wmJ0fhAzvH2a2jgjooKTp05\n09Jjx8S0PMbEtKzaGhmp6sKdTlVJdvKkWnd95zf8dY+su26GUN2fveJf75Cfn89jjz3G/fffT2Vl\nJfPmzfN5r+XDHuh6ZK+xsENH668ba6+Hh7dssODxqCWdWq+9jkfVI0fqNc8bjKx2n4DiW8/u8bh4\n+eWXmTNnDi+99BIAN910U/Ne7QbLhx1MqEf+/nsYNUqNvLfeYcXo0fv0ad4mqXmu3Vh3/W9/I+0/\n0sgdL+uuBwsr3SdgaF3PPjzxFA6Hg9tuu63T98g1e3dceaVawaX1kswdrbteX69uvHE48DocVL3/\nPj98/jmArLsuesS4T+Dxh1Koq6sDYMiQIZ2+R8LeHb/5DZw54xt0I+ytg97qcJ86xYBhwxgQFweW\nWgJJBLvU1FQAdu7c2enrQuI0PuCuvx4+/lidzqektCzS6Hb7bp9k/KypiQi3m4TEROrr6iA31+xv\nIEJIcnIyEyZMoLCwkMGDB5ORkcHKlSs5ePCgz+ukZ++up5+GY8dUz916gwVjF5U2PTvnzsGuXRzI\nyDC75SIEFRcXN1+3p6enU1VVRUqK77JV0rN3V1QUFBdDXp7aMDEtTZ3St91oweWCH39U5a8LF3K2\nzXSIEL1h1KhRVFZWUllZSVJSEvHx8e1eI2HvidhYtUJrcTF89ZXq2fv3V8fPo++Auq327bchIaHD\neV4hekN4eDgjO1k7QcLeG2bMgKNHoaoKDhxQPXlyshrIy8hQhxAmk7D3ltRUdchIe1DrST241ckA\nndBCb9SDW52EXWjhUvXgBQUFrF69OsCtCiwJu9DCperBN23aRHl5eQBbFHgSdiE0IWEXWmpbDw7Q\n0NDAs88+y/Dhw8nKyqKoqAiPx2NiK3uXjMaLHrNKLThcvB4cYNmyZTz88MMUFhZy9OhR5s2bR05O\nDjk5OSa3undI2EWPhML+5mHqgZkzZ7J06dLmV5WWlrJx40YJu1l0nicNBqG6v/liICoqyuc1SUlJ\nHD582JwG+oFlrtllntRcofbv37oe/GJsNlsAW+R/lunZZX9tc8m/v/VZpme32+1Mys42uxnakv3N\nrc8aPfvWrbByJSQkMO3gQbUiTJvrKyFE54I/7Nu3Q0kJPPAADB9O/KlT8MQTLBr/bEjtr20Vobq/\n+ebNm9v97N133zWhJf4TAcE9T3rz5s3YUlNp/Pe/4fx5on/3Oz5/fTn/u/1ffGtLDol1162ybrnu\n+5tbXQQQ/Gtmu1xqMYjaWmhsxHvhAh6POyT217bauuU6729udREQ5Otm33EH/PnPcMUVkJDA8RUr\nSB53Jf81fqIl99fuSFD/+7eh8/7mVhf8o/E5OfD73+MtK+P4Cy9wbMsWePVVy+2vbXWyv7n1Bf8A\nHcCvf03j1Vfj2bqV4ytWyEi8CYz9zc9adJ69K/cJfPnll0RHR5OZmRnAlgVO8PfsP4uy27l86lSz\nm6EtHfY3z83Nbd4rLRRZo2cXIgB27txJmFEVE4Is07ML0Zta17OXlZWRk5PDjTfe6FP1Fmos17PL\nuuvmCsX9zUueHcT8+fN54YUXQvref8uFXQQXq90noPjWsw8YOJg7bshi2bJl5jbLzyTsosesdJ+A\noW09uw4k7EJLcx5JYc4jF69lD0UyQCeEJiTsQmhCwi7Ez6w1yPjLSdiF+FldXR1Dhw41uxl+IwN0\nQmvl5eXU1NRQVVXF7t27Wbt2rdlN8hsJu9BaTU0N8+fPJy0tjVWrVjFs2DCzm+Q3EvYA07meOhhN\nnjyZyZMnm92MgJBr9gCRemphNunZA+RS9dSFhYUkJyczbdq0ALdM6EJ69gC5VD31li1b+OKLLwLY\nIqEbCbsQmpCwm6Sj/cGF8KeQuGa3ej21sT+4EP5k+bCHQj11CK+EJIKI5cMOoVFPXWR2g0TIC4mw\nW5GO9dTCXDJAJ4QmJOxCaELCLoQmwrzWGsYWQnST9OxCaCK8rKyMbdu2md2OHguV7yGEv4Rv27Yt\nJEISKt9DCH+R03ghNCFhF0ITEnYhNBFRXV1NdXW1ZarGLmbbtm2kpaWZ3QwhglYEYMGqsY6FwncQ\nwl8i0tLSGDZsmCUrx1rLz8+XsAvRCblmF0ITEnYhNCFhF0ITERMnTvS51nW5XLjdbvr27XvJN58/\nf57IyEhsNpvPz3/JZ/SWtt9DCOGrXdXbjBkzACguLu70jU6nk4yMDFauXMktt9zi81xXP0MIETjN\np/Effvghubm5LF++vNM3NDY2snjxYu68806OHz/u81xXP0MIEXjNYT969ChxcXFkZmZ2+oYLFy5Q\nW1vLmDFjiIyM9Hmuq58hhAi8dqfx06dPx2azdekUPC4ujg0bNrQ7jf8lnyGECAwZjRdCExJ2ITQh\nYRdCE10K+5o1a9i3b5+/2yKE8KMuhf2pp55i8+bN/m6LEMKPZClpITQh1+xCaELCLoQmJOxCaELC\nLoQmJOxCaELCLoQmJOxCaELCLoQmJOxCaELCLoQmJOxCaELCLoQmJOxCaELCLoQmJOxCaELCLoQm\nJOxCaELCLoQmJOxCaELCLoQmJOxCaELCLoQmJOxCaELCLoQmJOxCaELCLoQmJOxCaELCLoQmJOxC\naOL/AXBka0Tq6KmiAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tr, ts = balanceplot(pvalues, tree, mode='r') \n",
    "tr.render(file_name='%%inline', tree_style=ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is one thing to be wary about when using this module.  Since each balance is being calculated using log ratios,\n",
    "there can be no zeros in the original table.  This means that the zeros will need to be removed.  One approach is to add pseudocounts to\n",
    "all of the values.  The most common pseudocount to add is *1*.  There are many other approaches such as collapsing leaves.\n",
    "\n",
    "In conclusion, balance trees can be alternative, yet powerful approach at analyzing the community structure of a population."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}