[Homer-users] GLM analysis in NIRS: Partial response to baseline correction question

[Huppert, Ted]

I started to respond to Jared’s question on the HOMER-list, but got a bit carried away with my response.  Sorry to be long-winded on this.  Below I go through the math and an example comparing how I was suggesting Jared to analyze his data in comparison to the methods used in SPM and previous fMRI papers.

---------------------------
Let’s start with the standard time-series analysis model
 [cid:B05C8B63-AAEB-49B2-B4CC-349453149310]
 where Y is the data (vector) time course of size <number-time-points by 1> .  X is the so-called design matrix.  X is a matrix which is size <number time-points by  number-time-points in the HRF>.  Finally, beta is the vector of unknowns and models the “HRF” response.  Let’s use a specific example of 300s of data collected at 1Hz and we want to look at a 20s window for the HRF.  In this example, X is then size <300 x 20>.   The first column of X has a value of 1 for all the time points when the stimulus was presented and otherwise is 0.  The second column of X is the same as the first column but shifted by one row, etc.    The size of beta in this example is <20 by 1> and if I estimate beta and plot it, it should look like the “HRF” shape.    This is what we would call the FIR (finite impulse response) or deconvolution model.  It can be shown that in the limit that your stimuli events are far enough apart, this is the same as a block average and when the stimuli are closer together, this is a deconvolution.    This is the model used in HOMER-1 and mostly what people are familiar with in HOMER-2.

To solve this model, you would use the least-squares solution
 [cid:FA8F9A63-8AED-4753-AAE0-3790F2D232B9]

The multi-variate error in the estimate is given by
 [cid:8625AC90-8073-4F46-B272-870655A73EFD]

where sigma^2 is the variance of the model residual.  X^T means the transpose of X.

To do a statistical test (t-test) you would create a contrast vector (let’s call it “c”).  Let’s say we want to compute the contrast from 4-12s of our 20sec HRF window.  Thus, c is a vector of size <1x 20>  with a value of 1/8 at positions 4-12 and zeros elsewhere.    This is the same as you do in SPM.  Its 1/8 if you want the mean over the 8 points, otherwise you could use “1” for the sum of the 8 points.
 [cid:3C42A669-7D2E-4BAC-99F4-DBA5F747FB5E]

In this example, the T-value above is the statistical test on if the window 4-12s is statistically different from zero.  Likewise, if you wanted to test if the window (e.g) 4-12s was different from the window -5-0s, then you would create the design matrix by shifting all the points 5sec earlier and now you have 25 columns and 25 terms in beta (e.g. the first 5 points are the pre-time and the 6-25th point is the same as the HRF).  The “c” contrast is now a “-1/5” for the first 5 points and a “1/8” for points 9-17 (4-12sec) .    The T-contrast is now computed the same as before and represents the mean of the window 4-12s minus the mean of the baseline window from -5-0s.

Likewise if you have two conditions, you would have (e.g.) 25 x 2 =50 points in beta and the contrast is defined the same sort of way.

So taking Jared’s example:

There are two conditions, so the design matrix is (e.g.) <300x 50> and there are 50betas.  The first 25betas are the EXP condition and the second 25 are the NEUT condition.

Option 1.  Take the average activation during 10s of EXP across all presentations.  In this example, this is to set the contrast vector
“C” to unity (or 1/20) for the first block (the EXP condition) and zeros else.

Option 2 & 3.  Set C  to 1 for the first block and -1 for the last block.  Mathematically, since averaging commutes it doesn’t matter if you subtract the local NEUT response from each trial and then average or subtract the average NEUT response.  * At least under the assumption of random noise.  Drifts in baseline and non-stationarities will effect this however, but lets only deal with the theoretical version here.  In the presence of uncorrected drifts and non-stationarities, all the statistics are wrong anyway (see Barker et al 2013; http://www.ncbi.nlm.nih.gov/pubmed/24009999)

How does this relate to the canonical fMRI-GLM model?

Let’s now take a look at how SPM and the fMRI people do it.  In fMRI, its more common to assume a specific shape to the HRF.   Let’s call that shape GAM.  Sticking to the same example with a -5 to20sec response window, GAM is a vector of size <25x1> and if you plot it it looks like the idealized HRF response.

 With this canonical response, our model becomes
[cid:13246808-5023-4666-BE9F-FAE838A4EF4F]
where  now there is one model unknown (Beta-tilda) and you can see that Beta = GAM*Beta-tilda.   This is a new reduced basis set for the model.  We now have only one unknown when we used to have 25.  Thus, we now only have one statistical test to preform and we avoid multiple comparison corrections.

If you propagate this all the way through, you can see that
 [cid:E69EC2DF-9294-4855-83B4-E08477D9FB83]
 [cid:439C4A0C-5D85-423B-BC00-505304859A05]

Likewise the T-stat contrast is still defined as above.

A little algebra will show the T-stat obtained by assuming a canonical HRF (Gam) is the same as using a deconvolution model and setting the contrast vector (c) to be
C = GAM’ *  inv(GAM’*X’*X*GAM)*(X’*X)

Thus, we can see that there is a relationship between the canonical GLM model where the hypothesis of a response shape is provided upfront and the deconvolution model where the shape of the response is estimated but then a weighted average (given by the contrast vector) is used to compute the statistical test.

More specifically, in the limit that the stimulus events are not overlapping (e.g. a block average can be used), then solving for the deconvolution and then using C = GAM /|GAM^2| would produce the identical result as assuming the GAM canonical function in a GLM model (this is true no matter what the shape of the GAM function).    When you have overlapping events, the contrast vector is slightly corrected by the off-diagonal elements in X’*X so the contrast vector looks like a slightly filtered version of the canonical response depending on how well you jittered your stimulus timing.   In a well-designed event related experiment, the relationship of C = GAM /|GAM^2| will again be true as the off-diaginals in X’*X disappear.

In other words, (for a block average model or a well-designed event related design at least and in the limit of white noise) the canonical GLM model is the same as solving the deconvolution or block averaging model and then using a contrast vector that matches the canonical model to calculate a weighted average.

One important note however is that because of this relationship between the canonical model and computing your statistics from an equivalently weighted average, if your canonical response is wrong or equivalently the window you use for contrast is wrong, you introduce type-II errors.  It makes it harder to reject the null hypothesis and you have more have more false negatives (e.g. You fail to report areas of brain activity).    Its a trade off.  If you try only one window, then you might not see all the activity if you pick the wrong window or canonical response.  But if you try multiple contrast windows or develop your statistical test after looking at the data, now you introduce type-I (false discovery) errors and run the risk of reporting false results and need to correct for multiple comparisons.

So, hopefully this clarifies a bit about the math in HOMER (or for NIRS in general) and gives so guidance about how you can compare your NIRS analysis to previous fMRI designs.   So, my suggestion to Jared would be to not reinvent the wheel and model your analysis on previous results from fMRI.


----------------
Theodore Huppert, PhD
Associate Professor
Departments of Radiology and Bioengineering
Center for Neural Basis of Cognition
University of Pittsburgh
Email: huppertt at upmc.edu<mailto:huppertt at upmc.edu>
Phone: (412) 647-8459



"Insanity: doing the same thing over and over and expecting different results"-  Einstein
-------------- next part --------------
A non-text attachment was scrubbed...
Name: B05C8B63-AAEB-49B2-B4CC-349453149310.png
Type: image/png
Size: 1478 bytes
Desc: B05C8B63-AAEB-49B2-B4CC-349453149310.png
Url : http://mail.nmr.mgh.harvard.edu/pipermail/homer-users/attachments/20150723/263a8cb4/attachment.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: FA8F9A63-8AED-4753-AAE0-3790F2D232B9.png
Type: image/png
Size: 2075 bytes
Desc: FA8F9A63-8AED-4753-AAE0-3790F2D232B9.png
Url : http://mail.nmr.mgh.harvard.edu/pipermail/homer-users/attachments/20150723/263a8cb4/attachment-0001.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 8625AC90-8073-4F46-B272-870655A73EFD.png
Type: image/png
Size: 2420 bytes
Desc: 8625AC90-8073-4F46-B272-870655A73EFD.png
Url : http://mail.nmr.mgh.harvard.edu/pipermail/homer-users/attachments/20150723/263a8cb4/attachment-0002.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 3C42A669-7D2E-4BAC-99F4-DBA5F747FB5E.png
Type: image/png
Size: 2353 bytes
Desc: 3C42A669-7D2E-4BAC-99F4-DBA5F747FB5E.png
Url : http://mail.nmr.mgh.harvard.edu/pipermail/homer-users/attachments/20150723/263a8cb4/attachment-0003.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 13246808-5023-4666-BE9F-FAE838A4EF4F.png
Type: image/png
Size: 1826 bytes
Desc: 13246808-5023-4666-BE9F-FAE838A4EF4F.png
Url : http://mail.nmr.mgh.harvard.edu/pipermail/homer-users/attachments/20150723/263a8cb4/attachment-0004.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: E69EC2DF-9294-4855-83B4-E08477D9FB83.png
Type: image/png
Size: 2894 bytes
Desc: E69EC2DF-9294-4855-83B4-E08477D9FB83.png
Url : http://mail.nmr.mgh.harvard.edu/pipermail/homer-users/attachments/20150723/263a8cb4/attachment-0005.png 
-------------- next part --------------
A non-text attachment was scrubbed...
Name: 439C4A0C-5D85-423B-BC00-505304859A05.png
Type: image/png
Size: 2599 bytes
Desc: 439C4A0C-5D85-423B-BC00-505304859A05.png
Url : http://mail.nmr.mgh.harvard.edu/pipermail/homer-users/attachments/20150723/263a8cb4/attachment-0006.png