.. only:: html
.. note::
:class: sphx-glr-download-link-note
Click :ref:`here ` to download the full example code or to run this example in your browser via Binder
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_auto_examples_03_connectivity_plot_multi_subject_connectome.py:
Group Sparse inverse covariance for multi-subject connectome
=============================================================
This example shows how to estimate a connectome on a group of subjects
using the group sparse inverse covariance estimate.
.. code-block:: default
import numpy as np
from nilearn import plotting
n_subjects = 4 # subjects to consider for group-sparse covariance (max: 40)
def plot_matrices(cov, prec, title, labels):
"""Plot covariance and precision matrices, for a given processing. """
prec = prec.copy() # avoid side effects
# Put zeros on the diagonal, for graph clarity.
size = prec.shape[0]
prec[list(range(size)), list(range(size))] = 0
span = max(abs(prec.min()), abs(prec.max()))
# Display covariance matrix
plotting.plot_matrix(cov, cmap=plotting.cm.bwr,
vmin=-1, vmax=1, title="%s / covariance" % title,
labels=labels)
# Display precision matrix
plotting.plot_matrix(prec, cmap=plotting.cm.bwr,
vmin=-span, vmax=span, title="%s / precision" % title,
labels=labels)
Fetching datasets
------------------
.. code-block:: default
from nilearn import datasets
msdl_atlas_dataset = datasets.fetch_atlas_msdl()
rest_dataset = datasets.fetch_development_fmri(n_subjects=n_subjects)
# print basic information on the dataset
print('First subject functional nifti image (4D) is at: %s' %
rest_dataset.func[0]) # 4D data
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
/usr/lib/python3/dist-packages/numpy/lib/npyio.py:2358: VisibleDeprecationWarning: Reading unicode strings without specifying the encoding argument is deprecated. Set the encoding, use None for the system default.
output = genfromtxt(fname, **kwargs)
First subject functional nifti image (4D) is at: /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar123_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
Extracting region signals
--------------------------
.. code-block:: default
from nilearn import image
from nilearn import input_data
# A "memory" to avoid recomputation
from joblib import Memory
mem = Memory('nilearn_cache')
masker = input_data.NiftiMapsMasker(
msdl_atlas_dataset.maps, resampling_target="maps", detrend=True,
low_pass=None, high_pass=0.01, t_r=2, standardize=True,
memory='nilearn_cache', memory_level=1, verbose=2)
masker.fit()
subject_time_series = []
func_filenames = rest_dataset.func
confound_filenames = rest_dataset.confounds
for func_filename, confound_filename in zip(func_filenames,
confound_filenames):
print("Processing file %s" % func_filename)
# Computing some confounds
hv_confounds = mem.cache(image.high_variance_confounds)(
func_filename)
region_ts = masker.transform(func_filename,
confounds=[hv_confounds, confound_filename])
subject_time_series.append(region_ts)
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
[NiftiMapsMasker.fit] loading regions from /home/varoquau/nilearn_data/msdl_atlas/MSDL_rois/msdl_rois.nii
Processing file /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar123_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
________________________________________________________________________________
[Memory] Calling nilearn.image.image.high_variance_confounds...
high_variance_confounds('/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar123_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz')
__________________________________________high_variance_confounds - 0.9s, 0.0min
________________________________________________________________________________
[Memory] Calling nilearn.input_data.base_masker.filter_and_extract...
filter_and_extract('/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar123_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz',
,
{ 'allow_overlap': True,
'detrend': True,
'dtype': None,
'high_pass': 0.01,
'low_pass': None,
'maps_img': '/home/varoquau/nilearn_data/msdl_atlas/MSDL_rois/msdl_rois.nii',
'mask_img': None,
'smoothing_fwhm': None,
'standardize': True,
't_r': 2,
'target_affine': array([[ 4., 0., 0., -78.],
[ 0., 4., 0., -111.],
[ 0., 0., 4., -51.],
[ 0., 0., 0., 1.]]),
'target_shape': (40, 48, 35)}, confounds=[ array([[-0.174325, ..., -0.048779],
...,
[-0.044073, ..., 0.155444]]),
'/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar123_task-pixar_desc-reducedConfounds_regressors.tsv'], dtype=None, memory=Memory(location=nilearn_cache/joblib), memory_level=1, verbose=2)
[NiftiMapsMasker.transform_single_imgs] Loading data from /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar123_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
[NiftiMapsMasker.transform_single_imgs] Resampling images
[NiftiMapsMasker.transform_single_imgs] Extracting region signals
[NiftiMapsMasker.transform_single_imgs] Cleaning extracted signals
_______________________________________________filter_and_extract - 6.3s, 0.1min
Processing file /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar001_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
________________________________________________________________________________
[Memory] Calling nilearn.image.image.high_variance_confounds...
high_variance_confounds('/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar001_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz')
__________________________________________high_variance_confounds - 0.9s, 0.0min
________________________________________________________________________________
[Memory] Calling nilearn.input_data.base_masker.filter_and_extract...
filter_and_extract('/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar001_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz',
,
{ 'allow_overlap': True,
'detrend': True,
'dtype': None,
'high_pass': 0.01,
'low_pass': None,
'maps_img': '/home/varoquau/nilearn_data/msdl_atlas/MSDL_rois/msdl_rois.nii',
'mask_img': None,
'smoothing_fwhm': None,
'standardize': True,
't_r': 2,
'target_affine': array([[ 4., 0., 0., -78.],
[ 0., 4., 0., -111.],
[ 0., 0., 4., -51.],
[ 0., 0., 0., 1.]]),
'target_shape': (40, 48, 35)}, confounds=[ array([[-0.151677, ..., -0.057023],
...,
[-0.206928, ..., 0.102714]]),
'/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar001_task-pixar_desc-reducedConfounds_regressors.tsv'], dtype=None, memory=Memory(location=nilearn_cache/joblib), memory_level=1, verbose=2)
[NiftiMapsMasker.transform_single_imgs] Loading data from /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar001_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
[NiftiMapsMasker.transform_single_imgs] Resampling images
[NiftiMapsMasker.transform_single_imgs] Extracting region signals
[NiftiMapsMasker.transform_single_imgs] Cleaning extracted signals
_______________________________________________filter_and_extract - 6.3s, 0.1min
Processing file /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar002_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
________________________________________________________________________________
[Memory] Calling nilearn.image.image.high_variance_confounds...
high_variance_confounds('/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar002_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz')
__________________________________________high_variance_confounds - 0.9s, 0.0min
________________________________________________________________________________
[Memory] Calling nilearn.input_data.base_masker.filter_and_extract...
filter_and_extract('/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar002_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz',
,
{ 'allow_overlap': True,
'detrend': True,
'dtype': None,
'high_pass': 0.01,
'low_pass': None,
'maps_img': '/home/varoquau/nilearn_data/msdl_atlas/MSDL_rois/msdl_rois.nii',
'mask_img': None,
'smoothing_fwhm': None,
'standardize': True,
't_r': 2,
'target_affine': array([[ 4., 0., 0., -78.],
[ 0., 4., 0., -111.],
[ 0., 0., 4., -51.],
[ 0., 0., 0., 1.]]),
'target_shape': (40, 48, 35)}, confounds=[ array([[ 0.127944, ..., -0.087084],
...,
[-0.015679, ..., -0.02587 ]]),
'/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar002_task-pixar_desc-reducedConfounds_regressors.tsv'], dtype=None, memory=Memory(location=nilearn_cache/joblib), memory_level=1, verbose=2)
[NiftiMapsMasker.transform_single_imgs] Loading data from /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar002_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
[NiftiMapsMasker.transform_single_imgs] Resampling images
[NiftiMapsMasker.transform_single_imgs] Extracting region signals
[NiftiMapsMasker.transform_single_imgs] Cleaning extracted signals
_______________________________________________filter_and_extract - 6.3s, 0.1min
Processing file /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar003_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
________________________________________________________________________________
[Memory] Calling nilearn.image.image.high_variance_confounds...
high_variance_confounds('/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar003_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz')
__________________________________________high_variance_confounds - 0.9s, 0.0min
________________________________________________________________________________
[Memory] Calling nilearn.input_data.base_masker.filter_and_extract...
filter_and_extract('/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar003_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz',
,
{ 'allow_overlap': True,
'detrend': True,
'dtype': None,
'high_pass': 0.01,
'low_pass': None,
'maps_img': '/home/varoquau/nilearn_data/msdl_atlas/MSDL_rois/msdl_rois.nii',
'mask_img': None,
'smoothing_fwhm': None,
'standardize': True,
't_r': 2,
'target_affine': array([[ 4., 0., 0., -78.],
[ 0., 4., 0., -111.],
[ 0., 0., 4., -51.],
[ 0., 0., 0., 1.]]),
'target_shape': (40, 48, 35)}, confounds=[ array([[-0.089762, ..., -0.062316],
...,
[-0.065223, ..., -0.022868]]),
'/home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar003_task-pixar_desc-reducedConfounds_regressors.tsv'], dtype=None, memory=Memory(location=nilearn_cache/joblib), memory_level=1, verbose=2)
[NiftiMapsMasker.transform_single_imgs] Loading data from /home/varoquau/nilearn_data/development_fmri/development_fmri/sub-pixar003_task-pixar_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz
[NiftiMapsMasker.transform_single_imgs] Resampling images
[NiftiMapsMasker.transform_single_imgs] Extracting region signals
[NiftiMapsMasker.transform_single_imgs] Cleaning extracted signals
_______________________________________________filter_and_extract - 6.4s, 0.1min
Computing group-sparse precision matrices
------------------------------------------
.. code-block:: default
from nilearn.connectome import GroupSparseCovarianceCV
gsc = GroupSparseCovarianceCV(verbose=2)
gsc.fit(subject_time_series)
try:
from sklearn.covariance import GraphicalLassoCV
except ImportError:
# for Scitkit-Learn < v0.20.0
from sklearn.covariance import GraphLassoCV as GraphicalLassoCV
gl = GraphicalLassoCV(verbose=2)
gl.fit(np.concatenate(subject_time_series))
.. rst-class:: sphx-glr-script-out
Out:
.. code-block:: none
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 2
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 7
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 1.9s remaining: 0.0s
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 2
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 2
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 1
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 2
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 7
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 10.3s finished
[GroupSparseCovarianceCV.fit] [GroupSparseCovarianceCV] Done refinement 1 out of 4
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 3
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 6
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 1
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 2.0s remaining: 0.0s
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 4
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 3
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 3
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 6
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 17.7s finished
[GroupSparseCovarianceCV.fit] [GroupSparseCovarianceCV] Done refinement 2 out of 4
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 5
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 1
[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 3.2s remaining: 0.0s
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 9
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 10
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 5
[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 22.3s finished
[GroupSparseCovarianceCV.fit] [GroupSparseCovarianceCV] Done refinement 3 out of 4
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 6
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 1
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 1
[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 2.5s remaining: 0.0s
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 11
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 11
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 8
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 5
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[GroupSparseCovarianceCV.fit] Log-likelihood on test set is decreasing. Stopping at iteration 0
[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 16.2s finished
[GroupSparseCovarianceCV.fit] [GroupSparseCovarianceCV] Done refinement 4 out of 4
[GroupSparseCovarianceCV.fit] Final optimization
[GroupSparseCovarianceCV.fit] tolerance reached at iteration number 23: 6.583e-04
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
....[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 0.1s remaining: 0.0s
................[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 0.4s finished
[GraphicalLassoCV] Done refinement 1 out of 4: 0s
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
....[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 0.1s remaining: 0.0s
................[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 0.5s finished
[GraphicalLassoCV] Done refinement 2 out of 4: 0s
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
....[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 0.1s remaining: 0.0s
................[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 0.5s finished
[GraphicalLassoCV] Done refinement 3 out of 4: 1s
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
....[Parallel(n_jobs=1)]: Done 1 out of 1 | elapsed: 0.1s remaining: 0.0s
................[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 0.5s finished
[GraphicalLassoCV] Done refinement 4 out of 4: 1s
[Parallel(n_jobs=1)]: Using backend SequentialBackend with 1 concurrent workers.
[Parallel(n_jobs=1)]: Done 5 out of 5 | elapsed: 0.0s finished
[graphical_lasso] Iteration 0, cost 1.68e+02, dual gap 1.153e+00
[graphical_lasso] Iteration 1, cost 1.68e+02, dual gap -3.137e-03
[graphical_lasso] Iteration 2, cost 1.68e+02, dual gap 5.307e-04
[graphical_lasso] Iteration 3, cost 1.68e+02, dual gap 2.034e-04
[graphical_lasso] Iteration 4, cost 1.68e+02, dual gap 3.499e-05
GraphicalLassoCV(verbose=2)
Displaying results
-------------------
.. code-block:: default
atlas_img = msdl_atlas_dataset.maps
atlas_region_coords = plotting.find_probabilistic_atlas_cut_coords(atlas_img)
labels = msdl_atlas_dataset.labels
plotting.plot_connectome(gl.covariance_,
atlas_region_coords, edge_threshold='90%',
title="Covariance",
display_mode="lzr")
plotting.plot_connectome(-gl.precision_, atlas_region_coords,
edge_threshold='90%',
title="Sparse inverse covariance (GraphicalLasso)",
display_mode="lzr",
edge_vmax=.5, edge_vmin=-.5)
plot_matrices(gl.covariance_, gl.precision_, "GraphicalLasso", labels)
title = "GroupSparseCovariance"
plotting.plot_connectome(-gsc.precisions_[..., 0],
atlas_region_coords, edge_threshold='90%',
title=title,
display_mode="lzr",
edge_vmax=.5, edge_vmin=-.5)
plot_matrices(gsc.covariances_[..., 0],
gsc.precisions_[..., 0], title, labels)
plotting.show()
.. rst-class:: sphx-glr-horizontal
*
.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_multi_subject_connectome_001.png
:alt: plot multi subject connectome
:class: sphx-glr-multi-img
*
.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_multi_subject_connectome_002.png
:alt: plot multi subject connectome
:class: sphx-glr-multi-img
*
.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_multi_subject_connectome_003.png
:alt: plot multi subject connectome
:class: sphx-glr-multi-img
*
.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_multi_subject_connectome_004.png
:alt: plot multi subject connectome
:class: sphx-glr-multi-img
*
.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_multi_subject_connectome_005.png
:alt: plot multi subject connectome
:class: sphx-glr-multi-img
*
.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_multi_subject_connectome_006.png
:alt: plot multi subject connectome
:class: sphx-glr-multi-img
*
.. image:: /auto_examples/03_connectivity/images/sphx_glr_plot_multi_subject_connectome_007.png
:alt: plot multi subject connectome
:class: sphx-glr-multi-img
.. rst-class:: sphx-glr-timing
**Total running time of the script:** ( 1 minutes 46.088 seconds)
.. _sphx_glr_download_auto_examples_03_connectivity_plot_multi_subject_connectome.py:
.. only :: html
.. container:: sphx-glr-footer
:class: sphx-glr-footer-example
.. container:: binder-badge
.. image:: https://mybinder.org/badge_logo.svg
:target: https://mybinder.org/v2/gh/nilearn/nilearn.github.io/master?filepath=examples/auto_examples/03_connectivity/plot_multi_subject_connectome.ipynb
:width: 150 px
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_multi_subject_connectome.py `
.. container:: sphx-glr-download sphx-glr-download-jupyter
:download:`Download Jupyter notebook: plot_multi_subject_connectome.ipynb `
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery `_