Home > marsbar > @mardo > event_fitted_fir.m

event_fitted_fir

PURPOSE ^

method to compute fitted event time courses using FIR

SYNOPSIS ^

function [tc, dt] = event_fitted_fir(D, e_spec, bin_length, bin_no, opts)

DESCRIPTION ^

 method to compute fitted event time courses using FIR
 FORMAT [tc, dt] = event_fitted_fir(D, e_spec, bin_length, bin_no, opts)
 
 (defaults are in [])
 D          - design
 e_spec     - 2 by N array specifying events to combine
                 with row 1 giving session number
                 and row 2 giving event number in session
                 This may in due course become an object type
 bin_length  - duration of time bin for FIR in seconds [TR]
 bin_no      - number of time bins [24 seconds / TR]
 opts       - structure, containing fields with options
                'single' - if field present, gives single FIR 
                  This option removes any duration information, and
                  returns a simple per onset FIR model, where 1s in the
                  design matrix corresponds to 1 event at the given
                  offset.  
                'percent' - if field present, gives results as percent
                  of block means
 
 Returns
 tc         - fitted event time course, averaged over events
 dt         - time units (seconds per row in tc = bin_length)

 Here, just some notes to explain 'single' and 'stacked' FIR models.  Imagine
 you have an event of duration 10 seconds, and you want an FIR model.  To
 make things simple, let's say the TR is 1 second, and that a standard
 haemodynamic response function (HRF) lasts 24 seconds.
  
 In order to do the FIR model, there are two ways to go.  The first is to
 make an FIR model which estimates the signal (say) at every second (TR)
 after event onset, where your model (Impulse Response) lasts long enough
 to capture the event and its HRF response - say 10+24 = 24 seconds.  This
 is what I will call a 'single' FIR model.  Another approach - and this is
 what SPM does by default - is to think of the 10 second event as a (say)
 10 events one after the other, each starting 1 second after the last.
 Here the FIR model estimates the effect of one of these 1 second events,
 and the length of your FIR model (Impulse response) is just the length of
 the HRF (24 seconds).  This second approach I will call a 'stacked' FIR
 model, because the events are stacking up one upon another.
 
 The single and stacked models are the same thing, if you specify a
 duration of 0 for all your events.  If your events have different
 durations, it is very difficult to model the event response sensibly with
 a single FIR, because, for the later FIR time bins, some events will have
 stopped, and activity will be dropping to baseline, whereas other events
 will still be continuing.  In this case the stacked model can make sense,
 as it just models longer events as having more (say) 1 second events.
 However, if your events have non-zero durations, but each duration is the
 same, then it may be that you do not want the stacked model, because your
 interest is in the event time course across the whole event, rather than
 some average response which pools together responses in the start middle
 and end of your actual event response, as the stacked model does.  In such
 a case, you may want to switch to a single FIR model.

 There is an added problem for the stacked models, which is what to do
 about the actual height of the regressors.  That problem also requires
 a bit of exposition which I hope to get down to in due course.
  
 $Id$

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [tc, dt] = event_fitted_fir(D, e_spec, bin_length, bin_no, opts)
0002 % method to compute fitted event time courses using FIR
0003 % FORMAT [tc, dt] = event_fitted_fir(D, e_spec, bin_length, bin_no, opts)
0004 %
0005 % (defaults are in [])
0006 % D          - design
0007 % e_spec     - 2 by N array specifying events to combine
0008 %                 with row 1 giving session number
0009 %                 and row 2 giving event number in session
0010 %                 This may in due course become an object type
0011 % bin_length  - duration of time bin for FIR in seconds [TR]
0012 % bin_no      - number of time bins [24 seconds / TR]
0013 % opts       - structure, containing fields with options
0014 %                'single' - if field present, gives single FIR
0015 %                  This option removes any duration information, and
0016 %                  returns a simple per onset FIR model, where 1s in the
0017 %                  design matrix corresponds to 1 event at the given
0018 %                  offset.
0019 %                'percent' - if field present, gives results as percent
0020 %                  of block means
0021 %
0022 % Returns
0023 % tc         - fitted event time course, averaged over events
0024 % dt         - time units (seconds per row in tc = bin_length)
0025 %
0026 % Here, just some notes to explain 'single' and 'stacked' FIR models.  Imagine
0027 % you have an event of duration 10 seconds, and you want an FIR model.  To
0028 % make things simple, let's say the TR is 1 second, and that a standard
0029 % haemodynamic response function (HRF) lasts 24 seconds.
0030 %
0031 % In order to do the FIR model, there are two ways to go.  The first is to
0032 % make an FIR model which estimates the signal (say) at every second (TR)
0033 % after event onset, where your model (Impulse Response) lasts long enough
0034 % to capture the event and its HRF response - say 10+24 = 24 seconds.  This
0035 % is what I will call a 'single' FIR model.  Another approach - and this is
0036 % what SPM does by default - is to think of the 10 second event as a (say)
0037 % 10 events one after the other, each starting 1 second after the last.
0038 % Here the FIR model estimates the effect of one of these 1 second events,
0039 % and the length of your FIR model (Impulse response) is just the length of
0040 % the HRF (24 seconds).  This second approach I will call a 'stacked' FIR
0041 % model, because the events are stacking up one upon another.
0042 %
0043 % The single and stacked models are the same thing, if you specify a
0044 % duration of 0 for all your events.  If your events have different
0045 % durations, it is very difficult to model the event response sensibly with
0046 % a single FIR, because, for the later FIR time bins, some events will have
0047 % stopped, and activity will be dropping to baseline, whereas other events
0048 % will still be continuing.  In this case the stacked model can make sense,
0049 % as it just models longer events as having more (say) 1 second events.
0050 % However, if your events have non-zero durations, but each duration is the
0051 % same, then it may be that you do not want the stacked model, because your
0052 % interest is in the event time course across the whole event, rather than
0053 % some average response which pools together responses in the start middle
0054 % and end of your actual event response, as the stacked model does.  In such
0055 % a case, you may want to switch to a single FIR model.
0056 %
0057 % There is an added problem for the stacked models, which is what to do
0058 % about the actual height of the regressors.  That problem also requires
0059 % a bit of exposition which I hope to get down to in due course.
0060 %
0061 % $Id$
0062 
0063 if nargin < 2
0064   error('Need event specification');
0065 end
0066 if nargin < 3
0067   bin_length = [];
0068 end
0069 if nargin < 4
0070   bin_no = [];
0071 end
0072 if nargin < 5
0073   opts = [];
0074 end
0075 
0076 if ~is_fmri(D) | isempty(e_spec)
0077   tc = []; dt = [];
0078   return
0079 end
0080 if ~is_mars_estimated(D)
0081   error('Need a MarsBaR estimated design');
0082 end
0083 
0084 if size(e_spec, 1) == 1, e_spec = e_spec'; end
0085 [SN EN] = deal(1, 2);
0086 e_s_l = size(e_spec, 2);
0087 
0088 if isempty(bin_length)
0089   bin_length = tr(D);
0090 end
0091 if isempty(bin_no)
0092   bin_no = 25/bin_length;
0093 end
0094 bin_no = round(bin_no);
0095 
0096 % build a simple FIR model subpartition (X)
0097 %------------------------------------------
0098 dt          = bf_dt(D);
0099 blk_rows    = block_rows(D);
0100 oX          = design_matrix(D);
0101 [n_t_p n_eff] = size(oX);
0102 y           = summary_data(data(D));
0103 y           = apply_filter(D, y);
0104 n_rois      = size(y, 2);
0105 tc          = zeros(bin_no, n_rois);
0106 blk_mns     = block_means(D);
0107 
0108 % for each session
0109 for s = 1:length(blk_rows)
0110   sess_events = e_spec(EN, e_spec(SN, :) == s);
0111   brX         = blk_rows{s};
0112   iX_out      = [];
0113   X           = [];
0114   n_s_e       = length(sess_events);
0115   if isempty(n_s_e), break, end
0116   
0117   for ei = 1:n_s_e
0118     e           = sess_events(ei);
0119     
0120     % New design bit for FIR model for this trial type
0121     Xn          = event_x_fir(D, [s e]', bin_length, bin_no, opts);
0122     
0123     % Columns from original design that need to be removed
0124     iX_out      = [iX_out event_cols(D, [s e])];
0125     
0126     % Columns in new design matrix for basic FIR model
0127     iX_in(ei,:) = size(X, 2) + [1:size(Xn,2)];
0128     
0129     X           = [X Xn];
0130   end
0131 
0132   % put into previous design for this session, and filter
0133   %------------------------------------------------------
0134   iX0         = [1:n_eff];
0135   iX0(iX_out) = [];
0136   aX          = [X oX(brX,iX0)];
0137   KX          = apply_filter(D, aX, struct('sessions', s));
0138   
0139   % Reestimate to get FIR time courses
0140   %------------------------------------------------------
0141   xX          = spm_sp('Set',KX);
0142   pX          = spm_sp('x-',xX);
0143   betas       = pX*y(brX,:);
0144   tc_s        = betas(1:size(X,2), :);
0145   
0146   % Sum over events
0147   tc_s        = reshape(tc_s, bin_no, n_s_e, n_rois);
0148   tc_s        = squeeze(sum(tc_s, 2));  
0149   
0150   % Do percent if necessary
0151   if isfield(opts, 'percent'), tc_s = tc_s / blk_mns(s) * 100; end
0152   
0153   % Sum over sessions
0154   tc            = tc + tc_s;
0155   
0156 end
0157 tc = tc / e_s_l;
0158 dt = bin_length;

Generated on Thu 17-Mar-2016 01:14:46 by m2html © 2005