\n",
" **Question 1**

\n", " The `functional_data` is also in a numpy array. What is the shape of the functional_data array? How many timepoints are in the functional data?\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " The `functional_data` is also in a numpy array. What is the shape of the functional_data array? How many timepoints are in the functional data?\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here.\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Indexing numpy arrays\n",
"\n",
"The value of individual elements in a numpy array can be obtained by giving the appropriate indices. In both our structural and functional MRI data every voxel has a position in the $x$, $y$, $z$ space defined by the three spatial axes.\n",
"\n",
"\n",
"Let's look at the value of a voxel from the structural image:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x_coordinate = 30\n",
"y_coordinate = 40\n",
"z_coordinate = 25\n",
"one_voxel_value = structural_data[x_coordinate, y_coordinate, z_coordinate]\n",
"print(one_voxel_value)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ok... so the value of that voxel is 2. What can we do with that number? On its own we can't do much, but by looking at the values of many voxels at once we can begin to look at images of our data.\n",
"\n",
"The data are 3D, so we cannot directly view an image of the whole 3D array, so instead we will look at ***slices*** through the volume.\n",
"\n",
"In the illustration of structural and functional data that we've been looking there is an example of a slice shown in orange.\n",
"\n",
"Slicing the data like this means taking all of the values in two of the dimensions for a single location in the remaining dimension. In the illustration the orange slice is all of the elements of the $y$ and $z$ dimensions for a fixed value of $x$.\n",
"\n",
"In the indexing example above we gave a single coordinate or index value for each of the x,y, and z dimensions, but we can also easily take all of the values from a dimension using `:` as shown in the cell below.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# make a numpy array of sequential numbers\n",
"my_array = np.array(range(0,12))\n",
"\n",
"# check the shape\n",
"print(my_array.shape)\n",
"\n",
"# ok there is only one dimension and it has 12 elements in it\n",
"# as is usual in python the elements in each dimension are indexed starting at zero\n",
"# let's take a look at the first element of my_array\n",
"first_element = my_array[0]\n",
"print('the first element is ' + str(first_element))\n",
"\n",
"# how about the fifth element?\n",
"fifth_element = my_array[4]\n",
"print('the first element is ' + str(fifth_element))\n",
"\n",
"# and now let's take all of the elements using the `:` indexer\n",
"all_elements = my_array[:]\n",
"print('all the elements in the array:')\n",
"print(all_elements)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's try viewing 2D slices of our structural data. \n",
"\n",
"Take a look at the code in the next cell at the line starting with `ax.imshow` and see if you can figure out what's going on with the structural_data array.\n",
"\n",
"In that line we are taking all of the elements in the first two dimensions for a single position in the third dimension specified in the variable `slice`.\n",
"\n",
"This little widget allows you to drag the slider to change the slice that we are taking, letting us move back and forth through the brain viewing a series of 2D views.\n",
"\n",
"Give it a try."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@widgets.interact(slice=widgets.IntSlider(min=0, max=53, step=1, value=0))\n",
"def plot_struct(slice):\n",
" fig, ax = plt.subplots(1, 1, figsize=[10,10])\n",
" ax.imshow(structural_data[:,:,slice],'gray')\n",
" fig.subplots_adjust(wspace=0, hspace=0)\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Alright, that looks like a brain in a skull. For your orientation we are looking at the person's brain from the top with slice 0 being the lowest one and slice 53 at the top of the head. This is known as an \"axial\" or \"transverse\" view of the brain.\n",
"\n",
"If you want to know more about this concept check out the respective Wikipedia page (https://en.wikipedia.org/wiki/Anatomical_plane)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Delete this text and put your answer here.\n", "

\n",
" **Question 2**

\n", " Use the slider to move through slices in the structural image. Can you find the eyes? What is the value of the slices where you can see the eye?\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Use the slider to move through slices in the structural image. Can you find the eyes? What is the value of the slices where you can see the eye?\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here.\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's turn back to the functional data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"functional_data.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If we want to view 2D views of the fuctional data we need to select a plane in which we \"cut\" or slice the brain like in the structural image, but we also have a fourth dimension of time that is 96 elements long. This means we have a version of our 64x64x64 box at each of 96 timepoints, representing the measured BOLD signal in each voxel at each of 96 timepoints.\n",
"\n",
"So if we want to plot we need to choose one of the spatial dimensions to slice through and also a single timepoint for which we'd like to see that plane.\n",
"\n",
"Here is a version of the slider image we looked at before using the structural data, but now using the coronal data."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@widgets.interact(slice=widgets.IntSlider(min=0, max=53, step=1, value=0))\n",
"def plot_functional(slice):\n",
" tp=3\n",
" fig, ax = plt.subplots(1, 1, figsize=[10,10])\n",
" ax.imshow(functional_data[:, :, slice, tp],'gray')\n",
" fig.subplots_adjust(wspace=0, hspace=0)\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Can you find the eyes in this image? How does it look compared to the structural image?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Delete this text and put your answer here.\n", "

\n",
" **Question 3**

\n", " Take a look at where we index the functional_data array in the `plot_functional` function. Can you tell which timepoint we are using for plotting these data?\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Take a look at where we index the functional_data array in the `plot_functional` function. Can you tell which timepoint we are using for plotting these data?\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here.\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Timecourse of activation for a single voxel\n",
"\n",
"The images we just created show the BOLD activation value in all the voxels in a slice at a single point in time. But one of the things we'd like to use fMRI data for is to evaluate whether an individual voxel's signal change over time is related to a cognitive or perceptual task.\n",
"\n",
"So let's take a single voxel and look at its' timecourse of activation. To do this we will use the same logic as for plotting a 2d slice, but now we will give a single index for the x, y, and z dimensions, corresponding to a single voxel's coordinates in 3D space, and then take all the values from the fourth dimension (time) using ':'."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"x_coord = 30\n",
"y_coord = 30\n",
"z_coord = 35\n",
"\n",
"# get all of the values over time for a single voxel\n",
"one_voxel_timecourse = functional_data[x_coord, y_coord, z_coord, :]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Delete this text and put your answer here.\n", "

\n",
" **Question 4**

\n", " What is the shape of our one_voxel_timecourse variable? What does that number correspond to?\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " What is the shape of our one_voxel_timecourse variable? What does that number correspond to?\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here.\n", "

"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Answer\n",
"one_voxel_timecourse.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's plot the timecourse of activation for that voxel."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create an empty plot with defined aspect ratio\n",
"fig, ax = plt.subplots(1, 1, figsize=[18, 5])\n",
"\n",
"x_coord = 30\n",
"y_coord = 30\n",
"z_coord = 35\n",
"\n",
"# get all of the values over time for a single voxel\n",
"one_voxel_timecourse = functional_data[x_coord, y_coord, z_coord, :]\n",
"\n",
"# Plot the timecourse of a random voxel\n",
"ax.plot(one_voxel_timecourse, lw=3)\n",
"ax.set_xlim([0, one_voxel_timecourse.shape[0]])\n",
"ax.set_xlabel('time [s]', fontsize=20)\n",
"ax.set_ylabel('signal strength', fontsize=20)\n",
"ax.set_title('voxel time course', fontsize=25)\n",
"ax.tick_params(labelsize=12)\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hmm. Looks like data. What if we plot some the signal timecourse from some more voxels?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Create an empty plot with defined aspect ratio\n",
"fig, ax = plt.subplots(1, 1, figsize=[18, 5])\n",
"\n",
"# voxel indices for one voxel\n",
"x_coord = 30\n",
"y_coord = 30\n",
"z_coord = 35\n",
"\n",
"# Plot the timecourse of a random voxel from the transveral array\n",
"ax.plot(functional_data[x_coord, y_coord, z_coord, :], lw=3)\n",
"\n",
"# choose some more voxels at random and plot 100 more voxel timecourses\n",
"# COMPREHENSION CHECK: why do we choose a random index value between 0 and 64? Why not 0 and 100?\n",
"n_voxels_to_plot=100\n",
"for _ in range(0, n_voxels_to_plot):\n",
" x_coord = np.random.randint(64)\n",
" y_coord = np.random.randint(64)\n",
" z_coord = np.random.randint(64)\n",
" ax.plot(functional_data[x_coord, y_coord, z_coord, :], lw=3)\n",
"\n",
"ax.set_xlim([0, functional_data.shape[0]])\n",
"ax.set_xlabel('time [s]', fontsize=20)\n",
"ax.set_ylabel('signal strength', fontsize=20)\n",
"ax.set_title('voxel time course', fontsize=25)\n",
"ax.tick_params(labelsize=12)\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the timecourse for 100 voxels in our participants brain. How can we make sense of all this data? \n",
"\n",
"Even with the small number of voxels plotted above we can see that the activation timecourses differ in their average level of signal as well as their dynamics over time. In order to get more insights out of this kind of data we need to put more effort into analyzing the temporal dimension of the data."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Voxel timecourses and correlation\n",
"\n",
"One of the most direct ways to understand fMRI data is to see whether different voxel timecourses correlate with some predictor variable of interest. We will unpack this throughout the rest of this lab, but it should seem familiar: we have already seen examples of computing the correlation between two sets of numbers to understand whether there is a relationship between them. \n",
"\n",
"We also saw that in some ways regression operates in a related way where we'd like to find the closest match between predictor variable (sleep) and the values of some measured outcomes (grumpiness). Later we will explore this idea further by seeing if we can relate the stimuli that a participant was perceiving to the responses in their brain but for now let's warm up by revisiting correlation and using some of our new knowledge about indexing numpy arrays."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Delete this text and put your answer here.\n", "

\n",
" **Question 5**

\n", " Get the full timecourse for two voxels in the functional data (you'll need a different set of x,y,z coordinates for each voxel and all of the values in the time dimension of the functional_data array). Take a look back at the \"Timecourse of activation for a single voxel\" section above if you need a reminder on how to get the timeseries for one voxel.\n", " \n", "Once you have the activation values for each voxel you should (A) plot them on the same plot using the sample code below; (B) compute the correlation coefficient (r value) between the two timeseries.\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Get the full timecourse for two voxels in the functional data (you'll need a different set of x,y,z coordinates for each voxel and all of the values in the time dimension of the functional_data array). Take a look back at the \"Timecourse of activation for a single voxel\" section above if you need a reminder on how to get the timeseries for one voxel.\n", " \n", "Once you have the activation values for each voxel you should (A) plot them on the same plot using the sample code below; (B) compute the correlation coefficient (r value) between the two timeseries.\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here. You can adapt the example plotting code below and add your correlation code.\n", "

"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# example for plotting two timeseries on the same plot\n",
"\n",
"# Create an empty plot with defined aspect ratio\n",
"# ax is the 'handle' for the plot and we use it to tell python where to put the data we're plotting\n",
"fig, ax = plt.subplots(1, 1, figsize=[18, 5])\n",
"\n",
"\n",
"# You'll need to get the timeseries of activation for two voxels and put them in these two variables:\n",
"voxel_1_tc = \n",
"voxel_2_tc = \n",
"\n",
"# Plot the timecourse of each voxel to the current axes\n",
"ax.plot(voxel_1_tc, lw=3, label='voxel_1_tc')\n",
"ax.plot(voxel_2_tc, lw=3, label='voxel_2_tc')\n",
"\n",
"ax.set_xlim([0, voxel_1_tc.shape[0]])\n",
"ax.set_xlabel('time [s]', fontsize=20)\n",
"ax.set_ylabel('signal strength', fontsize=20)\n",
"ax.set_title('voxel time course', fontsize=25)\n",
"ax.tick_params(labelsize=12)\n",
"ax.legend()\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Relating brain data to the experimental protocol\n",
"\n",
"Congratulations! You're on your way to analyzing fMRI data. In the example above we computed the correlation between time voxel timeseries.\n",
"\n",
"The general approach to analyzing fMRI data builds on this approach but instead of correlating a voxel timeseries with another voxel timeseries (although we will return to this later) we see whether there is a relationship between a predictor variable representing the cognitive or perceptual context over the coures of the experiment and any individual voxels. If the predictor variable composed from the cognitive or perceptual considerations of the experiment can ***explain*** the variation in a voxel timeseries we would have some evidence on which to make inferences about what that part of the brain is doing.\n",
"\n",
"The data we have been playing with in this notebook come from an experiment in which the participant alternated between listening to spoken words or hearing nothing.\n",
"\n",
"So let's think: \n",
"\n",
"How would you imagine the brain signal response to look in a study like this? If a part of the brain was involved in processing spoken verbal stimuli what would we predict its activation to look?\n",
"\n",
"For now, lets keep things simple and just assume that without a stimulus (the playing of \"bi-syllabic words\") the signal in a sound-sensitive part of the brain is at a baseline level and during the auditory stimulation the signal goes up.\n",
"\n",
"Let's make a simple \"model\" of what that would look like. \n",
"\n",
"In order to do it we will need some details about the experiment and MRI scanning. We'll come back to this later, but for now let's use them to make a timecourse of what we think it would like for a brain region to be \"active\" during auditory stimulation and \"silent\" during rest (no stimulation)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# These are the main parameters of the fMRI scan and experimental design\n",
"block_design = ['rest', 'stim']\n",
"block_size = 6\n",
"block_RT = 7\n",
"block_total = 16\n",
"block_length = block_size*block_RT\n",
"\n",
"acq_num = block_size*block_total\n",
"data_time = block_length*block_total\n",
"data_time_vol = np.arange(acq_num)*block_RT\n",
"\n",
"x_size = 64\n",
"y_size = 64"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# we can use the experimental parameters to make timecourse plot showing when the auditory stimulation happened\n",
"rest = np.zeros(block_size)\n",
"stim = np.ones(block_size)\n",
"block = np.concatenate((rest, stim), axis=0)\n",
"predicted_response = np.tile(block, int(block_total/2))\n",
"\n",
"fig, ax = plt.subplots(1, 1, figsize=[18, 5])\n",
"\n",
"ax.plot(predicted_response, lw=3)\n",
"ax.set_xlim(0, acq_num-1)\n",
"ax.set_ylim(0, 1.5)\n",
"ax.set_title('predicted response', fontsize=25)\n",
"ax.set_yticks([0,1])\n",
"ax.set_xlabel('time [volumes]', fontsize=20)\n",
"ax.tick_params(labelsize=12)\n",
"ax.tick_params(labelsize=12)\n",
"\n",
"# fig.subplots_adjust(wspace=0, hspace=0.5)\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a moment to look at this plot. This represents the timecourse of auditory stimulation in the our experiment. The x-axis represents time and the y-axis shows whether sound was on ($y=1$) or off ($y=0$).\n",
"\n",
"This is the set of values, one for each point in time that we will use to try to understand the fMRI data collected in this experiment. We will search through the brain, looking to see whether any voxel timeseries are well explained by this predicted response and if they are we might conclude that they are involved in hearing."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Before we conduct a full brain regression let's try a correlation approach.\n",
"\n",
"We will use `stats.pearsonr` to compute the correlation between between the predicted response and all the voxel timeseries.\n",
"\n",
"Actually we won't do it for all the voxel timeseries because there are alot of them and it will either take a long time or crash our computer."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Delete this text and put your answer here. You can adapt the example plotting code below and add your correlation code.\n", "

\n",
" **Question 6**

\n", " How many total voxels are in our functional dataset? In other words, how many elements are there in the 3D grid contained in the functional data array?\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " How many total voxels are in our functional dataset? In other words, how many elements are there in the 3D grid contained in the functional data array?\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here.\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Yeah.. that's too many voxels to run all at once right now. \n",
"\n",
"Let's take a slice of our data:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"one_slice_data = functional_data[:,:,36,:]\n",
"one_slice_data.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we need to turn this 3D array (x,y, and t) into a 2D array to use in np.corrcoef. We can use the numpy reshape function:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n_slice_voxels = one_slice_data.shape[0]*one_slice_data.shape[1]\n",
"one_slice_data_2d = np.reshape(one_slice_data, (n_slice_voxels, one_slice_data.shape[2]))\n",
"one_slice_data_2d.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"OK great, now we have voxels by time. Let's see if any of the voxel timeseries correlates with our predicted response"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# compute the correlation betweeen each voxel timeseries and the predicted response\n",
"# store the pearsons r value in r_values and the corresponding significance level in p_values\n",
"\n",
"r_values = []\n",
"p_values = []\n",
"\n",
"# loop over every timeseries in one_slice_data_2d\n",
"for idx, val in enumerate(range(0,one_slice_data_2d.shape[0])):\n",
" r,p = stats.pearsonr(predicted_response, one_slice_data_2d[idx,:])\n",
"# r_values.append(stats.pearsonr(predicted_response, one_slice_data_2d[idx,:])[1])\n",
"# c_2.append(stats.pearsonr(predicted_response, one_slice_data_2d[idx,:])[1])\n",
" r_values.append(r)\n",
" p_values.append(p)\n",
" \n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now we have 4096 correlation values, one for each voxel in our slice, telling us how much that voxel timeseries correlated with the predictor timeseries. Let's reshape that so we can plot the correlation values as an image.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# we'll use some numpy convenience functions to find the indices of the voxel with the highest correlation value\n",
"\n",
"# make a 2d representation of the correlation values\n",
"r_values_as_slice = np.reshape(r_values, (one_slice_data.shape[0],one_slice_data.shape[1]))\n",
"p_values_as_slice = np.reshape(p_values, (one_slice_data.shape[0],one_slice_data.shape[1]))\n",
"\n",
"# get the x,y indices of the max correlation value\n",
"max_xy = np.unravel_index(r_values_as_slice.argmax(), r_values_as_slice.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's take a look and see how the timecourse from the most correlated voxel looks compared to the predictor timecourse:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# use the max_xy values to get the voxel timecourse that has the highest correlation\n",
"# max_xy is the x and y coordinate of the voxel whose activation is most correlated with the predictor timecourse\n",
"max_corr_timecourse = one_slice_data[max_xy[0], max_xy[1],:]\n",
"\n",
"# the raw timecourse is on a different scale that the predictor response so we will scale the values so \n",
"# they are all between zero and one\n",
"\n",
"# Define the min-max scaling function\n",
"def scale(data):\n",
" return (data - data.min()) / (data.max() - data.min())\n",
"\n",
"\n",
"# Create the plots\n",
"fig, ax = plt.subplots(1,1,figsize=(15, 5))\n",
"ax.plot(scale(max_corr_timecourse), lw=3, label='scaled voxel tc')\n",
"ax.plot(predicted_response, lw=3, label='predictor tc')\n",
"ax.set_xlim(0, acq_num-1)\n",
"ax.set_xlabel('time [volumes]', fontsize=20)\n",
"ax.tick_params(labelsize=12)\n",
"ax.legend()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Hey that looks pretty good!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"OK now we can see that the timecourse for this voxel looks similar to our assumed response if a voxel was to respond to sound like in the predicted timecourse. \n",
"\n",
"But how are voxels that correlate with our idealized timecourse distributed across the brain? Lets check.\n", " Delete this text and put your answer here.\n", "

\n", "\n", "We will make a figure that has three panels:\n", "\n", "One will show the image of the brain like before, the next will show an image where each voxel is colored according to the pearsons $r$ value calculated in the correlation above, and the third image will show a \"thresholded\" map of the correlation values, showing only those voxels whose correlation value is above a desired threshold." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# use some numpy magic to get an array that is nan everywhere that the correlation value is less than .2\n", "\n", "# first make a copy of the correlation values\n", "r_thresholded = r_values_as_slice.copy()\n", "\n", "# now set any values less than .2 to be equal to nan\n", "r_thresh = .2\n", "r_thresholded[r_thresholded < r_thresh] = np.nan" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Ok now lets visualize all the maps\n", "fig, ax = plt.subplots(1,3,figsize=(18, 6))\n", "\n", "# Create the plots\n", "ax[0].imshow(functional_data[:,:,36,10], cmap='gray')\n", "ax[0].set_title('brain slice', fontsize=25)\n", "ax[0].set_xticks([])\n", "ax[0].set_yticks([])\n", "\n", "\n", "ax[1].imshow(r_values_as_slice, cmap='afmhot')\n", "ax[1].set_title('un-thresholded map', fontsize=25)\n", "ax[1].set_xticks([])\n", "ax[1].set_yticks([])\n", "\n", "\n", "ax[2].imshow(functional_data[:,:,36,10], cmap='gray')\n", "ax[2].imshow(r_thresholded, cmap='afmhot')\n", "ax[2].set_title('thresholded map (r > ' + str(r_thresh) + ')', fontsize=25)\n", "ax[2].set_yticks([])\n", "ax[2].set_xticks([])\n", "ax[2].set_yticks([])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Look at that, we made some brain maps! \n", "\n", "And it looks like there is some structure to the data. \n", "\n", "The data we're looking at is oriented so that the front of the participant head is pointing up and the left/right of their head is in the left/right of the figure.\n", "\n", "The un-thresholded map in the middle is showing correlation value for each voxel.\n", "\n", "The thresholded map is made by finding the voxels in the untresholded map whose correlation values are above our threshold of 0.2 and overlaying them on the brain slice image. The color of the voxels indicates the strength of the correlation (the $r$ value) with higher $r$ values in brighter colors.\n", "\n", "The analysis we did was to measure whether any voxels correlated with our the predicted response (the timeseries of auditory activation) and it seems like the voxels with the highest values are on the side of the head, possibly near the ears. Very cool!\n", "\n", "But you also might notice that there are other voxels that are colored here. This is the issue we will return to in the next lab when we learn about dealing with multiple comparisons and false positives when you have lots of data." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a final exercise here, let's play with the significance threshold for plotting brain maps.\n", "\n", "This next cell will define a convenience function to wrap up the thresholded plotting code we used to make the figure above, except we will switch and adjust the ***p*** value associated with the calculated correlation coeffecients ($r$).\n", "\n", "The ***p*** value corresponds to the \"signifiance\" of a single correlation test and represents the chances of observing an $r$ value that high by chance. This is the value that we use for the classic $p < 0.05$ threshold for deciding whether a statistical result is significant.\n", "\n", "In the next cell we can use the slider to change the p value threshold, and only those voxels whose correlation had a p value ***less than*** the slide value will be plotted.\n", "\n", "The voxels will be colored according to their p value, with more significant voxels (lower p values) plotted in brighter colors." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "@widgets.interact(p_threshold=widgets.FloatSlider(min=.0000000000001, max=.3, step=.00001, value=0.3, height=40))\n", "def plot_thresholded_brains(p_threshold):\n", " assert p_threshold >= 0 and p_threshold <= 1,'p_reshold should be between 0 and 1'\n", " # first make a copy of the correlation values\n", " p_thresh = p_values_as_slice.copy()\n", "\n", " # now set any voxels with p greater than p_threshold to be equal to nan\n", " p_thresh[p_thresh > p_threshold] = np.nan\n", "\n", " fig, ax = plt.subplots(1, 1, figsize=[10,10])\n", " ax.imshow(functional_data[:,:,36,10], cmap='gray')\n", " ax.imshow(p_thresh, cmap='afmhot_r')\n", " ax.set_title('thresholded map (p < ' + str(p_threshold) + ')', fontsize=25)\n", " fig.subplots_adjust(wspace=0, hspace=0)\n", " ax.set_yticks([])\n", " ax.set_xticks([])\n", " ax.set_yticks([])\n", " plt.show()\n", " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "

\n",
" **Question 7**

\n", " Play with the slider, trying out different significance thresholds. What happens if you make the threshold very lenient? What happens if you make it very strict? Pay particular attention to the spatial distribution of the voxels in the thresholded map. What do you notice about the spatial arrangement of the voxels who have the most significant correlation with the predictor timecourse? Where are those voxels compared to each other? What about the voxels that are significant at p < .29? Do you notice anything funny about where some of these voxels are relative to the underlying gray brain image?\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Play with the slider, trying out different significance thresholds. What happens if you make the threshold very lenient? What happens if you make it very strict? Pay particular attention to the spatial distribution of the voxels in the thresholded map. What do you notice about the spatial arrangement of the voxels who have the most significant correlation with the predictor timecourse? Where are those voxels compared to each other? What about the voxels that are significant at p < .29? Do you notice anything funny about where some of these voxels are relative to the underlying gray brain image?\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here.\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Multiple comparison corrections with Bonferroni\n",
"\n",
"The traditional p < .05 threshold indicates that we could expect the statistical result that is significant at p < .05 to occur by chance less than 5% of the time. If we observe a correlation that has a significance value of p < .05 we would reject the null hypothesis that the two sets of values are unrelated with the caveat that we might be wrong up to 5% of the time.\n",
"\n",
"This works pretty well when we run a single test like the sleepy and grumpy example.\n",
"\n",
"But in our MRI analysis we have run ***many*** individual correlations, each with their own r and p value.\n",
"\n",
"Adding up the 5% chance of being wrong on each test across thousands of voxels means that at p < .05 we almost certainly are finding voxel correlations that are just due to chance.\n",
"\n",
"One approach to the problem of identifying significant results when running multiple statistica tests is the *Bonferroni correction*.\n",
"\n",
"The Bonferroni correction adjusts the p value (or the alpha level corresponding to the number of Type I errors we will accept in our results) to take account of the number of tests being run:\n",
"\n",
"$p_{bonferroni} = \\frac{\\alpha}{m}$\n",
"\n",
"$\\alpha$: desired overall false positive rate\n",
"\n",
"$m$: the number of comparisons (e.g., the number of voxels)\n",
"\n",
"To compute the appropriate p level using the Bonferroni correction we divide our desired false positive rate (e.g., $.05$) by the number of separate tests being conducted. In our fMRI data the number of tests is equal to the number of voxels since we are running one correlation for every voxel timeseries.\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Delete this text and put your answer here.\n", "

\n",
" **Question 8**

\n", " In order to do an appropriate Bonferroni correction to the significance level of our MRI data we need to know how many voxels there are in our functional_data. Using the shape() function that you learned about in this notebook, calculate the number of voxels in the functional data (hint: the volume of a 3D box is the length x width x height; the same answer your computed in Question 6). Then, calculate the Bonferroni corrected p value if our desired alpha level is .05.\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " In order to do an appropriate Bonferroni correction to the significance level of our MRI data we need to know how many voxels there are in our functional_data. Using the shape() function that you learned about in this notebook, calculate the number of voxels in the functional data (hint: the volume of a 3D box is the length x width x height; the same answer your computed in Question 6). Then, calculate the Bonferroni corrected p value if our desired alpha level is .05.\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here.\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If you want to see which voxels in the previous brain slice are significant at a particular p value you can play with the next cell. First set the p value you want using the `p_val_cutoff` variable and then run the cell."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"p_val_cutoff = 0.05\n",
"plot_thresholded_brains(p_val_cutoff)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Bonferroni correction is a valid statistical procedure, but it is often overly strict and can make it difficult to observe interesting true effects.\n",
"\n",
"Another approach to dealing with the problem of multiple comparisons is to reduce the number of comparisons!\n",
"\n",
"That sounds kind of simple and when stated like that it is, but this approach is where knowledge of the domain you are working in becomes crucial. It is often the case that there are things that are measured (or could be measured) that we don't really think are relevant..\n",
"\n",
"One thing you might have already noticed is that 3D box that contains each volume of the functional data includes voxels that are not inside the brain. This is what shows up as black when we plot slices of the structural and functional data.\n",
"\n",
"The MRI scanner doesn't know what it's making an image of. It just measures a 3D grid of voxels and saves those numbers. It's up to us as researchers to focus on the interesting parts of the image.\n",
"\n",
"Thus, we can reduce the number of comparisons by not even looking at voxels that aren't inside the brain. Our question with these fMRI data is how the person's brain responded to auditory stimulation so it doesn't really make sense to look outside of the person's brain."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### \"Masking\" and MRI analysis\n",
"\n",
"In order to restrict our MRI analysis to voxels of interest (like voxels that are actually in the brain) we need a way of identifying those voxels.\n",
"\n",
"In MRI this is commonly done by using \"masks\".\n",
"\n",
"Masks are arrays of the same size and shape as the MRI data, but each voxel is labeled with a 1 or a 0 indicating whether that voxel should be included or not. The voxels with a 1 in them are part of the \"mask\".\n",
"\n",
"Let's take a look at what a \"brain mask\" looks like for our data.\n",
"\n",
"The brain mask was created by making a volume where all of the voxels inside the brain were set to 1 and all the voxels outside were set to 0.\n",
"\n",
"Let's load the mask data and take a look:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# load the mask_data and check the shape\n",
"mask_path = './'\n",
"mask_data = nibabel.load(mask_path + 'fbrain_mask.nii').get_fdata() \n",
"mask_data = np.rot90(mask_data)\n",
"mask_data[mask_data<.9]=0\n",
"mask_data[mask_data>=.9]=1\n",
"mask_data.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That shape should look familiar: it's the same as the 3D shape of our functional data.\n",
"\n",
"Here's what the mask looks like side by side with our functional_data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"@widgets.interact(slice=widgets.IntSlider(min=0, max=53, step=1, value=0))\n",
"def plot_struct(slice):\n",
" #setup the figure:\n",
" fig, ax = plt.subplots(1,2,figsize=(18, 6))\n",
" \n",
" ax[0].imshow(functional_data[:,:,slice,10], 'gray')\n",
" ax[0].set_title('brain slice (functional_data)', fontsize=25)\n",
" ax[0].set_xticks([])\n",
" ax[0].set_yticks([])\n",
" \n",
" # index the structural_data:\n",
" ax[1].imshow(mask_data[:,:,slice],'gray')\n",
" ax[1].set_title('brain slice (brain_mask)', fontsize=25)\n",
" ax[1].set_xticks([])\n",
" ax[1].set_yticks([])\n",
" fig.subplots_adjust(wspace=0, hspace=0)\n",
" plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The brain mask looks sort of like the functional_data but there is no internal structure visible. Instead, all the brain voxels are white (you'll notice occasional holes in the mask -- this comes from the automated masking procedure which uses an activation-based heuristic for identifying voxels as brain or not; we can ignore those holes for now as they are likely to indicate places where the MRI signal was compromised).\n",
"\n",
"The mask lets us focus our analysis on only those voxels that are in the brain, and this will dramatically decrease the number of voxels we are considering.\n",
"\n",
"To get the number of voxels included in the mask we can use `np.sum()`\n",
"\n",
"np.sum gives the sum of values in a data array. Our brain_mask has a 1 in any voxel that we should consider and a zero everywhere else, so the sum of the brain_mask array is the number of voxels included in the mask"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"n_mask_voxels = np.sum(mask_data)\n",
"print(n_mask_voxels)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " Delete this text and put your answer here.\n", "

\n",
" **Question 9**

\n", " What is the Bonferroni corrected p < .05 threshold if we only consider the voxels in the mask?\n", "

"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n", " What is the Bonferroni corrected p < .05 threshold if we only consider the voxels in the mask?\n", "

\n",
" **Your Answer**

\n", " Delete this text and put your answer here.\n", "

"
]
}
],
"metadata": {
"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.7.6"
},
"toc": {
"base_numbering": 1,
"nav_menu": {},
"number_sections": true,
"sideBar": true,
"skip_h1_title": false,
"title_cell": "Table of Contents",
"title_sidebar": "Contents",
"toc_cell": false,
"toc_position": {},
"toc_section_display": true,
"toc_window_display": false
}
},
"nbformat": 4,
"nbformat_minor": 2
}
\n", " Delete this text and put your answer here.\n", "