<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistical testing for ensemble simulations\n",
    "\n",
    "In traditional COBRA simulations with a single model, most simulations result in a single quantity of interest, thus statistical comparisons usually don't make sense. For example, when simulating growth in two different media conditions, a single model can only output a single predicted growth rate for each condition.\n",
    "\n",
    "When accounting for uncertainty in model structure using an ensemble, these simulations generate a distribution rather than a single value. Because we are no longer comparing two individual values, proper interpretation requires statistical assessment of the distributions our ensemble simulations generate. In this example, we demonstrate this concept and one statistical option for univariate comparisons (e.g. comparisons between two conditions).\n",
    "\n",
    "First, let's load an ensemble for *Staphylococcus aureus* and the recipe for biolog growth media, which we'll use to simulate growth in single carbon source media."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import medusa\n",
    "from medusa.test import create_test_ensemble\n",
    "\n",
    "ensemble = create_test_ensemble(\"Staphylococcus aureus\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "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>Name</th>\n",
       "      <th>ID</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>H2O</td>\n",
       "      <td>cpd00001_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O2</td>\n",
       "      <td>cpd00007_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Phosphate</td>\n",
       "      <td>cpd00009_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CO2</td>\n",
       "      <td>cpd00011_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NH3</td>\n",
       "      <td>cpd00013_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Mn2+</td>\n",
       "      <td>cpd00030_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Zn2+</td>\n",
       "      <td>cpd00034_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Sulfate</td>\n",
       "      <td>cpd00048_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Cu2+</td>\n",
       "      <td>cpd00058_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Ca2+</td>\n",
       "      <td>cpd00063_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>H+</td>\n",
       "      <td>cpd00067_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>Cl-</td>\n",
       "      <td>cpd00099_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Co2+</td>\n",
       "      <td>cpd00149_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>K+</td>\n",
       "      <td>cpd00205_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Mg</td>\n",
       "      <td>cpd00254_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>Na+</td>\n",
       "      <td>cpd00971_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>Fe2+</td>\n",
       "      <td>cpd10515_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>fe3</td>\n",
       "      <td>cpd10516_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>Heme</td>\n",
       "      <td>cpd00028_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>H2S2O3</td>\n",
       "      <td>cpd00268_e</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Name          ID\n",
       "0         H2O  cpd00001_e\n",
       "1          O2  cpd00007_e\n",
       "2   Phosphate  cpd00009_e\n",
       "3         CO2  cpd00011_e\n",
       "4         NH3  cpd00013_e\n",
       "5        Mn2+  cpd00030_e\n",
       "6        Zn2+  cpd00034_e\n",
       "7     Sulfate  cpd00048_e\n",
       "8        Cu2+  cpd00058_e\n",
       "9        Ca2+  cpd00063_e\n",
       "10         H+  cpd00067_e\n",
       "11        Cl-  cpd00099_e\n",
       "12       Co2+  cpd00149_e\n",
       "13         K+  cpd00205_e\n",
       "14         Mg  cpd00254_e\n",
       "15        Na+  cpd00971_e\n",
       "16       Fe2+  cpd10515_e\n",
       "17        fe3  cpd10516_e\n",
       "18       Heme  cpd00028_e\n",
       "19     H2S2O3  cpd00268_e"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "biolog_base = pd.read_csv(\"../medusa/test/data/biolog_base_composition.csv\", sep=\",\")\n",
    "biolog_base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert the biolog base to a dictionary, which we can use to set ensemble.base_model.medium directly.\n",
    "biolog_base = {'EX_'+component:1000 for component in biolog_base['ID']}\n",
    "\n",
    "# Double check that the objective is set to the biomass reaction.\n",
    "# For this model, 'bio1' is the id of the biomass reaction.\n",
    "ensemble.base_model.objective = 'bio1'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's simulate growth on two different carbon sources, D-glucose (metabolite id: cpd00027) and maltose (metabolite id: cpd00179)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from medusa.flux_analysis import flux_balance\n",
    "\n",
    "carbon_sources = [\"EX_cpd00027_e\",\"EX_cpd00179_e\"]\n",
    "\n",
    "fluxes = {}\n",
    "for carbon_source in carbon_sources:\n",
    "    biolog_base[carbon_source] = 10\n",
    "    ensemble.base_model.medium = biolog_base\n",
    "    fluxes[carbon_source] = flux_balance.optimize_ensemble(ensemble,return_flux='bio1')\n",
    "    biolog_base[carbon_source] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's visualize the distributions of predicted flux through biomass using matplotlib. We'll generate a histogram for each condition, and plot the mean using a vertical line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "<KEY>//<KEY>",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "bins=np.histogram(np.hstack((fluxes[carbon_sources[0]]['bio1'],fluxes[carbon_sources[1]]['bio1'])), bins=20)[1]\n",
    "plt.hist(fluxes[carbon_sources[0]]['bio1'],\n",
    "         bins = bins,\n",
    "         label=carbon_sources[0],\n",
    "         color = \"red\",\n",
    "         alpha = 0.5)\n",
    "plt.hist(fluxes[carbon_sources[1]]['bio1'],\n",
    "         bins = bins,\n",
    "         label=carbon_sources[1], \n",
    "         color = \"blue\", \n",
    "         alpha = 0.5)\n",
    "plt.axvline(x=fluxes[carbon_sources[0]]['bio1'].mean(), c = 'red')\n",
    "plt.axvline(x=fluxes[carbon_sources[1]]['bio1'].mean(), c = 'blue')\n",
    "ax.set_ylabel('# ensemble members')\n",
    "ax.set_xlabel('Flux through biomass reaction')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually, we can see the mean for D-glucose (cpd00027) is slightly lower than for maltose (cpd00179). To evaluate this statistically, we'll use the Mann-Whitney U test (implemented in SciPy), which tests the null hypothesis that random sample drawn from a population (e.g. growth in D-glucose) is equally likely to be greater than or less than a random sample drawn from another population (e.g. growth in maltose). The Mann-Whitney U test is suitable for univariate comparisons regardless of the distribution of data (e.g. when data are non-normally distributed, replace a *t*-test with the Mann-Whitney U test)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "MannwhitneyuResult(statistic=362758.0, pvalue=1.103684856884951e-26)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import mannwhitneyu\n",
    "mannwhitneyu(fluxes[carbon_sources[0]],fluxes[carbon_sources[1]])"
   ]<commit_msg>docs: update stats compare with paired stats test<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Statistical testing for ensemble simulations\n",
    "\n",
    "In traditional COBRA simulations with a single model, most simulations result in a single quantity of interest, thus statistical comparisons usually don't make sense. For example, when simulating growth in two different media conditions, a single model can only output a single predicted growth rate for each condition.\n",
    "\n",
    "When accounting for uncertainty in model structure using an ensemble, these simulations generate a distribution rather than a single value. Because we are no longer comparing two individual values, proper interpretation requires statistical assessment of the distributions our ensemble simulations generate. In this example, we demonstrate this concept and one statistical option for univariate comparisons (e.g. comparisons between two conditions).\n",
    "\n",
    "First, let's load an ensemble for *Staphylococcus aureus* and the recipe for biolog growth media, which we'll use to simulate growth in single carbon source media."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import medusa\n",
    "from medusa.test import create_test_ensemble\n",
    "\n",
    "ensemble = create_test_ensemble(\"Staphylococcus aureus\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "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>Name</th>\n",
       "      <th>ID</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>H2O</td>\n",
       "      <td>cpd00001_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>O2</td>\n",
       "      <td>cpd00007_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Phosphate</td>\n",
       "      <td>cpd00009_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>CO2</td>\n",
       "      <td>cpd00011_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NH3</td>\n",
       "      <td>cpd00013_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Mn2+</td>\n",
       "      <td>cpd00030_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>Zn2+</td>\n",
       "      <td>cpd00034_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>Sulfate</td>\n",
       "      <td>cpd00048_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>Cu2+</td>\n",
       "      <td>cpd00058_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>Ca2+</td>\n",
       "      <td>cpd00063_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>H+</td>\n",
       "      <td>cpd00067_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>Cl-</td>\n",
       "      <td>cpd00099_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>Co2+</td>\n",
       "      <td>cpd00149_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>K+</td>\n",
       "      <td>cpd00205_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>Mg</td>\n",
       "      <td>cpd00254_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>Na+</td>\n",
       "      <td>cpd00971_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>Fe2+</td>\n",
       "      <td>cpd10515_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>fe3</td>\n",
       "      <td>cpd10516_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>Heme</td>\n",
       "      <td>cpd00028_e</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>H2S2O3</td>\n",
       "      <td>cpd00268_e</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "         Name          ID\n",
       "0         H2O  cpd00001_e\n",
       "1          O2  cpd00007_e\n",
       "2   Phosphate  cpd00009_e\n",
       "3         CO2  cpd00011_e\n",
       "4         NH3  cpd00013_e\n",
       "5        Mn2+  cpd00030_e\n",
       "6        Zn2+  cpd00034_e\n",
       "7     Sulfate  cpd00048_e\n",
       "8        Cu2+  cpd00058_e\n",
       "9        Ca2+  cpd00063_e\n",
       "10         H+  cpd00067_e\n",
       "11        Cl-  cpd00099_e\n",
       "12       Co2+  cpd00149_e\n",
       "13         K+  cpd00205_e\n",
       "14         Mg  cpd00254_e\n",
       "15        Na+  cpd00971_e\n",
       "16       Fe2+  cpd10515_e\n",
       "17        fe3  cpd10516_e\n",
       "18       Heme  cpd00028_e\n",
       "19     H2S2O3  cpd00268_e"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "biolog_base = pd.read_csv(\"../medusa/test/data/biolog_base_composition.csv\", sep=\",\")\n",
    "biolog_base"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert the biolog base to a dictionary, which we can use to set ensemble.base_model.medium directly.\n",
    "biolog_base = {'EX_'+component:1000 for component in biolog_base['ID']}\n",
    "\n",
    "# Double check that the objective is set to the biomass reaction.\n",
    "# For this model, 'bio1' is the id of the biomass reaction.\n",
    "ensemble.base_model.objective = 'bio1'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's simulate growth on two different carbon sources, D-glucose (metabolite id: cpd00027) and maltose (metabolite id: cpd00179)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from medusa.flux_analysis import flux_balance\n",
    "\n",
    "carbon_sources = [\"EX_cpd00027_e\",\"EX_cpd00179_e\"]\n",
    "\n",
    "fluxes = {}\n",
    "for carbon_source in carbon_sources:\n",
    "    biolog_base[carbon_source] = 10\n",
    "    ensemble.base_model.medium = biolog_base\n",
    "    fluxes[carbon_source] = flux_balance.optimize_ensemble(ensemble,return_flux='bio1', num_processes = 4)\n",
    "    biolog_base[carbon_source] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's visualize the distributions of predicted flux through biomass using matplotlib. We'll generate a histogram for each condition, and plot the mean using a vertical line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pylab as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "<KEY>",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "bins=np.histogram(np.hstack((fluxes[carbon_sources[0]]['bio1'],fluxes[carbon_sources[1]]['bio1'])), bins=20)[1]\n",
    "plt.hist(fluxes[carbon_sources[0]]['bio1'],\n",
    "         bins = bins,\n",
    "         label=carbon_sources[0],\n",
    "         color = \"red\",\n",
    "         alpha = 0.5)\n",
    "plt.hist(fluxes[carbon_sources[1]]['bio1'],\n",
    "         bins = bins,\n",
    "         label=carbon_sources[1], \n",
    "         color = \"blue\", \n",
    "         alpha = 0.5)\n",
    "plt.axvline(x=fluxes[carbon_sources[0]]['bio1'].mean(), c = 'red')\n",
    "plt.axvline(x=fluxes[carbon_sources[1]]['bio1'].mean(), c = 'blue')\n",
    "ax.set_ylabel('# ensemble members')\n",
    "ax.set_xlabel('Flux through biomass reaction')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Visually, we can see the mean for D-glucose (cpd00027) is slightly lower than for maltose (cpd00179). To evaluate this statistically, we'll use the Wilcoxon signed-rank test (implemented in SciPy), which tests the null hypothesis that the difference between paired samples (e.g. growth in D-glucose minus growth in maltose for each ensemble member) is symmetrically distributed around zero. Here, we choose a statistical test meant for paired data because each simulation result in one media condition has a related simulation result in the other condition which was generated using the same ensemble member. The Wilcoxon signed-rank test is suitable for paired univariate comparisons regardless of the distribution of data (e.g. when data are non-normally distributed, replace a paired *t*-test with the Wolcoxon signed-rank test)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "WilcoxonResult(statistic=0.0, pvalue=3.3257599356529824e-165)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import wilcoxon\n",
    "cond1 = fluxes[carbon_sources[0]].copy()\n",
    "cond2 = fluxes[carbon_sources[1]].copy()\n",
    "cond1.columns = [carbon_sources[0]]\n",
    "cond2.columns = [carbon_sources[1]]\n",
    "both_conditions = pd.concat([cond1,cond2], axis = 1, join_axes = [cond1.index])\n",
    "\n",
    "wilcoxon(x=both_conditions[carbon_sources[0]],y=both_conditions[carbon_sources[1]])"
   ]
{'fraction_non_alphanumeric': 0.19440124416796267, 'fraction_numerical': 0.037513549177623826, 'mean_word_length': 1.8105960264900662, 'pattern_counts': {'":': 130, '<': 364, '<?xml version=': 0, '>': 364, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6578403', 'n_tokens_mistral': 8484, 'n_tokens_neox': 7826, 'n_words': 1500}