Example 2: Inspect the model fits obtained by analyzePRF

Contents

Download dataset (if necessary) and add analyzePRF to the MATLAB path

setup;

Load and analyze the data (this follows the example1.m script)

% Start parallel MATLAB to speed up execution
if isempty(gcp)
  parpool;
end

% Load in the data
load('exampledataset.mat');

% Upsample the data to match the stimulus rate
data = tseriesinterp(data,2,1,2);

% Analyze the data
results = analyzePRF(stimulus,data,1,struct('seedmode',[0 1],'display','off'));

*** analyzePRF: started at 28-Mar-2023 08:43:26. ***
converting stimulus (run 1) to sparse.
converting stimulus (run 2) to sparse.
converting stimulus (run 3) to sparse.
converting stimulus (run 4) to sparse.
using the following maximum polynomial degrees: [3 3 3 3]
*** fitnonlinearmodel: started at 28-Mar-2023 08:43:26. ***
*** fitnonlinearmodel: loading data. ***
*** fitnonlinearmodel: outputdir = , chunksize = 10, chunknum = 1
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [52.578 74.662 6.263 91.286 0.500 ].
      for model 2 of 2, the seed is [52.578 74.662 6.263 91.286 0.500 ].
      the estimated parameters are [52.479 74.169 8.451 103.298 0.690 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [55.614 80.765 12.825 -5.985 0.500 ].
      for model 2 of 2, the seed is [55.614 80.765 12.825 -5.985 0.500 ].
      the estimated parameters are [55.614 80.765 12.825 -5.985 0.500 ].
    seed 1 was best. final estimated parameters are [52.479 74.169 8.451 103.298 0.690 ].
    trainperformance is 68.24. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [49.930 57.225 8.953 21.124 0.500 ].
      for model 2 of 2, the seed is [49.930 57.225 8.953 21.124 0.500 ].
      the estimated parameters are [49.734 61.859 6.828 16.422 0.143 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [49.931 57.221 8.945 21.123 0.500 ].
      for model 2 of 2, the seed is [49.931 57.221 8.945 21.123 0.500 ].
      the estimated parameters are [49.734 61.859 6.828 16.422 0.143 ].
    seed 2 was best. final estimated parameters are [49.734 61.859 6.828 16.422 0.143 ].
    trainperformance is 32.53. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [49.283 59.081 12.871 22.731 0.500 ].
      for model 2 of 2, the seed is [49.283 59.081 12.871 22.731 0.500 ].
      the estimated parameters are [48.086 67.773 -1.208 17.512 -0.120 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [49.283 59.079 12.868 22.731 0.500 ].
      for model 2 of 2, the seed is [49.283 59.079 12.868 22.731 0.500 ].
      the estimated parameters are [48.085 67.775 -1.125 17.511 -0.120 ].
    seed 1 was best. final estimated parameters are [48.086 67.773 -1.208 17.512 -0.120 ].
    trainperformance is 27.54. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [54.763 65.409 7.189 64.523 0.500 ].
      for model 2 of 2, the seed is [54.763 65.409 7.189 64.523 0.500 ].
      the estimated parameters are [55.453 67.443 2.085 50.435 0.057 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [54.760 65.406 -7.197 64.537 0.500 ].
      for model 2 of 2, the seed is [54.760 65.406 -7.197 64.537 0.500 ].
      the estimated parameters are [55.453 67.443 -2.085 50.435 0.057 ].
    seed 1 was best. final estimated parameters are [55.453 67.443 2.085 50.435 0.057 ].
    trainperformance is 66.62. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [285.313 193.906 38.321 718243.602 0.500 ].
      for model 2 of 2, the seed is [285.313 193.906 38.321 718243.602 0.500 ].
      the estimated parameters are [289.652 196.516 19.484 718243.602 0.143 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [285.313 193.906 -38.321 718244.396 0.500 ].
      for model 2 of 2, the seed is [285.313 193.906 -38.321 718244.396 0.500 ].
      the estimated parameters are [289.652 196.516 -19.484 718244.396 0.143 ].
    seed 2 was best. final estimated parameters are [289.652 196.516 -19.484 718244.396 0.143 ].
    trainperformance is 25.72. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [8.229 99.015 -7.631 458.063 0.500 ].
      for model 2 of 2, the seed is [8.229 99.015 -7.631 458.063 0.500 ].
      the estimated parameters are [18.547 87.142 -2.304 73.010 0.184 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [35.481 65.449 -0.750 -33.759 0.500 ].
      for model 2 of 2, the seed is [35.481 65.449 -0.750 -33.759 0.500 ].
      the estimated parameters are [35.481 65.449 -0.750 -33.759 0.500 ].
    seed 1 was best. final estimated parameters are [18.547 87.142 -2.304 73.010 0.184 ].
    trainperformance is 67.28. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [43.926 79.551 26.698 97.164 0.500 ].
      for model 2 of 2, the seed is [43.926 79.551 26.698 97.164 0.500 ].
      the estimated parameters are [46.075 134.227 18.850 80.423 0.049 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [43.926 79.551 26.692 97.152 0.500 ].
      for model 2 of 2, the seed is [43.926 79.551 26.692 97.152 0.500 ].
      the estimated parameters are [46.075 134.227 18.850 80.423 0.049 ].
    seed 1 was best. final estimated parameters are [46.075 134.227 18.850 80.423 0.049 ].
    trainperformance is 52.42. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [55.814 62.798 11.298 25.208 0.500 ].
      for model 2 of 2, the seed is [55.814 62.798 11.298 25.208 0.500 ].
      the estimated parameters are [69.940 83.701 -1.553 24.823 -0.099 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [52.800 52.801 -0.149 23.507 0.500 ].
      for model 2 of 2, the seed is [52.800 52.801 -0.149 23.507 0.500 ].
      the estimated parameters are [52.805 52.804 -0.135 23.507 5.051 ].
    seed 1 was best. final estimated parameters are [69.940 83.701 -1.553 24.823 -0.099 ].
    trainperformance is 9.90. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [49.970 50.119 -0.251 23.055 0.500 ].
      for model 2 of 2, the seed is [49.970 50.119 -0.251 23.055 0.500 ].
      the estimated parameters are [49.995 50.006 -0.027 23.130 0.065 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [49.060 50.978 -0.133 14.787 0.500 ].
      for model 2 of 2, the seed is [49.060 50.978 -0.133 14.787 0.500 ].
      the estimated parameters are [49.031 50.989 -0.043 14.788 0.230 ].
    seed 2 was best. final estimated parameters are [49.031 50.989 -0.043 14.788 0.230 ].
    trainperformance is 22.19. testperformance is NaN.
  aggregatedtestperformance is NaN.
  starting resampling case 1 of 1.
    trying seed 1 of 2.
      for model 1 of 2, the seed is [50.500 50.500 17.678 10.000 0.500 ].
      the estimated parameters are [50.427 50.551 -0.440 12.972 0.500 ].
      for model 2 of 2, the seed is [50.427 50.551 -0.440 12.972 0.500 ].
      the estimated parameters are [50.455 50.535 -0.248 18.342 0.191 ].
    trying seed 2 of 2.
      for model 1 of 2, the seed is [50.500 50.500 1.768 10.000 0.500 ].
      the estimated parameters are [50.427 50.551 -0.443 13.030 0.500 ].
      for model 2 of 2, the seed is [50.427 50.551 -0.443 13.030 0.500 ].
      the estimated parameters are [50.455 50.535 -0.248 18.343 0.191 ].
    seed 2 was best. final estimated parameters are [50.455 50.535 -0.248 18.343 0.191 ].
    trainperformance is 14.22. testperformance is NaN.
  aggregatedtestperformance is NaN.
*** fitnonlinearmodel: ended at 28-Mar-2023 08:43:46 (0.3 minutes). ***
saving results to /private/var/folders/bh/zdz8wgzs3cs6lwkkljzydb_00000gn/T/tpe16b46ba_c006_4929_9d21_d6e8bcc8397a.mat (just in case).
*** analyzePRF: ended at 28-Mar-2023 08:43:46 (0.3 minutes). ***

Perform some setup

% Define some variables
res = [100 100];                    % row x column resolution of the stimuli
resmx = 100;                        % maximum resolution (along any dimension)
hrf = results.options.hrf;          % HRF that was used in the model
degs = results.options.maxpolydeg;  % vector of maximum polynomial degrees used in the model

% Pre-compute cache for faster execution
[d,xx,yy] = makegaussian2d(resmx,2,2,2,2);

% Prepare the stimuli for use in the model
stimulusPP = {};
for p=1:length(stimulus)
  stimulusPP{p} = squish(stimulus{p},2)';  % this flattens the image so that the dimensionality is now frames x pixels
  stimulusPP{p} = [stimulusPP{p} p*ones(size(stimulusPP{p},1),1)];  % this adds a dummy column to indicate run breaks
end

% Define the model function.  This function takes parameters and stimuli as input and
% returns a predicted time-series as output.  Specifically, the variable <pp> is a vector
% of parameter values (1 x 5) and the variable <dd> is a matrix with the stimuli (frames x pixels).
% Although it looks complex, what the function does is pretty straightforward: construct a
% 2D Gaussian, crop it to <res>, compute the dot-product between the stimuli and the
% Gaussian, raise the result to an exponent, and then convolve the result with the HRF,
% taking care to not bleed over run boundaries.
modelfun = @(pp,dd) conv2run(posrect(pp(4)) * (dd*[vflatten(placematrix(zeros(res),makegaussian2d(resmx,pp(1),pp(2),abs(pp(3)),abs(pp(3)),xx,yy,0,0) / (2*pi*abs(pp(3))^2))); 0]) .^ posrect(pp(5)),hrf,dd(:,prod(res)+1));

% %%%%% BRIEF ASIDE:
%
% Note that the above line of code is highly compressed and is optimized for speed,
% not necessarily readability. Here we show how the line can be broken up, omitting
% some of the less critical steps:
%
% % here, we have a set of parameters, like those present in the
% % output from analyzePRF. pp is just a 5-element vector with
% % [ROW COL SIGMA GAIN EXPONENT].
% pp = results.params(1,:,1);
%
% % make a Gaussian and give it a specific scale
% mm = makegaussian2d(resmx,pp(1),pp(2),abs(pp(3)),abs(pp(3)),xx,yy,0,0);
% mm = mm / (2*pi*abs(pp(3))^2);
%
% % prepare the stimulus (here, we take just the first run of stimulus)
% dd = stimulus{1};  % e.g. 100 x 100 x 300
% dd = squish(dd,2);  % e.g. 10000 x 300
%
% % compute the neural drive (dot product of stimulus with Gaussian,
% % raise to an exponent, apply gain factor)
% neural = posrect(pp(4))  *  (  (dd' * mm(:)) .^ posrect(pp(5)) );  % time x 1
%
% % convolve with HRF
% ts = conv(neural,hrf);
% ts = ts(1:length(neural));  % crop to original length
%
% %%%%%

% Construct projection matrices that fit and remove the polynomials.
% Note that a separate projection matrix is constructed for each run.
polymatrix = {};
for p=1:length(degs)
  polymatrix{p} = projectionmatrix(constructpolynomialmatrix(size(data{p},2),0:degs(p)));
end

Inspect the data and the model fit

% Which voxel should we inspect?  Let's inspect the second voxel.
vx = 2;

% For each run, collect the data and the model fit.  We project out polynomials
% from both the data and the model fit.  This deals with the problem of
% slow trends in the data.
datats = {};
modelts = {};
for p=1:length(data)
  datats{p} =  polymatrix{p}*data{p}(vx,:)';
  modelts{p} = polymatrix{p}*modelfun(results.params(1,:,vx),stimulusPP{p});
end

% Visualize the results
figure; hold on;
set(gcf,'Units','points','Position',[100 100 1000 100]);
plot(cat(1,datats{:}),'r-');
plot(cat(1,modelts{:}),'b-');
straightline(300*(1:4)+.5,'v','g-');
xlabel('Time (s)');
ylabel('BOLD signal');
ax = axis;
axis([.5 1200+.5 ax(3:4)]);
title('Time-series data');

Published with MATLAB® R2022a