Added model probability estimation to model selection

edsgard · May 24, 2017 · febeb87 · febeb87
1 parent 278394c
commit febeb87
Show file tree

Hide file tree

Showing 5 changed files with 356 additions and 142 deletions.
diff --git a/Analysis/MouseOB/MOB_MS_results.csv b/Analysis/MouseOB/MOB_MS_results.csv
diff --git a/Analysis/MouseOB/MOB_final_results.csv b/Analysis/MouseOB/MOB_final_results.csv
diff --git a/Analysis/MouseOB/Mouse Olfactory Bulb Analysis.ipynb b/Analysis/MouseOB/Mouse Olfactory Bulb Analysis.ipynb
diff --git a/Analysis/SeqFISH/SeqFISH Results.ipynb b/Analysis/SeqFISH/SeqFISH Results.ipynb
@@ -1964,7 +1964,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python [default]",
+   "display_name": "Python 3",
    "language": "python",
    "name": "python3"
   },
@@ -1978,7 +1978,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.5.2"
+   "version": "3.6.0"
   }
  },
  "nbformat": 4,

diff --git a/Python-module/SpatialDE/base.py b/Python-module/SpatialDE/base.py
@@ -7,6 +7,7 @@
 from scipy import linalg
 from scipy import stats
 from scipy.misc import derivative
+from scipy.misc import logsumexp
 from tqdm import tqdm
 import pandas as pd
 
@@ -438,9 +439,20 @@ def model_search(X, exp_tab, DE_mll_results, kernel_space=None):
     logging.info('Performing model search')
     results = dyn_de(X, de_exp_tab, kernel_space)
     new_and_old_results = pd.concat((results, DE_mll_results))
+
+    # Calculate model probabilities
+    mask = new_and_old_results.groupby(['g', 'model'])['BIC'].transform(min) == new_and_old_results['BIC']
+    log_p_data_Hi = -new_and_old_results[mask].pivot_table(values='BIC', index='g', columns='model')
+    log_Z = logsumexp(log_p_data_Hi, 1)
+    log_p_Hi_data = (log_p_data_Hi.T - log_Z).T
+    p_Hi_data = np.exp(log_p_Hi_data).add_suffix('_prob')
+
+    # Select most likely model
     mask = new_and_old_results.groupby('g')['BIC'].transform(min) == new_and_old_results['BIC']
     ms_results = new_and_old_results[mask]
 
+    ms_results = ms_results.join(p_Hi_data, on='g')
+
     # Retain information from significance testing in the new table
     transfer_columns = ['pval', 'qval', 'max_ll_null']
     ms_results = ms_results.drop(transfer_columns, 1) \