second level RSA analysis

[Lukas Van Oudenhove]

Dear Cosmo team,

First of all thanks a ton for a great toolbox!

With the help of the documentation and earlier posts here, I managed to run first level whole-brain searchlight rsa without problems (at least so it seems) on 30 subjects, in a similar way than the example in the paper/docs.

The input and output cosmo dataset structures are attached in a .mat file FYI.

Then, I tried to run second-level analysis (simple one-sample t-test in every voxel) using cosmo_montecarlo_cluster stats (see code below), but I get the below warning, and all z-values in the output dataset z_ds.group.samples are zero, which seems very unlikely to me?

Could you please help clarifying the warning, and give advice on what may be causing the weird result with all zeros?

Thanks very much in advance!

Best wishes,

Lukas

Warning

Warning: A value of class "double" was indexed with no subscripts specified. Currently the result of this operation is the indexed
value itself, but in a future release, it will be an error.
> In cosmo_parallel_get_nproc_available>matlab_get_max_nproc_available_ge2013b (line 192)
In cosmo_parallel_get_nproc_available (line 60)
In cosmo_montecarlo_cluster_stat (line 261)
Warning: A value of class "double" was indexed with no subscripts specified. Currently the result of this operation is the indexed
value itself, but in a future release, it will be an error.
> In cosmo_parallel_get_nproc_available>matlab_get_max_nproc_available_ge2013b (line 192)
In cosmo_parallel_get_nproc_available (line 60)
In cosmo_parcellfun (line 37)
In cosmo_montecarlo_cluster_stat (line 344) 

Code

%% DO GROUP LEVEL STATS
% -------------------------------------------------------------------------

% CREATE STACKED COSMO DATASET
% ----------------------------
ds_group = cosmo_stack(cosmo_rsa_datasets);
ds_group.sa.chunks = (1:size(cosmo_rsa_datasets,2))';
ds_group.sa.targets = ones(size(cosmo_rsa_datasets,2),1);

% CREATE NEIGHBORHOOD
% -------------------
nh_group = cosmo_cluster_neighborhood(ds_group);

% DEFINE OPTIONS
% --------------
opt_group = struct();
opt_group.niter = 5000;
opt_group.h0_mean = 0;
opt_group.nproc = 12;

% RUN PERMUTATION TEST
% --------------------
z_ds_group = cosmo_montecarlo_cluster_stat(ds_group,nh_group,opt_group);

[Nick Oosterhof]

Greetings,

> On 26 Jun 2023, at 12:19, Lukas Van Oudenhove <lukas.vano...@gmail.com> wrote:
>
> First of all thanks a ton for a great toolbox!

Happy to hear that :-)

>
> With the help of the documentation and earlier posts here, I managed to run first level whole-brain searchlight rsa without problems (at least so it seems) on 30 subjects, in a similar way than the example in the paper/docs.
>
> The input and output cosmo dataset structures are attached in a .mat file FYI.
>
> Then, I tried to run second-level analysis (simple one-sample t-test in every voxel) using cosmo_montecarlo_cluster stats (see code below), but I get the below warning, and all z-values in the output dataset z_ds.group.samples are zero, which seems very unlikely to me?
>
> Could you please help clarifying the warning, and give advice on what may be causing the weird result with all zeros?

I wonder if it may have to do with the parallel threads / toolbox.
- Are you using a very recent version of Matlab?
- What is the output of cosmo_wtf?
- What happens if you set nproc=1?

[Lukas Van Oudenhove]

Thanks a lot!

I am using Matlab R2021a on a server running on Ubuntu 20.04.

Below is the output of cosmo_wtf (bar the Matlab path)

Running with nproc=1 does not make any difference in results (all zeroes), but does not produce the warning.

Maybe relevant to the problem, I reran the same analysis (with nproc = 12) on larger matrices including the individual runs rather than averaged over runs, in which case I do get some (negative) z-values.

There is indeed quite some intra-individual variability in behavioral ratings, which may be lost when averaging over runs, leading to true null results in the averaged analysis? 

Just thinking out aloud here - would love to hear your thoughts!

Cheers,

Lukas

Op maandag 26 juni 2023 om 22:21:55 UTC+2 schreef n.n.oosterhof:

[Nick Oosterhof]

Greetings,

> On 27 Jun 2023, at 12:50, Lukas Van Oudenhove <lukas.vano...@gmail.com> wrote:
>
> I am using Matlab R2021a on a server running on Ubuntu 20.04.
>
> Below is the output of cosmo_wtf (bar the Matlab path)

Actually I did not see it, but maybe it’s not so relevant for solving the issue.

>
> Running with nproc=1 does not make any difference in results (all zeroes), but does not produce the warning.

That’s good to know, thanks.

>
> Maybe relevant to the problem, I reran the same analysis (with nproc = 12) on larger matrices including the individual runs rather than averaged over runs, in which case I do get some (negative) z-values.

Is this a strong effect (let’s say z < -1.96)?

> There is indeed quite some intra-individual variability in behavioral ratings, which may be lost when averaging over runs, leading to true null results in the averaged analysis?

That definitely sounds like a possibility.

>
> Just thinking out aloud here - would love to hear your thoughts!

It might be the case that the effects (signal) is just not very strong. Or that is is highly variable across participants.

If you compute uncorrected group maps (using cosmo_stat, with z score output), what do you see? The presence of only weak / small clusters would be consistent with finding few or no features with a non-zero z-score.

[Lukas Van Oudenhove]

Hi Nick,

Many thanks!

I agree with your interpretation that the signal is not very strong due to high variability across participants, which is interesting on its own.

There are some negative z-values surviving the < -1.96 corrected threshold, and i explored the uncorrected group maps using cosmo_stat (thanks for your suggestion). 

I have a question related to this though: does this differ from cosmo_monte_carlo_cluster_stats in two ways being 1) non-parametric and 2) TFCE corrected? 

In the output of the latter, z-values range from -2.23 to 0 (no single positive value, but spanning the entire negative range), while in the output of the former, the range of z-values is -4.8587 to 2.5530 (which fits with stronger negative than positive effects).

Finally, I read in the helpful documentation of cosmo_monte_carlo_cluster_stats that inputting null values has more power, hence I would like to implement this (using cosmo_randomize_targets as suggested).

I see that I need a 1xP cell array (with P being number of subjects, right?) with approx 100 null datasets per subject.

My question is: at which level do I generate these null datasets? For each subject separately during my first-level loop over subjects running cosmo_searchlight, I presume? I could then set up a loop within each subject invoking cosmo_randomize_targets 100 times, assign the output to ds.targets, and run cosmo_searchlight on it 100 times)? Should I then stack those 100 datasets and put that stacked ds in the cell array for that subject? This will then result in a null dataset for each subject which has the same number of columns in ds.samples as the true dataset, but a different number of rows (30 in case of the true dataset, as I have 30 subjects, but 100 in case of each null dataset).

Any help on how to optimally implement this would be highly appreciated!

Best wishes,

Lukas

Op zondag 2 juli 2023 om 11:42:22 UTC+2 schreef n.n.oosterhof:

[Nick Oosterhof]

> On 3 Jul 2023, at 11:36, Lukas Van Oudenhove <lukas.vano...@gmail.com> wrote:
>
> I agree with your interpretation that the signal is not very strong due to high variability across participants, which is interesting on its own.
>
> There are some negative z-values surviving the < -1.96 corrected threshold, and i explored the uncorrected group maps using cosmo_stat (thanks for your suggestion).
> I have a question related to this though: does this differ from cosmo_monte_carlo_cluster_stats in two ways being 1) non-parametric and 2) TFCE corrected?
> In the output of the latter, z-values range from -2.23 to 0 (no single positive value, but spanning the entire negative range), while in the output of the former, the range of z-values is -4.8587 to 2.5530 (which fits with stronger negative than positive effects).

Indeed, when using cosmo_stats the z-scores are not corrected for multiple comparisons, so the range of z-scores is expected to be wider.

>
> Finally, I read in the helpful documentation of cosmo_monte_carlo_cluster_stats that inputting null values has more power, hence I would like to implement this (using cosmo_randomize_targets as suggested).
> I see that I need a 1xP cell array (with P being number of subjects, right?) with approx 100 null datasets per subject.
> My question is: at which level do I generate these null datasets? For each subject separately during my first-level loop over subjects running cosmo_searchlight, I presume?

Yes indeed.

> I could then set up a loop within each subject invoking cosmo_randomize_targets 100 times, assign the output to ds.targets,

I think you mean ds.sa.targets.

> and run cosmo_searchlight on it 100 times)?

Yes, 100 times for each participant. I would advise to store these results on disk. With 30 participants, one way to do so is to store 3000 (=30 * 100) files with the output from cosmo_searchlight.

> Should I then stack those 100 datasets and put that stacked ds in the cell array for that subject? This will then result in a null dataset for each subject which has the same number of columns in ds.samples as the true dataset, but a different number of rows (30 in case of the true dataset, as I have 30 subjects, but 100 in case of each null dataset).

No, the idea is to have 100 null datasets, each with with 30 rows in .samples. Also: cosmo_montecarlo_cluster_stat should complain (raise an error) if an incorrect shape is used.

[Lukas Van Oudenhove]

Thanks a ton Nick!

If I understand it correctly, in the 1xP cell array, P is 100 rather than the number of subjects, and each cell contains a dataset with 30 rows in ds.samples, one row with null data for every subject?

Should be pretty straightforward to set it up!

My code will be on Github soon so feel free to share or refer to it once it is done?!

Cheers,

Lukas

Op dinsdag 4 juli 2023 om 19:18:22 UTC+2 schreef n.n.oosterhof:

[Lukas Van Oudenhove]

Hi Nick,

Quick update here: I set up things as discussed before, but contrary to running the analysis without null data, all my values in .samples are zero now - any idea why that may be?

Thanks a lot again!

Best wishes,

Lukas

Op woensdag 5 juli 2023 om 09:36:46 UTC+2 schreef Lukas Van Oudenhove:

[Nick Oosterhof]

Hi Lukas,

> On 6 Jul 2023, at 16:45, Lukas Van Oudenhove <lukas.vano...@gmail.com> wrote:
>
> Quick update here: I set up things as discussed before, but contrary to running the analysis without null data, all my values in .samples are zero now - any idea why that may be?

Sorry to hear that. The only reason I can think of is that the signal in the data is quite weak, unfortunately.

Best,
Nick

[Lukas Van Oudenhove]

Hi Nick,

No worries!

Could you please confirm that my method for the null data (described above) is correct?

Thanks again!

Best wishes,

Lukas

Op vrijdag 14 juli 2023 om 19:52:11 UTC+2 schreef n.n.oosterhof:

[Lukas Van Oudenhove]

Hi Nick,

Also, I was wondering about the interpretation of the negative sign of the correlations I find between the behavioral and neural matrices at the group level.

I put in a behavioral dissimilarity matrix (Euclidean distance as default in cosmo_pdist) as input to cosmo_target_dsm_corr_measure (with the default metric which is correlation), which, if I understand things correctly, will use correlations to calculate the neural (dis)similarity matrices. 

My questions are 

1) are the resulting neural matrices similarity matrices, or does the 'correlation' option, as cosmo_pdist does, calculate 1 - corr hence yielding a dissimilarity matrix (in the former case the interpretation of the negative sign would be 'more neural similarity ~ more behavioral similarity)? 

2) do I need the set the center_data option to true when using the default correlation metric?

3) would it be better to change the metric to Euclidean distance to be consistent with the metric used for the behavioral dissimilarity matrix, and in this case, do I also need to set the center_data option to true?

Thanks a lot again!

Best wishes,

Lukas  

Op maandag 17 juli 2023 om 09:18:35 UTC+2 schreef Lukas Van Oudenhove:

[Nick Oosterhof]

Hi Lukas,

> On 19 Jul 2023, at 19:06, Lukas Van Oudenhove <lukas.vano...@gmail.com> wrote:
>
> Also, I was wondering about the interpretation of the negative sign of the correlations I find between the behavioral and neural matrices at the group level.
>
> I put in a behavioral dissimilarity matrix (Euclidean distance as default in cosmo_pdist) as input to cosmo_target_dsm_corr_measure (with the default metric which is correlation), which, if I understand things correctly, will use correlations to calculate the neural (dis)similarity matrices.
>
> My questions are
> 1) are the resulting neural matrices similarity matrices, or does the 'correlation' option, as cosmo_pdist does, calculate 1 - corr hence yielding a dissimilarity matrix (in the former case the interpretation of the negative sign would be 'more neural similarity ~ more behavioral similarity)?

A positive sign means more similarity. A negative sign means that pairs of conditions with higher neural similarity tend to have lower behavioural similarity and vice versa.

> 2) do I need the set the center_data option to true when using the default correlation metric?

I think the current opinion in the literature is that it is recommended. The main reason it is not the default is to avoid introducing breaking changes, since the first versions of the function did not do centering.

> 3) would it be better to change the metric to Euclidean distance to be consistent with the metric used for the behavioral dissimilarity matrix, and in this case, do I also need to set the center_data option to true?

I’m sorry, I don’t have a good intuition about which metric most be most suitable (and under which conditions). Have you tried looking into the literature?

Nevertheless, I hope this helps.

Best,
Nick