multi_model_solution.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "93fb24bf",
   "metadata": {},
   "source": [
    "# Basic linear modeling\n",
    "\n",
    "In this exercise we will run a simple regression on all voxels in a 4D\n",
    "FMRI image `ds114_sub009_t2r1.nii`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e100abd6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:17.783181Z",
     "iopub.status.busy": "2022-05-26T14:47:17.782489Z",
     "iopub.status.idle": "2022-05-26T14:47:18.627092Z",
     "shell.execute_reply": "2022-05-26T14:47:18.627848Z"
    }
   },
   "outputs": [],
   "source": [
    "#: Import some standard librares\n",
    "import numpy as np\n",
    "# Print to 4 DP\n",
    "np.set_printoptions(precision=4)\n",
    "import numpy.linalg as npl\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Set default imshow parameters\n",
    "plt.rcParams['image.cmap'] = 'gray'\n",
    "\n",
    "import nibabel as nib\n",
    "\n",
    "import nipraxis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "8980cdcf",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:18.632379Z",
     "iopub.status.busy": "2022-05-26T14:47:18.631767Z",
     "iopub.status.idle": "2022-05-26T14:47:18.819752Z",
     "shell.execute_reply": "2022-05-26T14:47:18.820237Z"
    }
   },
   "outputs": [],
   "source": [
    "# Fetch and load the data file.\n",
    "data_fname = nipraxis.fetch_file('ds114_sub009_t2r1.nii')\n",
    "img = nib.load(data_fname)\n",
    "data = img.get_fdata()\n",
    "# Knock off first four volumes.\n",
    "data = data[..., 4:]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a06faaaf",
   "metadata": {},
   "source": [
    "We make the design matrix from the convolved regressor from [Convolving with\n",
    "the hemodyamic response\n",
    "function](https://textbook.nipraxis.org/convolution_background.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "f0fd8db6",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:18.824264Z",
     "iopub.status.busy": "2022-05-26T14:47:18.823662Z",
     "iopub.status.idle": "2022-05-26T14:47:18.962058Z",
     "shell.execute_reply": "2022-05-26T14:47:18.962449Z"
    }
   },
   "outputs": [],
   "source": [
    "# Load the pre-written convolved time course\n",
    "conv_fname = nipraxis.fetch_file('ds114_sub009_t2r1_conv.txt')\n",
    "# Knock off the first four elements\n",
    "convolved = np.loadtxt(conv_fname)[4:]\n",
    "plt.plot(convolved)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "176fdcb5",
   "metadata": {},
   "source": [
    "Now compile the design matrix for a simple regression, where the regressor (independent variable) is the convolved time-course."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7d3c121d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:18.983841Z",
     "iopub.status.busy": "2022-05-26T14:47:18.983109Z",
     "iopub.status.idle": "2022-05-26T14:47:19.061474Z",
     "shell.execute_reply": "2022-05-26T14:47:19.061897Z"
    }
   },
   "outputs": [],
   "source": [
    "#- Compile the design matrix\n",
    "#- First column is convolved regressor\n",
    "#- Second column all ones\n",
    "#- Hint: investigate \"aspect\" keyword to ``plt.imshow`` for a nice\n",
    "#- looking image.\n",
    "design = np.ones((len(convolved), 2))\n",
    "design[:, 0] = convolved\n",
    "# Show an image of the resulting design matrix.\n",
    "plt.imshow(design, aspect=0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "51caa462",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.066264Z",
     "iopub.status.busy": "2022-05-26T14:47:19.065667Z",
     "iopub.status.idle": "2022-05-26T14:47:19.068037Z",
     "shell.execute_reply": "2022-05-26T14:47:19.067524Z"
    }
   },
   "outputs": [],
   "source": [
    "assert design.shape == (len(convolved), 2)\n",
    "assert np.allclose(np.mean(design, axis=0), [np.mean(convolved), 1])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "70f4d880",
   "metadata": {},
   "source": [
    "Next reshape the voxel data to two dimensions, where the first dimension is\n",
    "voxels, and the second dimension is time (volume).  Therefore there should be\n",
    "one column for each volume (time-point), and one row for each voxel.  Then transpose to get a time by voxels array.\n",
    "\n",
    "You may want to have another look at [reshaping 4D to 2D](https://textbook.nipraxis.org/reshape_and_4d.html).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "eeb48448",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.071522Z",
     "iopub.status.busy": "2022-05-26T14:47:19.070931Z",
     "iopub.status.idle": "2022-05-26T14:47:19.419300Z",
     "shell.execute_reply": "2022-05-26T14:47:19.419716Z"
    }
   },
   "outputs": [],
   "source": [
    "#- Reshape the 4D data to voxel by time 2D, and transpose to time by voxels.\n",
    "data_2d = np.reshape(data, (-1, data.shape[-1])).T\n",
    "# Show the resulting shape\n",
    "data_2d.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "3af455a0",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.424949Z",
     "iopub.status.busy": "2022-05-26T14:47:19.424319Z",
     "iopub.status.idle": "2022-05-26T14:47:19.426032Z",
     "shell.execute_reply": "2022-05-26T14:47:19.426446Z"
    }
   },
   "outputs": [],
   "source": [
    "assert len(data_2d.shape) == 2, 'Did you reshape to 2D?'\n",
    "assert data_2d.shape[0] == data.shape[-1], \\\n",
    "    'Did you transpose to time by voxels?'\n",
    "assert data_2d.shape[1] == np.prod(data.shape[:3])\n",
    "assert np.all(data_2d[:, 0] == data[0, 0, 0, :]), 'First col is first voxel'\n",
    "assert np.all(data_2d[:, -1] == data[-1, -1, -1, :]), 'Last col is last voxel'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab9854f5",
   "metadata": {},
   "source": [
    "Next estimate the betas for the design, by matrix-multiplying on the left with the pseudo-inverse of the design matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "75e7c721",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.441596Z",
     "iopub.status.busy": "2022-05-26T14:47:19.440854Z",
     "iopub.status.idle": "2022-05-26T14:47:19.449987Z",
     "shell.execute_reply": "2022-05-26T14:47:19.450428Z"
    }
   },
   "outputs": [],
   "source": [
    "#- Calculate the pseudoinverse of the design\n",
    "#- Apply to time by voxel array to get betas\n",
    "betas_2d = npl.pinv(design) @ data_2d\n",
    "# Show the resulting shape\n",
    "betas_2d.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4411199c",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.454628Z",
     "iopub.status.busy": "2022-05-26T14:47:19.453963Z",
     "iopub.status.idle": "2022-05-26T14:47:19.456014Z",
     "shell.execute_reply": "2022-05-26T14:47:19.456458Z"
    }
   },
   "outputs": [],
   "source": [
    "assert betas_2d.shape == (2, np.prod(data.shape[:3]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1f33ff4",
   "metadata": {},
   "source": [
    "You now want to reshape the betas into a four dimensional array, with the\n",
    "first three axes being voxel axes, and the last, of length 2, being the\n",
    "parameters.  So the resulting 4D array will be a stack of two 3D arrays, where the first 3D array has the parameters from the first column of the design matrix, and the second has the parameters from the second column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "60012dca",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.461108Z",
     "iopub.status.busy": "2022-05-26T14:47:19.460386Z",
     "iopub.status.idle": "2022-05-26T14:47:19.463342Z",
     "shell.execute_reply": "2022-05-26T14:47:19.463821Z"
    }
   },
   "outputs": [],
   "source": [
    "#- Transpose betas to give voxels by 2 array\n",
    "#- Reshape into 4D array, with same 3D shape as original data,\n",
    "#- last dimension length 2\n",
    "betas_4d = np.reshape(betas_2d.T, data.shape[:-1] + (-1,))\n",
    "# Show the resulting shape\n",
    "betas_4d.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dcd75a1e",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.468627Z",
     "iopub.status.busy": "2022-05-26T14:47:19.467947Z",
     "iopub.status.idle": "2022-05-26T14:47:19.469836Z",
     "shell.execute_reply": "2022-05-26T14:47:19.470371Z"
    },
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "assert betas_4d.shape[-1] == 2, 'Did you do any necessary transpose?'\n",
    "assert betas_4d.shape == data.shape[:3] + (2,)\n",
    "assert np.all(betas_4d[0, 0, 0, :] == betas_2d[:, 0]), 'First voxel is first col'\n",
    "assert np.all(betas_4d[-1, -1, -1, :] == betas_2d[:, -1]), 'Last voxel is last col'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "47a75099",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.521427Z",
     "iopub.status.busy": "2022-05-26T14:47:19.520370Z",
     "iopub.status.idle": "2022-05-26T14:47:19.611039Z",
     "shell.execute_reply": "2022-05-26T14:47:19.611479Z"
    }
   },
   "outputs": [],
   "source": [
    "#- Show the middle slice from the first beta volume\n",
    "plt.imshow(betas_4d[:, :, 14, 0], interpolation='nearest', cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9f5f007b",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2022-05-26T14:47:19.631486Z",
     "iopub.status.busy": "2022-05-26T14:47:19.630459Z",
     "iopub.status.idle": "2022-05-26T14:47:19.720362Z",
     "shell.execute_reply": "2022-05-26T14:47:19.720777Z"
    }
   },
   "outputs": [],
   "source": [
    "#- Show the middle slice from the second beta volume\n",
    "plt.imshow(betas_4d[:, :, 14, 1], interpolation='nearest', cmap='gray')"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "notebook_metadata_filter": "all,-language_info",
   "split_at_heading": true,
   "text_representation": {
    "extension": ".Rmd",
    "format_name": "rmarkdown",
    "format_version": "1.2",
    "jupytext_version": "1.13.7"
   }
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}