Example 3: Run GLMdenoise using finite impulse response (FIR) model

Contents

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

setup;

Load in the data

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

% Check the workspace
whos

  Name         Size                Bytes  Class     Attributes

  data         1x10            173671520  cell                
  design       1x10                 9600  cell                
  stimdur      1x1                     8  double              
  tr           1x1                     8  double              

Call GLMdenoisedata using FIR model

% It is typical to assume that the same fixed hemodynamic response function (HRF)
% applies across all experimental conditions and all voxels.  Indeed, this is the
% default behavior of GLMdenoisedata.
%
% In contrast, the finite impulse response (FIR) model makes no assumption about the
% shape of the hemodynamic response; it does this by estimating a separate parameter
% for each time point in the response to each condition.
%
% Applying the FIR model can be useful for datasets where differences in timecourses
% exist across conditions or voxels.  Moreover, even if timecourse differences do not
% necessarily exist, it may still be useful to apply the FIR model as a way of checking
% that this is indeed the case.  Also, applying the FIR model can be used as a means
% of troubleshooting (e.g., it may reveal some error in the handling or preparation
% of the data).
%
% GLMdenoisedata includes a flag that allows the FIR model to be used.  Denoising
% proceeds as usual because the FIR model simply changes the model of the experiment-
% related effects in the data.

% Here we apply the FIR model to the example dataset.  The 'fir' option indicates to
% use the FIR model and the 20 indicates to estimate timecourses up to 20 time
% points after condition onset.  Notice that the design matrix passed to GLMdenoisedata
% is in the usual format (a matrix that indicates condition onsets); there is no
% need to construct the FIR basis functions as that is done internally by the code.
% Also, notice that here we turn off bootstrapping to save computational time and
% memory requirements.
results = GLMdenoisedata(design,data,stimdur,tr, ...
            'fir',20,struct('numboots',0), ...
            'example3figures');

*** GLMdenoisedata: generating sanity-check figures. ***
*** GLMdenoisedata: performing cross-validation to determine R^2 values. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: determining noise pool. ***
*** GLMdenoisedata: calculating noise regressors. ***
*** GLMdenoisedata: performing cross-validation with 1 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 2 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 3 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 4 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 5 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 6 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 7 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 8 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 9 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 10 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 11 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 12 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 13 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 14 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 15 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 16 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 17 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 18 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 19 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: performing cross-validation with 20 PCs. ***
cross-validating model..........done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: selected number of PCs is 3. ***
*** GLMdenoisedata: fitting final model (no denoising, for comparison purposes). ***
fitting model...done.
preparing output...done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: fitting final model (with denoising). ***
fitting model...done.
preparing output...done.
computing model fits...done.
computing R^2...done.
computing SNR...done.
*** GLMdenoisedata: calculating denoised data and PC weights. ***
*** GLMdenoisedata: converting to percent BOLD change. ***
*** GLMdenoisedata: generating figures. ***

% Let's inspect the results for an example voxel.  By eye, it appears that the
% FIR timecourses for this voxel are reasonably approximated by a single HRF.
% Also, notice that the canonical HRF (red line) is substantially different from
% the FIR timecourses, suggesting that there may be value in fitting an HRF to the data.
xx = 53; yy = 27; zz = 1;
tcs = squeeze(results.modelmd(xx,yy,zz,:,:));  % 35 conditions x 21 time points
figure; hold on;
plot(0:tr:20*tr,tcs');  % plot FIR timecourses
plot(0:tr:20*tr,mean(tcs,1),'k-','LineWidth',3);  % plot the mean timecourse
hrf = getcanonicalhrf(stimdur,tr);  % compute a canonical HRF
  % scale canonical HRF to best match the mean timecourse
hrf = hrf * (pinv(hrf(1:21)')*mean(tcs,1)');
plot(0:tr:(length(hrf)-1)*tr,hrf,'r-','LineWidth',3);
straightline(0,'h','k-');
ax = axis;
axis([0 30 ax(3:4)]);
xlabel('Time from condition onset (s)');
ylabel('BOLD signal (% change)');
title('FIR timecourses (black = mean, red = canonical HRF)');

Published with MATLAB® R2012b