jupyter_notebooks/HypSelector_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# HypSelector Example\n",
    "\n",
    "This is an example use of the HypSelector utility with a grouped reactions set of model hypotheses\n",
    "as generated by HypBuilder. \n",
    "\n",
    "In this example, we do the selection\n",
    "with Nested Sampling via Gleipnir's MultiNest wrapper class, so you will need to download and build PyMultiNest and MultiNest to run this example. You also need to have HypBuilder installed in your PYTHONPATH to run this example. See the Gleipnir [README](https://github.com/LoLab-VU/Gleipnir/blob/master/README.md) for links to these dependencies.\n",
    "\n",
    "This example was adapted from the grouped_reactions_example from HypBuilder:\n",
    "https://github.com/LoLab-VU/HypBuilder/blob/master/grouped_reactions_example.csv\n",
    "There is also a script version of this example in the Gleipnir examples: \n",
    "https://github.com/LoLab-VU/Gleipnir/tree/master/examples/HypSelector/grouped_reactions\n",
    "\n",
    "The data used in this example is synthetic data generated from model_0 with the default\n",
    "parameters defined in the csv file; they are the last 10 timepoints. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pysb\n",
    "from pysb.simulator import ScipyOdeSimulator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The HypSelector class is imported from the pysb_utilities module:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gleipnir.pysb_utilities import HypSelector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can set the HypBuilder syntax model csv file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The HypBuilder format model csv file.\n",
    "model_csv = 'grouped_reactions.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we need to load the data that we want to use in the likelihood function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the data.\n",
    "data = np.load(\"model_0_AB_complex_data.npy\")\n",
    "# The timespan of the simulations.\n",
    "tspan = np.linspace(0, 5, 20)\n",
    "# Define the fancy indexer or mask for the time points that the data\n",
    "# corresponds to. -- In this case it is the last ten (out of 20) time points.\n",
    "data_time_idxs = np.array(list(range(len(tspan))))[10:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the data is a synthetic dataset composed of the timeseries output of the AB complex from model_0. Here is what all the model outputs look like (with the points used for calibration in red):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Model outputs](model_outputs_fitting-data.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We also need to prepare the observable data for HypSelector. The data should be packaged in a dictionary keyed to the name of the observable (as defined, or to be defined in our case, in the pysb model files). The value should be 3 element tuple of the form: (data, data_sd, time_idxs). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the observable data tuple for this observable: (data, data_sd, data_time_idxs)\n",
    "obs_data_t = tuple((data,None,data_time_idxs))\n",
    "# Generate the dictionary of observable data that is to be used in\n",
    "# computing the likelihood. -- Here we are just using the AB_complex\n",
    "# observable, which is the amount of A(B=1)%B(A=1).\n",
    "observable_data = dict()\n",
    "observable_data['AB_complex'] = obs_data_t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can initialize the HypSelector with model csv file. Internally, the HypSelector instance will pass the model to HypBuilder which builds all the model variants. The created models will be in a new directory/folder named 'hb_models' and each model will be named in the format 'model_idx'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['./output/grouped_reactions/model_1.py', './output/grouped_reactions/model_2.py', './output/grouped_reactions/model_0.py']\n",
      "['hb_models/model_1.py', 'hb_models/model_2.py', 'hb_models/model_0.py']\n"
     ]
    }
   ],
   "source": [
    "# Build the HypSelector.\n",
    "selector = HypSelector(model_csv)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check the number of models that were generated. In this case, we should have 3 model variants. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generated 3 models from input csv\n"
     ]
    }
   ],
   "source": [
    "# Check the number of models that were generated.\n",
    "n_models = selector.number_of_models()\n",
    "print(\"Generated {} models from input csv\".format(n_models))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to add the following line to each of the models so that the appropriate observable is defined. This can also be done in the via the HypBuilder syntax in the csv file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Append the needed observable to the model files\n",
    "obs_line = \"Observable(\\'AB_complex\\',A(B=1)%B(A=1))\"\n",
    "selector.append_to_models(obs_line)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we need to call the gen_nested_samplers function to have the HypSelector instance construct the Nested Sampler instances for each of the models. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "multinest\n"
     ]
    }
   ],
   "source": [
    "selector.gen_nested_samplers(tspan, observable_data,\n",
    "                             solver=ScipyOdeSimulator,\n",
    "                             ns_version='multinest',\n",
    "                             ns_population_size=100,\n",
    "                             log_likelihood_type='mse')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now ready to run. The output is a sorted pandas DataFrame with the natural logarithm of the evidences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  analysing data from multinest_run_model_0_.txt\n",
      "  analysing data from multinest_run_model_1_.txt\n",
      "  analysing data from multinest_run_model_2_.txt\n"
     ]
    }
   ],
   "source": [
    "sorted_log_evidences = selector.run_nested_sampling()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ouput is a sorted pandas DataFrame of the log_evidence, its error estimate and the model names; the output is sorted by the log_evidence values so that the frame indices are the rankings (lowest index => highest ranked), and the first row in the frame has the largest evidence value and the highest ranked model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The models and their log_evidence values-sorted:\n",
      "   log_evidence  log_evidence_error    model\n",
      "0     -9.142524            0.279983  model_0\n",
      "1   -155.750248            0.399444  model_2\n",
      "2   -960.328540            0.238591  model_1\n"
     ]
    }
   ],
   "source": [
    "print(\"The models and their log_evidence values-sorted:\")\n",
    "print(sorted_log_evidences)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Other model selection criteria\n",
    "\n",
    "In addition to the evidence values, there are a variety of [criteria](https://en.wikipedia.org/wiki/Model_selection#Criteria) used for model selection. HypSelector provides functions to estimate four additional model selection criteria from the Nested Sampling output, including Bayes Factors and three information criteria: Akaike Information Criterion, Bayesian Information Criterion, and Deviance Information Criterion. The following subsections briefly cover each of these criteria.  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayes Factors\n",
    "\n",
    "It's also a common practice in model selection applications to compare two models by computing the [Bayes factor](https://en.wikipedia.org/wiki/Bayes_factor) for one model over the other. From Nested Sampling, we can directly estimate the Bayes factor for two models by computing the quotient of their respective evidences; note that this assumes there is no a priori reason to prefer one model to another (i.e., each model's prior probability is the same). From the HypSelector instance we can call the bayes_factor function: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "bayes_factors = selector.bayes_factors()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This function returns a cross correlation-like pandas DataFrame of the Bayes for each model comparison. Each element is given by: (Evidence of the model along the column)/(Evidence of the model along the row/index)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DataFrame of the Bayes factor matrix:\n",
      "              model_0  model_1       model_2\n",
      "model_0  1.000000e+00      0.0  2.133410e-64\n",
      "model_1           inf      1.0           inf\n",
      "model_2  4.687330e+63      0.0  1.000000e+00\n"
     ]
    }
   ],
   "source": [
    "print(\"DataFrame of the Bayes factor matrix:\")\n",
    "print(bayes_factors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bayes factors for the model_0 numerator ratios:\n",
      "model_0    1.000000e+00\n",
      "model_1             inf\n",
      "model_2    4.687330e+63\n",
      "Name: model_0, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "print(\"Bayes factors for the model_0 numerator ratios:\")\n",
    "print(bayes_factors['model_0'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Akaike Information Criterion\n",
    "\n",
    "The [Akaike Information Criterion](https://en.wikipedia.org/wiki/Akaike_information_criterion) (AIC) is another commonly used metric for model selection. HypSelector provides a function that returns estimates of the AIC for each model, which can used if you are interested in comparing the model ranking by evidence to that by AIC value. From the HypSelector instance call the akaike_ic function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted_aic = selector.akaike_ic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function returns a sorted pandas DataFrame that is similar to the sorted_log_likelihoods returned from the run_nested_sampling function, but is sorted in ascending order by AIC value; with AIC, the model with lowest AIC value is the highest ranked. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The models and their AIC values-sorted:\n",
      "           AIC    model\n",
      "0     4.001650  model_0\n",
      "1   285.991843  model_2\n",
      "2  1919.185540  model_1\n"
     ]
    }
   ],
   "source": [
    "print(\"The models and their AIC values-sorted:\")\n",
    "print(sorted_aic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the AIC values rank the models in the same order as the log-evidence values. \n",
    "\n",
    "Note: It's probably worth noting that in this case the AIC value for a model is estimated from the highest likelihood sample found during the Nested Sampling run. In most cases this should be fine because the highest likelihood samples collected during the Nested Sampling will typically approach the maximum likelihood value (assuming the likelihood is actually bound from above), but it is worth keeping in mind that the output values may not actually correspond to the true global maximum of the likelihood function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bayesian Information Criterion\n",
    "\n",
    "The [Bayesian Information Criterion](https://en.wikipedia.org/wiki/Bayesian_information_criterion) (BIC) is another commonly used metric for model selection; the BIC is quite similar to the AIC but tends to have a stronger penalty for increasing the number of parameters in a model. HypSelector provides a function that returns estimates of the BIC for each model, which can used if you are interested in also comparing the model ranking by BIC value. From the HypSelector instance call the bayesian_ic function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted_bic = selector.bayesian_ic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function returns a sorted pandas DataFrame that is similar to the sorted_log_likelihoods returned from the run_nested_sampling function, but is sorted in ascending order by AIC value; with BIC; like with AIC, the model with lowest BIC value is the highest ranked. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The models and their BIC values-sorted:\n",
      "           BIC    model\n",
      "0     4.606820  model_0\n",
      "1   287.807354  model_2\n",
      "2  1921.001051  model_1\n"
     ]
    }
   ],
   "source": [
    "print(\"The models and their BIC values-sorted:\")\n",
    "print(sorted_bic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the BIC values rank the models in the same order as the log-evidence values and the AIC values. You can also see that in this case the BIC values are quite similar to the AIC values. \n",
    "\n",
    "The same note as with the AIC: It's probably worth noting that in this case the BIC value for a model is estimated from the highest likelihood sample found during the Nested Sampling run. In most cases this should be fine because the highest likelihood samples collected during the Nested Sampling will typically approach the maximum likelihood value (assuming the likelihood is actually bound from above), but it is worth keeping in mind that the output values may not actually correspond to the true global maximum of the likelihood function. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Deviance Information Criterion\n",
    "\n",
    "The [Deviance Information Criterion](https://en.wikipedia.org/wiki/Deviance_information_criterion) (DIC) is another metric used for model selection; unlike the AIC and BIC, DIC uses the Bayesian posterior distribution to make its estimate. HypSelector provides a function that returns estimates of the DIC for each model, which can used if you are interested in also comparing the model ranking by DIC value. From the HypSelector instance call the deviance_ic function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "sorted_dic = selector.deviance_ic()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function returns a sorted pandas DataFrame that is similar to the sorted_log_likelihoods returned from the run_nested_sampling function, but is sorted in ascending order by DIC value; with DIC; like with AIC and BIC, the model with lowest DIC value is the highest ranked. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The models and their DIC values-sorted:\n",
      "          DIC    model\n",
      "0    2.061977  model_0\n",
      "1  277.887056  model_2\n",
      "2         inf  model_1\n"
     ]
    }
   ],
   "source": [
    "print(\"The models and their DIC values-sorted:\")\n",
    "print(sorted_dic)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, the DIC values rank the models in the same order as the log-evidence values, the AIC values, and the BIC values. You can also see that for the most part the DIC values are similar to the AIC values and BIC values. In the case of model_1, the value of inf here implies that the model evidence evaluated to zero within the numerical precision.    "
   ]
  }
 ],
 "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.6.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}