Lesson 2: Linear encoding models#

A linear encoding approach#

In the previous lesson, we saw how we can define (non-linear) encoding models and fit their parameter \(\theta\) to predict voxel responses \(x\) using non-linear descent.

It is important to point out that in most of the literature, the parameters \(\theta\) are asssumed to be fixed. In such work, researchers assume a limited set of \(m\) neural populations, each with their own set of parameters \(\theta\).

The responses in different \(n\) neural measures (e.g., fMRI voxels) are then to assume to be a linear combination of these fixed neural populations. How much each neural population contributes to each voxel is then defined by a weight matrix \(W\) of size \(m \times n\).

\[x_j = \sum W_{i, j} \cdot f_j(\theta_j)\]

The big advantage of this approach is that it allows us fit the weight matrix \(W\) using linear regression. This is a much, much faster approach than fitting the parameters \(\theta\) using non-linear gradient descent.

In this lesson, we will see how we can use the EncodingModel class to fit linear encoding models.

Setting up a linear encoding model in braincoder#

In braincoder, we can set up a linear encoding model by defining a fixed number of neural encoding populations, each with their own parameters set \(\theta_j\).

Here we use a Von Mises tuning curve to define the neural populations that are senstive to the orientation of a grating stimulus. Note that orientations are given in radians, so lie between -pi and pi.

Set up a von Mises model#

from braincoder.models import VonMisesPRF
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# Set up six evenly spaced von Mises PRFs
centers = np.linspace(0.0, 2*np.pi, 6, endpoint=False)
parameters = pd.DataFrame({'mu':centers, 'kappa':1., 'amplitude':1.0, 'baseline':0.0},
                          index=pd.Index([f'Voxel {i+1}' for i in range(6)], name='voxel'))

# We have 3 voxels, each with a linear combination of the 6 von Mises functions:
weights = np.array([[1, 0, 1],
                    [1, .5, 1],
                    [0, 1, 0],
                    [0, .5, 0],
                    [0, 0, 1],
                    [0, 0, 1]]).astype(np.float32)

model = VonMisesPRF(parameters=parameters, weights=weights)

Once we have set up the model, we can first have a look at the predictions for the 6 different basis functions:

# Note that the function `basis_functions` returns a `tensorflow` `Tensor`,
# which has to be converted to a numpy array:
orientations = np.linspace(0, np.pi*2, 100)
basis_responses = model.basis_predictions(orientations, parameters).numpy()

_ = plt.plot(orientations, basis_responses)

basis_functions

Because the model also has a \(n{\times}m\) weight matrix (number of voxels x number of neural populations), we can also use the model to predict the responses of different voxels to the same orientation stimuli:

# Each voxel timeseries is a weighted sum of the six basis functions
pred = model.predict(paradigm=orientations)
_ = plt.plot(orientations, pred)

voxel_predictions1

Fit a linear encoding model using (regularised) OLS#

To fit linear encoding models we can use the braincoder.optimize.WeightFitter. This fits weights to the model using linear regression. Note that one can also provide an alpha-parameter to the WeightFitter to regularize the weights (pull them to 0; equivalent to putting a Gaussian prior on the weights). This is often a very good idea in real data!

from braincoder.optimize import WeightFitter
from braincoder.utils import get_rsq

# Simulate data
data = model.simulate(paradigm=orientations, noise=0.1)

# Fit the weights
weight_fitter = WeightFitter(model, parameters, data, orientations)
estimated_weights = weight_fitter.fit(alpha=0.1)

# Get predictions for the fitted weights
pred = model.predict(paradigm=orientations, weights=estimated_weights)
r2 = get_rsq(data, pred)

# Plot the data and the predictions
plt.figure()
plt.plot(orientations, data, c='k')
plt.plot(orientations, pred.values, c='k', ls='--')
plt.title(f'R2 = {r2.mean():.2f}')

voxel_predictions2

Note

The complete Python script and its output can be found here.

Summary#

In this lesson, we had a look at linear encoding models. These models

  • Assume a fixed number of neural populations, each with their own parameters \(\theta_j\)

  • Every voxel then is assumed to be a linear combination of these neural populations

  • The weights of this linear combination can be fit using linear regression

  • This is much faster than fitting the parameters \(\theta_j\) using non-linear gradient descent

In the next lesson, we will see how we can add a noise model to the encoding models, which yields a likelihood function which we can invert in a principled manner.