function [results,cache] = GLMestimatemodel(design,data,stimdur,tr,hrfmodel,hrfknobs,resampling,opt,cache)
<design> is the experimental design. There are three possible cases:
1. A where A is a matrix with dimensions time x conditions.
Each column should be zeros except for ones indicating condition onsets.
(Fractional values in the design matrix are also allowed.)
2. {A1 A2 A3 ...} where each of the A's are like the previous case.
The different A's correspond to different runs, and different runs
can have different numbers of time points.
3. {{C1_1 C2_1 C3_1 ...} {C1_2 C2_2 C3_2 ...} ...} where Ca_b is a vector of
onset times (in seconds) for condition a in run b. Time starts at 0
and is coincident with the acquisition of the first volume. This case
is compatible only with <hrfmodel> set to 'assume'.
<data> is the time-series data with dimensions X x Y x Z x time or a cell vector of
elements that are each X x Y x Z x time. XYZ can be collapsed such that the data
are given as a 2D matrix (XYZ x time); however, if you do this, then several of the
figures that are written out by this function will not be useful to look at. The
dimensions of <data> should mirror that of <design>. (For example, <design> and
<data> should have the same number of runs, the same number of time points, etc.
Thus, <data> should have at least two runs since <design> must have at least two
runs.) <data> should not contain any NaNs. We automatically convert <data> to
single format (if necessary).
<stimdur> is the duration of a trial in seconds
<tr> is the sampling rate in seconds
<hrfmodel> indicates the type of model to use for the HRF:
'fir' indicates a finite impulse response model (a separate timecourse
is estimated for every voxel and every condition)
'assume' indicates that the HRF is provided (see <hrfknobs>)
'optimize' indicates that we should estimate a global HRF from the data
<hrfknobs> (optional) is as follows:
if <hrfmodel> is 'fir', then <hrfknobs> should be the number of
time points in the future to model (N >= 0). For example, if N is 10,
then timecourses will consist of 11 points, with the first point
coinciding with condition onset.
if <hrfmodel> is 'assume', then <hrfknobs> should be time x 1 with
the HRF to assume.
if <hrfmodel> is 'optimize', then <hrfknobs> should be time x 1 with the
initial seed for the HRF. The length of this vector indicates the
number of time points that we will attempt to estimate in the HRF.
Note on normalization: In the case that <hrfmodel> is 'assume' or
'optimize', we automatically divide <hrfknobs> by the maximum value
so that the peak is equal to 1. And if <hrfmodel> is 'optimize',
then after fitting the HRF, we again normalize the HRF to peak at 1
(and adjust amplitudes accordingly). Default in the case of 'fir' is
20. Default in the case of 'assume' and 'optimize' is to use a
canonical HRF that is calculated based on <stimdur> and <tr>.
<resampling> specifies the resampling scheme:
0 means fit fully (don't bootstrap or cross-validate)
A means bootstrap A times (A >= 1)
-1 means perform leave-one-run-out cross-validation (in this case, there
must be at least two runs)
<opt> (optional) is a struct with the following fields:
<extraregressors> (optional) is time x regressors or a cell vector
of elements that are each time x regressors. The dimensions of
<extraregressors> should mirror that of <design> (i.e. same number of
runs, same number of time points). The number of extra regressors
does not have to be the same across runs, and each run can have zero
or more extra regressors. If [] or not supplied, we do
not use extra regressors in the model.
<maxpolydeg> (optional) is a non-negative integer with the maximum
polynomial degree to use for polynomial nuisance functions, which
are used to capture low-frequency noise fluctuations in each run.
Can be a vector with length equal to the number of runs (this
allows you to specify different degrees for different runs).
Default is to use round(L/2) for each run where L is the
duration in minutes of a given run.
<seed> (optional) is the random number seed to use (this affects
the selection of bootstrap samples). Default: sum(100*clock).
<bootgroups> (optional) is a vector of positive integers indicating
the grouping of runs to use when bootstrapping. For example,
a grouping of [1 1 1 2 2 2] means that of the six samples that are
drawn, three samples will be drawn (with replacement) from the first
three runs and three samples will be drawn (with replacement) from
the second three runs. This functionality is useful in situations
where different runs involve different conditions. Default: ones(1,D)
where D is the number of runs.
<numforhrf> (optional) is a positive integer indicating the number
of voxels (with the best R^2 values) to consider in fitting the
global HRF. This input matters only when <hrfmodel> is 'optimize'.
Default: 50. (If there are fewer than that number of voxels
available, we just use the voxels that are available.)
<hrffitmask> (optional) is X x Y x Z with 1s indicating all possible
voxels to consider for fitting the global HRF. This input matters
only when <hrfmodel> is 'optimize'. Special case is 1 which means
all voxels can be potentially chosen. Default: 1.
<wantpercentbold> (optional) is whether to convert the amplitude estimates
in 'models', 'modelmd', and 'modelse' to percent BOLD change. This is
done as the very last step, and is accomplished by dividing by the
absolute value of 'meanvol' and multiplying by 100. (The absolute
value prevents negative values in 'meanvol' from flipping the sign.)
Default: 1.
<hrfthresh> (optional) is an R^2 threshold. If the R^2 between the estimated
HRF and the initial HRF is less than <hrfthresh>, we decide to just use
the initial HRF. Set <hrfthresh> to -Inf if you never want to reject
the estimated HRF. Default: 50.
<suppressoutput> (optional) is whether to suppress fprintf statements. Default: 0.
<cache> (optional) is used for speeding up execution. If you are calling this
function with identical inputs except potentially for different <data>, then
if you can take the <cache> returned by the first call and re-use it for
subsequent calls.
Based on the experimental design (<design>, <stimdur>, <tr>) and the model
specification (<hrfmodel>, <hrfknobs>), fit a GLM model to the data (<data>)
using a certain resampling scheme (<resampling>).
Return <results> as a struct containing the following fields:
<models> contains the full set of model estimates (e.g. all bootstrap results)
<modelmd> contains the final model estimate (median of the estimates in <models>)
<modelse> contains the standard error of the final model estimate (half of the
68% range of the estimates in <models>). Note that <modelse> will be
computed in all resampling schemes (full-fit, bootstrapping, and
cross-validation) but can be interpreted as an estimate of standard
error only in the bootstrapping scheme.
<R2> is X x Y x Z with model accuracy expressed in terms of R^2 (percentage).
In the full-fit and bootstrap cases, <R2> is an R^2 value indicating how
well the final model estimate (<modelmd>) fits the data.
In the cross-validation case, <R2> is an R^2 value indicating how well
the cross-validated predictions of the model match the data. (The
predictions and the data are each aggregated across runs before
the computation of R^2.)
<R2run> is X x Y x Z x runs with R^2 values calculated on a per-run basis.
<signal> is X x Y x Z with the maximum absolute amplitude in <modelmd>
(this is computed over all conditions and time points in the case of 'fir'
and over all conditions in the case of 'assume' and 'optimize').
<noise> is X x Y x Z with the average amplitude error in <modelse>.
<SNR> is X x Y x Z with <signal> divided by <noise>.
<hrffitvoxels> is X x Y x Z with 1s indicating the voxels used for fitting
the global HRF. This input is returned as [] if <hrfmodel> is not 'optimize'.
In the bootstrap and cross-validation cases, <hrffitvoxels> indicates the
voxels corresponding to the last iteration.
<meanvol> is X x Y x Z with the mean of all volumes
<inputs> is a struct containing all inputs used in the call to this
function, excluding <data>. We additionally include a field called
'datasize' which contains the size of each element of <data>.
Additional details on the format of <models>, <modelmd>, and <modelse>:
- If <hrfmodel> is 'fir', then model estimates consist of timecourses:
<models> is X x Y x Z x conditions x time x resamples
<modelmd> is X x Y x Z x conditions x time
<modelse> is X x Y x Z x conditions x time
- If <hrfmodel> is 'assume' or 'optimize', then model estimates consist
of HRF estimates and amplitude estimates:
<models> is {A B} where A is time x resamples (HRF estimates)
and B is X x Y x Z x conditions x resamples (amplitude estimates)
<modelmd> is {A B} where A is time x 1 and B is X x Y x Z x conditions
<modelse> is {A B} where A is time x 1 and B is X x Y x Z x conditions
Notes on model accuracy (R^2):
- We quantify the accuracy of the GLM model as the amount of variance in the
time-series data that is explained by the deterministic portion of the model,
that is, the hemodynamic responses evoked by the various experimental conditions.
Note that this does not include the nuisance components of the model, that is,
the polynomial regressors and any extra regressors provided by the user
(see opt.extraregressors).
- The metric that we use for accuracy is R^2. Specifically:
R^2 = 100 * (1-sum((data-model)^2)/sum(data^2))
- Before computing R^2 between the model and the data, we project out
polynomial regressors from both the model and the data. The purpose of
this is to reduce the influence of low-frequency fluctuations (which
can be quite large in fMRI data) on the model accuracy metric.
Notes on bootstrapping:
- Bootstrap samples are drawn from entire runs. (Bootstrapping individual
data points would be inappropriate due to temporal correlations in fMRI noise.)
For example, if there are 10 runs, each bootstrap sample consists of 10 runs
drawn with replacement from the 10 runs.
- In cases of unbalanced designs, it is possible that a bootstrap sample contains
no occurrences of a given condition; in this case, a warning is reported and
the beta weight estimated for that condition is set to zero.
Notes on the estimation of a global HRF:
- When <hrfmodel> is 'optimize', we estimate a global HRF from the data.
This is achieved using an iterative fitting strategy: First, the HRF is fixed
to the initial seed provided by <hrfknobs>, and we estimate the amplitudes
using OLS. Then, the amplitudes are fixed (to the estimates obtained in
the previous step), and we estimate the HRF using OLS. Next, the HRF is fixed
(to the estimate obtained in the previous step), and we re-estimate the
amplitudes using OLS. This process is repeated until convergence.
- The reason for the iterative fitting strategy is that the entire model
cannot be estimated at once using linear fitting techniques (and nonlinear
techniques would be too costly).
- At the HRF-estimation steps of the fitting process, the entire dataset can
in theory be fit. However, this is undesirable for two reasons. One,
fitting the entire dataset may have exorbitant memory requirements.
Two, assuming that most voxels are unrelated to the experimental paradigm
(as is typically the case in an fMRI experiment), fitting the entire dataset
will result in a poor-quality (noisy) HRF. To resolve these issues, we use
a strategy in which we determine the best voxels in terms of R^2 at a given
amplitude-estimation step and fit only these voxels in the subsequent
HRF-estimation step. The number of voxels that are chosen is controlled
by opt.numforhrf, and the pool of chosen voxels is updated at each
amplitude-estimation step.
- In some cases, the fitted global HRF may diverge wildly from the initial
seed. This may indicate extremely low SNR and/or a problem with the coding
of the experimental design and/or a poor initial seed for the HRF. If the
R^2 between the initial seed and the fitted global HRF is less than opt.hrfthresh,
we issue a warning and simply use the initial seed as the HRF (instead of
relying on the fitted global HRF). These cases should be inspected and
troubleshooted on a case-by-case basis. (In GLMdenoisedata.m, a figure
named "HRF.png" is created --- if the initial and estimated HRF are
exactly overlapping on the figure, this indicates that the exception
case occured.)
Additional information:
- In some circumstances (e.g. using a FIR model with insufficient data),
the design matrix may be singular and there is no unique solution. Our
strategy for these cases is as follows: If MATLAB issues a warning during
the inversion of the autocorrelation matrix (i.e. X'*X), then program
execution halts.
History:
- 2014/07/31: return rawdesign in cache; change cubic to pchip to avoid warnings
- 2013/12/11: now, after we are done using opt.seed, we reset the random number seed
to something random (specifically, sum(100*clock)).
- 2013/11/18: add cache input/output; update documentation; new default for
maxpolydeg (it used to always be 2); add opt.hrfthresh;
add opt.suppressoutput; some speed-ups
- 2013/05/12: allow <design> to specify onset times
- 2013/05/12: update to indicate fractional values in design matrix are allowed.
- 2013/05/12 - regressors that are all zero now receive a 0 weight (instead of crashing)
- 2013/05/12 - fixed a bug regarding how the extraregressors were being handled.
previously, the extraregressors and the polynomial regressors were being regressed
out sequentially, which is improper. now, the two regressors are being fit
simultaneously, which is the correct way to do it.
- 2012/12/06: automatically convert data to single format
- 2012/12/03: *** Tag: Version 1.02 ***. Use faster OLS computation (less
error-checking; program execution will halt if design matrix is singular);
implement various speed-ups; minor bug fixes.
- 2012/11/24:
- INPUTS: add stimdur and tr; hrfknobs is optional now; add opt.hrffitmask; add opt.wantpercentbold
- OUTPUTS: add signal,noise,SNR; add hrffitvoxels; add meanvol; add inputs
- add a speed-up (design2pre)
- 2012/11/02 - Initial version.
- 2012/10/30 - Automatic division of HRF. Ensure one complete round of fitting in optimize case.
Add sanity check on HRF.0001 function [results,cache] = GLMestimatemodel(design,data,stimdur,tr,hrfmodel,hrfknobs,resampling,opt,cache,mode) 0002 0003 % function [results,cache] = GLMestimatemodel(design,data,stimdur,tr,hrfmodel,hrfknobs,resampling,opt,cache) 0004 % 0005 % <design> is the experimental design. There are three possible cases: 0006 % 1. A where A is a matrix with dimensions time x conditions. 0007 % Each column should be zeros except for ones indicating condition onsets. 0008 % (Fractional values in the design matrix are also allowed.) 0009 % 2. {A1 A2 A3 ...} where each of the A's are like the previous case. 0010 % The different A's correspond to different runs, and different runs 0011 % can have different numbers of time points. 0012 % 3. {{C1_1 C2_1 C3_1 ...} {C1_2 C2_2 C3_2 ...} ...} where Ca_b is a vector of 0013 % onset times (in seconds) for condition a in run b. Time starts at 0 0014 % and is coincident with the acquisition of the first volume. This case 0015 % is compatible only with <hrfmodel> set to 'assume'. 0016 % <data> is the time-series data with dimensions X x Y x Z x time or a cell vector of 0017 % elements that are each X x Y x Z x time. XYZ can be collapsed such that the data 0018 % are given as a 2D matrix (XYZ x time); however, if you do this, then several of the 0019 % figures that are written out by this function will not be useful to look at. The 0020 % dimensions of <data> should mirror that of <design>. (For example, <design> and 0021 % <data> should have the same number of runs, the same number of time points, etc. 0022 % Thus, <data> should have at least two runs since <design> must have at least two 0023 % runs.) <data> should not contain any NaNs. We automatically convert <data> to 0024 % single format (if necessary). 0025 % <stimdur> is the duration of a trial in seconds 0026 % <tr> is the sampling rate in seconds 0027 % <hrfmodel> indicates the type of model to use for the HRF: 0028 % 'fir' indicates a finite impulse response model (a separate timecourse 0029 % is estimated for every voxel and every condition) 0030 % 'assume' indicates that the HRF is provided (see <hrfknobs>) 0031 % 'optimize' indicates that we should estimate a global HRF from the data 0032 % <hrfknobs> (optional) is as follows: 0033 % if <hrfmodel> is 'fir', then <hrfknobs> should be the number of 0034 % time points in the future to model (N >= 0). For example, if N is 10, 0035 % then timecourses will consist of 11 points, with the first point 0036 % coinciding with condition onset. 0037 % if <hrfmodel> is 'assume', then <hrfknobs> should be time x 1 with 0038 % the HRF to assume. 0039 % if <hrfmodel> is 'optimize', then <hrfknobs> should be time x 1 with the 0040 % initial seed for the HRF. The length of this vector indicates the 0041 % number of time points that we will attempt to estimate in the HRF. 0042 % Note on normalization: In the case that <hrfmodel> is 'assume' or 0043 % 'optimize', we automatically divide <hrfknobs> by the maximum value 0044 % so that the peak is equal to 1. And if <hrfmodel> is 'optimize', 0045 % then after fitting the HRF, we again normalize the HRF to peak at 1 0046 % (and adjust amplitudes accordingly). Default in the case of 'fir' is 0047 % 20. Default in the case of 'assume' and 'optimize' is to use a 0048 % canonical HRF that is calculated based on <stimdur> and <tr>. 0049 % <resampling> specifies the resampling scheme: 0050 % 0 means fit fully (don't bootstrap or cross-validate) 0051 % A means bootstrap A times (A >= 1) 0052 % -1 means perform leave-one-run-out cross-validation (in this case, there 0053 % must be at least two runs) 0054 % <opt> (optional) is a struct with the following fields: 0055 % <extraregressors> (optional) is time x regressors or a cell vector 0056 % of elements that are each time x regressors. The dimensions of 0057 % <extraregressors> should mirror that of <design> (i.e. same number of 0058 % runs, same number of time points). The number of extra regressors 0059 % does not have to be the same across runs, and each run can have zero 0060 % or more extra regressors. If [] or not supplied, we do 0061 % not use extra regressors in the model. 0062 % <maxpolydeg> (optional) is a non-negative integer with the maximum 0063 % polynomial degree to use for polynomial nuisance functions, which 0064 % are used to capture low-frequency noise fluctuations in each run. 0065 % Can be a vector with length equal to the number of runs (this 0066 % allows you to specify different degrees for different runs). 0067 % Default is to use round(L/2) for each run where L is the 0068 % duration in minutes of a given run. 0069 % <seed> (optional) is the random number seed to use (this affects 0070 % the selection of bootstrap samples). Default: sum(100*clock). 0071 % <bootgroups> (optional) is a vector of positive integers indicating 0072 % the grouping of runs to use when bootstrapping. For example, 0073 % a grouping of [1 1 1 2 2 2] means that of the six samples that are 0074 % drawn, three samples will be drawn (with replacement) from the first 0075 % three runs and three samples will be drawn (with replacement) from 0076 % the second three runs. This functionality is useful in situations 0077 % where different runs involve different conditions. Default: ones(1,D) 0078 % where D is the number of runs. 0079 % <numforhrf> (optional) is a positive integer indicating the number 0080 % of voxels (with the best R^2 values) to consider in fitting the 0081 % global HRF. This input matters only when <hrfmodel> is 'optimize'. 0082 % Default: 50. (If there are fewer than that number of voxels 0083 % available, we just use the voxels that are available.) 0084 % <hrffitmask> (optional) is X x Y x Z with 1s indicating all possible 0085 % voxels to consider for fitting the global HRF. This input matters 0086 % only when <hrfmodel> is 'optimize'. Special case is 1 which means 0087 % all voxels can be potentially chosen. Default: 1. 0088 % <wantpercentbold> (optional) is whether to convert the amplitude estimates 0089 % in 'models', 'modelmd', and 'modelse' to percent BOLD change. This is 0090 % done as the very last step, and is accomplished by dividing by the 0091 % absolute value of 'meanvol' and multiplying by 100. (The absolute 0092 % value prevents negative values in 'meanvol' from flipping the sign.) 0093 % Default: 1. 0094 % <hrfthresh> (optional) is an R^2 threshold. If the R^2 between the estimated 0095 % HRF and the initial HRF is less than <hrfthresh>, we decide to just use 0096 % the initial HRF. Set <hrfthresh> to -Inf if you never want to reject 0097 % the estimated HRF. Default: 50. 0098 % <suppressoutput> (optional) is whether to suppress fprintf statements. Default: 0. 0099 % <cache> (optional) is used for speeding up execution. If you are calling this 0100 % function with identical inputs except potentially for different <data>, then 0101 % if you can take the <cache> returned by the first call and re-use it for 0102 % subsequent calls. 0103 % 0104 % Based on the experimental design (<design>, <stimdur>, <tr>) and the model 0105 % specification (<hrfmodel>, <hrfknobs>), fit a GLM model to the data (<data>) 0106 % using a certain resampling scheme (<resampling>). 0107 % 0108 % Return <results> as a struct containing the following fields: 0109 % <models> contains the full set of model estimates (e.g. all bootstrap results) 0110 % <modelmd> contains the final model estimate (median of the estimates in <models>) 0111 % <modelse> contains the standard error of the final model estimate (half of the 0112 % 68% range of the estimates in <models>). Note that <modelse> will be 0113 % computed in all resampling schemes (full-fit, bootstrapping, and 0114 % cross-validation) but can be interpreted as an estimate of standard 0115 % error only in the bootstrapping scheme. 0116 % <R2> is X x Y x Z with model accuracy expressed in terms of R^2 (percentage). 0117 % In the full-fit and bootstrap cases, <R2> is an R^2 value indicating how 0118 % well the final model estimate (<modelmd>) fits the data. 0119 % In the cross-validation case, <R2> is an R^2 value indicating how well 0120 % the cross-validated predictions of the model match the data. (The 0121 % predictions and the data are each aggregated across runs before 0122 % the computation of R^2.) 0123 % <R2run> is X x Y x Z x runs with R^2 values calculated on a per-run basis. 0124 % <signal> is X x Y x Z with the maximum absolute amplitude in <modelmd> 0125 % (this is computed over all conditions and time points in the case of 'fir' 0126 % and over all conditions in the case of 'assume' and 'optimize'). 0127 % <noise> is X x Y x Z with the average amplitude error in <modelse>. 0128 % <SNR> is X x Y x Z with <signal> divided by <noise>. 0129 % <hrffitvoxels> is X x Y x Z with 1s indicating the voxels used for fitting 0130 % the global HRF. This input is returned as [] if <hrfmodel> is not 'optimize'. 0131 % In the bootstrap and cross-validation cases, <hrffitvoxels> indicates the 0132 % voxels corresponding to the last iteration. 0133 % <meanvol> is X x Y x Z with the mean of all volumes 0134 % <inputs> is a struct containing all inputs used in the call to this 0135 % function, excluding <data>. We additionally include a field called 0136 % 'datasize' which contains the size of each element of <data>. 0137 % 0138 % Additional details on the format of <models>, <modelmd>, and <modelse>: 0139 % - If <hrfmodel> is 'fir', then model estimates consist of timecourses: 0140 % <models> is X x Y x Z x conditions x time x resamples 0141 % <modelmd> is X x Y x Z x conditions x time 0142 % <modelse> is X x Y x Z x conditions x time 0143 % - If <hrfmodel> is 'assume' or 'optimize', then model estimates consist 0144 % of HRF estimates and amplitude estimates: 0145 % <models> is {A B} where A is time x resamples (HRF estimates) 0146 % and B is X x Y x Z x conditions x resamples (amplitude estimates) 0147 % <modelmd> is {A B} where A is time x 1 and B is X x Y x Z x conditions 0148 % <modelse> is {A B} where A is time x 1 and B is X x Y x Z x conditions 0149 % 0150 % Notes on model accuracy (R^2): 0151 % - We quantify the accuracy of the GLM model as the amount of variance in the 0152 % time-series data that is explained by the deterministic portion of the model, 0153 % that is, the hemodynamic responses evoked by the various experimental conditions. 0154 % Note that this does not include the nuisance components of the model, that is, 0155 % the polynomial regressors and any extra regressors provided by the user 0156 % (see opt.extraregressors). 0157 % - The metric that we use for accuracy is R^2. Specifically: 0158 % R^2 = 100 * (1-sum((data-model)^2)/sum(data^2)) 0159 % - Before computing R^2 between the model and the data, we project out 0160 % polynomial regressors from both the model and the data. The purpose of 0161 % this is to reduce the influence of low-frequency fluctuations (which 0162 % can be quite large in fMRI data) on the model accuracy metric. 0163 % 0164 % Notes on bootstrapping: 0165 % - Bootstrap samples are drawn from entire runs. (Bootstrapping individual 0166 % data points would be inappropriate due to temporal correlations in fMRI noise.) 0167 % For example, if there are 10 runs, each bootstrap sample consists of 10 runs 0168 % drawn with replacement from the 10 runs. 0169 % - In cases of unbalanced designs, it is possible that a bootstrap sample contains 0170 % no occurrences of a given condition; in this case, a warning is reported and 0171 % the beta weight estimated for that condition is set to zero. 0172 % 0173 % Notes on the estimation of a global HRF: 0174 % - When <hrfmodel> is 'optimize', we estimate a global HRF from the data. 0175 % This is achieved using an iterative fitting strategy: First, the HRF is fixed 0176 % to the initial seed provided by <hrfknobs>, and we estimate the amplitudes 0177 % using OLS. Then, the amplitudes are fixed (to the estimates obtained in 0178 % the previous step), and we estimate the HRF using OLS. Next, the HRF is fixed 0179 % (to the estimate obtained in the previous step), and we re-estimate the 0180 % amplitudes using OLS. This process is repeated until convergence. 0181 % - The reason for the iterative fitting strategy is that the entire model 0182 % cannot be estimated at once using linear fitting techniques (and nonlinear 0183 % techniques would be too costly). 0184 % - At the HRF-estimation steps of the fitting process, the entire dataset can 0185 % in theory be fit. However, this is undesirable for two reasons. One, 0186 % fitting the entire dataset may have exorbitant memory requirements. 0187 % Two, assuming that most voxels are unrelated to the experimental paradigm 0188 % (as is typically the case in an fMRI experiment), fitting the entire dataset 0189 % will result in a poor-quality (noisy) HRF. To resolve these issues, we use 0190 % a strategy in which we determine the best voxels in terms of R^2 at a given 0191 % amplitude-estimation step and fit only these voxels in the subsequent 0192 % HRF-estimation step. The number of voxels that are chosen is controlled 0193 % by opt.numforhrf, and the pool of chosen voxels is updated at each 0194 % amplitude-estimation step. 0195 % - In some cases, the fitted global HRF may diverge wildly from the initial 0196 % seed. This may indicate extremely low SNR and/or a problem with the coding 0197 % of the experimental design and/or a poor initial seed for the HRF. If the 0198 % R^2 between the initial seed and the fitted global HRF is less than opt.hrfthresh, 0199 % we issue a warning and simply use the initial seed as the HRF (instead of 0200 % relying on the fitted global HRF). These cases should be inspected and 0201 % troubleshooted on a case-by-case basis. (In GLMdenoisedata.m, a figure 0202 % named "HRF.png" is created --- if the initial and estimated HRF are 0203 % exactly overlapping on the figure, this indicates that the exception 0204 % case occured.) 0205 % 0206 % Additional information: 0207 % - In some circumstances (e.g. using a FIR model with insufficient data), 0208 % the design matrix may be singular and there is no unique solution. Our 0209 % strategy for these cases is as follows: If MATLAB issues a warning during 0210 % the inversion of the autocorrelation matrix (i.e. X'*X), then program 0211 % execution halts. 0212 % 0213 % History: 0214 % - 2014/07/31: return rawdesign in cache; change cubic to pchip to avoid warnings 0215 % - 2013/12/11: now, after we are done using opt.seed, we reset the random number seed 0216 % to something random (specifically, sum(100*clock)). 0217 % - 2013/11/18: add cache input/output; update documentation; new default for 0218 % maxpolydeg (it used to always be 2); add opt.hrfthresh; 0219 % add opt.suppressoutput; some speed-ups 0220 % - 2013/05/12: allow <design> to specify onset times 0221 % - 2013/05/12: update to indicate fractional values in design matrix are allowed. 0222 % - 2013/05/12 - regressors that are all zero now receive a 0 weight (instead of crashing) 0223 % - 2013/05/12 - fixed a bug regarding how the extraregressors were being handled. 0224 % previously, the extraregressors and the polynomial regressors were being regressed 0225 % out sequentially, which is improper. now, the two regressors are being fit 0226 % simultaneously, which is the correct way to do it. 0227 % - 2012/12/06: automatically convert data to single format 0228 % - 2012/12/03: *** Tag: Version 1.02 ***. Use faster OLS computation (less 0229 % error-checking; program execution will halt if design matrix is singular); 0230 % implement various speed-ups; minor bug fixes. 0231 % - 2012/11/24: 0232 % - INPUTS: add stimdur and tr; hrfknobs is optional now; add opt.hrffitmask; add opt.wantpercentbold 0233 % - OUTPUTS: add signal,noise,SNR; add hrffitvoxels; add meanvol; add inputs 0234 % - add a speed-up (design2pre) 0235 % - 2012/11/02 - Initial version. 0236 % - 2012/10/30 - Automatic division of HRF. Ensure one complete round of fitting in optimize case. 0237 % Add sanity check on HRF. 0238 0239 % Internal input: 0240 % <mode> (optional) is 0241 % 1 means that only the 'R2' output is desired (to save computation time) 0242 % 2 means that hrfmodel is 'optimize', resampling is 0, and we only care 0243 % about the hrf and hrffitvoxels outputs (to save time and memory) 0244 % Default: 0. 0245 0246 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DEAL WITH INPUTS, ETC. 0247 0248 % input 0249 if ~exist('hrfknobs','var') || isempty(hrfknobs) 0250 if isequal(hrfmodel,'fir') 0251 hrfknobs = 20; 0252 else 0253 hrfknobs = normalizemax(getcanonicalhrf(stimdur,tr)'); 0254 end 0255 end 0256 if ~exist('opt','var') || isempty(opt) 0257 opt = struct(); 0258 end 0259 if ~exist('cache','var') || isempty(cache) 0260 cache = []; 0261 end 0262 if ~exist('mode','var') || isempty(mode) 0263 mode = 0; 0264 end 0265 0266 % massage input 0267 if ~iscell(design) 0268 design = {design}; 0269 end 0270 if ~iscell(data) 0271 data = {data}; 0272 end 0273 for p=1:length(data) 0274 if ~isa(data{p},'single') 0275 data{p} = single(data{p}); 0276 end 0277 end 0278 0279 % calc 0280 numruns = length(design); 0281 if resampling == 0 0282 resamplecase = 'full'; 0283 elseif resampling >= 1 0284 resamplecase = 'boot'; 0285 else 0286 resamplecase = 'xval'; 0287 end 0288 is3d = size(data{1},4) > 1; 0289 if is3d 0290 dimdata = 3; 0291 dimtime = 4; 0292 xyzsize = sizefull(data{1},3); 0293 else 0294 dimdata = 1; 0295 dimtime = 2; 0296 xyzsize = size(data{1},1); 0297 end 0298 numvoxels = prod(xyzsize); 0299 0300 % deal with defaults 0301 if ~isfield(opt,'extraregressors') || isempty(opt.extraregressors) 0302 opt.extraregressors = cell(1,numruns); 0303 end 0304 if ~isfield(opt,'maxpolydeg') || isempty(opt.maxpolydeg) 0305 opt.maxpolydeg = zeros(1,numruns); 0306 for p=1:numruns 0307 opt.maxpolydeg(p) = round(((size(data{p},dimtime)*tr)/60)/2); 0308 end 0309 end 0310 if ~isfield(opt,'seed') || isempty(opt.seed) 0311 opt.seed = sum(100*clock); 0312 end 0313 if ~isfield(opt,'bootgroups') || isempty(opt.bootgroups) 0314 opt.bootgroups = ones(1,numruns); 0315 end 0316 if ~isfield(opt,'numforhrf') || isempty(opt.numforhrf) 0317 opt.numforhrf = 50; 0318 end 0319 if ~isfield(opt,'hrffitmask') || isempty(opt.hrffitmask) 0320 opt.hrffitmask = 1; 0321 end 0322 if ~isfield(opt,'wantpercentbold') || isempty(opt.wantpercentbold) 0323 opt.wantpercentbold = 1; 0324 end 0325 if ~isfield(opt,'hrfthresh') || isempty(opt.hrfthresh) 0326 opt.hrfthresh = 50; 0327 end 0328 if ~isfield(opt,'suppressoutput') || isempty(opt.suppressoutput) 0329 opt.suppressoutput = 0; 0330 end 0331 if isequal(hrfmodel,'assume') || isequal(hrfmodel,'optimize') 0332 hrfknobs = normalizemax(hrfknobs); 0333 end 0334 if length(opt.maxpolydeg) == 1 0335 opt.maxpolydeg = repmat(opt.maxpolydeg,[1 numruns]); 0336 end 0337 0338 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CALCULATE MEAN VOLUME 0339 0340 volcnt = cellfun(@(x) size(x,dimtime),data); 0341 meanvol = reshape(catcell(2,cellfun(@(x) squish(mean(x,dimtime),dimdata),data,'UniformOutput',0)) ... 0342 * (volcnt' / sum(volcnt)),[xyzsize 1]); 0343 0344 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% DEAL WITH NUISANCE COMPONENTS 0345 0346 % construct projection matrices for the nuisance components 0347 polymatrix = {}; 0348 exmatrix = {}; 0349 combinedmatrix = {}; 0350 for p=1:numruns 0351 0352 % this projects out polynomials 0353 pmatrix = constructpolynomialmatrix(size(data{p},dimtime),0:opt.maxpolydeg(p)); 0354 polymatrix{p} = projectionmatrix(pmatrix); 0355 0356 % this projects out the extra regressors 0357 if isempty(opt.extraregressors{p}) 0358 exmatrix{p} = 1; 0359 else 0360 exmatrix{p} = projectionmatrix(opt.extraregressors{p}); 0361 end 0362 0363 % this projects out both of them 0364 combinedmatrix{p} = projectionmatrix(cat(2,pmatrix,opt.extraregressors{p})); 0365 0366 end 0367 0368 % project out nuisance components from the data. 0369 % after this step, data will have polynomials removed, 0370 % and data2 will have both polynomials and extra regressors removed. 0371 data2 = {}; % NOTE: data and data2 are big --- be careful of MEMORY usage. 0372 for p=1:numruns 0373 data{p} = squish(data{p},dimdata)'; 0374 data2{p} = combinedmatrix{p}*data{p}; 0375 data{p} = polymatrix{p}*data{p}; 0376 end 0377 % note that data and data2 are now in flattened format (time x voxels)!! 0378 0379 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIT MODELS 0380 0381 % note: cache.rawdesign will exist for 'fir' and 'assume' but not 'optimize' and 0382 % for 'full' and 'boot' but not 'xval'. 0383 0384 switch resamplecase 0385 0386 case 'full' 0387 0388 % this is the full-fit case 0389 0390 % fit the model to the entire dataset. we obtain just one analysis result. 0391 if ~opt.suppressoutput, fprintf('fitting model...');, end 0392 results = {}; 0393 [results{1},hrffitvoxels,cache] = ... 0394 fitmodel_helper(design,data2,tr,hrfmodel,hrfknobs, ... 0395 opt,combinedmatrix,dimdata,dimtime,cache); 0396 if ~opt.suppressoutput, fprintf('done.\n');, end 0397 0398 case 'boot' 0399 0400 % this is the bootstrap case 0401 0402 % set random seed 0403 setrandstate({opt.seed}); 0404 0405 % there are two reasons to call this line. one is that in the case of 0406 % (bootstrap + optimize), we have to do a pre-call to get some cache. 0407 % another is that we may need the cache.rawdesign output. so, let's just call it. 0408 [d,d,cache] = ... 0409 fitmodel_helper(design,data2,tr,hrfmodel,hrfknobs, ... 0410 opt,combinedmatrix,dimdata,dimtime,cache); 0411 0412 % loop over bootstraps and collect up the analysis results. 0413 results = {}; 0414 if ~opt.suppressoutput, fprintf('bootstrapping model');, end 0415 for p=1:resampling 0416 statusdots(p,resampling); 0417 0418 % figure out bootstrap sample 0419 ix = []; 0420 for q=1:max(opt.bootgroups) 0421 num = sum(opt.bootgroups==q); % number in this group 0422 ix = [ix subscript(find(opt.bootgroups==q),ceil(rand(1,num)*num))]; 0423 end 0424 0425 % fit the model to the bootstrap sample 0426 if isequal(hrfmodel,'optimize') 0427 cache2 = struct('design2pre',{cache.design2pre(ix)}); 0428 else 0429 cache2 = []; 0430 end 0431 [results{p},hrffitvoxels] = fitmodel_helper(design(ix),data2(ix), ... 0432 tr,hrfmodel,hrfknobs,opt,combinedmatrix(ix),dimdata,dimtime,cache2); 0433 0434 end 0435 if ~opt.suppressoutput, fprintf('done.\n');, end 0436 0437 % randomize the random seed 0438 setrandstate(1); 0439 0440 case 'xval' 0441 0442 % this is the cross-validation case 0443 0444 % loop over cross-validation iterations. in each iteration, we record 0445 % the analysis result and also record the time-series predictions. 0446 modelfit = {}; 0447 results = {}; 0448 if ~opt.suppressoutput, fprintf('cross-validating model');, end 0449 for p=1:numruns 0450 statusdots(p,numruns); 0451 0452 % figure out resampling scheme 0453 ix = setdiff(1:numruns,p); 0454 0455 % fit the model 0456 [results{p},hrffitvoxels] = fitmodel_helper(design(ix),data2(ix), ... 0457 tr,hrfmodel,hrfknobs,opt,combinedmatrix(ix),dimdata,dimtime,[]); % NOTE: no cache 0458 0459 % compute the prediction 0460 modelfit(p) = GLMpredictresponses(results{p},{design{p}},tr,size(data2{p},1),1); % 1 because results{p} is in flattened format 0461 0462 % massage format 0463 modelfit{p} = reshape(modelfit{p},[xyzsize size(modelfit{p},2)]); 0464 0465 end 0466 if ~opt.suppressoutput, fprintf('done.\n');, end 0467 0468 end 0469 0470 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE MODEL ESTIMATES FOR OUTPUT 0471 0472 % in this special case, we do not have to perform this section, 0473 % so let's skip it to save computational time. 0474 if isequal(resamplecase,'xval') && mode==1 0475 results = struct(); 0476 0477 % otherwise, do it as usual 0478 else 0479 0480 if ~opt.suppressoutput, fprintf('preparing output...');, end 0481 switch hrfmodel 0482 0483 case 'fir' 0484 0485 results = struct('models',cat(4,results{:})); % voxels x conditions x time x resamples 0486 if size(results.models,4) == 1 0487 results.modelmd = results.models; 0488 results.modelse = zeros(size(results.models),class(results.models)); 0489 else 0490 temp = zeros([sizefull(results.models,3) 3],class(results.models)); 0491 for p=1:size(results.models,3) % ugly to avoid memory usage 0492 temp(:,:,p,:) = prctile(results.models(:,:,p,:),[16 50 84],4); 0493 end 0494 results.modelmd = temp(:,:,:,2); 0495 results.modelse = diff(temp(:,:,:,[1 3]),1,4)/2; 0496 clear temp; 0497 end 0498 0499 % massage format 0500 sz = sizefull(results.models,4); 0501 results.models = reshape(results.models,[xyzsize sz(2:4)]); 0502 results.modelmd = reshape(results.modelmd,[xyzsize sz(2:3)]); 0503 results.modelse = reshape(results.modelse,[xyzsize sz(2:3)]); 0504 0505 case {'assume' 'optimize'} 0506 0507 temp = catcell(2,cellfun(@(x) x(1),results)); 0508 results = catcell(3,cellfun(@(x) x(2),results)); % beware of MEMORY here 0509 results = struct('models',{{temp results}}); % ugly to avoid memory usage 0510 0511 % deal with {1} 0512 if size(results.models{1},2) == 1 0513 results.modelmd{1} = results.models{1}; 0514 results.modelse{1} = zeros(size(results.models{1}),class(results.models{1})); 0515 else 0516 temp = prctile(results.models{1},[16 50 84],2); 0517 results.modelmd{1} = temp(:,2); 0518 results.modelse{1} = diff(temp(:,[1 3]),1,2)/2; 0519 end 0520 0521 % deal with {2} 0522 if size(results.models{2},3) == 1 0523 results.modelmd{2} = results.models{2}; 0524 results.modelse{2} = zeros(size(results.models{2}),class(results.models{2})); 0525 else 0526 temp = zeros([sizefull(results.models{2},2) 3],class(results.models{2})); 0527 for p=1:size(results.models{2},2) % ugly to avoid memory usage 0528 temp(:,p,:) = prctile(results.models{2}(:,p,:),[16 50 84],3); 0529 end 0530 results.modelmd{2} = temp(:,:,2); 0531 results.modelse{2} = diff(temp(:,:,[1 3]),1,3)/2; 0532 clear temp; 0533 end 0534 0535 % massage format 0536 sz = sizefull(results.models{2},3); 0537 results.models{2} = reshape(results.models{2},[xyzsize sz(2:3)]); 0538 results.modelmd{2} = reshape(results.modelmd{2},[xyzsize sz(2)]); 0539 results.modelse{2} = reshape(results.modelse{2},[xyzsize sz(2)]); 0540 0541 end 0542 if ~opt.suppressoutput, fprintf('done.\n');, end 0543 0544 end 0545 0546 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMPUTE MODEL FITS (IF NECESSARY) 0547 0548 if ~(mode==2) 0549 0550 if ~opt.suppressoutput, fprintf('computing model fits...');, end 0551 switch resamplecase 0552 0553 case {'full' 'boot'} 0554 0555 % compute the time-series fit corresponding to the final model estimate 0556 modelfit = GLMpredictresponses(results.modelmd,design,tr,cellfun(@(x) size(x,1),data),dimdata); 0557 0558 case 'xval' 0559 0560 % in the cross-validation case, we have already computed the cross-validated 0561 % predictions of the model and stored them in the variable 'modelfit'. 0562 0563 end 0564 if ~opt.suppressoutput, fprintf('done.\n');, end 0565 0566 end 0567 0568 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMPUTE R^2 0569 0570 if ~(mode==2 || mode==3) 0571 0572 if ~opt.suppressoutput, fprintf('computing R^2...');, end 0573 0574 % remove polynomials from the model fits (or predictions) 0575 modelfit = cellfun(@(a,b) a*squish(b,dimdata)',polymatrix,modelfit,'UniformOutput',0); % result is time x voxels 0576 0577 % calculate overall R^2 [beware: MEMORY] 0578 results.R2 = reshape(calccodcell(modelfit,data,1)',[xyzsize 1]); % notice that we use 'data' not 'data2' 0579 0580 % calculate R^2 on a per-run basis [beware: MEMORY] 0581 results.R2run = catcell(dimdata+1,cellfun(@(x,y) reshape(calccod(x,y,1,0,0)',[xyzsize 1]),modelfit,data,'UniformOutput',0)); 0582 0583 % clear 0584 clear modelfit; % big memory usage 0585 0586 if ~opt.suppressoutput, fprintf('done.\n');, end 0587 0588 end 0589 0590 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% COMPUTE SNR 0591 0592 if ~opt.suppressoutput, fprintf('computing SNR...');, end 0593 0594 if ~(isequal(resamplecase,'xval') && mode==1) && ~(mode==2) 0595 0596 switch hrfmodel 0597 0598 case 'fir' 0599 results.signal = max(max(abs(results.modelmd),[],dimdata+1),[],dimdata+2); 0600 results.noise = mean(mean(results.modelse,dimdata+1),dimdata+2); 0601 results.SNR = results.signal ./ results.noise; 0602 0603 case {'assume' 'optimize'} 0604 results.signal = max(abs(results.modelmd{2}),[],dimdata+1); 0605 results.noise = mean(results.modelse{2},dimdata+1); 0606 results.SNR = results.signal ./ results.noise; 0607 0608 end 0609 0610 end 0611 0612 if ~opt.suppressoutput, fprintf('done.\n');, end 0613 0614 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE ADDITIONAL OUTPUTS 0615 0616 % this is a special case 0617 if isempty(hrffitvoxels) 0618 results.hrffitvoxels = []; 0619 else 0620 results.hrffitvoxels = copymatrix(zeros([xyzsize 1]),hrffitvoxels,1); 0621 end 0622 results.meanvol = meanvol; 0623 0624 % return all the inputs (except for the data) in the output. 0625 % also, include a new field 'datasize'. 0626 results.inputs.design = design; 0627 results.inputs.datasize = cellfun(@(x) size(x),data,'UniformOutput',0); 0628 results.inputs.stimdur = stimdur; 0629 results.inputs.tr = tr; 0630 results.inputs.hrfmodel = hrfmodel; 0631 results.inputs.hrfknobs = hrfknobs; 0632 results.inputs.resampling = resampling; 0633 results.inputs.opt = opt; 0634 0635 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CONVERT TO % BOLD CHANGE 0636 0637 if opt.wantpercentbold && ~(isequal(resamplecase,'xval') && mode==1) 0638 con = 1./abs(results.meanvol) * 100; 0639 switch hrfmodel 0640 case 'fir' 0641 results.models = bsxfun(@times,results.models,con); 0642 results.modelmd = bsxfun(@times,results.modelmd,con); 0643 results.modelse = bsxfun(@times,results.modelse,con); 0644 case {'assume' 'optimize'} 0645 results.models{2} = bsxfun(@times,results.models{2},con); 0646 results.modelmd{2} = bsxfun(@times,results.modelmd{2},con); 0647 results.modelse{2} = bsxfun(@times,results.modelse{2},con); 0648 end 0649 end 0650 0651 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0652 0653 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HELPER FUNCTION: 0654 0655 function [f,hrffitvoxels,cache] = ... 0656 fitmodel_helper(design,data2,tr,hrfmodel,hrfknobs,opt,combinedmatrix,dimdata,dimtime,cache) 0657 0658 % if hrfmodel is 'fir', then <f> will be voxels x conditions x time (flattened format) 0659 % if hrfmodel is 'assume' or 'optimize', then <f> will be {A B} 0660 % where A is time x 1 and B is voxels x conditions (flattened format). 0661 % <hrffitvoxels> is [] unless hrfmodel is 'optimize', in which case it will be 0662 % a column vector of voxel indices. 0663 % 0664 % note: cache.rawdesign will exist for 'fir' and 'assume' but not 'optimize'. 0665 0666 % internal constants 0667 minR2 = 99; % in 'optimize' mode, if R^2 between previous HRF and new HRF 0668 % is above this threshold (and we have at least gone through 0669 % one complete round of fitting (just so that we can change 0670 % a little from the initial seed)), then we stop fitting. 0671 0672 % init 0673 hrffitvoxels = []; 0674 0675 switch hrfmodel 0676 case 'fir' 0677 0678 % since 'fir', we can assume design is not the onset case, but check it 0679 assert(~iscell(design{1})); 0680 0681 % calc 0682 numconditions = size(design{1},2); 0683 0684 % prepare design matrix 0685 for p=1:length(design) 0686 0687 % expand original design matrix using delta basis functions. 0688 % the length of each timecourse is L. 0689 design{p} = constructstimulusmatrices(full(design{p})',0,hrfknobs,0); % time x L*conditions 0690 0691 % save a record of the raw design matrix 0692 cache.rawdesign{p} = design{p}; 0693 0694 % remove polynomials and extra regressors 0695 design{p} = combinedmatrix{p}*design{p}; % time x L*conditions 0696 0697 end 0698 0699 % fit model 0700 f = mtimescell(olsmatrix2(cat(1,design{:})),data2); % L*conditions x voxels 0701 f = permute(reshape(f,hrfknobs+1,numconditions,[]),[3 2 1]); % voxels x conditions x L 0702 0703 case 'assume' 0704 0705 % prepare design matrix 0706 for p=1:length(design) 0707 0708 % if onset-time case 0709 if iscell(design{p}) 0710 0711 % calc 0712 alltimes = linspacefixeddiff(0,tr,size(data2{p},1)); 0713 hrftimes = linspacefixeddiff(0,tr,length(hrfknobs)); 0714 0715 % loop over conditions 0716 temp = zeros(size(data2{p},1),length(design{p})); % this will be time x conditions 0717 for q=1:length(design{p}) 0718 0719 % onset times for qth condition in run p 0720 otimes = design{p}{q}; 0721 0722 % intialize 0723 yvals = 0; 0724 0725 % loop over onset times 0726 for r=1:length(otimes) 0727 0728 % interpolate to find values at the data sampling time points 0729 yvals = yvals + interp1(otimes(r) + hrftimes,hrfknobs',alltimes,'pchip',0); 0730 0731 end 0732 0733 % record 0734 temp(:,q) = yvals; 0735 0736 end 0737 0738 % save a record of the raw design matrix 0739 cache.rawdesign{p} = temp; 0740 0741 % remove polynomials and extra regressors 0742 design{p} = combinedmatrix{p}*temp; % time x conditions 0743 0744 % if regular matrix case 0745 else 0746 0747 % convolve original design matrix with HRF 0748 ntime = size(design{p},1); % number of time points 0749 design{p} = conv2(full(design{p}),hrfknobs); % convolve 0750 design{p} = design{p}(1:ntime,:); % extract desired subset 0751 0752 % save a record of the raw design matrix 0753 cache.rawdesign{p} = design{p}; 0754 0755 % remove polynomials and extra regressors 0756 design{p} = combinedmatrix{p}*design{p}; % time x conditions 0757 0758 end 0759 0760 end 0761 0762 % fit model 0763 f = mtimescell(olsmatrix2(cat(1,design{:})),data2); % conditions x voxels 0764 f = {hrfknobs f'}; 0765 0766 case 'optimize' 0767 0768 % since 'optimize', we can assume design is not the onset case, but check it 0769 assert(~iscell(design{1})); 0770 0771 % calc 0772 numinhrf = length(hrfknobs); 0773 numcond = size(design{1},2); 0774 0775 % if cache is empty, fill it 0776 if isempty(cache) 0777 0778 % precompute for speed 0779 design2pre = {}; 0780 for p=1:length(design) 0781 0782 % expand design matrix using delta functions 0783 ntime = size(design{p},1); % number of time points 0784 design2pre{p} = constructstimulusmatrices(full(design{p})',0,numinhrf-1,0); % time x L*conditions 0785 design2pre{p} = reshape(design2pre{p},[],numcond); % time*L x conditions 0786 0787 end 0788 0789 % record it 0790 cache.design2pre = design2pre; 0791 0792 % otherwise, use the cache 0793 else 0794 design2pre = cache.design2pre; 0795 end 0796 0797 % loop until convergence 0798 currenthrf = hrfknobs; % initialize 0799 cnt = 1; 0800 while 1 0801 0802 % fix the HRF, estimate the amplitudes 0803 if mod(cnt,2)==1 0804 0805 % prepare design matrix 0806 design2 = {}; 0807 for p=1:length(design) 0808 0809 % convolve original design matrix with HRF 0810 ntime = size(design{p},1); % number of time points 0811 design2{p} = conv2(full(design{p}),currenthrf); % convolve 0812 design2{p} = design2{p}(1:ntime,:); % extract desired subset 0813 0814 % remove polynomials and extra regressors 0815 design2{p} = combinedmatrix{p}*design2{p}; % time x conditions 0816 0817 end 0818 0819 % estimate the amplitudes 0820 currentbeta = mtimescell(olsmatrix2(cat(1,design2{:})),data2); % conditions x voxels 0821 0822 % calculate R^2 0823 modelfit = cellfun(@(x) x*currentbeta,design2,'UniformOutput',0); 0824 R2 = calccodcell(modelfit,data2,1)'; 0825 clear modelfit; 0826 0827 % figure out indices of good voxels 0828 if isequal(opt.hrffitmask,1) 0829 temp = R2; 0830 else 0831 temp = copymatrix(R2,~opt.hrffitmask(:),-Inf); % shove -Inf in where invalid 0832 end 0833 temp = nanreplace(temp,-Inf); 0834 [dd,ii] = sort(temp); 0835 iichosen = ii(max(1,end-opt.numforhrf+1):end); 0836 iichosen = setdiff(iichosen,iichosen(temp(iichosen)==-Inf)); 0837 hrffitvoxels = iichosen; 0838 0839 % fix the amplitudes, estimate the HRF 0840 else 0841 0842 % prepare design matrix 0843 design2 = {}; 0844 for p=1:length(design) 0845 0846 % calc 0847 ntime = size(design{p},1); % number of time points 0848 0849 % weight and sum based on the current amplitude estimates. only include the good voxels. 0850 design2{p} = design2pre{p} * currentbeta(:,hrffitvoxels); % time*L x voxels 0851 0852 % remove polynomials and extra regressors 0853 design2{p} = reshape(design2{p},ntime,[]); % time x L*voxels 0854 design2{p} = combinedmatrix{p}*design2{p}; % time x L*voxels 0855 design2{p} = permute(reshape(design2{p},ntime,numinhrf,[]),[1 3 2]); % time x voxels x L 0856 0857 end 0858 0859 % estimate the HRF 0860 previoushrf = currenthrf; 0861 datasubset = cellfun(@(x) x(:,hrffitvoxels),data2,'UniformOutput',0); 0862 prev = warning('query'); warning off; 0863 currenthrf = olsmatrix2(squish(cat(1,design2{:}),2)) * vflatten(cat(1,datasubset{:})); 0864 warning(prev); 0865 0866 % if HRF is all zeros (this can happen when the data are all zeros), get out prematurely 0867 if all(currenthrf==0) 0868 break; 0869 end 0870 0871 % check for convergence 0872 if calccod(previoushrf,currenthrf,[],0,0) >= minR2 && cnt > 2 0873 break; 0874 end 0875 0876 end 0877 0878 cnt = cnt + 1; 0879 0880 end 0881 0882 % sanity check 0883 if calccod(hrfknobs,previoushrf,[],0,0) < opt.hrfthresh 0884 warning('Global HRF estimate is far from the initial seed, probably indicating low SNR. We are just going to use the initial seed as the HRF estimate.'); 0885 [f,hrffitvoxels] = fitmodel_helper(design,data2,tr,'assume',hrfknobs,opt,combinedmatrix,dimdata,dimtime,[]); 0886 return; 0887 end 0888 0889 % normalize results 0890 mx = max(previoushrf); 0891 previoushrf = previoushrf / mx; 0892 currentbeta = currentbeta * mx; 0893 0894 % return 0895 f = {previoushrf currentbeta'}; 0896 0897 end