BayesianGLM_Power.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9199e51a",
   "metadata": {},
   "source": [
    "# Bayesian analysis of emotion-mediated memory as a function of power"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "54a2808f",
   "metadata": {},
   "source": [
    "Here, I am going to integrate neural data into my mixed-effects logistic regression models. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f89f70c3-f2ac-49d9-93d0-6b592fea3dc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "import bambi as bmb\n",
    "import argparse\n",
    "from os.path import join\n",
    "import pandas as pd\n",
    "from pandas.api.types import CategoricalDtype\n",
    "import seaborn as sns \n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.backends.backend_pdf import PdfPages\n",
    "from Bayesian_model_utils import run_model, plot_res, print_latex_table, plot_predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e647e104-9954-47d4-a7e8-751dffa4df5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model 1: No interaction effect here: \n",
    "\n",
    "pow_df = pd.read_csv('/home1/salman.qasim/Salman_Project/FR_Emotion/final_NHB_power_data_revised.csv', memory_map=True)\n",
    "\n",
    "# rename the electrodes to properly do the elec random effects \n",
    "pow_df.elec = pow_df.apply(lambda x: f'{x.subj}_{x.elec}', axis=1)\n",
    "\n",
    "# De-mean the valence  \n",
    "pow_df['valence'] = pow_df['valence'] - 0.5\n",
    "\n",
    "# Set categorical orders \n",
    "\n",
    "cat_type = CategoricalDtype(categories=['right', 'left'], ordered=True)\n",
    "pow_df.hemi = pow_df.hemi.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['Amy', 'Hipp'], ordered=True)\n",
    "pow_df.region = pow_df.region.astype(cat_type)\n",
    "\n",
    "# Let's do separate models for band\n",
    "\n",
    "for band in ['theta']:\n",
    "        \n",
    "    model_df = pow_df[pow_df.band==band]\n",
    "\n",
    "    y = 'memory'\n",
    "    X = ['arousal', 'valence', 'region', 'hemi', 'power'] \n",
    "    Intx = ['power:arousal', \n",
    "            'power:valence',\n",
    "            'power:hemi',\n",
    "            'power:region',\n",
    "            'power:region:hemi']\n",
    "\n",
    "    label = (f\"{y}\" + \"_{}\"*len(X)).format(*X) + f'_{band}' + 'no_intx'\n",
    "\n",
    "#             'power:hemi:region:arousal',\n",
    "#             'power:hemi:region:CV'\n",
    "    rand_effect = ['subj', 'elec']\n",
    "    categorical = ['hemi', 'region']\n",
    "\n",
    "    # Drop nan data\n",
    "    model_df = model_df.dropna(subset=X)\n",
    "\n",
    "    run_model(model_df, y, X, Intx, rand_effect, rand_slopes=False, categorical=categorical,\n",
    "             chains=4, cores=4, tune=500, draws=1000, label=label, save_model_res=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40bc447f-8391-4b03-8fab-da3de653deb7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model 2: Examine the effect of HFA power and word features on recall. \n",
    "\n",
    "pow_df = pd.read_csv('/home1/salman.qasim/Salman_Project/FR_Emotion/final_NHB_power_data_revised.csv', memory_map=True)\n",
    "\n",
    "# rename the electrodes to properly do the elec random effects \n",
    "pow_df.elec = pow_df.apply(lambda x: f'{x.subj}_{x.elec}', axis=1)\n",
    "\n",
    "# De-mean the valence  \n",
    "pow_df['valence'] = pow_df['valence'] - 0.5\n",
    "\n",
    "# Set categorical orders \n",
    "\n",
    "cat_type = CategoricalDtype(categories=['right', 'left'], ordered=True)\n",
    "pow_df.hemi = pow_df.hemi.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['Amy', 'Hipp'], ordered=True)\n",
    "pow_df.region = pow_df.region.astype(cat_type)\n",
    "\n",
    "# Let's do separate models for band\n",
    "\n",
    "for band in ['hfa', 'theta']:\n",
    "        \n",
    "    model_df = pow_df[pow_df.band==band]\n",
    "\n",
    "    y = 'memory'\n",
    "    X = ['arousal', 'valence', 'region', 'hemi', 'power'] \n",
    "    Intx = ['power:arousal', \n",
    "            'power:valence',\n",
    "            'power:hemi',\n",
    "            'power:region',\n",
    "            'power:region:hemi',\n",
    "            'power:arousal:valence']\n",
    "\n",
    "    label = (f\"{y}\" + \"_{}\"*len(X)).format(*X) + f'_{band}' + '_all_elecs' + 'no_soz'\n",
    "\n",
    "#             'power:hemi:region:arousal',\n",
    "#             'power:hemi:region:CV'\n",
    "    rand_effect = ['subj', 'elec']\n",
    "    categorical = ['hemi', 'region']\n",
    "\n",
    "    # Drop nan data\n",
    "    model_df = model_df.dropna(subset=X)\n",
    "\n",
    "    run_model(model_df, y, X, Intx, rand_effect, rand_slopes=False, categorical=categorical,\n",
    "             chains=4, cores=4, tune=500, draws=1000, label=label, save_model_res=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d94baae2-5ffb-4540-af5e-a4d589fdaedc",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model 3: Fullest model \n",
    "\n",
    "pow_df = pd.read_csv('/home1/salman.qasim/Salman_Project/FR_Emotion/final_NHB_power_data_revised.csv', memory_map=True)\n",
    "\n",
    "# rename the electrodes to properly do the elec random effects \n",
    "pow_df.elec = pow_df.apply(lambda x: f'{x.subj}_{x.elec}', axis=1)\n",
    "\n",
    "# De-mean the valence  \n",
    "pow_df['valence'] = pow_df['valence'] - 0.5\n",
    "\n",
    "# Set categorical orders \n",
    "\n",
    "cat_type = CategoricalDtype(categories=['right', 'left'], ordered=True)\n",
    "pow_df.hemi = pow_df.hemi.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['Amy', 'Hipp'], ordered=True)\n",
    "pow_df.region = pow_df.region.astype(cat_type)\n",
    "\n",
    "# Let's do separate models for band\n",
    "\n",
    "for band in ['hfa', 'theta']:\n",
    "        \n",
    "    model_df = pow_df[pow_df.band==band]\n",
    "\n",
    "    y = 'memory'\n",
    "    X = ['arousal', 'valence', 'region', 'hemi', 'power'] \n",
    "    Intx = ['power:arousal', \n",
    "            'power:valence',\n",
    "            'power:hemi',\n",
    "            'power:region',\n",
    "            'power:region:hemi',\n",
    "            'power:valence:hemi',\n",
    "            'power:valence:region', \n",
    "            'power:arousal:hemi',\n",
    "            'power:arousal:region',\n",
    "            'power:arousal:valence']\n",
    "\n",
    "    label = (f\"{y}\" + \"_{}\"*len(X)).format(*X) + f'_{band}' + '_all_elecs' + 'no_soz'\n",
    "\n",
    "#             'power:hemi:region:arousal',\n",
    "#             'power:hemi:region:CV'\n",
    "    rand_effect = ['subj', 'elec']\n",
    "    categorical = ['hemi', 'region']\n",
    "\n",
    "    # Drop nan data\n",
    "    model_df = model_df.dropna(subset=X)\n",
    "\n",
    "    run_model(model_df, y, X, Intx, rand_effect, rand_slopes=False, categorical=categorical,\n",
    "             chains=4, cores=4, tune=500, draws=1000, label=label, save_model_res=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11668434-e7b8-4f5c-98ab-5edced4ec279",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model (Control): Examine the effect of HFA power and word features on recall, plus hipp AP position (sontinuous)\n",
    "\n",
    "pow_df = pd.read_csv('/home1/salman.qasim/Salman_Project/FR_Emotion/final_NHB_power_data_revised.csv', memory_map=True)\n",
    "\n",
    "# rename the electrodes to properly do the elec random effects \n",
    "pow_df.elec = pow_df.apply(lambda x: f'{x.subj}_{x.elec}', axis=1)\n",
    "\n",
    "# De-mean the valence  \n",
    "pow_df['valence'] = pow_df['valence'] - 0.5\n",
    "\n",
    "# Set categorical orders \n",
    "\n",
    "cat_type = CategoricalDtype(categories=['right', 'left'], ordered=True)\n",
    "pow_df.hemi = pow_df.hemi.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['Amy', 'Hipp'], ordered=True)\n",
    "pow_df.region = pow_df.region.astype(cat_type)\n",
    "\n",
    "\n",
    "# Restrict this to gamma \n",
    "pow_df = pow_df[pow_df.band=='hfa']\n",
    "\n",
    "# Restrict this to hippocampus\n",
    "pow_df = pow_df[pow_df.region=='Hipp']\n",
    "\n",
    "# Bin longitudinal position\n",
    "APbins = [-50, -20, 10]\n",
    "pow_df['longitudinal_bin' ] = np.nan\n",
    "pow_df['longitudinal_bin'] = pd.cut(pow_df['y_coord'], \n",
    "                                            bins=APbins, \n",
    "                                            labels=['posterior', 'anterior'])\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['posterior', 'anterior'], ordered=True)\n",
    "pow_df.longitudinal_bin = pow_df.longitudinal_bin.astype(cat_type)\n",
    "\n",
    "\n",
    "pow_df = pow_df.dropna(subset=['arousal', 'valence', 'region', 'power', 'hemi', 'y_coord'])\n",
    "\n",
    "\n",
    "y = 'memory'\n",
    "X = ['arousal', 'valence', 'y_coord', 'power', 'hemi'] \n",
    "\n",
    "Intx = ['power:arousal', \n",
    "        'power:valence',\n",
    "        'power:hemi',\n",
    "        'power:y_coord',\n",
    "        'power:y_coord:hemi',\n",
    "        'power:arousal:valence',\n",
    "        'power:arousal:y_coord', \n",
    "        'power:valence:y_coord']\n",
    "\n",
    "rand_effect = ['subj', 'elec']\n",
    "categorical = ['hemi']\n",
    "\n",
    "# Drop nan data\n",
    "pow_df = pow_df.dropna(subset=X)\n",
    "\n",
    "band='hfa'\n",
    "label = (f\"{y}\" + \"_{}\"*len(X)).format(*X) + f'_{band}' + '_all_elecs' + 'no_soz'\n",
    "\n",
    "run_model(pow_df, y, X, Intx, rand_effect, rand_slopes=False, categorical=categorical,\n",
    "         chains=4, cores=4, tune=1000, draws=1000, label=label, save_model_res=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a601410-c36a-45be-ac18-0f5ac6875537",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model (Control): Examine the effect of HFA power and word features on recall, plus hipp AP position (binned)\n",
    "\n",
    "pow_df = pd.read_csv('/home1/salman.qasim/Salman_Project/FR_Emotion/final_NHB_power_data_revised.csv', memory_map=True)\n",
    "\n",
    "# rename the electrodes to properly do the elec random effects \n",
    "pow_df.elec = pow_df.apply(lambda x: f'{x.subj}_{x.elec}', axis=1)\n",
    "\n",
    "# De-mean the valence  \n",
    "pow_df['valence'] = pow_df['valence'] - 0.5\n",
    "\n",
    "# Set categorical orders \n",
    "\n",
    "cat_type = CategoricalDtype(categories=['right', 'left'], ordered=True)\n",
    "pow_df.hemi = pow_df.hemi.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['Amy', 'Hipp'], ordered=True)\n",
    "pow_df.region = pow_df.region.astype(cat_type)\n",
    "\n",
    "\n",
    "# Restrict this to gamma \n",
    "pow_df = pow_df[pow_df.band=='hfa']\n",
    "\n",
    "# Restrict this to hippocampus\n",
    "pow_df = pow_df[pow_df.region=='Hipp']\n",
    "\n",
    "# Bin longitudinal position\n",
    "APbins = [-50, -20, 10]\n",
    "pow_df['longitudinal_bin' ] = np.nan\n",
    "pow_df['longitudinal_bin'] = pd.cut(pow_df['y_coord'], \n",
    "                                            bins=APbins, \n",
    "                                            labels=['posterior', 'anterior'])\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['posterior', 'anterior'], ordered=True)\n",
    "pow_df.longitudinal_bin = pow_df.longitudinal_bin.astype(cat_type)\n",
    "\n",
    "\n",
    "pow_df = pow_df.dropna(subset=['arousal', 'valence', 'region', 'power', 'hemi', 'longitudinal_bin'])\n",
    "\n",
    "\n",
    "y = 'memory'\n",
    "X = ['arousal', 'valence', 'longitudinal_bin', 'power', 'hemi'] \n",
    "\n",
    "Intx = ['power:arousal', \n",
    "        'power:valence',\n",
    "        'power:hemi',\n",
    "        'power:longitudinal_bin',\n",
    "        'power:longitudinal_bin:hemi',\n",
    "        'power:arousal:valence',\n",
    "        'power:arousal:longitudinal_bin', \n",
    "        'power:valence:longitudinal_bin']\n",
    "\n",
    "rand_effect = ['subj', 'elec']\n",
    "categorical = ['hemi', 'longitudinal_bin']\n",
    "\n",
    "# Drop nan data\n",
    "pow_df = pow_df.dropna(subset=X)\n",
    "\n",
    "band='hfa'\n",
    "label = (f\"{y}\" + \"_{}\"*len(X)).format(*X) + f'_{band}' + '_all_elecs' + 'no_soz'\n",
    "\n",
    "run_model(pow_df, y, X, Intx, rand_effect, rand_slopes=False, categorical=categorical,\n",
    "         chains=4, cores=4, tune=1000, draws=1000, label=label, save_model_res=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2423c9b5-98e3-4175-9ceb-d39c2a356d17",
   "metadata": {},
   "source": [
    "Below you will find our one and only Bayesian linear mixed-effects model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4f73259-5077-48d9-9c38-e8064d98aacb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Model 5: effect of stimulation on post-pre HFA and also every other power band. \n",
    "\n",
    "# run separate models for stimulation \n",
    "\n",
    "stim_t_test_data = pd.read_csv('/home1/salman.qasim/Salman_Project/FR_Emotion/FR2_stim_allbands_diff_newfilt_600ms_50buff_hfa.csv')\n",
    "\n",
    "# replace the non-stim data with the word 2 stim data only\n",
    "non_stim_path = '/home1/salman.qasim/Salman_Project/FR_Emotion/FR2_nostim_allbands_diff_newfilt_600ms_50buff_stimOFf_only.csv'\n",
    "\n",
    "# Set categorical orders \n",
    "\n",
    "cat_type = CategoricalDtype(categories=['MTL', 'Amy', 'Hipp'], ordered=True)\n",
    "stim_t_test_data.stim_reg = stim_t_test_data.stim_reg.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['MTL', 'Amy', 'Hipp'], ordered=True)\n",
    "stim_t_test_data.rec_reg = stim_t_test_data.rec_reg.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['MTL', 'Amy', 'Hipp'], ordered=True)\n",
    "stim_t_test_data.rec_reg = stim_t_test_data.rec_reg.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['right', 'left'], ordered=True)\n",
    "stim_t_test_data.stim_hemi = stim_t_test_data.stim_hemi.astype(cat_type)\n",
    "\n",
    "cat_type = CategoricalDtype(categories=['right', 'left'], ordered=True)\n",
    "stim_t_test_data.rec_hemi = stim_t_test_data.rec_hemi.astype(cat_type)\n",
    "\n",
    "\n",
    "\n",
    "for sr in ['MTL', 'Hipp']:\n",
    "    for rr in ['Hipp']:\n",
    "        for y in ['theta_diff', 'alpha_diff', 'beta_diff', 'hfa_diff']:\n",
    "            model_df = stim_t_test_data[(stim_t_test_data.stim_reg==sr)]\n",
    "            model_df = model_df[(model_df.rec_reg==rr)]\n",
    "\n",
    "            X = ['stim', 'stim_hemi'] \n",
    "            Intx = ['stim',\n",
    "                    'stim:stim_hemi']\n",
    "\n",
    "            label = (f\"{y}\" + \"_{}\"*len(X)).format(*X)+f'stim_{sr}_rec{rr}'\n",
    "\n",
    "            rand_effect = 'sub'\n",
    "            categorical = X\n",
    "\n",
    "            # Drop nan data\n",
    "            model_df = model_df.dropna()\n",
    "\n",
    "            model, results = run_model(model_df, y, X, Intx, rand_effect, rand_slopes=False, categorical=categorical,\n",
    "                     chains=4, cores=4, tune=1000, draws=1500, target_accept=0.95, model_fam='gaussian', return_model=True, label = label)\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "bambi_env",
   "language": "python",
   "name": "bambi_env"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}