% Compute voxel response without noise
nTrials = nStimuli * nRepeats;
for iTrial = 1:nTrials
% get the neural response to this stimulus, by indexing the correct column of the neuralResponse matrix
thisNeuralResponse = neuralResponse(:,stimuli(iTrial));
% multiply this by the neuronToVoxelWeights to get the voxel response on this trial. Note that you need
% to get the matrix dimensions right, so transpose is needed on thisNeuralResponse
voxelResponseNoisy(iTrial,:) = thisNeuralResponse' * neuronToVoxelWeights;
end
% add noise
noiseStandardDeviation = 0.5;
% normalize response
voxelResponseNoisy = voxelResponseNoisy / mean(voxelResponseNoisy(:));
% add gaussian noise
voxelResponseNoisy = voxelResponseNoisy + noiseStandardDeviation*randn(size(voxelResponseNoisy));
% split into train and tes
trainVoxelResponseNoisy = voxelResponseNoisy(firstHalf,:);
testVoxelResponseNoisy = voxelResponseNoisy(secondHalf:end,:);
% compute weights on train data
estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponseNoisy;
% compute model prediction on test data
modelPrediction = channelResponse(secondHalf,:) * estimatedWeights;
% compute residual
residualResponse = testVoxelResponseNoisy-modelPrediction;
% compute r2
r2 = 1-var(residualResponse(:))/var(testVoxelResponseNoisy(:))
% invert model and compute channel response
estimatedChannelResponse = testVoxelResponseNoisy * pinv(estimatedWeights);
% plot estimated channel profiles
figure;colors = hsv(nStimuli);
for iStimuli = 1:nStimuli
plot(prefOrientation,mean(estimatedChannelResponse(iStimuli:nStimuli:end,:),1),'-','Color',colors(iStimuli,:));
hold on
end
xlabel('Channel orientation preference (deg)');
ylabel('Estimated channel response (percentile of max)');
title(sprintf('r2 = %0.4f',r2));
++++
This should make estimated channel response profiles that are more realistic
++++Answer|
{{:tutorials:channelresponseprofileshinoise.png|}}
++++
===== Break here =====
For UCSB Brain Camp, if you get to this point, take a break and let us know. We'll explain the stimulus likelihood part with a few slides.
===== Stimulus likelihood =====
In the above, an important distinction is that we are inverting to estimate the model responses and not the stimulus itself. There are some cases where this may be exactly what we want, for example, when we model responses as [[https://www.physiology.org/doi/full/10.1152/jn.00540.2011|linear combinations of target and mask]] or when we want to know about [[http://serenceslab.ucsd.edu/Files/Publications/Scolari2012.pdf|off-target gain]]. But, often what we really want to know is what the population (given our model) tells us about the stimulus. To be clear, there is distinction between [[https://www.eneuro.org/content/6/2/ENEURO.0363-18.2019|reconstructing the model and reconstructing the stimulus]].
{{:tutorials:channelreconstructmodel.png?600x220|}}
We can do this using a [[http://www.jeheelab.org/bergenmaprattejehee2015_nn.pdf|Bayesian analysis]] that creates a stimulus likelihood function - that is a function that tells you the probability of any stimulus given a particular response . For this, we need to model not just the mean response which the encoding model above does, but also the variance around that mean response.
The basic idea is to look at the residuals after you have fit the encoding model. That is, one removes the prediction of the channel encoding model by subtraction and looks at the residuals:
{{:tutorials:channelresidual.png?600x360|}}
In particular, one wants to fit a multi-variate gaussian noise model to these residuals. In [[http://www.jeheelab.org/bergenmaprattejehee2015_nn.pdf|Van Bergen & Jehee, 2015]] they use a noise model which has individual variances on each voxel, covariance across voxels and variance on each one of the channels:
{{:tutorials:channelnoisemodel.png?600x350|}}
Once you know this you can figure out the probability of any stimulus given the measured response, by seeing where that response lives in the multivariate gaussian who's mean is the channel encoding model prediction for that orientation and whose covariance is modeled with the above equations:
{{:tutorials:channellikelihoodbayesian.png?600x425|}}
For the purpose of this tutorial, we will just model the noise as independent, identically distributed gaussian noise on each voxel. Of course, this is probably wrong, as voxels may have different amount of variation and will have covariation with each other - both due to spatially local effects and to common sources of signal. All of this has been modeled elegantly and validated in [[http://www.jeheelab.org/bergenjehee2017_ni.pdf|Van Bergen & Jehee, 2017]].
To compute our simple noise model, we just get the residual variance on each voxel after fitting the encoding model. So, look at the residual, voxel-by-voxel and compute the variance.
++++Matlab cheat sheet|
% split half into train and test
firstHalf = 1:round(nTrials/2);
secondHalf = round(nTrials/2)+1:nTrials;
trainVoxelResponse = voxelResponse(firstHalf,:);
testVoxelResponse = voxelResponse(secondHalf:end,:);
% compute weights on train data
estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse;
% compute weights on train data
estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse;
% compute model prediction on test data
modelPrediction = channelResponse(secondHalf,:) * estimatedWeights;
% compute residual
residualResponse = testVoxelResponseNoisy-modelPrediction;
% compute residual variance, note that this is a scalar
residualVariance = var(residualResponse(:));
% make this into a covariance matrix in which the diagonal contains the variance for each voxel
% and off diagonals (in this case all 0) contain covariance between voxels
modelCovar = eye(nVoxels)*residualVariance;
++++
Ok, now we have our full model. The weights give you the mean response expected for any orientation
stimulus and the covariance matrix tells us the variation around that mean. So, now if we want
to figure out the probability of seeing any test response given any orientation, we simply compute
the probability as a multivariate gaussian probability distribution
++++Matlab cheat sheet|
% cycle over each trial
nTestTrials = size(testVoxelResponse,1);
for iTrial = 1:nTestTrials
% now cycle over all possible orientation
for iOrientation = 1:179
% compute the mean voxel response predicted by the channel encoding model
predictedResponse = channelBasis(iOrientation,:)*estimatedWeights;
% now use that mean response and the model covariance to estimate the probability
% of seeing this orientation given the response on this trial
likelihood(iTrial,iOrientation) = mvnpdf(testVoxelResponse(iTrial,:),predictedResponse,modelCovar);
end
end
% now plot the likelihood function averaged over repeats
figure
for iStimuli = 1:nStimuli
plot(1:179,mean(likelihood(iStimuli:nStimuli:end,:),1));
hold on
end
xlabel('stimulus orientation (deg)');
ylabel('probability given trial response');
++++
Should look like something like this:
++++Answer|
{{:tutorials:channellikelihood.png|}}
++++
===== Inverted Encoding model with different channel basis functions =====
Now, let's do a little trick. We're going to show that the channel basis functions in the simple inverted encoding model above are only constrained up to a linear transform, which means that you can recreate channel response profiles of arbitrary shapes (as long as they are linearly transformable from the original unimodal exponentiated half-sinusoids we used).
First, we will design a transformation matrix that will convert are original channel basis functions into bimodal ones. Can you figure out an 8x8 transform that will do that?
++++Matlab cheat sheet|
% reweight the channels
channelReweighting = [0 0.8 0.4 0 0 0 0.4 0.8]';
% make into a full matrix xform to transform the original channels
for iChannel = 1:nChannels
xform(iChannel,:) = circshift(channelReweighting,iChannel-1);
end
% and get new bimodal channels
bimodalChannelBasis = channelBasis * xform;
++++
Now, take a look at those channels
++++Matlab cheat sheet|
% display a figure with one of the channels
figure
plot(orientations,bimodalChannelBasis(:,5));
xlabel('orientation (deg)');
ylabel('Channel response (normalized to 1)');
++++
Should look like this
++++Answer|
{{:tutorials:channelbimodal.png|}}
++++
Ok, now fit this model to the data and invert - remember, we are not changing anything about the data, we are just changing the model we are using to analyze the data.
++++Matlab cheat sheet|
% compute the channelResponse for each trial
for iTrial = 1:nTrials
channelResponse(iTrial,:) = bimodalChannelBasis(stimuli(iTrial),:);
end
% compute estimated weights
estimatedWeights = pinv(channelResponse) * voxelResponse;
% compute model prediction
modelPrediction = channelResponse * estimatedWeights;
% compute residual
residualResponse = voxelResponse-modelPrediction;
% compute r2
r2 = 1-var(residualResponse(:))/var(voxelResponse(:))
% compute estimated channel response profiles
estimatedChannelResponse = testVoxelResponse * pinv(estimatedWeights);
% and plot one of the channels averaged across all trials
figure;
plot(prefOrientation,mean(estimatedChannelResponse(5:nStimuli:end,:),1));
xlabel('Channel preferred orientation (deg)');
ylabel('Estimated channel response (percentile of full)');
title(sprintf('r2 = %0.4f',r2));
++++
Should look like the following
++++Answer|
{{:tutorials:channelbimodelestimatedchannelresponse.png|}}
++++
Remember, there was nothing in our simulated data that was bimodal. The bimodal channel response functions that we get out are simply because we modeled the responses in this way. You get out what you put in. And in this case, since the channel basis functions are linearly related to the original ones we used, they account for exactly the same amount of variance - since they span the same subspace.
Ok. So, we need to be careful about what to interpret when we invert to get the model back. That is, when we do that, we get back the model that we put in. If we use the Bayesian model from above, we can avoid the ambiguity with the channel basis functions. Compute the Bayesian model with this bimodal basis set and you will see that you get back a unimodal likelihood function.
To compute our simple noise model, we just get the residual variance on each voxel after fitting the encoding model. So,
++++Matlab cheat sheet|
% compute weights on train data
estimatedWeights = pinv(channelResponse(firstHalf,:))*trainVoxelResponse;
% compute model prediction on test data
modelPrediction = channelResponse(secondHalf,:) * estimatedWeights;
% compute residual
residualResponse = testVoxelResponseNoisy-modelPrediction;
% compute residual variance, note that this is a scalar
residualVariance = var(residualResponse(:));
% make this into a covariance matrix in which the diagonal contains the variance for each voxel
% and off diagonals (in this case all 0) contain covariance between voxels
modelCovar = eye(nVoxels)*residualVariance;
% cycle over each trial
nTestTrials = size(testVoxelResponse,1);
for iTrial = 1:nTestTrials
% now cycle over all possible orientation
for iOrientation = 1:179
% compute the mean voxel response predicted by the channel encoding model
predictedResponse = bimodalChannelBasis(iOrientation,:)*estimatedWeights;
% now use that mean response and the model covariance to estimate the probability
% of seeing this orientation given the response on this trial
likelihood(iTrial,iOrientation) = mvnpdf(testVoxelResponse(iTrial,:),predictedResponse,modelCovar);
end
end
% now plot the likelihood function averaged over repeats
figure
for iStimuli = 1:nStimuli
plot(1:179,mean(likelihood(iStimuli:nStimuli:end,:),1));
hold on
end
xlabel('stimulus orientation (deg)');
ylabel('probability given trial response');
++++
Should look like something like this:
++++Answer|
{{:tutorials:channellikelihood.png|}}
++++