Dear MNE experts,
I have a question about creating a mask for
spatio_temporal_cluster_1samp_test. I'm trying to restrict the cluster
analysis on the occipital and temporal region for both left and right
hemispheres. The way I tried is as following:
1. ########### Exclude_indicies except for occipital and temporal ######
2. label_lh = mne.read_labels_from_annot('fsaverage',
parc='PALS_B12_Lobes', subjects_dir=subjects_dir,hemi='lh')
3. label_rh = mne.read_labels_from_annot('fsaverage',
parc='PALS_B12_Lobes', subjects_dir=subjects_dir,hemi='rh')
4.
5. stc_Ocp_lh = stc.in_label(label_lh[3])
6. stc_Tmp_lh = stc.in_label(label_lh[5])
7. stc_Ocp_rh = stc.in_label(label_rh[6])
8. stc_Tmp_rh = stc.in_label(label_rh[8])
9.
10. Ocp_vertices_lh = stc_Ocp_lh.vertices[0]
11. Tmp_vertices_lh = stc_Tmp_lh.vertices[0]
12. All_vertices_lh = stc.vertices[0]
13. Ocp_vertices_rh = stc_Ocp_rh.vertices[1]
14. Tmp_vertices_rh = stc_Tmp_rh.vertices[1]
15. All_vertices_rh = stc.vertices[1]
16.
17. Ocp_indices_lh = np.searchsorted(All_vertices_lh, Ocp_vertices_lh)
18. Tmp_indices_lh = np.searchsorted(All_vertices_lh, Tmp_vertices_lh)
19. Ocp_indices_rh = np.searchsorted(All_vertices_rh, Ocp_vertices_rh) +
2562
20. Tmp_indices_rh = np.searchsorted(All_vertices_rh, Tmp_vertices_rh) +
2562
21.
22. temp_all_indices = range(5124)
23. Exclude_indices_1 = np.setdiff1d(temp_all_indices, Ocp_indices_lh)
24. Exclude_indices_2 = np.setdiff1d(Exclude_indices_1, Tmp_indices_lh)
25. Exclude_indices_3 = np.setdiff1d(Exclude_indices_2, Ocp_indices_rh)
26. *Exclude_indices* = np.setdiff1d(Exclude_indices_3, Tmp_indices_rh)
27.
28. ########### Compute statistic ###########
29. print('Computing connectivity.')
30. connectivity = spatial_tris_connectivity(grade_to_tris(4))
31.
32. # Note that X needs to be a multi-dimensional array of shape
33. # samples (subjects) x time x space, so we permute dimensions
34. X = np.transpose(X, [2, 1, 0]) # X.shape was [5124, 81, 9],
vertices * times * subjects
35.
36. # Now let's actually do the clustering. This can take a long
time...
37. p_threshold = 0.05
38. t_threshold = -stats.distributions.t.ppf(p_threshold / 2.,
n_subjects - 1)
39. print('Clustering.')
40. T_obs, clusters, cluster_p_values, H0 = clu = \
41. spatio_temporal_cluster_1samp_test(X, connectivity=connectivity,
n_jobs=2,
42. spatial_exclude =
*Exclude_indices*,
43. threshold=t_threshold)
44. # Now select the clusters that are sig. at p < 0.05 (note that
this value
45. # is multiple-comparisons corrected).
46. good_cluster_inds = np.where(cluster_p_values < 0.05)[0]
47.
48. ########### Visualize the clusters ###########
49. print('Visualizing clusters.')
50. fsave_vertices = [np.arange(2562), np.arange(2562)]
51.
52. # Now let's build a convenient representation of each cluster,
where each
53. # cluster becomes a "time point" in the SourceEstimate
54. stc_all_cluster_vis = summarize_clusters_stc(clu, tstep=tstep,
55. vertices=fsave_vertices,
56. subject='fsaverage')
57.
58. # Let's actually plot the first "time point" in the
SourceEstimate, which
59. # shows all the clusters, weighted by duration
60.
61. # blue blobs are for condition A < condition B, red for A > B
62. stc_all_cluster_vis.plot(subject='fsaverage', hemi='both',
views='lateral',
63. subjects_dir=subjects_dir,
time_viewer = False,
64. time_label='Duration significant
(ms)', colormap = 'mne')
-----------End of the code---------
I thought I selected the vertices of occipital and temporal lobes
using *spatial_exclude
= Exclude_indices. * However, when I visualize the clusters, there are
still clusters show up in other regions e.g. frontal lobe, although overall
there are less clusters as expected when compared to when using
*spatial_exclude
= None.*
It's likely that the way I define Exclude_indices is wrong. It would be
greatly appreciated if anyone can help me correct this.
Thanks in advance,
Emma