<jupyter_start><jupyter_text>Dimensionality and neural data?Due to some combination of simple experimental paradigms and more fundamental unknown principles of cortical processing, the neural data we recorded is often much lower dimensional than the recording dimensionality. This tutorial explores a few different ways of thinking about neural dimensionality. It contrasts DCA with the commonly used PCA. Run once to download the M1 reaching dataThe file is 1.1 GB, so it may take a few minutes to download. This data is from the Sabes lab and is recordings from M1 while a monkey is reaching to different locations on a grid. This is the same data used in the DCA paper. In this tutorial, we won't cross validate the results as was done in the paper to keep things simple.More information and data can be found here:<NAME>., <NAME>., <NAME>., & <NAME>. (2017). Nonhuman Primate Reaching with Multichannel Sensorimotor Cortex Electrophysiology [Data set]. Zenodo. http://doi.org/10.5281/zenodo.583331<jupyter_code>fname = 'indy_20160627_01.mat'
if not os.path.isfile(fname): # check if file was already downloaded
    tmp = f"{fname}_tmp"
    urllib.request.urlretrieve('https://zenodo.org/record/583331/files/indy_20160627_01.mat?download=1', tmp)
    os.rename(tmp, fname)<jupyter_output><empty_output><jupyter_text>Let's load and visualize some of the dataWe'll use 50ms bins and preprocess the data by removing neurons with very low firing rates, square-root transforming the spike counts, and high-pass filtering the data to remove slow nonstationarity (30s timescale).This should load a dictionary that contains the (preprocessed) spike counts for 109 neurons along with the cursor location sampled at the same rate. We'll visualize the spike raster and cursor location for 1 minute of data.<jupyter_code>data = load_sabes_data(fname, bin_width_s=.05, preprocess=True)
keys = data.keys()
print(data.keys())
print(*[(key, data[key].shape) for key in keys])
X = data['M1']
Xn = X / X.std(axis=0, keepdims=True) # normalized version will be used later
Y = data['cursor']
plt.figure(figsize=(10, 5))
plt.imshow(X[:1200].T, extent=[0, 1199*.05, 0, 108], cmap='gray_r', aspect='auto')
plt.xlabel('Time (s)')
plt.ylabel('Neuron')

plt.figure(figsize=(5, 5))
plt.plot(*Y[:1200].T, c='k')
plt.xlabel('cursor x')
plt.ylabel('cursor y')<jupyter_output><empty_output><jupyter_text>What is the dimensionality of the neural data?There are many ways of defining dimensionality. PCA organizes the data along projections of decreasing variance explained. If the variance of your dataset is confined to a lower dimensional subspace in the original measurement space, PCA will find this manifold. DCA, on the other hand, looks for subspaces with highly predictive dynamics, as measured by Predictive Information (PI).We'll focus on the first 30 dimensions, but you could extend the analysis out to 109. We'll look at the objective of PCA (variance explained) and DCA (PI) as a function of projection dimensionality. This is a purely unsupervised analysis of dimensionality. We can also ask how well the projections (found in an unsupervised manner) can be used to predict behavior for each method.One weakness of PCA, which motivated the development of DCA, is that it cannot distinguish high variance dynamics from high variance noise. Let's first look at the variance explained by PCA projections and their $R^2$ in predicting the behavioral data. We'll also plot the max $R^2$ for a fully supervised linear method at each dimensionality.<jupyter_code>max_dim = 30
lag = 4 # 200ms lag for the neural data for predicting behavior
ds = np.arange(1, max_dim+1)
pca_model = PCA()
pca_model.fit(X)
var = np.sum(pca_model.explained_variance_)

pca1_scores = np.zeros(ds.size)
max_scores = np.zeros(ds.size)
for ii, d in enumerate(ds):
    Xd = pca_model.transform(X)[:, :d]
    rr_model = RR(alpha=1e-6)
    rr_model.fit(Xd[:-lag], Y[lag:])
    pca1_scores[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))
rr_model = RR(alpha=1e-6)
rr_model.fit(X[:-lag], Y[lag:])
max_scores[:] = r2_score(Y[lag:], rr_model.predict(X[:-lag]))
u, s, v = np.linalg.svd(rr_model.coef_)
for ii, d in enumerate(range(1, Y.shape[1])):
    rr_model.coef_ = (u[:, :d] * s[:d]) @ v[:d]
    max_scores[ii] = r2_score(Y[lag:], rr_model.predict(X[:-lag]))


plt.plot(ds, np.cumsum(pca_model.explained_variance_)[:ds.size] / var, label='Var. explained', c='C0')
plt.ylim(0, 1.01)
plt.plot(ds, pca1_scores / max_scores[-1], label=r'PCA $R^2$', c='C1')
plt.plot(ds, max_scores / max_scores[-1], label=r'Max $R^2$', c='C3')
plt.legend(loc='best')
plt.xlabel('Projected dimensionality')
plt.ylabel('0-1 normalized metric')<jupyter_output><empty_output><jupyter_text>To make visualization easier, we've normalized the y-axis so that 1 is the value of the metric at the full dimensionality of the dataset (109). So, all plots would eventually go to 1 at $d=109$. Questions:- Does the variance explained look low dimensional?- Does the $R^2$ look low dimensional?We can run the equivalent analysis for DCA. Rather than explained variance, we'll look at PI. $R^2$ is computed in the same way across methods. This analysis will be somewhat slower since the projections need to be re-fit for each dimensionality.<jupyter_code>pi = np.zeros(ds.size)
dca_scores = np.zeros(ds.size)
dca_model = DCA(T=10, d=109)
dca_model.estimate_data_statistics(X) # only need to estimate this once
max_pi = dca_model.score() # PI of data with no dimensionality reduction
for ii, d in enumerate(ds):
    print(d)
    dca_model.fit_projection(d=d)
    pi[ii] = dca_model.score()
    Xd = dca_model.transform(X)
    rr_model = RR(alpha=1e-6)
    rr_model.fit(Xd[:-lag], Y[lag:])
    dca_scores[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))
plt.plot(ds, pi / max_pi, label='PI', c='C0')
plt.plot(ds, pca1_scores / max_scores[-1], label=r'PCA $R^2$', c='C1')
plt.plot(ds, dca_scores / max_scores[-1], label=r'DCA $R^2$', c='C2')
plt.plot(ds, max_scores / max_scores[-1], label=r'Max $R^2$', c='C3')

plt.ylim(0, 1.01)
plt.legend(loc='lower right')
plt.xlabel('Projected dimensionality')
plt.ylabel('0-1 normalized metric')<jupyter_output><empty_output><jupyter_text>Questions:- Does the PI look low dimensional?- Does the $R^2$ look low dimensional? PreprocessingData analysis methods are often sensitive to certain types of preprocessing. Spiking data is sometimes variance normalized per neuron before PCA is performed. Does this change the picture of the data PCA gives?What about DCA? What does this say about the invariances encoded in the choice of metrics: PI versus variance?<jupyter_code>pca_model = PCA()
pca_model.fit(Xn) # Xn rather than X
var = np.sum(pca_model.explained_variance_)

pca2_scores = np.zeros(ds.size)
max_scores2 = np.zeros(ds.size)
for ii, d in enumerate(ds):
    Xd = pca_model.transform(Xn)[:, :d]
    rr_model = RR(alpha=1e-6)
    rr_model.fit(Xd[:-lag], Y[lag:])
    pca2_scores[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))
rr_model = RR(alpha=1e-6)
rr_model.fit(Xd[:-lag], Y[lag:])
max_scores2[:] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))
u, s, v = np.linalg.svd(rr_model.coef_)
for ii, d in enumerate(range(1, Y.shape[1])):
    rr_model.coef_ = (u[:, :d] * s[:d]) @ v[:d]
    max_scores2[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))

plt.plot(ds, np.cumsum(pca_model.explained_variance_)[:ds.size] / var, label='Var. explained', c='C0')
plt.ylim(0, 1.01)
plt.plot(ds, pca2_scores / max_scores2[-1], label=r'PCA $R^2$', c='C1')
plt.plot(ds, max_scores2 / max_scores2[-1], label=r'Max $R^2$', c='C3')
plt.legend(loc='lower right')
plt.xlabel('Projected dimensionality')
plt.ylabel('0-1 normalized metric')
pi2 = np.zeros(ds.size)
dca_scores2 = np.zeros(ds.size)
dca_model = DCA(T=10, d=109)
dca_model.estimate_data_statistics(Xn) # only need to estimate this once
for ii, d in enumerate(ds):
    print(d)
    dca_model.fit_projection(d=d)
    pi2[ii] = dca_model.score()
    Xd = dca_model.transform(Xn)
    rr_model = RR(alpha=1e-6)
    rr_model.fit(Xd[:-lag], Y[lag:])
    dca_scores[ii] = r2_score(Y[lag:], rr_model.predict(Xd[:-lag]))
plt.plot(ds, pi / max_pi, label='PI')
plt.plot(ds, pca2_scores / max_scores2[-1], label=r'PCA $R^2$', c='C1')
plt.plot(ds, dca_scores / max_scores2[-1], label=r'DCA $R^2$', c='C2')
plt.plot(ds, max_scores2 / max_scores2[-1], label=r'Max $R^2$', c='C3')
plt.ylim(0, 1.01)
plt.legend(loc='lower right')
plt.xlabel('Projected dimensionality')
plt.ylabel('0-1 normalized metric')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.13359896870971522, 'fraction_numerical': 0.02589944919723427, 'mean_word_length': 5.31680236861584, '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': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '367168', 'license': 'BSD-3-Clause-LBNL', 'n_tokens_mistral': 3238, 'n_tokens_neox': 2984, 'n_words': 922, 'path': 'notebooks/tutorial.ipynb', 'repo_name': 'JesseLivezey/DynamicalComponentsAnalysis', 'size': '8579'}