Single-event GLM design matrix construction and inference

Dependencies and setup

Stats resources

• Basic statistics terminology This video gently introduces some of the key concepts that provide the foundation for statistics.
• Simple Linear Regression This video explains how a regression works using a single variable.
• Matrix Algebra Basics This video provides the background linear algebra needed for understanding the GLM.
• Multiple Linear Regression This video explains how multiple regression works using linear algebra.
• Hypothesis Testing This video covers the basics of hypothesis testing.
• Contrasts in Linear Models This video provides an overview of how to test hypotheses using contrasts in the context of the GLM.
• Intepreting Regression Parameters This video covers how to interpret the results from a regression analysis.
• Mean Centering Regressors This video covers a more subtle detail of why you might consider mean centering your continuour regression variables.

Conceptual foundation: A simple experimental design

We're going to construct a model (X), which is the predicted fMRI signal in a brain area that responds to faces, and fit it to the data. We can use the same model to estimate the amplitude of the response to faces () in any voxel or region of a single person's brain (y). Fitting the model involves solving the equation , where y is a brain time series (a vector), Xis the model or "design matrix", and is the unknown vector of regression slopes or "parameters". or "hat" means something that is estimated from the data during model fitting.

In our examples here, the regression slopes are equivalent to amplitudes for different psychological conditions. In this first, simplest example, we only have one regressor, corresponding to the "face" condition. There are only two model parameters: one for the "face" regressor of interest () and an intercept (). Xis a [T x 2] matrix, where T is the number of time points (fMRI images) collected. It has two columns, an intercept (with constant values for every time point) and a "face" regressor whose values are the predicted rise and fall of a hypothetical "face processing" region across time. is a [2 x 1] vector consisting of .

Create a simple, one-event design

figure; imagesc(X); colorbar

drawnow, snapnow

Questions to answer

Create simulated data

noise_ac = noise_arp(size(X, 1), [.7 .3]);

y = X(:, 1) + .5 * noise_ac; % simulated data

create_figure('noise and simulated data', 3, 1);

plot(noise)

title('Random Gaussian noise')

subplot(3, 1, 2)

plot(noise_ac)

title('Autocorrelated noise')

subplot(3, 1, 3)

plot(y)

title('Simulated Data')

Questions to answer

Fit the model

The beautiful thing about the GLM in its simplest form is that you can solve for the unknown using one algebraic equation:

. This projects the data (y) onto the subspace spanned by the columns of X, and provides a unique, optimal solution for , as long as the columns of X are not perfectly colinear. The fitted response is the "reconstructed" data, as close as I can get to the data using a linear combination of model predictors. We can construct the fit .

bhat = inv(X' * X) * X' * y;

yhat = X * bhat;

% Now plot the coefficients:

figure; bar(bhat')

set(gca, 'XTickLabel', {'Intercept' 'B_1 (Faces)'});

title('Regression coefficients');

% Now reconstruct the fits and plot it against the "data"

figure; hold on;

plot(y)

plot(yhat)

legend({'data' 'fit'})

xlabel('time')

Questions to answer

2. Re-run the simulation several times. What do you notice about how the beta for the "faces" regressor changes across repeats? How variable is the parameter estimate? The variance of the parameter estimate across runs is called the standard error of the estimate, .

3. If you re-ran the simulation many (e.g., 100,000) times, what do you think the average value would be? The long-run average value is called the expected value of the estimate, . Will the observed ever equal exactly? Why/why not?

4. If you increase the noise level, what happens to ? What happens to ?

5. If you decrease the noise level, what happens to ? What happens to ?

A statistical test for face activation ()

To perform statistical inference--that is, conduct a hypothesis test--we want to know whether the value we have observed for is unlikely to be due to chance. To perform the test, we propose a null hypothesis. This is that there is no true activation, , and the observed deviation from 0 is due to chance. We then evaluate how (un)likely it would be to observe the value under the null hypothesis. This value is the P-value, where P stands for probability. If P is low (e.g., P < 0.05 by convention, or often P < 0.001 in neuroimaging studies) we might reject the null hypothesis in favor of the alternative hypothesis, here that . We'll learn more about thresholding, or more principled methods of choosing P-values, later.

How do we evaluate this likelihood? The standard way is to make a few assumptions and use parametric statistics, which gives us theoretically informed values based on a test statistic, which is generally a ratio between an estimate of signal strength and error, or signal / noise, with a known probability distribution. This gives us an easy-to-calculate and usually robust estimate of the likelihood of observing (and, more generally, the data I observe) under the null hypothesis.

The most widely used test statistic is the t statistic, which is a estimate divided by its standard error, or . We already know how to estimate using regression. Now we have to estimate its standard error.

Statistical inference: A simulation-based example

If we could repeat the test many times (e.g., 10,000 or more) given the same ground truth () and noise level, we could estimate directly. It would be across j repetitions. is the variance, a measure of spread, and is calculated as , where j indexes each of J repetitions and the "bar" denotes the mean across J repetitions. Let's try it here:

% Set the noise level (formally, noise standard deviation) here:

% -------------------------------------------------------

Noise_std =10

Noise_std = 10

% Observed effect in sample

% -------------------------------------------------------

noise_ac = noise_arp(size(X, 1), [.7 .3]);

y = X(:, 1) + Noise_std * noise_ac; % simulated data

bhat_obs = inv(X' * X) * X' * y; % parameter estimates

bhat_obs = bhat_obs(1);

% Null hypothesis simulation - Monte Carlo

% -------------------------------------------------------

bhat_null = zeros(1000, 1); % initialize a vector to store 1,000 repeats, each on a new sample

for j = 1:1000

noise_ac_j = noise_arp(size(X, 1), [.7 .3]); % new noise sample

y = 0 * X(:, 1) + Noise_std * noise_ac_j; % simulated data - no true effect

bhat = inv(X' * X) * X' * y; % parameter estimates

bhat_null(j, 1) = bhat(1);

end

create_figure('null hypothesis sim');

histogram(bhat_null, 50)

h = plot_vertical_line(bhat_obs);

set(h, 'LineWidth', 3, 'Color', [1 .5 .2])

Questions to answer

Statistical inference: A parametric t-test

It turns out that we can estimate from a single sample, if we know two quantities: The sample error variance, which is related to how much error variability there is among your individual obervations, and the error degrees of freedom, which is how many independent observations there are in your sample. We can estimate the sample variance from the residuals (e), the errors that I cannot explain with my model: . Because of the central limit theorem, is the variance of e divided by the square root of the degrees of freedom, or .

Now for the degrees of freedom. If all observations are independent random draws from a population, this is generally the number of samples minus the number of parameters we have estimated. Here, . Note that in an fMRI time series, the observations are not independent, because they are autocorrelated -- and we have fewer degrees of freedom than we think. Estimating them requires specialized time series statistics that we'll discuss later, so for now we'll just use the naive estimate of df. Thus:

% Set the noise level (formally, noise standard deviation) here:

% -------------------------------------------------------

Noise_std =10

Noise_std = 10

% Observed effect in sample

% -------------------------------------------------------

noise_ac = noise_arp(size(X, 1), [.7 .3]);

y = X(:, 1) + Noise_std * noise_ac; % simulated data

bhat_obs = inv(X' * X) * X' * y; % parameter estimates

e = y - X * bhat_obs;

% Parametric t-test

% -------------------------------------------------------

df = length(y) - size(X, 2); % Assumed degrees of freedom. Correct if no autocorrelation

t = bhat_obs ./ sqrt(var(e) ./ df);

P = 2 .* min(tcdf(t, df), 1-tcdf(t, df));

P = P(1); % we are interested in "face" regressor (param 1), not the intercept

disp(P)

0.8149

if P < 0.05, disp('Significant!'), else, disp('Not significant :('), end

Not significant :(

% Visualize the t-distribution and observed t-value

% -------------------------------------------------------

tprob = tpdf([-3.5:.01:3.5], df);

figure; hold on;

plot([-3.5:.01:3.5], tprob);

h = plot_vertical_line(t(1));

set(h, 'LineWidth', 3, 'Color', [1 .5 .2])

title('t distribution and observed t value')

Questions to answer

Statistics in Matlab

Alternatives to parametric tests