3_ISC_ECG.ipynb

{
    "cells": [
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Script for conducting ISC on emotion ECG data"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "###  Section 1: Basic setting & Data Import"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "%matplotlib inline\n",
                "import os\n",
                "import neurokit2 as nk\n",
                "import numpy as np\n",
                "import pandas as pd\n",
                "import seaborn as sns\n",
                "from scipy.stats import rankdata\n",
                "from sklearn.metrics import pairwise_distances\n",
                "from nltools.stats import isc\n",
                "from nltools.data import Adjacency\n",
                "from Function import import_ecg, rain_ecg, rain_hrv, mean_df2_phy, mean_df3_phy\n",
                "import matplotlib.pyplot as plt\n",
                "from jupyterthemes import jtplot\n",
                "jtplot.style(theme='grade3') \n",
                "\n",
                "### Subject list\n",
                "sub_ag = [1, 2] + list(range(4,12)) + [13,14,16,17] + list(range(19,24)) + [25,26,28,29,30,31] + list(range(33,43)) + [46] # 36\n",
                "sub_ax = list(range(1,5)) + list(range(7,17)) + [19,20,21,22,23,25,26,28,29,30,31] + list(range(33,42)) + [43,44,45,46] # 38\n",
                "sub_fe = list(range(1,13)) + [14,17,19,20,21,22,23] + [25,26,28,29,30,31] + list(range(33,41)) + [42,44,45,46] # 37\n",
                "sub_hl = [1,2] + list(range(4,18)) + [19,20,21,22,23,25,26,28,29,30,31] + list(range(33,42)) + [43,44,46] # 39\n",
                "sub_ha = list(range(1,18)) + [19,20,21,22,23,25,26,28,29,31] + list(range(33,47)) # 41"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Define the directory\n",
                "ag_dir = 'F:/1_Emotion_Data/Data/Physiology/Angry/'\n",
                "ax_dir = 'F:/1_Emotion_Data/Data/Physiology/Anxiety/'\n",
                "fe_dir = 'F:/1_Emotion_Data/Data/Physiology/Fear/'\n",
                "hl_dir = 'F:/1_Emotion_Data/Data/Physiology/Helpless/'\n",
                "ha_dir = 'F:/1_Emotion_Data/Data/Physiology/Happy/'\n",
                "results_dir = 'F:/1_Emotion_Data/Results/0_ISC/Physiology/'"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "## Import the flattened 1-D array\n",
                "ag1_meta = import_ecg(sub_ag, 'angry', '1', ag_dir)\n",
                "ag2_meta = import_ecg(sub_ag, 'angry', '2', ag_dir)\n",
                "ag3_meta = import_ecg(sub_ag, 'angry', '3', ag_dir)\n",
                "ax1_meta = import_ecg(sub_ax, 'anxiety', '1', ax_dir)\n",
                "ax2_meta = import_ecg(sub_ax, 'anxiety', '2', ax_dir)\n",
                "fe1_meta = import_ecg(sub_fe, 'fear', '1', fe_dir)\n",
                "fe2_meta = import_ecg(sub_fe, 'fear', '2', fe_dir)\n",
                "fe3_meta = import_ecg(sub_fe, 'fear', '3', fe_dir)\n",
                "hl2_meta = import_ecg(sub_hl, 'helpless', '2', hl_dir)\n",
                "hl3_meta = import_ecg(sub_hl, 'helpless', '3', hl_dir)\n",
                "ha1_meta = import_ecg(sub_ha, 'happy', '1', ha_dir)\n",
                "ha3_meta = import_ecg(sub_ha, 'happy', '3', ha_dir)"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "### Section 2: ISC for *Angry* ECG"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "* ag1: (1) 220-240, the headmaster mangae to rap the young girl\n",
                "* ag2: (1) 125 - 145, Japanese army kill armless ordinary Chinese people\n",
                "* ag3: (1) 65-85, the man beat the boy rudely  "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                " Preprocess the raw ECG and aggregate results\n",
                "* Features: Cleaned signal, Heart Rate, ECG derived Respiration(EDR), HVR_MeanNN, HRV_SDNN"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Construct a meta dict for angry results\n",
                "ag_ecg, ag_hrv = {},{}\n",
                "# Loop for all subjects\n",
                "for sub in sub_ag:\n",
                "    ecg_ag1, ecg_ag2, ecg_ag3 = ag1_meta[sub][220*2000:240*2000], ag2_meta[sub][125*2000:145*2000], ag3_meta[sub][65*2000:85*2000]\n",
                "    ag1_results, ag2_results, ag3_results = rain_ecg(ecg_ag1, 2000, 250), rain_ecg(ecg_ag2, 2000, 250), rain_ecg(ecg_ag3, 2000, 250)\n",
                "    ag1 = ag1_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    ag2 = ag2_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    ag3 = ag3_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    ag_ecg[sub] = mean_df3_phy(ag1, ag2, ag3)\n",
                "    ag_hrv[sub] = mean_df3_phy(rain_hrv(ecg_ag1, 2000, 250), rain_hrv(ecg_ag2, 2000, 250), rain_hrv(ecg_ag3, 2000, 250))\n",
                "\n",
                "## Save the angry result to .npy file\n",
                "np.save(os.path.join(results_dir, 'ag_ecg.npy'), ag_ecg) \n",
                "np.save(os.path.join(results_dir, 'ag_hrv.npy'), ag_hrv) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Extract features to conduct ISC"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "ag_ecg = np.load(os.path.join(results_dir, 'ag_ecg.npy'), allow_pickle=True).item()\n",
                "ag_hrv = np.load(os.path.join(results_dir, 'ag_hrv.npy'), allow_pickle=True).item()\n",
                "ag_hr, ag_edr, ag_hrvm, ag_hrvsd = np.zeros((5000, 36)), np.zeros((5000, 36)), np.zeros((36)), np.zeros((36))\n",
                "for sub in sub_ag:\n",
                "    ## For the index 1 for 'ECG_Rate', 2 for 'EDR'\n",
                "    hr = ag_ecg[sub].iloc[:,1].to_numpy()\n",
                "    edr = ag_ecg[sub].iloc[:,2].to_numpy()\n",
                "    mean = ag_hrv[sub].iloc[:,0]\n",
                "    sd = ag_hrv[sub].iloc[:,1]\n",
                "    hr = np.ndarray.flatten(hr)\n",
                "    edr = np.ndarray.flatten(edr)\n",
                "    idx = sub_ag.index(sub)\n",
                "    ag_hr[:,idx] = hr\n",
                "    ag_edr[:,idx] = edr\n",
                "    ag_hrvm[idx] = mean\n",
                "    ag_hrvsd[idx] = sd"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Calculate ISC on Angry ECG features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Calculate ISC on angry heart rate\n",
                "ag_hr_similarity, ag_hr_stats = isc(ag_hr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "ag_hr_similarity = ag_hr_similarity.squareform()\n",
                "ag_edr_similarity, ag_edr_stats = isc(ag_edr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "ag_edr_similarity = ag_edr_similarity.squareform()\n",
                "\n",
                "# Plot\n",
                "mask =np.zeros_like(ag_hr_similarity)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "sns.heatmap(ag_hr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax1.set_title(\"Angry HR ISC\", fontsize=25, fontweight='bold')\n",
                "sns.heatmap(ag_edr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax2.set_title(\"Angry EDR ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "# Save the result\n",
                "np.save(os.path.join(results_dir, 'ag_hr_isc.npy'), ag_hr_similarity) \n",
                "np.save(os.path.join(results_dir, 'ag_edr_isc.npy'), ag_edr_similarity) \n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Calculate ISC on Angry HRV features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Set the paprameters and import data\n",
                "n_subs = 36\n",
                "hrvm_rank, hrvsd_rank = rankdata(ag_hrvm), rankdata(ag_hrvsd) # explicity convert the raw scores to ranks\n",
                "hrvm_simi, hrvsd_simi = np.zeros((n_subs, n_subs)), np.zeros((n_subs, n_subs))\n",
                "\n",
                "# Calculate the ann isc\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvm_rank[i], hrvm_rank[j]])/n_subs\n",
                "            hrvm_simi[i,j] = sim_ij\n",
                "            hrvm_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvm_simi[i,j] = 1\n",
                "\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvsd_rank[i], hrvsd_rank[j]])/n_subs\n",
                "            hrvsd_simi[i,j] = sim_ij\n",
                "            hrvsd_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvsd_simi[i,j] = 1\n",
                "\n",
                "# Plot the isc matrix\n",
                "mask =np.zeros_like(hrvm_simi)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "# fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "# sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax1.set_title(\"Angry HRV mean ISC\", fontsize=25, fontweight='bold')\n",
                "# sns.heatmap(hrvsd_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax2.set_title(\"Angry HRV SD ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "plt.figure(figsize=(20,10))\n",
                "sns.heatmap(hrvm_simi, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "#plt.title('Angry HRV Intersubject Similarity', fontsize=25, fontweight='bold')\n",
                "# Save the result\n",
                "# np.save(os.path.join(results_dir, 'ag_hrvm_isc.npy'), hrvm_simi) \n",
                "# np.save(os.path.join(results_dir, 'ag_hrvsd_isc.npy'), hrvsd_simi) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Section 3: ISC for *Anxiety* EEG"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "* ax1: 135 - 155, The bomb is about to explode\n",
                "* ax2: 135 - 155, The girl ask for the man to open the door\n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                " Preprocess the raw ECG and aggregate results\n",
                "* Features: Cleaned signal, Heart Rate, ECG derived Respiration(EDR), HVR_MeanNN, HRV_SDNN"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Construct a meta dict for anxiety results\n",
                "ax_ecg, ax_hrv = {},{}\n",
                "# Loop for all subjects\n",
                "for sub in sub_ax:\n",
                "    ecg_ax1, ecg_ax2 = ax1_meta[sub][135*2000:155*2000], ax2_meta[sub][135*2000:155*2000]\n",
                "    ax1_results, ax2_results = rain_ecg(ecg_ax1, 2000, 250), rain_ecg(ecg_ax2, 2000, 250)\n",
                "    ax1 = ax1_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    ax2 = ax2_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    ax_ecg[sub] = mean_df2_phy(ax1, ax2)\n",
                "    ax_hrv[sub] = mean_df2_phy(rain_hrv(ecg_ax1, 2000, 250), rain_hrv(ecg_ax2, 2000, 250))\n",
                "\n",
                "## Save the anxiety result to .npy file\n",
                "np.save(os.path.join(results_dir, 'ax_ecg.npy'), ax_ecg) \n",
                "np.save(os.path.join(results_dir, 'ax_hrv.npy'), ax_hrv) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Extract features to conduct ISC"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "ax_ecg = np.load(os.path.join(results_dir, 'ax_ecg.npy'), allow_pickle=True).item()\n",
                "ax_hrv = np.load(os.path.join(results_dir, 'ax_hrv.npy'), allow_pickle=True).item()\n",
                "ax_hr, ax_edr, ax_hrvm, ax_hrvsd = np.zeros((5000, 38)), np.zeros((5000, 38)), np.zeros((38)), np.zeros((38))\n",
                "for sub in sub_ax:\n",
                "    ## For the index 1 for 'ECG_Rate', 2 for 'EDR'\n",
                "    hr = ax_ecg[sub].iloc[:,1].to_numpy()\n",
                "    edr = ax_ecg[sub].iloc[:,2].to_numpy()\n",
                "    mean = ax_hrv[sub].iloc[:,0]\n",
                "    sd = ax_hrv[sub].iloc[:,1]\n",
                "    hr = np.ndarray.flatten(hr)\n",
                "    edr = np.ndarray.flatten(edr)\n",
                "    idx = sub_ax.index(sub)\n",
                "    ax_hr[:,idx] = hr\n",
                "    ax_edr[:,idx] = edr\n",
                "    ax_hrvm[idx] = mean\n",
                "    ax_hrvsd[idx] = sd"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                " Calculate the ISC on anxiety ECG features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Calculate ISC on angry heart rate\n",
                "ax_hr_similarity, ax_hr_stats = isc(ax_hr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "ax_hr_similarity = ax_hr_similarity.squareform()\n",
                "ax_edr_similarity, ax_edr_stats = isc(ax_edr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "ax_edr_similarity = ax_edr_similarity.squareform()\n",
                "\n",
                "# Plot\n",
                "mask =np.zeros_like(ax_hr_similarity)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "sns.heatmap(ax_hr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax1.set_title(\"Anxiety HR ISC\", fontsize=25, fontweight='bold')\n",
                "sns.heatmap(ax_edr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax2.set_title(\"Anxiety EDR ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "# Save the result\n",
                "np.save(os.path.join(results_dir, 'ax_hr_isc.npy'), ax_hr_similarity) \n",
                "np.save(os.path.join(results_dir, 'ax_edr_isc.npy'), ax_edr_similarity) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Calculate ISC on Anxiety HRV features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Set the paprameters and import data\n",
                "n_subs = 38\n",
                "hrvm_rank, hrvsd_rank = rankdata(ax_hrvm), rankdata(ax_hrvsd) # explicity convert the raw scores to ranks\n",
                "hrvm_simi, hrvsd_simi = np.zeros((n_subs, n_subs)), np.zeros((n_subs, n_subs))\n",
                "\n",
                "# Calculate the ann isc\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvm_rank[i], hrvm_rank[j]])/n_subs\n",
                "            hrvm_simi[i,j] = sim_ij\n",
                "            hrvm_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvm_simi[i,j] = 1\n",
                "\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvsd_rank[i], hrvsd_rank[j]])/n_subs\n",
                "            hrvsd_simi[i,j] = sim_ij\n",
                "            hrvsd_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvsd_simi[i,j] = 1\n",
                "\n",
                "# Plot the isc matrix\n",
                "mask =np.zeros_like(hrvm_simi)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "# fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "# sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax1.set_title(\"Anxiety HRV mean ISC\", fontsize=25, fontweight='bold')\n",
                "# sns.heatmap(hrvsd_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax2.set_title(\"Anxiety HRV SD ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "plt.figure(figsize=(20,10))\n",
                "sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "plt.title('Anxiety HRV Intersubject Similarity', fontsize=25, fontweight='bold')\n",
                "# Save the result\n",
                "# np.save(os.path.join(results_dir, 'ax_hrvm_isc.npy'), hrvm_simi) \n",
                "# np.save(os.path.join(results_dir, 'ax_hrvsd_isc.npy'), hrvsd_simi) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Section 4: ISC for *Fear* ECG"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "* fe1: 110 - 130\n",
                "* fe2: 30 - 50\n",
                "* fe3: 45 - 65"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Preprocess the raw ECG and aggregate results\n",
                "* Features: Cleaned signal, Heart Rate, ECG derived Respiration(EDR), HVR_MeanNN, HRV_SDNN"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Construct a meta dict for angry results\n",
                "fe_ecg, fe_hrv = {},{}\n",
                "# Loop for all subjects\n",
                "for sub in sub_fe:\n",
                "    ecg_fe1, ecg_fe2, ecg_fe3 = fe1_meta[sub][110*2000:130*2000], fe2_meta[sub][30*2000:50*2000], fe3_meta[sub][45*2000:65*2000]\n",
                "    fe1_results, fe2_results, fe3_results = rain_ecg(ecg_fe1, 2000, 250), rain_ecg(ecg_fe2, 2000, 250), rain_ecg(ecg_fe3, 2000, 250)\n",
                "    fe1 = fe1_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    fe2 = fe2_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    fe3 = fe3_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    fe_ecg[sub] = mean_df3_phy(fe1, fe2, fe3)\n",
                "    fe_hrv[sub] = mean_df3_phy(rain_hrv(ecg_fe1, 2000, 250), rain_hrv(ecg_fe2, 2000, 250), rain_hrv(ecg_fe3, 2000, 250))\n",
                "\n",
                "## Save the fear result to .npy file\n",
                "np.save(os.path.join(results_dir, 'fe_ecg.npy'), fe_ecg) \n",
                "np.save(os.path.join(results_dir, 'fe_hrv.npy'), fe_hrv) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Extract features to conduct ISC"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "fe_ecg = np.load(os.path.join(results_dir, 'fe_ecg.npy'), allow_pickle=True).item()\n",
                "fe_hrv = np.load(os.path.join(results_dir, 'fe_hrv.npy'), allow_pickle=True).item()\n",
                "fe_hr, fe_edr, fe_hrvm, fe_hrvsd = np.zeros((5000, 37)), np.zeros((5000, 37)), np.zeros((37)), np.zeros((37))\n",
                "for sub in sub_fe:\n",
                "    hr = fe_ecg[sub].iloc[:,1].to_numpy()\n",
                "    edr = fe_ecg[sub].iloc[:,2].to_numpy()\n",
                "    mean = fe_hrv[sub].iloc[:,0]\n",
                "    sd = fe_hrv[sub].iloc[:,1]\n",
                "    hr = np.ndarray.flatten(hr)\n",
                "    edr = np.ndarray.flatten(edr)\n",
                "    idx = sub_fe.index(sub)\n",
                "    fe_hr[:,idx] = hr\n",
                "    fe_edr[:,idx] = edr\n",
                "    fe_hrvm[idx] = mean\n",
                "    fe_hrvsd[idx] = sd"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                " Calculate the ISC on fear ECG features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Calculate ISC on angry heart rate\n",
                "fe_hr_similarity, fe_hr_stats = isc(fe_hr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "fe_hr_similarity = fe_hr_similarity.squareform()\n",
                "fe_edr_similarity, fe_edr_stats = isc(fe_edr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "fe_edr_similarity = fe_edr_similarity.squareform()\n",
                "\n",
                "# Plot\n",
                "mask =np.zeros_like(fe_hr_similarity)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "sns.heatmap(fe_hr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax1.set_title(\"Fear HR ISC\", fontsize=25, fontweight='bold')\n",
                "sns.heatmap(fe_edr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax2.set_title(\"Fear EDR ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "# Save the result\n",
                "np.save(os.path.join(results_dir, 'fe_hr_isc.npy'), fe_hr_similarity) \n",
                "np.save(os.path.join(results_dir, 'fe_edr_isc.npy'), fe_hr_similarity) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                " Conduct ISC on fear HRV features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Set the paprameters and import data\n",
                "n_subs = 37\n",
                "hrvm_rank, hrvsd_rank = rankdata(fe_hrvm), rankdata(fe_hrvsd) # explicity convert the raw scores to ranks\n",
                "hrvm_simi, hrvsd_simi = np.zeros((n_subs, n_subs)), np.zeros((n_subs, n_subs))\n",
                "\n",
                "# Calculate the ann isc\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvm_rank[i], hrvm_rank[j]])/n_subs\n",
                "            hrvm_simi[i,j] = sim_ij\n",
                "            hrvm_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvm_simi[i,j] = 1\n",
                "\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvsd_rank[i], hrvsd_rank[j]])/n_subs\n",
                "            hrvsd_simi[i,j] = sim_ij\n",
                "            hrvsd_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvsd_simi[i,j] = 1\n",
                "\n",
                "# Plot the isc matrix\n",
                "mask =np.zeros_like(hrvm_simi)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "# fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "# sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax1.set_title(\"Fear HRV mean ISC\", fontsize=25, fontweight='bold')\n",
                "# sns.heatmap(hrvsd_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax2.set_title(\"Fear HRV SD ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "plt.figure(figsize=(20,10))\n",
                "sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "plt.title('Fear HRV Intersubject Similarity', fontsize=25, fontweight='bold')\n",
                "# Save the result\n",
                "# np.save(os.path.join(results_dir, 'fe_hrvm_isc.npy'), hrvm_simi) \n",
                "# np.save(os.path.join(results_dir, 'fe_hrvsd_isc.npy'), hrvsd_simi) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Section 5: ISC for *Helpless* ECG"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "* hl2: 115 - 135, The astronaut has run out of oxygen and she managed to broke her arm, \n",
                "* hl3: 50 - 70, The woman ask for the man for medicine "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Preprocess the raw ECG and aggregate results\n",
                "* Features: Cleaned signal, Heart Rate, ECG derived Respiration(EDR), HVR_MeanNN, HRV_SDNN"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Construct a meta dict for anxiety results\n",
                "hl_ecg, hl_hrv = {},{}\n",
                "# Loop for all subjects\n",
                "for sub in sub_hl:\n",
                "    ecg_hl2, ecg_hl3 = hl2_meta[sub][115*2000:135*2000], hl3_meta[sub][50*2000:70*2000]\n",
                "    hl2_results, hl3_results = rain_ecg(ecg_hl2, 2000, 250), rain_ecg(ecg_hl3, 2000, 250)\n",
                "    hl2 = hl2_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    hl3 = hl3_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    hl_ecg[sub] = mean_df2_phy(hl2, hl3)\n",
                "    hl_hrv[sub] = mean_df2_phy(rain_hrv(ecg_hl2, 2000, 250), rain_hrv(ecg_hl3, 2000, 250))\n",
                "\n",
                "## Save the anxiety result to .npy file\n",
                "np.save(os.path.join(results_dir, 'hl_ecg.npy'), hl_ecg) \n",
                "np.save(os.path.join(results_dir, 'hl_hrv.npy'), hl_hrv) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Extract features to conduct ISC"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "hl_ecg = np.load(os.path.join(results_dir, 'hl_ecg.npy'), allow_pickle=True).item()\n",
                "hl_hrv = np.load(os.path.join(results_dir, 'hl_hrv.npy'), allow_pickle=True).item()\n",
                "hl_hr, hl_edr, hl_hrvm, hl_hrvsd = np.zeros((5000, 39)), np.zeros((5000, 39)), np.zeros((39)), np.zeros((39))\n",
                "for sub in sub_hl:\n",
                "    ## For the index 1 for 'ECG_Rate', 2 for 'EDR'\n",
                "    hr = hl_ecg[sub].iloc[:,1].to_numpy()\n",
                "    edr = hl_ecg[sub].iloc[:,2].to_numpy()\n",
                "    mean = hl_hrv[sub].iloc[:,0]\n",
                "    sd = hl_hrv[sub].iloc[:,1]\n",
                "    hr = np.ndarray.flatten(hr)\n",
                "    edr = np.ndarray.flatten(edr)\n",
                "    idx = sub_hl.index(sub)\n",
                "    hl_hr[:,idx] = hr\n",
                "    hl_edr[:,idx] = edr\n",
                "    hl_hrvm[idx] = mean\n",
                "    hl_hrvsd[idx] = sd"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Conduct ISC on helpless ECG features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Calculate ISC on angry heart rate\n",
                "hl_hr_similarity, hl_hr_stats = isc(hl_hr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "hl_hr_similarity = hl_hr_similarity.squareform()\n",
                "hl_edr_similarity, hl_edr_stats = isc(hl_edr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "hl_edr_similarity = hl_edr_similarity.squareform()\n",
                "\n",
                "# Plot\n",
                "mask =np.zeros_like(hl_hr_similarity)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "sns.heatmap(hl_hr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax1.set_title(\"Helpless HR Intersubject Similarity\", fontsize=20, fontweight='bold')\n",
                "sns.heatmap(hl_edr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax2.set_title(\"Helpless EDR Intersubject Similarity\", fontsize=20, fontweight='bold')\n",
                "\n",
                "# Save the result\n",
                "# np.save(os.path.join(results_dir, 'hl_hr_isc.npy'), hl_hr_similarity) \n",
                "# np.save(os.path.join(results_dir, 'hl_edr_isc.npy'), hl_edr_similarity) \n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Conduct ISC on helpless HRV features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Set the paprameters and import data\n",
                "n_subs = 39\n",
                "hrvm_rank, hrvsd_rank = rankdata(hl_hrvm), rankdata(hl_hrvsd) # explicity convert the raw scores to ranks\n",
                "hrvm_simi, hrvsd_simi = np.zeros((n_subs, n_subs)), np.zeros((n_subs, n_subs))\n",
                "\n",
                "# Calculate the ann isc\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvm_rank[i], hrvm_rank[j]])/n_subs\n",
                "            hrvm_simi[i,j] = sim_ij\n",
                "            hrvm_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvm_simi[i,j] = 1\n",
                "\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvsd_rank[i], hrvsd_rank[j]])/n_subs\n",
                "            hrvsd_simi[i,j] = sim_ij\n",
                "            hrvsd_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvsd_simi[i,j] = 1\n",
                "\n",
                "# Plot the isc matrix\n",
                "mask =np.zeros_like(hrvm_simi)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "# fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "# sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax1.set_title(\"Fear HRV mean ISC\", fontsize=25, fontweight='bold')\n",
                "# sns.heatmap(hrvsd_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax2.set_title(\"Fear HRV SD ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "plt.figure(figsize=(20,10))\n",
                "sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "plt.title('Helpless HRV Intersubject Similarity', fontsize=25, fontweight='bold')\n",
                "# Save the result\n",
                "# np.save(os.path.join(results_dir, 'hl_hrvm_isc.npy'), hrvm_simi) \n",
                "# np.save(os.path.join(results_dir, 'hl_hrvsd_isc.npy'), hrvsd_simi) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "## Section 6: ISC for *happy* ECG"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "* ha1: 180:200\n",
                "* ha3: 5-25"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Process raw ECG and aggregate results\n",
                "* Features: Cleaned signal, Heart Rate, ECG derived Respiration(EDR), HVR_MeanNN, HRV_SDNN"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Construct a meta dict for happy results\n",
                "ha_ecg, ha_hrv = {},{}\n",
                "# Loop for all subjects\n",
                "for sub in sub_ha:\n",
                "    ecg_ha1, ecg_ha3 = ha1_meta[sub][180*2000:200*2000], ha3_meta[sub][5*2000:25*2000]\n",
                "    ha1_results, ha3_results = rain_ecg(ecg_ha1, 2000, 250), rain_ecg(ecg_ha3, 2000, 250)\n",
                "    ha1 = ha1_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    ha3 = ha3_results.loc[:,['ECG_Clean', 'ECG_Rate', 'ECG_EDR']]\n",
                "    ha_ecg[sub] = mean_df2_phy(ha1, hl3)\n",
                "    ha_hrv[sub] = mean_df2_phy(rain_hrv(ecg_ha1, 2000, 250), rain_hrv(ecg_ha3, 2000, 250))\n",
                "\n",
                "## Save the happy result to .npy file\n",
                "np.save(os.path.join(results_dir, 'ha_ecg.npy'), ha_ecg) \n",
                "np.save(os.path.join(results_dir, 'ha_hrv.npy'), ha_hrv) "
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Extract features to conduct ISC"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "ha_ecg = np.load(os.path.join(results_dir, 'ha_ecg.npy'), allow_pickle=True).item()\n",
                "ha_hrv = np.load(os.path.join(results_dir, 'ha_hrv.npy'), allow_pickle=True).item()\n",
                "ha_hr, ha_edr, ha_hrvm, ha_hrvsd = np.zeros((5000, 41)), np.zeros((5000, 41)), np.zeros((41)), np.zeros((41))\n",
                "for sub in sub_ha:\n",
                "    ## For the index 1 for 'ECG_Rate', 2 for 'EDR'\n",
                "    hr = ha_ecg[sub].iloc[:,1].to_numpy()\n",
                "    edr = ha_ecg[sub].iloc[:,2].to_numpy()\n",
                "    mean = ha_hrv[sub].iloc[:,0]\n",
                "    sd = ha_hrv[sub].iloc[:,1]\n",
                "    hr = np.ndarray.flatten(hr)\n",
                "    edr = np.ndarray.flatten(edr)\n",
                "    idx = sub_ha.index(sub)\n",
                "    ha_hr[:,idx] = hr\n",
                "    ha_edr[:,idx] = edr\n",
                "    ha_hrvm[idx] = mean\n",
                "    ha_hrvsd[idx] = sd"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Conduct ISC on happy ECG features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Calculate ISC on angry heart rate\n",
                "ha_hr_similarity, ha_hr_stats = isc(ha_hr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "ha_hr_similarity = ha_hr_similarity.squareform()\n",
                "ha_edr_similarity, ha_edr_stats = isc(ha_edr, n_bootstraps=5000, metric='median', method='bootstrap', n_jobs=-1)\n",
                "ha_edr_similarity = ha_edr_similarity.squareform()\n",
                "\n",
                "# Plot\n",
                "mask =np.zeros_like(ha_hr_similarity)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "sns.heatmap(ha_hr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax1.set_title(\"Happy HR ISC\", fontsize=25, fontweight='bold')\n",
                "sns.heatmap(ha_edr_similarity, square=True, cmap='RdBu_r', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "ax2.set_title(\"Happy EDR ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "# Save the result\n",
                "np.save(os.path.join(results_dir, 'ha_hr_isc.npy'), ha_hr_similarity) \n",
                "np.save(os.path.join(results_dir, 'ha_edr_isc.npy'), ha_edr_similarity) \n"
            ]
        },
        {
            "cell_type": "markdown",
            "metadata": {},
            "source": [
                "Conduct ISC on happy HRV features"
            ]
        },
        {
            "cell_type": "code",
            "execution_count": null,
            "metadata": {},
            "outputs": [],
            "source": [
                "# Set the paprameters and import data\n",
                "n_subs = 41\n",
                "hrvm_rank, hrvsd_rank = rankdata(ha_hrvm), rankdata(ha_hrvsd) # explicity convert the raw scores to ranks\n",
                "hrvm_simi, hrvsd_simi = np.zeros((n_subs, n_subs)), np.zeros((n_subs, n_subs))\n",
                "\n",
                "# Calculate the ann isc\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvm_rank[i], hrvm_rank[j]])/n_subs\n",
                "            hrvm_simi[i,j] = sim_ij\n",
                "            hrvm_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvm_simi[i,j] = 1\n",
                "\n",
                "for i in range(n_subs):\n",
                "    for j in range(n_subs):\n",
                "        if i < j:\n",
                "            sim_ij = np.mean([hrvsd_rank[i], hrvsd_rank[j]])/n_subs\n",
                "            hrvsd_simi[i,j] = sim_ij\n",
                "            hrvsd_simi[j,i] = sim_ij\n",
                "        elif i==j:\n",
                "            hrvsd_simi[i,j] = 1\n",
                "\n",
                "# Plot the isc matrix\n",
                "mask =np.zeros_like(hrvm_simi)\n",
                "mask[np.triu_indices_from(mask)] = True\n",
                "# fig, (ax1, ax2) = plt.subplots(1,2, figsize=(20,8))\n",
                "# sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax1.set_title(\"Fear HRV mean ISC\", fontsize=25, fontweight='bold')\n",
                "# sns.heatmap(hrvsd_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, ax=ax2, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "# ax2.set_title(\"Fear HRV SD ISC\", fontsize=25, fontweight='bold')\n",
                "\n",
                "plt.figure(figsize=(20,10))\n",
                "sns.heatmap(hrvm_simi, square=True, cmap='Reds', mask=mask, linewidths=0.1, xticklabels=False, yticklabels=False, cbar_kws={'label': 'similarity', \"shrink\": 0.8})\n",
                "plt.title('Happy HRV Intersubject Similarity', fontsize=25, fontweight='bold')\n",
                "# Save the result\n",
                "# np.save(os.path.join(results_dir, 'ha_hrvm_isc.npy'), hrvm_simi) \n",
                "# np.save(os.path.join(results_dir, 'ha_hrvsd_isc.npy'), hrvsd_simi) "
            ]
        }
    ],
    "metadata": {
        "interpreter": {
            "hash": "0637cefd005ce32f53abe336a7c6e0206dbf0dc1ed7f0aa326c4767b89ff415c"
        },
        "kernelspec": {
            "display_name": "Python 3.8.11 64-bit ('EEG': conda)",
            "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.11"
        },
        "orig_nbformat": 4
    },
    "nbformat": 4,
    "nbformat_minor": 2
}