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 = 34 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 (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$
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 = 34 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 (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;