examples.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "571dfc7e-4d93-4f46-a31d-5ea54e8514bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d80c3343-a5a9-407d-a1bb-a4f53125038a",
   "metadata": {},
   "source": [
    "# BinBVS examples"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "354e3238-c911-4e9f-b8a4-383c56ff1650",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import binbvs\n",
    "import pathlib\n",
    "from scipy import sparse"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9157a06-6200-472a-b2eb-8714b4e77657",
   "metadata": {},
   "source": [
    "### Generate data\n",
    "Select the first `p0` features to be the 'true' features.\n",
    "Create the directory for MCMC output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7837571a-85fb-4842-9cd7-86d96c7f732c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "making examples_out\n"
     ]
    }
   ],
   "source": [
    "n = 100000 # number of samples\n",
    "p = 100 # number of simulated features\n",
    "p0 = 5 # number of 'true' features - i.e. features that are associated to the outcome\n",
    "\n",
    "# set the random seed\n",
    "np.random.seed(1)\n",
    "\n",
    "# Generate the sparse matrix\n",
    "X = sparse.random(n,p, density=0.01).tocsc()\n",
    "X_col_names = [f'gene_{i}' for i in range(p)]\n",
    "\n",
    "# initialise the outcome variable\n",
    "y = (X[:,:p0].sum(1)>0).A1.astype(int)\n",
    "\n",
    "# generate the output directory\n",
    "out_dir = pathlib.Path('examples_out')\n",
    "out_dir.mkdir(exist_ok=True)\n",
    "print('making', out_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de446f90-1a52-4b09-b2eb-0a01843434f7",
   "metadata": {},
   "source": [
    "### MCMC example\n",
    "1. Set the MCMC output (MCMC outputs can be large so we require that an output folder be specified)\n",
    "2. Generate the results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1921e624-4992-4e5e-a2ed-8072cefc4b44",
   "metadata": {},
   "source": [
    "##### Run MCMC using Metropolis-Hastings (rather than Gibbs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "53b80293-34f8-423a-bc8d-391cd53be63a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 100, acc=0.040000\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 200, acc=0.035000\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 300, acc=0.023333\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 400, acc=0.017500\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 500, acc=0.014000\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 600, acc=0.011667\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 700, acc=0.010000\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 800, acc=0.008750\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 900, acc=0.007778\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "iter: 1000, acc=0.007000\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "time taken for 1000 iterations\n",
      "0:00:13.297226\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PIP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>gene_0</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_2</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_3</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_4</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_1</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_64</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_73</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_72</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_71</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_70</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_69</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_68</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_67</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_66</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_65</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_61</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_63</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_62</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_75</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_60</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         PIP\n",
       "gene_0   1.0\n",
       "gene_2   1.0\n",
       "gene_3   1.0\n",
       "gene_4   1.0\n",
       "gene_1   1.0\n",
       "gene_64  0.0\n",
       "gene_73  0.0\n",
       "gene_72  0.0\n",
       "gene_71  0.0\n",
       "gene_70  0.0\n",
       "gene_69  0.0\n",
       "gene_68  0.0\n",
       "gene_67  0.0\n",
       "gene_66  0.0\n",
       "gene_65  0.0\n",
       "gene_61  0.0\n",
       "gene_63  0.0\n",
       "gene_62  0.0\n",
       "gene_75  0.0\n",
       "gene_60  0.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_rep = 10*p # Number of MCMC iterations\n",
    "\n",
    "\n",
    "mcmc = binbvs.MCMC(X, X_col_names, y, out_dir, n_rep, use_gibbs=False)\n",
    "\n",
    "# Run MCMC with initial condition given by a single randomly chosen feature\n",
    "mcmc.run(r_n_gamma0=1)\n",
    "display(mcmc.pip_df[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "397220af-6b81-4257-8622-94aab5322440",
   "metadata": {},
   "source": [
    "##### Run MCMC using Gibbs sampling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "84244dca-28c1-4285-9ffa-d2f678386973",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iter: 10, acc=0.100000\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n",
      "time taken for 10 iterations\n",
      "0:00:13.555676\n",
      "written to  examples_out/pip.csv\n",
      "written to  examples_out/out.pcl\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>PIP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>gene_0</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_2</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_3</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_4</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_1</th>\n",
       "      <td>1.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_64</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_73</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_72</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_71</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_70</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_69</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_68</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_67</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_66</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_65</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_61</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_63</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_62</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_75</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_60</th>\n",
       "      <td>0.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         PIP\n",
       "gene_0   1.0\n",
       "gene_2   1.0\n",
       "gene_3   1.0\n",
       "gene_4   1.0\n",
       "gene_1   1.0\n",
       "gene_64  0.0\n",
       "gene_73  0.0\n",
       "gene_72  0.0\n",
       "gene_71  0.0\n",
       "gene_70  0.0\n",
       "gene_69  0.0\n",
       "gene_68  0.0\n",
       "gene_67  0.0\n",
       "gene_66  0.0\n",
       "gene_65  0.0\n",
       "gene_61  0.0\n",
       "gene_63  0.0\n",
       "gene_62  0.0\n",
       "gene_75  0.0\n",
       "gene_60  0.0"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_rep = 10 # Number of MCMC iterations is reduced as when Gibbs sampling is used, a single iteration contains an entire sweep of all variables\n",
    "\n",
    "gibbs = binbvs.MCMC(X, X_col_names, y, out_dir, n_rep, use_gibbs=True)\n",
    "gibbs.run(r_n_gamma0=1)\n",
    "display(gibbs.pip_df[:20])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce3d2cba-bbb2-4e55-bc89-73a4087d54b8",
   "metadata": {},
   "source": [
    "### Example using Variational Bayes inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7a7e696e-9c82-453e-b0e4-082530cbd06b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0:00:04.610630\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>qPIP</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>gene_0</th>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_2</th>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_3</th>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_4</th>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_1</th>\n",
       "      <td>1.000000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_94</th>\n",
       "      <td>3.954042e-267</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_32</th>\n",
       "      <td>4.827377e-271</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_59</th>\n",
       "      <td>6.034413e-272</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_91</th>\n",
       "      <td>1.508635e-272</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_24</th>\n",
       "      <td>3.771668e-273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_7</th>\n",
       "      <td>3.771668e-273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_28</th>\n",
       "      <td>3.771668e-273</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_6</th>\n",
       "      <td>2.357393e-274</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_51</th>\n",
       "      <td>1.473433e-275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_50</th>\n",
       "      <td>1.473433e-275</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_31</th>\n",
       "      <td>9.209347e-277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_37</th>\n",
       "      <td>4.604723e-277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_20</th>\n",
       "      <td>4.604723e-277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_84</th>\n",
       "      <td>2.302386e-277</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>gene_87</th>\n",
       "      <td>2.878074e-278</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                  qPIP\n",
       "gene_0    1.000000e+00\n",
       "gene_2    1.000000e+00\n",
       "gene_3    1.000000e+00\n",
       "gene_4    1.000000e+00\n",
       "gene_1    1.000000e+00\n",
       "gene_94  3.954042e-267\n",
       "gene_32  4.827377e-271\n",
       "gene_59  6.034413e-272\n",
       "gene_91  1.508635e-272\n",
       "gene_24  3.771668e-273\n",
       "gene_7   3.771668e-273\n",
       "gene_28  3.771668e-273\n",
       "gene_6   2.357393e-274\n",
       "gene_51  1.473433e-275\n",
       "gene_50  1.473433e-275\n",
       "gene_31  9.209347e-277\n",
       "gene_37  4.604723e-277\n",
       "gene_20  4.604723e-277\n",
       "gene_84  2.302386e-277\n",
       "gene_87  2.878074e-278"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n_rep = 2 # Two iterations of variational inference attains good results with this simulated data\n",
    "\n",
    "vb = binbvs.VB(X, X_col_names, y)\n",
    "vb.run(n_rep)\n",
    "\n",
    "vb.qpip_df[:20] # note that the column name is named 'qPIP' to emphasise that this is an approximation to the true PIP"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "generic",
   "language": "python",
   "name": "generic"
  },
  "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.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}