brain-bridge2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BRAIN-BRIDGE\n",
    "## Let's explore individual structural differences among the Human Connectome Project dataset! How can we ~bridge~ one brain with the rest?\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "## We extracted all FreeSurfer parcellation metrics for the 74 regions of the Destrieux Atlas for 1042 subjects."
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "### For each region, we will examine:"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  1. Gray Matter Volume"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  2. Average Cortical Thickness"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  3. Surface Area"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  4. Number of Vertices"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  5. Thickness Std Dev"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  6. Mean Curvature"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  7. Gaussian Curvature"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  8. Folding Index"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/markdown": [
       "####  9. Intrinsic Curvature Index"
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "from nilearn import datasets, plotting \n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display, Markdown\n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from sklearn.pipeline import Pipeline\n",
    "from sklearn.feature_selection import SelectKBest, f_regression\n",
    "from sklearn.linear_model import LinearRegression, Ridge\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "struct_df = pd.read_csv('data/all_structural_metrics_HCP.csv')\n",
    "#display(Markdown(f'## There are {struct_df.shape[0]} subjects in the HCP dataset.'))\n",
    "display(Markdown(f'## We extracted all FreeSurfer parcellation metrics for the 74 regions of the Destrieux Atlas for {struct_df.shape[0]} subjects.'))\n",
    "display(Markdown('### For each region, we will examine:'))\n",
    "display(Markdown('####  1. Gray Matter Volume'))\n",
    "display(Markdown('####  2. Average Cortical Thickness'))\n",
    "display(Markdown('####  3. Surface Area'))\n",
    "display(Markdown('####  4. Number of Vertices'))\n",
    "display(Markdown('####  5. Thickness Std Dev'))\n",
    "display(Markdown('####  6. Mean Curvature'))\n",
    "display(Markdown('####  7. Gaussian Curvature'))\n",
    "display(Markdown('####  8. Folding Index'))\n",
    "display(Markdown('####  9. Intrinsic Curvature Index'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create individual dataframes for each measure\n",
    "thick_df = struct_df.filter(regex=('ThickAvg$'))\n",
    "thick_df.insert(0, 'id_number', struct_df['id_number'])\n",
    "vol_df = struct_df.filter(regex='GrayVol$')\n",
    "vol_df.insert(0, 'id_number', struct_df['id_number'])\n",
    "area_df = struct_df.filter(regex='SurfArea$')\n",
    "area_df.insert(0, 'id_number', struct_df['id_number'])\n",
    "numvert_df = struct_df.filter(regex='NumVert$')\n",
    "numvert_df.insert(0, 'id_number', struct_df['id_number'])\n",
    "meancurv_df = struct_df.filter(regex='MeanCurv$')\n",
    "meancurv_df.insert(0, 'id_number', struct_df['id_number']) \n",
    "curvind_df = struct_df.filter(regex='CurvInd$')\n",
    "curvind_df.insert(0, 'id_number', struct_df['id_number']) \n",
    "foldind_df = struct_df.filter(regex='FoldInd$')\n",
    "foldind_df.insert(0, 'id_number', struct_df['id_number']) \n",
    "gauscurv_df = struct_df.filter(regex='GausCurv$')\n",
    "gauscurv_df.insert(0, 'id_number', struct_df['id_number'])\n",
    "thickstd_df = struct_df.filter(regex='ThickStd$')\n",
    "thickstd_df.insert(0, 'id_number', struct_df['id_number'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_brain = widgets.Output()\n",
    "out_table = widgets.Output()\n",
    "out_bar = widgets.Output()\n",
    "out_plot = widgets.Output()\n",
    "\n",
    "\n",
    "def zscore(id_number, dtype):\n",
    "    if dtype == 'Gray Matter Volume':\n",
    "        df = vol_df\n",
    "    elif dtype == 'Surface Area':\n",
    "        df = area_df\n",
    "    elif dtype == 'Average Cortical Thickness':\n",
    "        df = thick_df\n",
    "    elif dtype == 'Mean Curvature':\n",
    "        df = meancurv_df\n",
    "    elif dtype == 'Number Of Vertices':\n",
    "        df = numvert_df\n",
    "    elif dtype == 'Curvature Index':\n",
    "        df = curvind_df\n",
    "    elif dtype == 'Folding Index':\n",
    "        df = foldind_df\n",
    "    elif dtype == 'Gaussian Curvature':\n",
    "        df = meancurv_df\n",
    "    elif dtype == 'Thickness Std Dev':\n",
    "        df = thickstd_df\n",
    "        \n",
    "    patient_df = df[df['id_number'] == id_number].drop(columns=['id_number'])\n",
    "    control_df = df[df['id_number'] != id_number].drop(columns=['id_number'])\n",
    "    control_means = np.mean(control_df, axis=0)\n",
    "    control_stds = np.std(control_df, axis=0)\n",
    "    return (patient_df - control_means) / control_stds\n",
    "\n",
    "def plot_brain(z_data, dtype):\n",
    "    # The left and right hemisphere views will be stored here\n",
    "    left_hemi_widget = None\n",
    "    right_hemi_widget = None\n",
    "\n",
    "    # 1. Determine the global min and max values across both hemispheres for consistent color scaling\n",
    "    global_vmin = min(z_data.min())\n",
    "    global_vmax = max(z_data.max())\n",
    "\n",
    "    for hemi in ['left', 'right']:\n",
    "        z_data_tmp = z_data.copy()\n",
    "\n",
    "        # Fetching Destrieux atlas and fsaverage\n",
    "        destrieux_atlas = datasets.fetch_atlas_surf_destrieux()\n",
    "        fsaverage = datasets.fetch_surf_fsaverage()\n",
    "\n",
    "        # Filter z_data to only include variables starting with lh_ or rh_\n",
    "        if hemi == 'left':\n",
    "            z_data_tmp = z_data_tmp.filter(regex='^lh_')\n",
    "            z_data_tmp.columns = z_data_tmp.columns.str.replace('lh_', '')\n",
    "        else:\n",
    "            z_data_tmp = z_data_tmp.filter(regex='^rh_')\n",
    "            z_data_tmp.columns = z_data_tmp.columns.str.replace('rh_', '')\n",
    "\n",
    "        # remove the ending of columns after last underscore\n",
    "        z_data_tmp.columns = z_data_tmp.columns.str.rsplit(pat='_', n=1).str[0]\n",
    "        \n",
    "        if 'Lat_Fis-post' not in z_data_tmp.columns:\n",
    "            # Handle the case when 'Lat_Fis-post' doesn't exist. Raise an error.\n",
    "            raise ValueError('Lat_Fis-post does not exist. I must place a Medial_wall column in the dataframe after this column to match up with the atlas labels.')\n",
    "\n",
    "        loc = z_data_tmp.columns.get_loc('Lat_Fis-post')\n",
    "        z_data_tmp.insert(loc+1, 'Medial_wall', 0)\n",
    "        z_data_tmp.insert(0, 'Unknown', 0)\n",
    "\n",
    "        # Zip the values of z_data to the labels of the Destrieux atlas\n",
    "        region_values = dict(zip(destrieux_atlas['labels'], z_data_tmp.values[0]))\n",
    "\n",
    "        # Initialize an array with zeros\n",
    "        mapped_values = np.zeros_like(destrieux_atlas[f'map_{hemi}'], dtype=float)\n",
    "\n",
    "        # Populate the mapped_values array using regi)on_values dictionary\n",
    "        for label, value in region_values.items():\n",
    "            region_idx = destrieux_atlas['labels'].index(label)\n",
    "            mapped_values[destrieux_atlas[f'map_{hemi}'] == region_idx] = value\n",
    "\n",
    "        fsavg = fsaverage.infl_left if hemi == 'left' else fsaverage.infl_right\n",
    "\n",
    "        # Plotting\n",
    "  \n",
    "        view = plotting.view_surf(fsavg, mapped_values,\n",
    "                                cmap='coolwarm', symmetric_cmap=True,\n",
    "                                vmax=np.max([np.abs(global_vmin), np.abs(global_vmax)]))\n",
    "                                #view=hemi_view)\n",
    "\n",
    "        widget = widgets.HTML(view.get_iframe())\n",
    "\n",
    "        if hemi == 'left':\n",
    "            left_hemi_widget = widget\n",
    "        else:\n",
    "            right_hemi_widget = widget\n",
    "\n",
    "    # Display the widgets side by side  \n",
    "    with out_brain:\n",
    "        display(widgets.HBox([left_hemi_widget, right_hemi_widget]))\n",
    "\n",
    "def plot_bar_for_thresholded_regions(z_data, dtype, thresh):\n",
    "\n",
    "    prominent_regions = [col for col in z_data.columns if z_data[col].abs().mean() > thresh]\n",
    "\n",
    "    # Sort in descending order\n",
    "    prominent_regions = sorted(prominent_regions, key=lambda x: z_data[x].mean(), reverse=True)\n",
    "\n",
    "    # set color to be red if z-score is positive, blue if negative\n",
    "    # color_map = lambda z: 'red' if z > 0 else 'blue'\n",
    "    fig = plt.figure(figsize=(10, 5))\n",
    "    plt.bar(prominent_regions, z_data[prominent_regions].values[0], color=(0.2, 0.4, 0.6, 0.6))\n",
    "    plt.ylabel('Z-Score')\n",
    "    plt.title(f'Regions with Abs(Z-Score) > {thresh}')\n",
    "    plt.xticks(rotation=90)\n",
    "    # add a horizontal line at z=0\n",
    "    plt.axhline(y=0, color='black', linestyle='--')\n",
    "    plt.tight_layout()\n",
    "    plt.close(fig)\n",
    "    # make the table into a widget object\n",
    "    # plot_widget = widgets.HTML(value=plt.gcf().get_iframe().decode())\n",
    "\n",
    "    # with out_table:\n",
    "    #     # plt.show()\n",
    "    #     # display in a horizontal box with the table\n",
    "    #     display(widgets.HBox([plot_widget]))\n",
    "    # convert figure into a widget object\n",
    "\n",
    "    barplot = plot_bar_for_thresholded_regions(z_data, dtype, thresh)\n",
    "\n",
    "    with out_bar:\n",
    "        out_bar.clear_output(wait=True)\n",
    "        display(barplot)\n",
    "\n",
    "    # return fig\n",
    "\n",
    "\n",
    "def create_plot(id_number, dtype, region):\n",
    "    if dtype == 'Gray Matter Volume':\n",
    "        df = vol_df\n",
    "    elif dtype == 'Surface Area':\n",
    "        df = area_df\n",
    "    elif dtype == 'Average Cortical Thickness':\n",
    "        df = thick_df\n",
    "    elif dtype == 'Mean Curvature':\n",
    "        df = meancurv_df\n",
    "    elif dtype == 'Number Of Vertices':\n",
    "        df = numvert_df\n",
    "    elif dtype == 'Curvature Index':\n",
    "        df = curvind_df\n",
    "    elif dtype == 'Folding Index':\n",
    "        df = foldind_df\n",
    "    elif dtype == 'Gaussian Curvature':\n",
    "        df = meancurv_df\n",
    "    elif dtype == 'Thickness Std Dev':\n",
    "        df = thickstd_df\n",
    "\n",
    "    region_data = df[region]\n",
    "    \n",
    "    subject_data = df[df['id_number'] == id_number][region].values[0]\n",
    "    \n",
    "    # initialize subplots but make the right plot wider than the left\n",
    "    fig, ax = plt.subplots(1, 2, figsize=(14, 7), gridspec_kw={'width_ratios': [2, 3]})\n",
    "\n",
    "    # Scatter plot with jitter\n",
    "    sns.boxplot(y=region_data, ax=ax[0], color='lightgray', showfliers=False)\n",
    "    sns.stripplot(y=region_data, jitter=0.3, size=3, ax=ax[0], alpha=0.6)\n",
    "    ax[0].scatter(x=0, y=subject_data, color='red', s=50, label=f'Subject {id_number}: Val={subject_data:.2f}')\n",
    "    ax[0].set_title(f'Distribution of {region}')\n",
    "    ax[0].set_ylabel(dtype)\n",
    "    ax[0].set_xticks([])  # Hide x-axis ticks as they are not necessary in this context.\n",
    "    ax[0].set_xlabel('Subjects')\n",
    "    ax[0].legend()\n",
    "\n",
    "    # Distribution with z-score\n",
    "    sns.kdeplot(region_data, ax=ax[1], shade=True)\n",
    "    z_val = (subject_data - region_data.mean()) / region_data.std()\n",
    "    ax[1].axvline(x=subject_data, color='r', linestyle='--', label=f'Subject {id_number}: Z={z_val:.2f}')\n",
    "    ax[1].set_title(f'Z-Score Distribution for {region}')\n",
    "    ax[1].set_xlabel(dtype)\n",
    "    ax[1].legend()\n",
    "    \n",
    "    with out_plot:\n",
    "        out_plot.clear_output(wait=True)\n",
    "        plt.tight_layout()\n",
    "        plt.show()\n",
    "\n",
    "def create_interactive_table(id_number, z_data, dtype, thresh):\n",
    "    prominent_regions = [col for col in z_data.columns if z_data[col].abs().mean() > thresh]\n",
    "    # Sort in descending order\n",
    "    # prominent_regions.sort(key=lambda x: z_data[x].abs().mean(), reverse=True)\n",
    "    prominent_regions = sorted(prominent_regions, key=lambda x: z_data[x].mean(), reverse=True)\n",
    "\n",
    "    region_selector = widgets.Select(options=prominent_regions, description='Region:', rows=25)\n",
    "    # make region_selector wider\n",
    "    region_selector.layout.width = '400px'\n",
    "  \n",
    "    def on_region_selected(change):\n",
    "        region = change['new']\n",
    "        create_plot(id_number, dtype, region)\n",
    "\n",
    "    region_selector.observe(on_region_selected, names='value')\n",
    "\n",
    "    # barplot = plot_bar_for_thresholded_regions(z_data, dtype, thresh)\n",
    "\n",
    "    # with out_bar:\n",
    "    #     out_bar.clear_output(wait=True)\n",
    "    #     display(barplot)\n",
    "\n",
    "    with out_table:\n",
    "        out_table.clear_output(wait=True)\n",
    "        display(region_selector)\n",
    "\n",
    "    # Trigger an initial plot when the table is first created.\n",
    "    create_plot(id_number, dtype, prominent_regions[0])\n",
    "\n",
    "def submit_id(b):\n",
    "    out_brain.clear_output(wait=True)\n",
    "    # Fetch the current selected value from the radio buttons inside the function\n",
    "    dtype = data_type.value\n",
    "\n",
    "    # get the id number from the input cell\n",
    "    id_number = id.value.strip()\n",
    "    # if the id number is not a number, raise an error\n",
    "    if not id_number.isdigit():\n",
    "        raise ValueError('Please enter a number for the ID number.')\n",
    "    \n",
    "    id_number = np.int64(id_number)\n",
    "\n",
    "    z_data = zscore(id_number, dtype)\n",
    "    plot_brain(z_data, dtype)\n",
    "\n",
    "def submit_thresh(b):\n",
    "    out_table.clear_output(wait=True)\n",
    "    out_bar.clear_output(wait=True)\n",
    "    out_plot.clear_output(wait=True)\n",
    "\n",
    "    # Fetch the current selected value from the radio buttons inside the function\n",
    "    dtype = data_type.value\n",
    "\n",
    "    # get the id number from the input cell\n",
    "    id_number = id.value.strip()\n",
    "    # if the id number is not a number, raise an error\n",
    "    if not id_number.isdigit():\n",
    "        raise ValueError('Please enter a number for the ID number.')\n",
    "    \n",
    "    id_number = np.int64(id_number)\n",
    "    thresh_value = thresh.value\n",
    "    thresh_value = np.float64(thresh_value)\n",
    "\n",
    "    z_data = zscore(id_number, dtype)\n",
    "\n",
    "    # plot_bar_for_thresholded_regions(z_data, dtype, thresh_value)\n",
    "    #interactive_plot_output = create_interactive_table(id_number, z_data, dtype, thresh_value)\n",
    "    #display(interactive_plot_output)\n",
    "    create_interactive_table(id_number, z_data, dtype, thresh_value)\n",
    "    plot_bar_for_thresholded_regions(z_data, dtype, thresh)\n",
    "\n",
    "# Add an input cell to enter an id_number\n",
    "id = widgets.Text(\n",
    "    value='',\n",
    "    placeholder='Enter ID number',\n",
    "    description='ID Number:',\n",
    "    disabled=False  \n",
    "    # make it larger\n",
    "\n",
    ")\n",
    "\n",
    "# Add an input cell for the threshold\n",
    "thresh = widgets.FloatText(\n",
    "    value=2,\n",
    "    description='Z-Score Threshold:',\n",
    "    disabled=False\n",
    ")\n",
    "\n",
    "# Add radio buttons to select the data type\n",
    "data_type = widgets.RadioButtons(\n",
    "    options=['Gray Matter Volume', 'Surface Area', 'Average Cortical Thickness', 'Mean Curvature', 'Number Of Vertices', 'Curvature Index', 'Folding Index', 'Gaussian Curvature', 'Thickness Std Dev'],\n",
    "    description='Structural Metric:',\n",
    "    disabled=False,\n",
    "    value='Gray Matter Volume'\n",
    ")\n",
    "\n",
    "# create a submit button for the threshold that will update the table and plot\n",
    "thresh_button = widgets.Button(description=\"Submit\")\n",
    "thresh_button.on_click(submit_thresh)  # Bind the button click to your function\n",
    "\n",
    "# Create the submit button\n",
    "submit_button = widgets.Button(description=\"Submit\")\n",
    "submit_button.on_click(submit_id)  # Bind the button click to your function\n",
    "\n",
    "\n",
    "id.layout.width = '300px'\n",
    "thresh.layout.width = '300x'\n",
    "# make more space for the text in thresh\n",
    "id.style.description_width = '70px'\n",
    "thresh.style.description_width = '120px'\n",
    "\n",
    "box_layout = widgets.Layout(display='flex',\n",
    "                flex_flow='column',\n",
    "                align_items='flex-start',\n",
    "                width='100%')\n",
    "\n",
    "# group_id = widgets.HBox([id, submit_button], layout=box_layout)\n",
    "group_thresh = widgets.HBox([thresh, thresh_button], layout=box_layout)\n",
    "\n",
    "# display(Markdown('## Enter a subject ID number and select a structural metric to plot the z-scores on the brain.'))\n",
    "# Display widgets for Z-score brains\n",
    "# display(widgets.VBox([id, data_type, submit_button, out_brain]))\n",
    "\n",
    "# display(Markdown('## Enter a z-score threshold to plot the regions with absolute value of z-scores above the threshold.'))\n",
    "# Display widgets for threshold/plots\n",
    "# display(widgets.VBox([group_thresh, out_table, out_plot]))\n",
    "\n",
    "# # Display widgets for Z-score brains\n",
    "# display(id)\n",
    "# display(data_type)\n",
    "# display(submit_button)\n",
    "# display(out_brain)\n",
    "\n",
    "# #Display widgets for threshold/plots\n",
    "# display(thresh)\n",
    "# display(thresh_button)\n",
    "# display(out_table)\n",
    "# display(out_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "## Enter a subject ID number and select a structural metric to plot the z-scores on the brain."
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9bc20d1b8535473087fd55a6b5fc22c7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Text(value='', description='ID Number:', layout=Layout(width='300px'), placeholder='Enter ID number', style=Te…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "56283af4cabc4484b60ff53a2ce7f4f0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "RadioButtons(description='Structural Metric:', options=('Gray Matter Volume', 'Surface Area', 'Average Cortica…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a001011eba9a414dae0cfc75a325fcdb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Submit', style=ButtonStyle())"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Markdown('## Enter a subject ID number and select a structural metric to plot the z-scores on the brain.'))\n",
    "display(id, data_type, submit_button)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "46053acd9ec14f2eb48d54142607751c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(out_brain)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/markdown": [
       "## Enter a z-score threshold to plot the regions with absolute value of z-scores above the threshold."
      ],
      "text/plain": [
       "<IPython.core.display.Markdown object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3ff6ca152f304669935bb3106a5d0a1b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatText(value=2.0, description='Z-Score Threshold:', layout=Layout(width='300x'), style=Descr…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(Markdown('## Enter a z-score threshold to plot the regions with absolute value of z-scores above the threshold.'))\n",
    "display(group_thresh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "beb9a8b51230451bbc7868f0d158bb69",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(out_table)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ee52941251f8475f8d0ffa1ee9452aba",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(out_plot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/37/67fqfg3s0nx798zjptbk3rl80000gn/T/ipykernel_94216/1928321273.py:37: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  grayvol_columns['id'] = grayvol_columns['id'].astype(int)\n",
      "/var/folders/37/67fqfg3s0nx798zjptbk3rl80000gn/T/ipykernel_94216/1928321273.py:44: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  surf_columns['id'] = surf_columns['id'].astype(int)\n",
      "/var/folders/37/67fqfg3s0nx798zjptbk3rl80000gn/T/ipykernel_94216/1928321273.py:51: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  thick_columns['id'] = thick_columns['id'].astype(int)\n"
     ]
    }
   ],
   "source": [
    "# Get the data for the regression model\n",
    "\n",
    "#import cognitive data metrics\n",
    "cogdata = pd.read_csv('data/cogtotals.csv')\n",
    "\n",
    "#rename column to match imaging data file\n",
    "cogdata.rename(columns={\"Subject\": \"id\"}, inplace = True)\n",
    "cogdata.head()\n",
    "\n",
    "# delete ids without cog score - assuming missing data is represented by NaN values\n",
    "missing_data_rows = cogdata[cogdata.isnull().any(axis=1)]\n",
    "clean_cog_data = cogdata.dropna()\n",
    "clean_cog_data.to_csv('data/clean_cog_data.csv', index=False)\n",
    "clean_cog_data = pd.read_csv('data/clean_cog_data.csv')\n",
    "clean_cog_data.head()\n",
    "\n",
    "#import ids for subjects with imaging metrics\n",
    "imagingdata = pd.read_csv('data/freesurfer_diffusion_HCP.csv')\n",
    "\n",
    "#find overlap in subjects\n",
    "set1 = (clean_cog_data['id']).to_list()\n",
    "set2 = (imagingdata['id']).to_list()\n",
    "new_list = [value for value in set1 if value in set2]\n",
    "df = pd.DataFrame(new_list)\n",
    "\n",
    "#pull cog data for those that have overlapping structural data\n",
    "selected_cog_data = clean_cog_data[clean_cog_data['id'].isin(new_list)]\n",
    "#pull structural data for those that have overlapping cog data\n",
    "structuraldata = pd.read_csv('data/all_structural_metrics_HCP.csv')\n",
    "structuraldata.rename(columns={'id_number': 'id'}, inplace=True)\n",
    "selected_struc_data = structuraldata[structuraldata['id'].isin(new_list)]\n",
    "selected_struc_data = selected_struc_data.dropna()\n",
    "\n",
    "#GRAY MATTER VOLUME\n",
    "grayvol_columns = selected_struc_data.filter(like='GrayVol', axis=1)\n",
    "grayvol_columns.insert(0, 'id', structuraldata['id'])\n",
    "grayvol_columns['id'] = grayvol_columns['id'].astype(int)\n",
    "grayvol_df = grayvol_columns.sort_values(by='id')\n",
    "cognitive_gray_vol_data = pd.merge(selected_cog_data, grayvol_df)\n",
    "\n",
    "#SURFACE AREA\n",
    "surf_columns = selected_struc_data.filter(like='SurfArea', axis=1)\n",
    "surf_columns.insert(0, 'id', structuraldata['id'])\n",
    "surf_columns['id'] = surf_columns['id'].astype(int)\n",
    "surf_df = surf_columns.sort_values(by='id')\n",
    "cognitive_surf_data = pd.merge(selected_cog_data, surf_df)\n",
    "\n",
    "#THICKNESS\n",
    "thick_columns = selected_struc_data.filter(like='ThickAvg', axis=1)\n",
    "thick_columns.insert(0, 'id', structuraldata['id'])\n",
    "thick_columns['id'] = thick_columns['id'].astype(int)\n",
    "thick_df = thick_columns.sort_values(by='id')\n",
    "\n",
    "#merge cognitive data and surface area data\n",
    "cognitive_thick_data = pd.merge(selected_cog_data, thick_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Let's use this brain data to predict cognitive performance. We will try to predict the NIH Toolbox Cognitive Functioning Battery composite score adjusted for age."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c5ffda5211734ed88c5d48ac7ca14bbf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "RadioButtons(description='Structural Metric:', options=('Volume', 'Area', 'Thickness'), value='Volume')"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8ece1a4f69ff498c9f13599ea702e8bc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Button(description='Submit', style=ButtonStyle())"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "359d5ac51516427798374b5112e43622",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "out_plot_regression = widgets.Output()\n",
    "\n",
    "def run_regression(dtype):\n",
    "    if dtype == 'Volume':\n",
    "        metrics_df = grayvol_df\n",
    "        cog_df = cognitive_gray_vol_data\n",
    "    elif dtype == 'Area':\n",
    "        metrics_df = surf_df\n",
    "        cog_df = cognitive_surf_data\n",
    "    elif dtype == 'Thickness':\n",
    "        metrics_df = thick_df\n",
    "        cog_df = cognitive_thick_data\n",
    "\n",
    "    X = metrics_df.drop(columns=['id'])\n",
    "    y = cog_df['CogTotalComp_AgeAdj']\n",
    "\n",
    "    (X_train, X_test, y_train, y_test)= train_test_split(X, y, test_size = 0.3, random_state = 42)\n",
    "\n",
    "    ## Standardize the features using StandardScaler\n",
    "    scaler = StandardScaler()\n",
    "    X_train_scaled = scaler.fit_transform(X_train)\n",
    "    X_test_scaled = scaler.transform(X_test)\n",
    "\n",
    "    # Initialize a regression model (e.g., Linear Regression)\n",
    "    model = LinearRegression()\n",
    "\n",
    "    # Train the model on the training data\n",
    "    model.fit(X_train_scaled, y_train)\n",
    "\n",
    "    # Predict cognitive scores on the test data\n",
    "    y_pred = model.predict(X_test_scaled)\n",
    "    #plot_regression(y_test, y_pred)\n",
    "    return y_test, y_pred\n",
    "\n",
    "def plot_regression(y_test, y_pred):\n",
    "    # initialize subplots but make the right plot wider than the left\n",
    "    fig, ax = plt.subplots(1, 2, figsize=(14, 7), gridspec_kw={'width_ratios': [1, 1]})\n",
    "\n",
    "    # Plot the actual vs. predicted cognitive scores\n",
    "    # plt.figure(figsize=(10,6))\n",
    "    ax[0].scatter(y_test, y_pred, alpha=0.6)\n",
    "    ax[0].set_title('Actual vs. Predicted Cognitive Scores')\n",
    "    ax[0].set_xlabel('Actual Scores')\n",
    "    ax[0].set_ylabel('Predicted Scores')\n",
    "\n",
    "    # Plot the line of best fit\n",
    "    slope, intercept = np.polyfit(y_test, y_pred, 1)\n",
    "    ax[0].plot(y_test, slope*y_test + intercept, color='red', linewidth=2)\n",
    "\n",
    "    # Annotate with MSE and R2 values\n",
    "    mse = mean_squared_error(y_test, y_pred)\n",
    "    r2 = r2_score(y_test, y_pred)\n",
    "    ax[0].annotate(f'MSE: {mse:.2f}', xy=(0.05, 0.95), xycoords='axes fraction')\n",
    "    ax[0].annotate(f'R^2: {r2:.2f}', xy=(0.05, 0.90), xycoords='axes fraction')\n",
    "\n",
    "    # Plot residuals vs. predicted values\n",
    "    residuals = y_test - y_pred\n",
    "    # plt.figure(figsize=(10,6))\n",
    "    ax[1].scatter(y_pred, residuals, alpha=0.6)\n",
    "    ax[1].axhline(y=0, color='red', linestyle='--')\n",
    "    ax[1].set_title('Residuals vs. Predicted Cognitive Scores')\n",
    "    ax[1].set_xlabel('Predicted Scores')\n",
    "    ax[1].set_ylabel('Residuals')\n",
    "    \n",
    "    with out_plot_regression:\n",
    "        # out_plot_regression.clear_output(wait=True)\n",
    "        # plt.tight_layout()\n",
    "        plt.show()\n",
    "        \n",
    "button = widgets.RadioButtons(\n",
    "    options=['Volume', 'Area', 'Thickness'],\n",
    "    disabled=False,\n",
    "    value='Volume',\n",
    "    description='Structural Metric:',\n",
    ")\n",
    "\n",
    "def submit_regression(b):\n",
    "    out_plot_regression.clear_output(wait=True)\n",
    "    dtype = button.value\n",
    "    y_test, y_pred = run_regression(dtype)\n",
    "    plot_regression(y_test, y_pred)\n",
    "\n",
    "# Create the submit button\n",
    "submit_button = widgets.Button(description=\"Submit\")\n",
    "# run submit_regression when radio button is clicked\n",
    "submit_button.on_click(submit_regression)\n",
    "\n",
    "display(button)\n",
    "display(submit_button)\n",
    "display(out_plot_regression)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.9.13 ('base')",
   "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.9.6"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "7d7010f9d425683ae41e5b0b8dc0d6eb35e44e03f18d381e022146c22a330439"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}