{
 "metadata": {
  "name": "",
  "signature": "sha256:fcc826c0d940ab56b2e6f8e8a4dd7cca2354afe2f8bdf37474ca5df2e4cbbfbc"
 },
 "name": "multivoxel",
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - compatibility with Python 3\n",
      "from __future__ import print_function  # print('me') instead of print 'me'\n",
      "from __future__ import division  # 1/2 == 0.5, not 0"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - show figures inside the notebook\n",
      "%matplotlib inline"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import common modules\n",
      "import numpy as np  # the Python array package\n",
      "import matplotlib.pyplot as plt  # the Python plotting package"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - set gray colormap and nearest neighbor interpolation by default\n",
      "plt.rcParams['image.cmap'] = 'gray'\n",
      "plt.rcParams['image.interpolation'] = 'nearest'"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import numpy.linalg with a shorter name\n",
      "import numpy.linalg as npl"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import nibabel\n",
      "import nibabel as nib"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - import stimuli, spm_funcs from pna_code\n",
      "# Call over your friendly instructor if you get an error importing\n",
      "# You might want to make sure you have the 'for-class-day9' branch checked out in pna_code\n",
      "import stimuli\n",
      "import spm_funcs"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Estimating for the whole brain"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Here we load the processed dataset.  The dataset has had the following transformations:\n",
      "\n",
      "* drop first 4 volumes;\n",
      "* slice time correction;\n",
      "* motion correction;\n",
      "* coregistration of structural to functional scans;\n",
      "* registration of structural scan to template;\n",
      "* resampling of functional and structural scans to match template.\n",
      "\n",
      "The processed functional image is `wafds114_sub009_t2r1.nii`.\n",
      "\n",
      "See `nipype_ds114_sub009_t2r1.py` for details."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - load 'wafds114_sub009_t2r1.nii'\n",
      "# Get 4D data\n",
      "img = nib.load('wafds114_sub009_t2r1.nii')\n",
      "data = img.get_data()\n",
      "data.shape"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - get number of volumes, 3D vol shape, number of voxels in one volume\n",
      "n_vols = data.shape[-1]\n",
      "vol_shape = data.shape[:3]\n",
      "n_voxels = np.prod(vol_shape)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "During processing, we dropped the first 4 volumes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - number of volumes dropped from the original 4D data\n",
      "n_dropped = 4"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are going to run our linear model on these data.  To run the model we want a two-dimensional array with one column for each voxel.  The voxel time courses are therefore in the columns.  The array will be shape (`n_vols`, `n_voxels`).\n",
      "\n",
      "Be careful - remember the difference between reshape and transpose."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Reshape data to n_voxels by n_vols, call this `data_2d`\n",
      "# Then make result into n_vols by n_voxels\n",
      "# Check shape is (n_vols, n_voxels)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we have our data ready for the estimation, we can compile our design matrix.\n",
      "\n",
      "To do this, we go back over the steps we covered in the [day6 regression](../day6/regression_solutions.ipynb) notebok."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The TR for this run was 2.5 seconds."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - TR in seconds\n",
      "TR = 2.5"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "There was only one event type for this run.  See the [ds114 data description](http://practical-neuroimaging.github.io/example_data.html#dataset-ds114) for details.\n",
      "\n",
      "The experiment was an on-off block design.  During \"on\" periods the subject saw a written noun every 3 seconds and had to think of an associated verb. During \"off\" periods subjects saw scrambled versions of the written nouns.  The run started with 10 seconds \"off\" followed by 6 repeats of (30 seconds \"on\" followed by 30 seconds \"off\").  The run finshed with 30 seconds \"on\" and 25 seconds \"off\".\n",
      "\n",
      "The event onsets in seconds and event durations are in the text file `ds114_sub009_t2r1_cond.txt`.\n",
      "\n",
      "We can load this text file and return the neural prediction using our `stimuli.events2neural` function."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use IPython to check the signature of stimuli.events2neural"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The event onsets in `ds114_sub009_t2r1_cond.txt` refer to seconds in terms of the original run.  During processing we dropped 4 scans at the beginning of the run.\n",
      "\n",
      "So, the number of volumes we need to specify to `stimuli.events2neural` is the number of volumes in the original run, before we dropped these volumes."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create neural time course prediction with \"stimuli.events2neural\"\n",
      "# Remember the extra volume!\n",
      "# Plot your neural time just to check it is what we expect."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We dropped the first four volumes, so we need to drop the first four elements of the predicted neural time course:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Remove the first four elements of the neural time course\n",
      "# Check the neural time course now has the same length as number of volumes"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Remember HRF convolution?  This is revision from [day5](day5/hrf_convolution_solutions.ipynb).  You might want to check back there to remember how we did this last time.\n",
      "\n",
      "We are going to use the hemodynamic response function in `spm_funcs.spm_hrf`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use IPython to check the signature of spm_funcs.spm_hrf"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First we need a time vector giving times at which we are going to sample our HRF function.  We want to sample from time 0, up to 35 seconds, sampling every 2.5 seconds (our TR):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make `times` vector at which to sample our HRF"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now sample the HRF at these times to get our *convolution kernel*:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Sample the HRF at given times to give convolution kernel\n",
      "# Plot the HRF values against times to make sure it looks OK"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Use `np.convolve` to convolve the predicted neural time course with the HRF samples.  This will give us a predicted hemodynamic time course.\n",
      "\n",
      "Don't forget to drop the last `length_of_kernel - 1` values that resulted from the convolution.  As you remember, these values that we are trimming are the predicted values for volumes after the end of our run."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Convolve the neural prediction with the HRF samples\n",
      "# Drop the last `length_of_kernel - 1` values\n",
      "# Check the result has the same length `n_vols`\n",
      "# Plot the result"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We have the hemodynamic prediction for the event.\n",
      "\n",
      "We now need something reasonable to model the drift of the signal over time.\n",
      "\n",
      "Start with a linear drift.  Use `np.linspace` to make a linear drift term that starts at -1 for the first volume and ends with 1 for the last volume:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Linear drift term using `np.linspace`\n",
      "# Starts at -1, ends at 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now compile the design matrix.  We want a design matrix `X` with three columns:\n",
      "\n",
      "* The hemodynamic prediction\n",
      "* The linear trend regressor\n",
      "* The vector of ones that will model the mean\n",
      "\n",
      "It should be shape `n_vols` by 3."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Compile design matrix X with columns:\n",
      "# * The hemodynamic prediction\n",
      "# * The linear trend regressor\n",
      "# * The vector of ones that will model the mean\n",
      "# Show with `plt.imshow`"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can estimate the betas on our reshaped `n_vols` by `n_voxels` data, using `npl.pinv`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Estimate betas for every voxel, using pinv\n",
      "# Check resulting betas array is 3 by n_voxels"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Define a contrast to select only the first beta.  Matrix multiply (`dot`) your contrast by the betas array to get a 1D array with length `n_voxels`."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Define contrast to select the first beta\n",
      "# Take dot of contrast with betas array to get `con_data_1d` array\n",
      "# (this is the top half of the t-statistic)\n",
      "# Check `con_data_1d` has shape (n_voxels,)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Reshape this array to the original shape of the 3D volumes of `wafds114_sub009_t2r1.nii`.\n",
      "\n",
      "Display the middle slice (over the third dimension):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Reshape contrast result to shape `vol_shape`\n",
      "# Show middle slice (over the third dimension)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next we want to get some idea of the error in our estimates.\n",
      "\n",
      "Remember our t statistic formula:\n",
      "\n",
      "$$\n",
      "t = \\frac{c^T \\hat\\beta}{\\sqrt{\\hat{\\sigma}^2 c^T (X^T X)^+ c}}\n",
      "$$\n",
      "\n",
      "We get $\\hat{\\sigma}^2$ from the data.  We will start with this estimate of error in our data.\n",
      "\n",
      "$\\hat{\\sigma}^2$ can also be called the mean residual sum of squares (MRSS)."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Our first task is to get the residual sum of squares (RSS):\n",
      "\n",
      "The residual sum of squares is :\n",
      "\n",
      "$$\n",
      "\\textrm{RSS} = \\sum_{i=1}^{n} (Y_i - \\widehat{Y_i})^2\n",
      "$$\n",
      "\n",
      "$\\widehat{Y}$ is the fitted data - that is - the data as predicted from the design and the estimated betas $\\hat\\beta$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the fitted data from the design and the betas\n",
      "# Check the fitted data is shape (n_vols, n_voxels)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The residuals are the real data minus the fitted data:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Subtract fitted data from real data to get residuals"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The RSS is sum of squared residuals.  The sum is over volumes:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the RSS\n",
      "# Check this is shape (n_voxels,)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$\\hat{\\sigma}^2$ is the mean residual sum of squares.   To get this we divide the RSS by the *degrees of freedom due to error*.\n",
      "\n",
      "This is turn is given by the number of observations (`n_vols`) minus the number of *independent* columns in the design `X`.\n",
      "\n",
      "Use `npl.matrix_rank` to calculate the degrees of freedom due to error."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the degrees of freedom due to error"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the MRSS"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we have an esimate of the error in the data, let us compare this to the contrast data.\n",
      "\n",
      "Convert the MRSS data to a 3D image, and show the middle slice over the third axis (or the same slice you showed above for the contrast data):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Reshape MRSS result to shape `vol_shape`\n",
      "# Show middle slice (over the third dimension)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Returning to our t statistic:\n",
      "\n",
      "$$\n",
      "t = \\frac{c^T \\hat\\beta}{\\sqrt{\\hat{\\sigma}^2 c^T (X^T X)^+ c}}\n",
      "$$\n",
      "\n",
      "We have $\\hat{\\sigma}^2$ (the MRSS).\n",
      "\n",
      "Now we need to calculate $c^T (X^T X)^+ c$ (the variance due to the design)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the pseudoinverse of (the dot product of X.T and X)\n",
      "# Do the left and right dot product with the c vector\n",
      "# Call the result `design_variance`"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now we can calculate the t statistic.\n",
      "\n",
      "But - before we do this - we have a small problem - which is that the data outside the brain have very small values.  This will waste processing time, and lead to division by zero errors, for voxels we are not interested in.\n",
      "\n",
      "So, first we want to make a mask to select voxels that have a reasonable signal.\n",
      "\n",
      "To find voxels with a reasonable signal, let's first take the mean over volumes of the `wafds114_sub009_t2r1.nii` 4D series:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get mean over original 4D functional image"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Next plot the histogram of the (raveled, 1D version of the) mean image.  Select a signal threshold that seems to separate the low signal from the high signal."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot a histogram of the raveled mean image.\n",
      "# Hint: remember the `bins` keyword"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Make a mask that is True for voxels that are above this threshold, False otherwise.\n",
      "\n",
      "Show the middle slice (over the third dimension).  Adjust the threshold if it doesn't look right."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make mask for mean image that is True for voxels over threshold\n",
      "# Show middle slice (over third axis)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We are going to calculate the t-statistics within this mask.\n",
      "\n",
      "* Find the contrast (top half of t-statistic) data within this mask;\n",
      "* Find the MRSS data within this mask; (use the 3D version of MRSS)\n",
      "* Use the MRSS data and `design_variance` to calculate t-statistics within the mask;\n",
      "* Make a 3D volume of zeros called `t_vol`;\n",
      "* Fill the values of `t_vol` within the mask, with the t statistic values.\n",
      "* Show the middle slice (over the third axis)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# * Find the contrast (top half of t-statistic) data within this mask;\n",
      "# * Find the MRSS data within this mask; (use the 3D version of MRSS)\n",
      "# * Use the MRSS data and `design_variance` to calculate t-statistics within the mask;\n",
      "# * Make a 3D volume of zeros called `t_vol`;\n",
      "# * Fill the values of `t_vol` within the mask, with the t statistic values.\n",
      "# * Show the middle slice (over the third axis)."
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now calculate the Bonferroni p-value threshold for a family-wise error rate of 0.05 and the number of voxels within the mask:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate the Bonferroni threshold for FWE alpha 0.05, number of voxels in mask"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Calculate corresponding t-value using the Scipy t-distribution, and the degrees of freedom due to error.  You want the inverse survival function method (`isf`) to get the t statistic value for a given probability value and degrees of freedom:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# - get t distribution code from scipy library\n",
      "from scipy.stats import t as t_dist"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Get t statistic corresponding to Bonferroni p-threshold"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Make a volume that has True in voxels where the t statistic value is greater than the Bonferroni t threshold and False otherwise.   Show the middle slice (over the third dimension)."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Make volume that is True where t statistic is greater than Bonferroni t thresh\n",
      "# Show middle slice over third dimension"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Save the t volume and contrast images to disk for us to play with later:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Save the t statistic 3D volume to disk as a NIfTI image\n",
      "# Save the contrast value 3D volume to disk as a NIfTI image"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "### Plotting with nilearn"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "See if you can follow along using the examples at:\n",
      "\n",
      "https://nilearn.github.io/auto_examples/manipulating_visualizing/plot_demo_plotting.html"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# import plotting package from nilearn"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot the t statistic image with the glass brain display"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot the t statistic image thresholded at the Bonferroni t threshold"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Use the `plot_stat_map` command to display the t image on the default template\n",
      "# Use the Bonferroni t statistic threshold"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}