Spatio-temporal custer permutation with linear regression

💬 Support & Discussions
,
1

Hi all,

In my MEG project, I want to run a linear regression across sensors and time and combine it with cluster-based permutation to control for multiple comparisons. Based on a few tutorials I managed to write some code that gives me reasonably-looking output, but I am not sure whether I did it correctly, that’s why I am here.

For some context, I have a within-subject design with 58 subjects that each did 4 conditions. The four conditions are on a continuous scale from zero to one, so I think linear regression would be best suited to analyse it even though an Anova with 1 Factor and 4 levels would also work.

To run the regression, I created a dataframe with the shape (4, 58, 150, 102) (space last), and defined a custom stat_fun like so:

def linreg_cluster(*args):
    labels = ["Intercept", "Coherence"]
    n_cond = len(args) # 4
    n_obs = args[0].shape[0] # 58
    dm = np.ones((n_cond * n_obs, 2))
    dm[:, 1] = np.repeat([0.03, 0.06, 0.12, 0.24], n_obs) # levels of the factor coherence
    data = np.concatenate(args)
    # use private function, because public API only works on Epochs or Sourceestimates
    beta, stderr, t_val, p_val, mlog10_p_val = mne.stats.regression._fit_lm(data, dm, labels)
    return t_val['Coherence']

Then I computed the cluster-forming t-statistics threshold like this:

pval = 0.001 
df = 58 - 4 - 1  # degrees of freedom for the test
thresh = scipy.stats.t.ppf(1 - pval / 2, df)  # two-tailed, t distribution

and finally ran the test:

    cluster_stats = mne.stats.spatio_temporal_cluster_test(
        data,
        n_permutations=500,
        threshold=thresh,
        stat_fun=linreg_cluster,
        buffer_size=None,
        adjacency=adjacency)

Does this approach look right to you? I also tried the same procedure but then with a one-way F-test instead of a linear regression as statistical function. This was a bit a easier and I am reasonably confident that the results check out. But while the linear regression results are somewhat similar, they look much worse than I expected (i.e. the clusters have essentially no temporal extent, see below). So, I was wondering whether I did do anything wrong.

Here an example of one cluster (ignore the F values in the labels, I forgot to change them to ts:

mag_1_reg_cluster
mag_1_reg_cluster1000×300 96.7 KB

Thanks already!
Eduard

2

Hello Eduard,

I would suggest to first check the linear regression output. I assume you run the regression on evoked data? I think regression is best suited for single trial analysis, so instead of running on evokeds I would run a regression model over epochs for each participant and then use the resulting evoked beta coefficients for further group level statistics. Maybe this tutorial helps?

Best,

Carina

1 Like
3

Hi Carina,
Thanks for your reply!

I would suggest to first check the linear regression output.

Could you explain for what I should check the output?

I would run a regression model over epochs for each participant and then use the resulting evoked beta coefficients for further group level statistics

Yeah, that sounds reasonable, I wasn’t sure though how I would run the 2nd level statistic then. Suppose I have the single-subject regression betas, would I then simply construct a matrix (n_subs, n_conds, n_time, n_sensors), and permute over subjects? Or do I need to specify some stat_fun for that level as well?

Thanks,
Eduard