second-level/3_run_second_level_rsa.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Run second-level RSA on the cross-validated Feedback History RSA clusters\n",
    "\n",
    "## Goal\n",
    "\n",
    "We first found the two important cortical ROIs encoding feedback histories, the MiOG and IFG clusters. We thought that these clusters might represent other task-relevant information. Therefore, we conducted second-level RSA on these clusters to identify whether the feedback history clusters could encode other task-related aspects, especially the spatial information."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "\n",
    "from pathlib import Path\n",
    "from typing import Literal\n",
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from nilearn.image import load_img\n",
    "\n",
    "from rsatoolbox.rdm import compare_rho_a\n",
    "from scipy.spatial.distance import pdist\n",
    "from scipy.stats import (\n",
    "    permutation_test,\n",
    "    ttest_1samp,\n",
    ")\n",
    "from sklearn.linear_model import LinearRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Paths\n",
    "\n",
    "# Cluster mask directory path\n",
    "GIT_ROOT_LINES = !git rev-parse --show-toplevel\n",
    "GIT_ROOT = Path(GIT_ROOT_LINES[0])\n",
    "feedback_history_cluster_mask_dir = (\n",
    "    GIT_ROOT\n",
    "    / \"second-level\"\n",
    "    / \"output\"\n",
    "    / \"cross_validated_cluster_mask\"\n",
    "    / \"recent_2_and_recent_3_trial\"\n",
    ")\n",
    "\n",
    "# Feedback history RSA cluster masks\n",
    "feedback_history_cluster_mask_list = [\n",
    "    (\n",
    "        \"1_MiOG\",\n",
    "        load_img(\n",
    "            feedback_history_cluster_mask_dir\n",
    "            / \"feedback_history_cluster_1_MiOG_mask.nii\"\n",
    "        )\n",
    "        .get_fdata()\n",
    "        .astype(bool),\n",
    "    ),\n",
    "    (\n",
    "        \"2_IFG\",\n",
    "        load_img(\n",
    "            feedback_history_cluster_mask_dir\n",
    "            / \"feedback_history_cluster_2_IFG_mask.nii\"\n",
    "        )\n",
    "        .get_fdata()\n",
    "        .astype(bool),\n",
    "    ),\n",
    "]\n",
    "\n",
    "# Data paths\n",
    "discovery_first_level_output_path = Path(\n",
    "    \"/users/jss/data/ISL_photographer/2022/bids/derivatives/first-level/\"  # Please update accordingly\n",
    ")\n",
    "validation_first_level_output_path = Path(\n",
    "    \"/users/jss/data/ISL_photographer/2023/bids/derivatives/first-level\"  # Please update accordingly\n",
    ")\n",
    "subgroup_first_level_output_path_dict = {\n",
    "    \"discovery\": discovery_first_level_output_path,\n",
    "    \"validation\": validation_first_level_output_path,\n",
    "}\n",
    "\n",
    "# Result dataframe output path\n",
    "rsa_result_output_path = GIT_ROOT / \"second-level\" / \"output\" / \"second_level_rsa\"\n",
    "\n",
    "# Subject ID lists\n",
    "discovery_subject_list = [\n",
    "    \"sub-09\",\n",
    "    \"sub-10\",\n",
    "    \"sub-16\",\n",
    "    \"sub-17\",\n",
    "    \"sub-19\",\n",
    "    \"sub-20\",\n",
    "    \"sub-21\",\n",
    "    \"sub-22\",\n",
    "    \"sub-23\",\n",
    "    \"sub-24\",\n",
    "    \"sub-26\",\n",
    "    \"sub-27\",\n",
    "    \"sub-28\",\n",
    "    \"sub-29\",\n",
    "    \"sub-30\",\n",
    "    \"sub-31\",\n",
    "]\n",
    "\n",
    "validation_subject_list = [\n",
    "    \"sub-32\",\n",
    "    \"sub-33\",\n",
    "    \"sub-34\",\n",
    "    \"sub-35\",\n",
    "    \"sub-36\",\n",
    "    \"sub-37\",\n",
    "    \"sub-38\",\n",
    "    \"sub-39\",\n",
    "    \"sub-41\",\n",
    "    \"sub-42\",\n",
    "    \"sub-45\",\n",
    "    \"sub-47\",\n",
    "    \"sub-48\",\n",
    "    \"sub-49\",\n",
    "    \"sub-50\",\n",
    "    \"sub-51\",\n",
    "]\n",
    "\n",
    "subgroup_subject_list_dict = {\n",
    "    \"discovery\": discovery_subject_list,\n",
    "    \"validation\": validation_subject_list,\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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>subgroup</th>\n",
       "      <th>subject_id</th>\n",
       "      <th>run_id</th>\n",
       "      <th>city</th>\n",
       "      <th>trial</th>\n",
       "      <th>exploration_time</th>\n",
       "      <th>exploration_distance</th>\n",
       "      <th>capture_action</th>\n",
       "      <th>capture_lat</th>\n",
       "      <th>capture_lon</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2022</td>\n",
       "      <td>sub-09</td>\n",
       "      <td>run-01</td>\n",
       "      <td>Boston</td>\n",
       "      <td>1</td>\n",
       "      <td>1.710</td>\n",
       "      <td>0.0048</td>\n",
       "      <td>capture</td>\n",
       "      <td>42.359721</td>\n",
       "      <td>-71.059837</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2022</td>\n",
       "      <td>sub-09</td>\n",
       "      <td>run-01</td>\n",
       "      <td>Boston</td>\n",
       "      <td>2</td>\n",
       "      <td>27.062</td>\n",
       "      <td>0.1169</td>\n",
       "      <td>capture</td>\n",
       "      <td>42.360616</td>\n",
       "      <td>-71.060584</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2022</td>\n",
       "      <td>sub-09</td>\n",
       "      <td>run-01</td>\n",
       "      <td>Boston</td>\n",
       "      <td>3</td>\n",
       "      <td>30.895</td>\n",
       "      <td>0.1266</td>\n",
       "      <td>capture</td>\n",
       "      <td>42.359631</td>\n",
       "      <td>-71.059813</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2022</td>\n",
       "      <td>sub-09</td>\n",
       "      <td>run-01</td>\n",
       "      <td>Boston</td>\n",
       "      <td>4</td>\n",
       "      <td>15.715</td>\n",
       "      <td>0.0477</td>\n",
       "      <td>capture</td>\n",
       "      <td>42.359269</td>\n",
       "      <td>-71.059501</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2022</td>\n",
       "      <td>sub-09</td>\n",
       "      <td>run-01</td>\n",
       "      <td>Boston</td>\n",
       "      <td>5</td>\n",
       "      <td>16.853</td>\n",
       "      <td>0.0103</td>\n",
       "      <td>capture</td>\n",
       "      <td>42.359315</td>\n",
       "      <td>-71.059610</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   subgroup subject_id  run_id    city  trial  exploration_time  \\\n",
       "0      2022     sub-09  run-01  Boston      1             1.710   \n",
       "1      2022     sub-09  run-01  Boston      2            27.062   \n",
       "2      2022     sub-09  run-01  Boston      3            30.895   \n",
       "3      2022     sub-09  run-01  Boston      4            15.715   \n",
       "4      2022     sub-09  run-01  Boston      5            16.853   \n",
       "\n",
       "   exploration_distance capture_action  capture_lat  capture_lon  \n",
       "0                0.0048        capture    42.359721   -71.059837  \n",
       "1                0.1169        capture    42.360616   -71.060584  \n",
       "2                0.1266        capture    42.359631   -71.059813  \n",
       "3                0.0477        capture    42.359269   -71.059501  \n",
       "4                0.0103        capture    42.359315   -71.059610  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Behavioral data\n",
    "# Preprocessed exploration-related behavioral data\n",
    "# Note: For subgroup IDs, we used \"2022\" for the Discovery group and \"2023\" for the Validation group\n",
    "# `exploration_time`: Elapsed time from the exploration start to capture\n",
    "# `exploration_distance`: Distance between the start and captue points\n",
    "# `capture_action`: Whether the participant captured a scene within 40 s time limit ('capture') or\n",
    "#                   failed to capture within the time limit ('capture_failed')\n",
    "# `capture_lat` / `capture_lon`: Latitude/Longitude for each capture location\n",
    "\n",
    "exploration_info_df = pd.read_csv(\n",
    "    GIT_ROOT / \"second-level\" / \"data\" / \"photographer_exploration_info.csv\"\n",
    ")\n",
    "exploration_info_df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_model_rdm_list = [\n",
    "    # feedback history model names\n",
    "    \"current_trial\",\n",
    "    \"one_back_trial\",\n",
    "    \"two_back_trial\",\n",
    "    \"recent_2_trial\",\n",
    "    \"recent_3_trial\",\n",
    "    \"previous_2_trial\",\n",
    "    # exploration model names\n",
    "    \"exploration_time\",\n",
    "    \"exploration_distance\",\n",
    "    \"capture_distance\",\n",
    "]\n",
    "\n",
    "\n",
    "# Haversine formula to compute a distance between two coordinates on the earth\n",
    "# REF: https://www.movable-type.co.uk/scripts/latlong.html\n",
    "def haversine(coord_1: tuple[float, float], coord_2: tuple[float, float]):\n",
    "    R = 6371.0088\n",
    "\n",
    "    lat_1 = np.radians(coord_1[0])\n",
    "    lon_1 = np.radians(coord_1[1])\n",
    "\n",
    "    lat_2 = np.radians(coord_2[0])\n",
    "    lon_2 = np.radians(coord_2[1])\n",
    "\n",
    "    d_lat = lat_2 - lat_1\n",
    "    d_lon = lon_2 - lon_1\n",
    "\n",
    "    a = (np.sin(d_lat / 2) * np.sin(d_lat / 2)) + np.cos(lat_1) * np.cos(lat_2) * (\n",
    "        np.sin(d_lon / 2) * np.sin(d_lon / 2)\n",
    "    )\n",
    "    c = 2 * np.arctan2(np.sqrt(a), np.sqrt(1 - a))\n",
    "\n",
    "    distance = R * c\n",
    "    return distance\n",
    "\n",
    "\n",
    "# A helper function to load feedback history/exploration model RDMs\n",
    "def load_model_rdm_vector(\n",
    "    rdm_name: str,\n",
    "    subgroup: Literal[\"discovery\", \"validation\"],\n",
    "    subject_id: str,\n",
    "    run_id: str,\n",
    "    exploration_info_df=exploration_info_df,\n",
    "):\n",
    "    if rdm_name not in all_model_rdm_list:\n",
    "        raise ValueError\n",
    "\n",
    "    if rdm_name in all_model_rdm_list[:6]:  # Feedback history model RDMs\n",
    "        subgroup_first_level_output_path = subgroup_first_level_output_path_dict[\n",
    "            subgroup\n",
    "        ]\n",
    "\n",
    "        subject_model_rdm_dir = (\n",
    "            subgroup_first_level_output_path\n",
    "            / subject_id\n",
    "            / \"rsa_model_rdm\"\n",
    "            / \"feedback_model\"\n",
    "        )\n",
    "\n",
    "        model_rdm_path = (\n",
    "            subject_model_rdm_dir\n",
    "            / f\"{subject_id}_{run_id}_task-photographer_{rdm_name}_vector.npy\"\n",
    "        )\n",
    "        return np.load(model_rdm_path)\n",
    "\n",
    "    elif rdm_name in all_model_rdm_list[6:]:  # Exploration model RDMs\n",
    "        run_exploration_df = exploration_info_df[\n",
    "            (exploration_info_df[\"subject_id\"] == subject_id)\n",
    "            & (exploration_info_df[\"run_id\"] == run_id)\n",
    "        ]\n",
    "\n",
    "        if (rdm_name == \"exploration_time\") or (rdm_name == \"exploration_distance\"):\n",
    "            run_exploration_data_list = run_exploration_df[rdm_name].tolist()\n",
    "            run_exploration_data_vector = [\n",
    "                [t] for t in run_exploration_data_list[2:]\n",
    "            ]  # Using Trial 3 to 8\n",
    "            return pdist(run_exploration_data_vector)  # Euclidean distance\n",
    "\n",
    "        elif rdm_name == \"capture_distance\":\n",
    "            run_exploration_coord_vector = [\n",
    "                (lat, lon)\n",
    "                for lat, lon in zip(\n",
    "                    run_exploration_df[\"capture_lat\"][2:],\n",
    "                    run_exploration_df[\"capture_lon\"][2:],\n",
    "                )\n",
    "            ]\n",
    "            return pdist(run_exploration_coord_vector, metric=haversine)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Second-level RSA using exploration models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A helper function to compute neural RDM vector\n",
    "def compute_neural_rdm_vector(neural_beta_array: np.ndarray):\n",
    "    neural_rdm_vector = pdist(neural_beta_array, \"correlation\")\n",
    "    neural_rdm_vector = np.nan_to_num(\n",
    "        neural_rdm_vector, nan=1.0\n",
    "    )  # correlation distances can be nan!\n",
    "\n",
    "    return neural_rdm_vector\n",
    "\n",
    "\n",
    "# Main function to create second-level RSA dataframe from a cluster mask and subgroup name\n",
    "def create_exploration_second_level_rsa_df(\n",
    "    cluster_mask: tuple[str, np.ndarray], subgroup: Literal[\"discovery\", \"validation\"]\n",
    "):\n",
    "    cluster_name, cluster_mask = cluster_mask\n",
    "\n",
    "    print(f\"Processing {cluster_name} (Subgroup: {subgroup})...\")\n",
    "\n",
    "    exploration_model_list = all_model_rdm_list[6:]\n",
    "    run_id_list = [\"run-01\", \"run-02\", \"run-03\", \"run-04\", \"run-05\"]\n",
    "    trial_id_list = [\n",
    "        \"trial3\",\n",
    "        \"trial4\",\n",
    "        \"trial5\",\n",
    "        \"trial6\",\n",
    "        \"trial7\",\n",
    "        \"trial8\",\n",
    "    ]  # Using Trial 3 to Trial 8\n",
    "\n",
    "    subject_list = subgroup_subject_list_dict[subgroup]\n",
    "    subgroup_first_level_output_dir = subgroup_first_level_output_path_dict[subgroup]\n",
    "\n",
    "    rsa_result_dict_list: list[dict[str, list]] = []\n",
    "\n",
    "    for subject_id in subject_list:\n",
    "        subject_feedback_beta_dir = (\n",
    "            subgroup_first_level_output_dir\n",
    "            / subject_id\n",
    "            / \"rsa_neural_data\"\n",
    "            / \"feedback_beta\"\n",
    "        )\n",
    "\n",
    "        subject_rsa_dict: dict[str, list] = {}\n",
    "\n",
    "        for model_rdm in exploration_model_list:\n",
    "            subject_rsa_dict[model_rdm] = []\n",
    "\n",
    "        for run_id in run_id_list:\n",
    "            # 1. Compute neural rdm vector\n",
    "            subject_run_feedback_beta_list = []\n",
    "\n",
    "            for trial_id in trial_id_list:\n",
    "                norm_beta_path = (\n",
    "                    subject_feedback_beta_dir\n",
    "                    / f\"{subject_id}_task-photographer_{run_id}_{trial_id}_feedback_beta.nii\"\n",
    "                )\n",
    "                norm_beta_array = load_img(norm_beta_path).get_fdata()\n",
    "                masked_beta_array = norm_beta_array[cluster_mask]\n",
    "                subject_run_feedback_beta_list.append(masked_beta_array)\n",
    "\n",
    "            subject_run_feedback_beta_array = np.vstack(subject_run_feedback_beta_list)\n",
    "            subject_run_feedback_neural_rdm_vector = compute_neural_rdm_vector(\n",
    "                subject_run_feedback_beta_array\n",
    "            )\n",
    "\n",
    "            # 2. Load model rdm vector and perform RSA\n",
    "            for model_rdm in exploration_model_list:\n",
    "                subject_run_model_rdm_vector = load_model_rdm_vector(\n",
    "                    model_rdm, subgroup, subject_id, run_id\n",
    "                )\n",
    "\n",
    "                subject_run_rsa_corr = compare_rho_a(\n",
    "                    subject_run_feedback_neural_rdm_vector, subject_run_model_rdm_vector\n",
    "                )[0][0]  # Spearman's rank correlation\n",
    "                subject_run_rsa_corr = np.arctanh(\n",
    "                    subject_run_rsa_corr\n",
    "                )  # Fisher's r-to-z transform\n",
    "                subject_rsa_dict[model_rdm].append(subject_run_rsa_corr)\n",
    "\n",
    "        subject_rsa_result_dict = {\n",
    "            \"subgroup\": subgroup,\n",
    "            \"subject_id\": subject_id,\n",
    "            \"cluster\": cluster_name,\n",
    "        }\n",
    "\n",
    "        # Average individual RSA correlations across runs\n",
    "        for model_rdm in exploration_model_list:\n",
    "            subject_run_rsa_corr_list = subject_rsa_dict[model_rdm]\n",
    "            assert len(subject_run_rsa_corr_list) == len(run_id_list)\n",
    "            subject_mean_rsa_corr = np.mean(subject_run_rsa_corr_list)\n",
    "\n",
    "            subject_rsa_result_dict[model_rdm] = subject_mean_rsa_corr\n",
    "\n",
    "        rsa_result_dict_list.append(subject_rsa_result_dict)\n",
    "\n",
    "    second_level_RSA_df = pd.DataFrame.from_records(rsa_result_dict_list)\n",
    "\n",
    "    return second_level_RSA_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper functions to perform permutation-based one-sample t-tests on a numpy array\n",
    "def t_statistic(x: np.ndarray, axis=1):\n",
    "    return ttest_1samp(x, popmean=0, axis=axis).statistic\n",
    "\n",
    "\n",
    "def permute_1samp(x: np.ndarray, alternative: Literal[\"two-sided\", \"less\", \"greater\"]):\n",
    "    rng = np.random.RandomState(605)\n",
    "    return permutation_test(\n",
    "        (x,),\n",
    "        t_statistic,\n",
    "        permutation_type=\"samples\",\n",
    "        n_resamples=10000,\n",
    "        random_state=rng,\n",
    "        alternative=alternative,\n",
    "    )\n",
    "\n",
    "\n",
    "# Performing permuation tests (one-sample t-tests, one-sided) for a given second-level RSA results df and model RDM candidates\n",
    "def perm_test_on_second_level_rsa(rsa_df: pd.DataFrame, model_rdm_list: list[str]):\n",
    "    for rdm_model in model_rdm_list:\n",
    "        print(f\"{rdm_model} RDM\")\n",
    "\n",
    "        print(f\"\\tmean (z-scored) r = {round(np.mean(rsa_df[rdm_model]), 3)}\")\n",
    "\n",
    "        res = permute_1samp(rsa_df[rdm_model], alternative=\"greater\")  # one-sided test\n",
    "        print(\n",
    "            f\"\\tt = {round(res.statistic, 3)} (P = {round(res.pvalue, 4)}) {'*' if res.pvalue < 0.05 else ''}\"\n",
    "        )\n",
    "\n",
    "        d = np.mean(rsa_df[rdm_model]) / np.std(\n",
    "            rsa_df[rdm_model]\n",
    "        )  # standardized effect size\n",
    "        print(f\"\\td = {round(d, 3)}\")\n",
    "\n",
    "        print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 1_MiOG (Subgroup: discovery)...\n",
      "\n",
      "exploration_time RDM\n",
      "\tmean (z-scored) r = -0.043\n",
      "\tt = -1.364 (P = 0.9044) \n",
      "\td = -0.352\n",
      "\n",
      "exploration_distance RDM\n",
      "\tmean (z-scored) r = -0.049\n",
      "\tt = -1.557 (P = 0.9286) \n",
      "\td = -0.402\n",
      "\n",
      "capture_distance RDM\n",
      "\tmean (z-scored) r = 0.282\n",
      "\tt = 9.45 (P = 0.0001) *\n",
      "\td = 2.44\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (1) MiOG cluster, Discovery group\n",
    "exploration_second_level_MiOG_disovery_df = create_exploration_second_level_rsa_df(\n",
    "    feedback_history_cluster_mask_list[0], \"discovery\"\n",
    ")\n",
    "\n",
    "exploration_second_level_rsa_output_dir = (\n",
    "    rsa_result_output_path / \"exploration_second_level\"\n",
    ")\n",
    "os.makedirs(exploration_second_level_rsa_output_dir, exist_ok=True)\n",
    "\n",
    "exploration_second_level_MiOG_disovery_df.to_csv(\n",
    "    exploration_second_level_rsa_output_dir / \"discovery_MiOG_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_second_level_rsa(\n",
    "    exploration_second_level_MiOG_disovery_df, all_model_rdm_list[6:]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 1_MiOG (Subgroup: validation)...\n",
      "\n",
      "exploration_time RDM\n",
      "\tmean (z-scored) r = 0.002\n",
      "\tt = 0.058 (P = 0.4806) \n",
      "\td = 0.015\n",
      "\n",
      "exploration_distance RDM\n",
      "\tmean (z-scored) r = 0.02\n",
      "\tt = 0.649 (P = 0.2602) \n",
      "\td = 0.168\n",
      "\n",
      "capture_distance RDM\n",
      "\tmean (z-scored) r = 0.309\n",
      "\tt = 6.917 (P = 0.0001) *\n",
      "\td = 1.786\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (2) MiOG cluster, Validation group\n",
    "exploration_second_level_MiOG_validation_df = create_exploration_second_level_rsa_df(\n",
    "    feedback_history_cluster_mask_list[0], \"validation\"\n",
    ")\n",
    "\n",
    "exploration_second_level_MiOG_validation_df.to_csv(\n",
    "    exploration_second_level_rsa_output_dir / \"validation_MiOG_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_second_level_rsa(\n",
    "    exploration_second_level_MiOG_validation_df, all_model_rdm_list[6:]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 2_IFG (Subgroup: discovery)...\n",
      "\n",
      "exploration_time RDM\n",
      "\tmean (z-scored) r = -0.068\n",
      "\tt = -2.295 (P = 0.98) \n",
      "\td = -0.593\n",
      "\n",
      "exploration_distance RDM\n",
      "\tmean (z-scored) r = -0.044\n",
      "\tt = -1.378 (P = 0.9049) \n",
      "\td = -0.356\n",
      "\n",
      "capture_distance RDM\n",
      "\tmean (z-scored) r = 0.209\n",
      "\tt = 5.981 (P = 0.0003) *\n",
      "\td = 1.544\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (3) IFG cluster, Discovery group\n",
    "exploration_second_level_IFG_disovery_df = create_exploration_second_level_rsa_df(\n",
    "    feedback_history_cluster_mask_list[1], \"discovery\"\n",
    ")\n",
    "\n",
    "exploration_second_level_IFG_disovery_df.to_csv(\n",
    "    exploration_second_level_rsa_output_dir / \"discovery_IFG_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_second_level_rsa(\n",
    "    exploration_second_level_IFG_disovery_df, all_model_rdm_list[6:]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 2_IFG (Subgroup: validation)...\n",
      "\n",
      "exploration_time RDM\n",
      "\tmean (z-scored) r = 0.009\n",
      "\tt = 0.324 (P = 0.3819) \n",
      "\td = 0.084\n",
      "\n",
      "exploration_distance RDM\n",
      "\tmean (z-scored) r = 0.037\n",
      "\tt = 1.481 (P = 0.0809) \n",
      "\td = 0.382\n",
      "\n",
      "capture_distance RDM\n",
      "\tmean (z-scored) r = 0.255\n",
      "\tt = 6.806 (P = 0.0001) *\n",
      "\td = 1.757\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (4) IFG cluster, Validation group\n",
    "exploration_second_level_IFG_validation_df = create_exploration_second_level_rsa_df(\n",
    "    feedback_history_cluster_mask_list[1], \"validation\"\n",
    ")\n",
    "\n",
    "exploration_second_level_IFG_validation_df.to_csv(\n",
    "    exploration_second_level_rsa_output_dir / \"validation_IFG_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_second_level_rsa(\n",
    "    exploration_second_level_IFG_validation_df, all_model_rdm_list[6:]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Partial correlation-based second-level RSA on feedback history models (controlling for the Capture Distance model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main function to create partial correlation-based second-level RSA dataframe\n",
    "# We regress out the Capture Distance model from both neural and feedback history model RDMs\n",
    "def create_partial_correlation_feedback_history_second_level_rsa_df(\n",
    "    cluster_mask: tuple[str, np.ndarray], subgroup: Literal[\"discovery\", \"validation\"]\n",
    "):\n",
    "    cluster_name, cluster_mask = cluster_mask\n",
    "\n",
    "    print(f\"Processing {cluster_name} (Subgroup: {subgroup})...\")\n",
    "\n",
    "    feedback_history_model_rdm_list = all_model_rdm_list[:6]\n",
    "    run_id_list = [\"run-01\", \"run-02\", \"run-03\", \"run-04\", \"run-05\"]\n",
    "    trial_id_list = [\"trial3\", \"trial4\", \"trial5\", \"trial6\", \"trial7\", \"trial8\"]\n",
    "\n",
    "    subject_list = subgroup_subject_list_dict[subgroup]\n",
    "    subgroup_first_level_output_dir = subgroup_first_level_output_path_dict[subgroup]\n",
    "\n",
    "    rsa_result_dict_list: list[dict[str, list]] = []\n",
    "\n",
    "    for subject_id in subject_list:\n",
    "        subject_feedback_beta_dir = (\n",
    "            subgroup_first_level_output_dir\n",
    "            / subject_id\n",
    "            / \"rsa_neural_data\"\n",
    "            / \"feedback_beta\"\n",
    "        )\n",
    "\n",
    "        subject_rsa_dict: dict[str, list] = {}\n",
    "\n",
    "        for rdm_model in feedback_history_model_rdm_list:\n",
    "            subject_rsa_dict[rdm_model] = []\n",
    "\n",
    "        for run_id in run_id_list:\n",
    "            # 1. Compute neural rdm vector\n",
    "            subject_run_feedback_beta_list = []\n",
    "\n",
    "            for trial_id in trial_id_list:\n",
    "                norm_beta_path = (\n",
    "                    subject_feedback_beta_dir\n",
    "                    / f\"{subject_id}_task-photographer_{run_id}_{trial_id}_feedback_beta.nii\"\n",
    "                )\n",
    "                norm_beta_array = load_img(norm_beta_path).get_fdata()\n",
    "                masked_beta_array = norm_beta_array[cluster_mask]\n",
    "                subject_run_feedback_beta_list.append(masked_beta_array)\n",
    "\n",
    "            subject_run_feedback_beta_array = np.vstack(subject_run_feedback_beta_list)\n",
    "            subject_run_feedback_neural_rdm_vector = compute_neural_rdm_vector(\n",
    "                subject_run_feedback_beta_array\n",
    "            )\n",
    "\n",
    "            # 2. Load the Capture Distance rdm vector\n",
    "            subject_run_capture_distance_vector = load_model_rdm_vector(\n",
    "                \"capture_distance\", subgroup, subject_id, run_id\n",
    "            )\n",
    "\n",
    "            # 3. load model rdm vector and perform partial correlation RSA\n",
    "            for rdm_model in feedback_history_model_rdm_list:\n",
    "                subject_run_model_rdm_vector = load_model_rdm_vector(\n",
    "                    rdm_model, subgroup, subject_id, run_id\n",
    "                )\n",
    "\n",
    "                fit_neural = LinearRegression(fit_intercept=True)\n",
    "                fit_model = LinearRegression(fit_intercept=True)\n",
    "\n",
    "                fit_neural.fit(\n",
    "                    subject_run_capture_distance_vector.reshape(-1, 1),\n",
    "                    subject_run_feedback_neural_rdm_vector,\n",
    "                )\n",
    "                fit_model.fit(\n",
    "                    subject_run_capture_distance_vector.reshape(-1, 1),\n",
    "                    subject_run_model_rdm_vector,\n",
    "                )\n",
    "\n",
    "                # Regress out information associated with the Capture Distance model\n",
    "                # from both model and neural RDM vectors\n",
    "                residual_neural_vector = (\n",
    "                    subject_run_feedback_neural_rdm_vector\n",
    "                    - fit_neural.predict(\n",
    "                        subject_run_capture_distance_vector.reshape(-1, 1)\n",
    "                    )\n",
    "                )\n",
    "                residual_model_vector = (\n",
    "                    subject_run_model_rdm_vector\n",
    "                    - fit_model.predict(\n",
    "                        subject_run_capture_distance_vector.reshape(-1, 1)\n",
    "                    )\n",
    "                )\n",
    "\n",
    "                # Perform RSA between residual vectors\n",
    "                subject_run_rsa_corr = compare_rho_a(\n",
    "                    residual_neural_vector, residual_model_vector\n",
    "                )[0][0]\n",
    "                subject_run_rsa_corr = np.arctanh(subject_run_rsa_corr)\n",
    "                subject_rsa_dict[rdm_model].append(subject_run_rsa_corr)\n",
    "\n",
    "        subject_rsa_result_dict = {\n",
    "            \"subgroup\": subgroup,\n",
    "            \"subject_id\": subject_id,\n",
    "            \"cluster\": cluster_name,\n",
    "        }\n",
    "\n",
    "        for rdm_model in feedback_history_model_rdm_list:\n",
    "            subject_run_rsa_corr_list = subject_rsa_dict[rdm_model]\n",
    "            assert len(subject_run_rsa_corr_list) == len(run_id_list)\n",
    "            subject_mean_rsa_corr = np.mean(subject_run_rsa_corr_list)\n",
    "\n",
    "            subject_rsa_result_dict[rdm_model] = subject_mean_rsa_corr\n",
    "\n",
    "        rsa_result_dict_list.append(subject_rsa_result_dict)\n",
    "\n",
    "    second_level_RSA_df = pd.DataFrame.from_records(rsa_result_dict_list)\n",
    "\n",
    "    return second_level_RSA_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 1_MiOG (Subgroup: discovery)...\n",
      "\n",
      "current_trial RDM\n",
      "\tmean (z-scored) r = 0.027\n",
      "\tt = 0.751 (P = 0.2393) \n",
      "\td = 0.194\n",
      "\n",
      "one_back_trial RDM\n",
      "\tmean (z-scored) r = 0.032\n",
      "\tt = 1.468 (P = 0.0779) \n",
      "\td = 0.379\n",
      "\n",
      "two_back_trial RDM\n",
      "\tmean (z-scored) r = -0.003\n",
      "\tt = -0.073 (P = 0.5324) \n",
      "\td = -0.019\n",
      "\n",
      "recent_2_trial RDM\n",
      "\tmean (z-scored) r = 0.067\n",
      "\tt = 2.407 (P = 0.0171) *\n",
      "\td = 0.621\n",
      "\n",
      "recent_3_trial RDM\n",
      "\tmean (z-scored) r = 0.048\n",
      "\tt = 1.263 (P = 0.1177) \n",
      "\td = 0.326\n",
      "\n",
      "previous_2_trial RDM\n",
      "\tmean (z-scored) r = 0.042\n",
      "\tt = 1.157 (P = 0.1291) \n",
      "\td = 0.299\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (1) MiOG cluster, Discovery group\n",
    "partial_feedback_history_second_level_MiOG_disovery_df = (\n",
    "    create_partial_correlation_feedback_history_second_level_rsa_df(\n",
    "        feedback_history_cluster_mask_list[0], \"discovery\"\n",
    "    )\n",
    ")\n",
    "\n",
    "partial_feedback_history_second_level_rsa_output_dir = (\n",
    "    rsa_result_output_path / \"partial_feedback_history_second_level\"\n",
    ")\n",
    "os.makedirs(partial_feedback_history_second_level_rsa_output_dir, exist_ok=True)\n",
    "\n",
    "partial_feedback_history_second_level_MiOG_disovery_df.to_csv(\n",
    "    partial_feedback_history_second_level_rsa_output_dir\n",
    "    / \"discovery_MiOG_partial_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_second_level_rsa(\n",
    "    partial_feedback_history_second_level_MiOG_disovery_df, all_model_rdm_list[:6]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 1_MiOG (Subgroup: validation)...\n",
      "\n",
      "current_trial RDM\n",
      "\tmean (z-scored) r = 0.038\n",
      "\tt = 1.225 (P = 0.1151) \n",
      "\td = 0.316\n",
      "\n",
      "one_back_trial RDM\n",
      "\tmean (z-scored) r = 0.082\n",
      "\tt = 2.008 (P = 0.0339) *\n",
      "\td = 0.518\n",
      "\n",
      "two_back_trial RDM\n",
      "\tmean (z-scored) r = -0.014\n",
      "\tt = -0.452 (P = 0.6667) \n",
      "\td = -0.117\n",
      "\n",
      "recent_2_trial RDM\n",
      "\tmean (z-scored) r = 0.115\n",
      "\tt = 3.354 (P = 0.0023) *\n",
      "\td = 0.866\n",
      "\n",
      "recent_3_trial RDM\n",
      "\tmean (z-scored) r = 0.107\n",
      "\tt = 3.111 (P = 0.0044) *\n",
      "\td = 0.803\n",
      "\n",
      "previous_2_trial RDM\n",
      "\tmean (z-scored) r = 0.078\n",
      "\tt = 2.142 (P = 0.0219) *\n",
      "\td = 0.553\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (2) MiOG cluster, Validation group\n",
    "partial_feedback_history_second_level_MiOG_validation_df = (\n",
    "    create_partial_correlation_feedback_history_second_level_rsa_df(\n",
    "        feedback_history_cluster_mask_list[0], \"validation\"\n",
    "    )\n",
    ")\n",
    "\n",
    "partial_feedback_history_second_level_MiOG_validation_df.to_csv(\n",
    "    partial_feedback_history_second_level_rsa_output_dir\n",
    "    / \"validation_MiOG_partial_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_second_level_rsa(\n",
    "    partial_feedback_history_second_level_MiOG_validation_df, all_model_rdm_list[:6]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 2_IFG (Subgroup: discovery)...\n",
      "\n",
      "current_trial RDM\n",
      "\tmean (z-scored) r = 0.042\n",
      "\tt = 1.431 (P = 0.0881) \n",
      "\td = 0.369\n",
      "\n",
      "one_back_trial RDM\n",
      "\tmean (z-scored) r = 0.021\n",
      "\tt = 1.008 (P = 0.1849) \n",
      "\td = 0.26\n",
      "\n",
      "two_back_trial RDM\n",
      "\tmean (z-scored) r = 0.019\n",
      "\tt = 0.571 (P = 0.2941) \n",
      "\td = 0.147\n",
      "\n",
      "recent_2_trial RDM\n",
      "\tmean (z-scored) r = 0.064\n",
      "\tt = 2.702 (P = 0.0095) *\n",
      "\td = 0.698\n",
      "\n",
      "recent_3_trial RDM\n",
      "\tmean (z-scored) r = 0.073\n",
      "\tt = 3.211 (P = 0.0031) *\n",
      "\td = 0.829\n",
      "\n",
      "previous_2_trial RDM\n",
      "\tmean (z-scored) r = 0.041\n",
      "\tt = 1.252 (P = 0.1147) \n",
      "\td = 0.323\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (3) IFG cluster, Discovery group\n",
    "partial_feedback_history_second_level_IFG_disovery_df = (\n",
    "    create_partial_correlation_feedback_history_second_level_rsa_df(\n",
    "        feedback_history_cluster_mask_list[1], \"discovery\"\n",
    "    )\n",
    ")\n",
    "\n",
    "partial_feedback_history_second_level_IFG_disovery_df.to_csv(\n",
    "    partial_feedback_history_second_level_rsa_output_dir\n",
    "    / \"discovery_IFG_partial_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_second_level_rsa(\n",
    "    partial_feedback_history_second_level_IFG_disovery_df, all_model_rdm_list[:6]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 2_IFG (Subgroup: validation)...\n",
      "\n",
      "current_trial RDM\n",
      "\tmean (z-scored) r = 0.026\n",
      "\tt = 1.256 (P = 0.1118) \n",
      "\td = 0.324\n",
      "\n",
      "one_back_trial RDM\n",
      "\tmean (z-scored) r = 0.102\n",
      "\tt = 2.388 (P = 0.018) *\n",
      "\td = 0.617\n",
      "\n",
      "two_back_trial RDM\n",
      "\tmean (z-scored) r = 0.018\n",
      "\tt = 0.439 (P = 0.33) \n",
      "\td = 0.113\n",
      "\n",
      "recent_2_trial RDM\n",
      "\tmean (z-scored) r = 0.115\n",
      "\tt = 4.315 (P = 0.0004) *\n",
      "\td = 1.114\n",
      "\n",
      "recent_3_trial RDM\n",
      "\tmean (z-scored) r = 0.131\n",
      "\tt = 3.798 (P = 0.0014) *\n",
      "\td = 0.981\n",
      "\n",
      "previous_2_trial RDM\n",
      "\tmean (z-scored) r = 0.128\n",
      "\tt = 3.672 (P = 0.0003) *\n",
      "\td = 0.948\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (4) IFG cluster, Validation group\n",
    "partial_feedback_history_second_level_IFG_validation_df = (\n",
    "    create_partial_correlation_feedback_history_second_level_rsa_df(\n",
    "        feedback_history_cluster_mask_list[1], \"validation\"\n",
    "    )\n",
    ")\n",
    "\n",
    "partial_feedback_history_second_level_IFG_validation_df.to_csv(\n",
    "    partial_feedback_history_second_level_rsa_output_dir\n",
    "    / \"validation_IFG_partial_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_second_level_rsa(\n",
    "    partial_feedback_history_second_level_IFG_validation_df, all_model_rdm_list[:6]\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Partial correlation-based second-level RSA on Capture Distance model (controlling for either Recent-2 or Recent-3 Trial model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main function to create partial correlation-based second-level RSA dataframe\n",
    "# We regress out either Recent-2 or Recent-3 Trial model (feedback history models) from both neural and Capture Distance model RDMs\n",
    "def create_partial_correlation_capture_distance_second_level_rsa_df(\n",
    "    cluster_mask: tuple[str, np.ndarray], subgroup: Literal[\"discovery\", \"validation\"]\n",
    "):\n",
    "    cluster_name, cluster_mask = cluster_mask\n",
    "\n",
    "    print(f\"Processing {cluster_name} (Subgroup {subgroup})...\")\n",
    "\n",
    "    feedback_history_model_rdm_list = [\n",
    "        \"recent_2_trial\",\n",
    "        \"recent_3_trial\",\n",
    "    ]\n",
    "\n",
    "    run_id_list = [\"run-01\", \"run-02\", \"run-03\", \"run-04\", \"run-05\"]\n",
    "    trial_id_list = [\"trial3\", \"trial4\", \"trial5\", \"trial6\", \"trial7\", \"trial8\"]\n",
    "\n",
    "    subject_list = subgroup_subject_list_dict[subgroup]\n",
    "    subgroup_first_level_output_dir = subgroup_first_level_output_path_dict[subgroup]\n",
    "\n",
    "    rsa_result_dict_list: list[dict[str, list]] = []\n",
    "\n",
    "    for subject_id in subject_list:\n",
    "        subject_feedback_beta_dir = (\n",
    "            subgroup_first_level_output_dir\n",
    "            / subject_id\n",
    "            / \"rsa_neural_data\"\n",
    "            / \"feedback_beta\"\n",
    "        )\n",
    "\n",
    "        subject_rsa_dict: dict[str, list] = {}\n",
    "\n",
    "        for rdm_model in feedback_history_model_rdm_list:\n",
    "            subject_rsa_dict[rdm_model] = []\n",
    "\n",
    "        for run_id in run_id_list:\n",
    "            # 1. Compute neural rdm vector\n",
    "            subject_run_feedback_beta_list = []\n",
    "\n",
    "            for trial_id in trial_id_list:\n",
    "                norm_beta_path = (\n",
    "                    subject_feedback_beta_dir\n",
    "                    / f\"{subject_id}_task-photographer_{run_id}_{trial_id}_feedback_beta.nii\"\n",
    "                )\n",
    "                norm_beta_array = load_img(norm_beta_path).get_fdata()\n",
    "                masked_beta_array = norm_beta_array[cluster_mask]\n",
    "                subject_run_feedback_beta_list.append(masked_beta_array)\n",
    "\n",
    "            subject_run_feedback_beta_array = np.vstack(subject_run_feedback_beta_list)\n",
    "            subject_run_feedback_neural_rdm_vector = compute_neural_rdm_vector(\n",
    "                subject_run_feedback_beta_array\n",
    "            )\n",
    "\n",
    "            # 2. Load the Capture Distance rdm vector\n",
    "            subject_run_capture_distance_vector = load_model_rdm_vector(\n",
    "                \"capture_distance\", subgroup, subject_id, run_id\n",
    "            )\n",
    "\n",
    "            # 3. Load model feedback history rdm vector and perform RSA\n",
    "            for rdm_model in feedback_history_model_rdm_list:\n",
    "                subject_run_model_rdm_vector = load_model_rdm_vector(\n",
    "                    rdm_model, subgroup, subject_id, run_id\n",
    "                )\n",
    "\n",
    "                fit_neural = LinearRegression(fit_intercept=True)\n",
    "                fit_capture = LinearRegression(fit_intercept=True)\n",
    "\n",
    "                fit_neural.fit(\n",
    "                    subject_run_model_rdm_vector.reshape(-1, 1),\n",
    "                    subject_run_feedback_neural_rdm_vector,\n",
    "                )\n",
    "                fit_capture.fit(\n",
    "                    subject_run_model_rdm_vector.reshape(-1, 1),\n",
    "                    subject_run_capture_distance_vector,\n",
    "                )\n",
    "\n",
    "                # In this case, regress out the feedback history model from\n",
    "                # neural and \"Capture Distance\" RDMs\n",
    "                residual_neural_vector = (\n",
    "                    subject_run_feedback_neural_rdm_vector\n",
    "                    - fit_neural.predict(subject_run_model_rdm_vector.reshape(-1, 1))\n",
    "                )\n",
    "                residual_capture_vector = (\n",
    "                    subject_run_capture_distance_vector\n",
    "                    - fit_capture.predict(subject_run_model_rdm_vector.reshape(-1, 1))\n",
    "                )\n",
    "\n",
    "                subject_run_rsa_corr = compare_rho_a(\n",
    "                    residual_neural_vector, residual_capture_vector\n",
    "                )[0][0]\n",
    "                subject_run_rsa_corr = np.arctanh(subject_run_rsa_corr)\n",
    "                subject_rsa_dict[rdm_model].append(subject_run_rsa_corr)\n",
    "\n",
    "        subject_rsa_result_dict = {\n",
    "            \"subgroup\": subgroup,\n",
    "            \"subject_id\": subject_id,\n",
    "            \"cluster\": cluster_name,\n",
    "        }\n",
    "\n",
    "        for rdm_model in feedback_history_model_rdm_list:\n",
    "            rsa_name = f\"capture_distance__{rdm_model}\"\n",
    "            subject_run_rsa_corr_list = subject_rsa_dict[rdm_model]\n",
    "            assert len(subject_run_rsa_corr_list) == len(run_id_list)\n",
    "            subject_mean_rsa_corr = np.mean(subject_run_rsa_corr_list)\n",
    "\n",
    "            subject_rsa_result_dict[rsa_name] = subject_mean_rsa_corr\n",
    "\n",
    "        rsa_result_dict_list.append(subject_rsa_result_dict)\n",
    "\n",
    "    second_level_RSA_df = pd.DataFrame.from_records(rsa_result_dict_list)\n",
    "\n",
    "    return second_level_RSA_df\n",
    "\n",
    "\n",
    "# A helper function to perform permutation tests on dataframes from above\n",
    "# since the model RDM column names are distinct from previous second-level RSA results\n",
    "def perm_test_on_capture_distance_partial_correlation(rsa_df: pd.DataFrame):\n",
    "    feedback_history_model_list = [\n",
    "        \"recent_2_trial\",\n",
    "        \"recent_3_trial\",\n",
    "    ]\n",
    "\n",
    "    for rdm_model in feedback_history_model_list:\n",
    "        rsa_name = f\"capture_distance__{rdm_model}\"\n",
    "        print(f\"Capture Distance <- {rdm_model} RDM (regressed out)\")\n",
    "\n",
    "        print(f\"\\tmean (z-scored) r = {round(np.mean(rsa_df[rsa_name]), 3)}\")\n",
    "\n",
    "        res = permute_1samp(rsa_df[rsa_name], alternative=\"greater\")\n",
    "        print(\n",
    "            f\"\\tt = {round(res.statistic, 3)} (P = {round(res.pvalue, 4)}) {'*' if res.pvalue < 0.05 else ''}\"\n",
    "        )\n",
    "\n",
    "        d = np.mean(rsa_df[rsa_name]) / np.std(rsa_df[rsa_name])\n",
    "        print(f\"\\td = {round(d, 3)}\")\n",
    "\n",
    "        print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 1_MiOG (Subgroup discovery)...\n",
      "\n",
      "Capture Distance <- recent_2_trial RDM (regressed out)\n",
      "\tmean (z-scored) r = 0.267\n",
      "\tt = 7.67 (P = 0.0002) *\n",
      "\td = 1.98\n",
      "\n",
      "Capture Distance <- recent_3_trial RDM (regressed out)\n",
      "\tmean (z-scored) r = 0.257\n",
      "\tt = 7.886 (P = 0.0001) *\n",
      "\td = 2.036\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (1) MiOG cluster, Discovery group\n",
    "partial_capture_distance_second_level_MiOG_disovery_df = (\n",
    "    create_partial_correlation_capture_distance_second_level_rsa_df(\n",
    "        feedback_history_cluster_mask_list[0], \"discovery\"\n",
    "    )\n",
    ")\n",
    "\n",
    "partial_capture_distance_second_level_rsa_output_dir = (\n",
    "    rsa_result_output_path / \"partial_capture_distance_second_level\"\n",
    ")\n",
    "os.makedirs(partial_capture_distance_second_level_rsa_output_dir, exist_ok=True)\n",
    "\n",
    "partial_capture_distance_second_level_MiOG_disovery_df.to_csv(\n",
    "    partial_capture_distance_second_level_rsa_output_dir\n",
    "    / \"discovery_MiOG_partial_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_capture_distance_partial_correlation(\n",
    "    partial_capture_distance_second_level_MiOG_disovery_df\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 1_MiOG (Subgroup validation)...\n",
      "\n",
      "Capture Distance <- recent_2_trial RDM (regressed out)\n",
      "\tmean (z-scored) r = 0.285\n",
      "\tt = 6.502 (P = 0.0001) *\n",
      "\td = 1.679\n",
      "\n",
      "Capture Distance <- recent_3_trial RDM (regressed out)\n",
      "\tmean (z-scored) r = 0.269\n",
      "\tt = 5.941 (P = 0.0001) *\n",
      "\td = 1.534\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (2) MiOG cluster, Validation group\n",
    "partial_capture_distance_second_level_MiOG_validation_df = (\n",
    "    create_partial_correlation_capture_distance_second_level_rsa_df(\n",
    "        feedback_history_cluster_mask_list[0], \"validation\"\n",
    "    )\n",
    ")\n",
    "\n",
    "partial_capture_distance_second_level_MiOG_validation_df.to_csv(\n",
    "    partial_capture_distance_second_level_rsa_output_dir\n",
    "    / \"validation_MiOG_partial_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_capture_distance_partial_correlation(\n",
    "    partial_capture_distance_second_level_MiOG_validation_df\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 2_IFG (Subgroup discovery)...\n",
      "\n",
      "Capture Distance <- recent_2_trial RDM (regressed out)\n",
      "\tmean (z-scored) r = 0.185\n",
      "\tt = 4.85 (P = 0.0004) *\n",
      "\td = 1.252\n",
      "\n",
      "Capture Distance <- recent_3_trial RDM (regressed out)\n",
      "\tmean (z-scored) r = 0.172\n",
      "\tt = 4.87 (P = 0.0003) *\n",
      "\td = 1.258\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (3) IFG cluster, Discovery group\n",
    "partial_capture_distance_second_level_IFG_disovery_df = (\n",
    "    create_partial_correlation_capture_distance_second_level_rsa_df(\n",
    "        feedback_history_cluster_mask_list[1], \"discovery\"\n",
    "    )\n",
    ")\n",
    "\n",
    "partial_capture_distance_second_level_IFG_disovery_df.to_csv(\n",
    "    partial_capture_distance_second_level_rsa_output_dir\n",
    "    / \"discovery_IFG_partial_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_capture_distance_partial_correlation(\n",
    "    partial_capture_distance_second_level_IFG_disovery_df\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing 2_IFG (Subgroup validation)...\n",
      "\n",
      "Capture Distance <- recent_2_trial RDM (regressed out)\n",
      "\tmean (z-scored) r = 0.227\n",
      "\tt = 5.689 (P = 0.0002) *\n",
      "\td = 1.469\n",
      "\n",
      "Capture Distance <- recent_3_trial RDM (regressed out)\n",
      "\tmean (z-scored) r = 0.229\n",
      "\tt = 5.456 (P = 0.0002) *\n",
      "\td = 1.409\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# (4) IFG cluster, Validation group\n",
    "partial_capture_distance_second_level_IFG_disovery_df = (\n",
    "    create_partial_correlation_capture_distance_second_level_rsa_df(\n",
    "        feedback_history_cluster_mask_list[1], \"validation\"\n",
    "    )\n",
    ")\n",
    "\n",
    "partial_capture_distance_second_level_IFG_disovery_df.to_csv(\n",
    "    partial_capture_distance_second_level_rsa_output_dir\n",
    "    / \"validation_IFG_partial_second_level.csv\",\n",
    "    index=False,\n",
    ")\n",
    "\n",
    "print()\n",
    "perm_test_on_capture_distance_partial_correlation(\n",
    "    partial_capture_distance_second_level_IFG_disovery_df\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.10.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}