<jupyter_start><jupyter_text>Linear Regression + TimeseriesFor this mini-lab, you will learn a bit about how to use linear regression to fit _filters_ that you can use to model timeseries. You'll be using linear regression techniques to model fMRI responses to video stimuli (the ones that we've talked about in class) based on the presence of two categories in the videos: people and buildings.<jupyter_code>%config InlineBackend.figure_format = 'retina'
import matplotlib.pyplot as plt
import numpy as np
from sklearn import datasets, linear_model
from sklearn.metrics import r2_score

# Load the fMRI dataset
datafile = np.load('data.npz')

# list all the variables in the file
print(datafile.files)
# Load all the data
# (you don't need to worry about subtracting the mean later because that's done here)

# the features (X) say whether people (column zero) or buildings (column one) are present in each video clip
X_trn = datafile['X_trn'] # time x features
X_trn -= X_trn.mean(0) # subtract the mean over time
X_test = datafile['X_test'] # time x features
X_test -= X_test.mean(0) # ditto

# y_trn and y_test have the fMRI response of one voxel
y_trn = datafile['y_trn'] # time
y_trn -= y_trn.mean(0)
y_test = datafile['y_test'] # time
y_test -= y_test.mean(0)

# ybig_trn and ybig_test have the response of 10,000 voxels
ybig_trn = datafile['ybig_trn']
ybig_trn -= ybig_trn.mean(0)
ybig_test = datafile['ybig_test']
ybig_test -= ybig_test.mean(0)
# print the size of each array<jupyter_output><empty_output><jupyter_text>Plot the regression target (y) and each of the featuresStep 1 in any analysis should _always_ be to **LOOK AT YOUR DATA**. Let's plot it and see what it looks like: what's the range, does it look uniform, etc.<jupyter_code># plot y_trn, X_trn[:,0], and X_trn[:,1] on the same axis. do you see anything?
# compute the correlation between X_trn[:,0] and y_trn (ditto for X_trn[:,1])
# plot a histogram of y_trn<jupyter_output><empty_output><jupyter_text>The Finite Impulse Response (FIR) Regression ModelOne of the key things you should know about fMRI is that it measures the Blood-Oxygen-Level Dependent (BOLD) signal, which tells you about blood flow in each area of the brain. This signal doesn't directly track neural activity! After there is a burst of neural activity in an area, nearby blood vessels slowly respond and recruit more blood over the next 2-8 seconds. So we can't directly model the fMRI response with our stimuli! We need to account for the slow [hemodynamic response](https://en.wikipedia.org/wiki/Haemodynamic_response).We can think of the hemodynamic response as (to first approximation) a _filter_ on the underlying neural activity. So if the stimulus features (here, the presence of people or buildings in a video) is correlated with the neural activity, we can try to find a filter that we can apply to the stimulus features in order to make them look like the response.We'll do this by creating new features that are _lagged_ versions of the stimulus features. If we use lags (1,2,3,4) then it's like we're modeling the response $y_t$ as a function of the stimulus features $x_{t-1}, x_{t-2}, x_{t-3}, x_{t-4}$. This is just like fitting a filter that is (effectively) convolved with the stimulus feature to get the response!<jupyter_code># let's start with just one feature, feature zero
# create a matrix that will hold 4 lagged feature vectors
X_trn_lagged = np.zeros((###SOMETHING###))
# now copy the feature X_trn[:,0] into each of the 4 columns of your new matrix
# we want the first column to be lagged by 1 timepoint relative to the original,
# the second is lagged 2, etc.

for ii,lag in enumerate([1,2,3,4]):
    # `frm` should be the timepoint that you start copying from in X_trn
    frm = ## SOMETHING ##
    # and `to` is the last timepoint
    to = ## SOMETHING
    X_trn_lagged[lag:,ii] = X_trn[frm:to,0]
# plot each of these 4 feature vectors on the same plot here to make sure they look right
# they should look identical, but shifted in time
# now do the same thing to create X_val_lagged<jupyter_output><empty_output><jupyter_text>Fit the model!Next, use `np.linalg.lstsq` to fit a linear regression model that predicts `y_trn` from `X_trn_lagged`. Save the `wt, rank, res, sing` as before.<jupyter_code># do regression here
# now plot the weights as dots connected by lines. what do these weights mean? why do they look like that?<jupyter_output><empty_output><jupyter_text>Test the model!Use the dot product between `wt` and `X_val_lagged` to predict responses in the test dataset, then compute the correlation between them. How well does our model work?<jupyter_code># create predicted test timecourse
y_test_pred = ## SOMETHING ##
# compute correlation
model1_corr = ## SOMETHING ##
print(model1_corr)<jupyter_output><empty_output><jupyter_text>Repeat the modeling procedure for every voxel in the `ybig_*` datasetThe model above was fit to just one voxel. Let's repeat this process for each of the 10,000 voxels in the bigger dataset. Ideally we would do this by reshaping the `ybig_*` datasets into a more sensible shape, but it might be easier to do using a for-loop.<jupyter_code>wt_big = np.zeros((4, 100, 100))
for ii in range(100):
    for jj in range(100):
        ## SOMETHING ##
        wt_big[:,ii,jj] = ## SOMETHING ##
# compute predictions & find the prediction correlation for every voxel in the big set
corr_big = np.zeros((100,100))
for ii in range(100):
    for jj in range(100):
        corr_big[ii,jj] = ## SOMETHING ##
# plot the correlations using plt.matshow. 
# this shows correlation across one axial slice through a subject's brain. what do you see?<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07830371282773183, 'fraction_numerical': 0.012317437972901637, 'mean_word_length': 4.841726618705036, 'pattern_counts': {'":': 0, '<': 25, '<?xml version=': 0, '>': 25, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '125472', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 1811, 'n_tokens_neox': 1715, 'n_words': 784, 'path': 'lecture38_linear_regression_lab_2/lecture38-lab_stretch.ipynb', 'repo_name': 'alexhuth/ndap-fa2018', 'size': '5683'}