perform_FDR_correction.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.io import loadmat\n",
    "from scipy import stats\n",
    "from IPython.core.debugger import set_trace\n",
    "import pickle\n",
    "import os\n",
    "from tqdm import tqdm\n",
    "import random\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "random.seed(97)\n",
    "np.random.seed(97)\n",
    "\n",
    "root_save_dir = \"predictions/\"\n",
    "sub_data = \"sub_space_data/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def FDR(vector, q, do_correction = False):\n",
    "    original_shape = vector.shape\n",
    "    vector = vector.flatten()\n",
    "    N = vector.shape[0]\n",
    "    sorted_vector = sorted(vector)\n",
    "    if do_correction:\n",
    "        C = np.sum([1.0/i for i in range(N)])\n",
    "    else:\n",
    "        C = 1.0\n",
    "    thresh = 0\n",
    "    #a=b\n",
    "    for i in range(N-1, 0, -1):\n",
    "        if sorted_vector[i]<= (i*1.0)/N*q/C:\n",
    "            thresh = sorted_vector[i]\n",
    "            break\n",
    "    thresh_vector = vector<=thresh\n",
    "    thresh_vector = thresh_vector.reshape(original_shape)\n",
    "    thresh_vector = thresh_vector*1.0\n",
    "    print(\"FDR threshold is : {}, {} voxels rejected\".format(thresh, thresh_vector.sum()))\n",
    "    return thresh_vector, thresh"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_subjects = [\"F\",\"G\",\"H\",\"I\",\"J\",\"K\",\"L\",\"M\",\"N\"]\n",
    "all_feature_pairs = [\n",
    "                     (\"node_count_punct\", \"punct_final\"), # (NC + PU) - (PU)\n",
    "                     (\"syntactic_surprisal_punct\", \"punct_final\"), # (SS + PU) - (PU)                 \n",
    "                     (\"word_frequency_punct\", \"punct_final\"), # (WF + PU) - (PU)\n",
    "                     (\"word_length_punct\", \"punct_final\"), # (WL + PU) - (PU)\n",
    "                     (\"all_complexity_metrics_punct\", \"punct_final\"), # (CM + PU) - (PU)\n",
    "                     (\"pos_dep_tags_all_complexity_metrics\", \"all_complexity_metrics_punct\"), # since PD already contains PU, this tests (PD + CM + PU) - (CM + PU)\n",
    "                     (\"aggregated_contrege_comp_pos_dep_tags_all_complexity_metrics\", \"pos_dep_tags_all_complexity_metrics\"), # (CC + PD + CM + PU) - (PD + CM + PU)\n",
    "                     (\"aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics\", \"pos_dep_tags_all_complexity_metrics\"), # (CI + PD + CM + PU) - (PD + CM + PU)\n",
    "                     (\"aggregated_incontrege_pos_dep_tags_all_complexity_metrics\", \"pos_dep_tags_all_complexity_metrics\"), # (INC + PD + CM + PU) - (PD + CM + PU)\n",
    "                     (\"aggregated_bert_PCA_dims_15_contrege_incomp_pos_dep_tags_all_complexity_metrics\", \"aggregated_contrege_incomp_pos_dep_tags_all_complexity_metrics\") # (BERT + CI + PD + CM + PU) - (CI + PD + CM + PU)\n",
    "                    ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_uncorrected_sig = []\n",
    "for sub in all_subjects:\n",
    "    print(sub)\n",
    "    punct_uncorrected_sig = np.load(root_save_dir + \"punct_final\" + \"/{}_sig.npy\".format(sub))\n",
    "    all_uncorrected_sig.append(punct_uncorrected_sig)\n",
    "    \n",
    "    for feat in all_feature_pairs:\n",
    "        uncorrected_sig = np.load(root_save_dir + \"{}_diff_{}\".format(feat[0],feat[1]) + \"/{}_sig_boot.npy\".format(sub))\n",
    "        all_uncorrected_sig.append(uncorrected_sig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "q = 0.05\n",
    "all_corrected_sig, _ = FDR(np.hstack(all_uncorrected_sig), q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "last_end = 0\n",
    "ind = 0\n",
    "for sub in all_subjects:\n",
    "    print(sub)\n",
    "    punct_corrected_sig = all_corrected_sig[last_end:last_end + all_uncorrected_sig[ind].shape[0]]\n",
    "    np.save(root_save_dir + \"punct_final\" + \"/{}_sig_group_corrected\".format(sub), punct_corrected_sig)\n",
    "    last_end += all_uncorrected_sig[ind].shape[0]\n",
    "    ind += 1\n",
    "    \n",
    "    print(\"punct_final: Num voxels rejected = {}\".format(punct_corrected_sig.sum()))\n",
    "    \n",
    "    for feat in all_feature_pairs:\n",
    "        corrected_sig = all_corrected_sig[last_end:last_end + all_uncorrected_sig[ind].shape[0]]    \n",
    "        np.save(root_save_dir + \"{}_diff_{}\".format(feat[0],feat[1]) + \"/{}_sig_bootstrap_group_corrected\".format(sub), corrected_sig)\n",
    "        last_end += all_uncorrected_sig[ind].shape[0]\n",
    "        ind += 1\n",
    "        \n",
    "        print(\"{}: Num voxels rejected = {}\".format(\"{}_diff_{}\".format(feat[0],feat[1]), corrected_sig.sum()))"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}