直接连接的逆协方差_逆协方差矩阵

直接连接的逆协方差 inverse covariance for direct connections
对于直接连接的逆协方差矩阵的解释：
给定一组时间序列（如前一节中提取的），功能性连接组是一组表示区域之间大脑相互作用的连接。这里我们展示了使用稀疏逆协方差来提取只关注区域间直接相互作用的功能性连接体。
这个非常有用。比如你在前面计算出了感兴趣ROI，然后你想作ROI的功能连接，如果是直接提取时间序列作相关性矩阵，虽然你去掉了混淆信号，但是这些连接还是不太精确，因为他们也指向了其他labels。
所以最正确的做法，就是用直接连接的逆协方差矩阵。

# 导入稀疏的逆协方差估计器
from sklearn.covariance import GraphicalLassoCV
estimator = GraphicalLassoCV()
# 提取时间序列
estimator.fit(time_series)  
# 协方差矩阵和逆协方差矩阵（精度矩阵）可以分别在估计器的covariance_和precision_属性中找到
estimator.covariance_  #协方差矩阵
estimator.precision_  #逆协方差矩阵，这就是我们想要的
# Display the covariance
from nilearn import plotting
# 这个是显示协方差矩阵
plotting.plot_matrix(estimator.covariance_, labels=labels,
                     figure=(9, 7), vmax=1, vmin=-1,
                     title='Covariance')
# 协方差矩阵可视化
from nilearn import datasets
atlas = datasets.fetch_atlas_msdl()
coords = atlas.region_coords
plotting.plot_connectome(estimator.covariance_, coords,
                         title='Covariance')
# 稀疏的逆协方差矩阵可视化
plotting.plot_matrix(-estimator.precision_, labels=labels,
                     figure=(9, 7), vmax=1, vmin=-1,
                     title='Sparse inverse covariance')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

这里的逆协方差矩阵，使用的是-。即：-estimator.precision_，因为负的话，得到的就是偏相关。偏相关就是把其他所有ROI当成协变量，然后去掉他们的影响，类似于协方差分析，很实用。

plotting.plot_connectome(-estimator.precision_, coords,
                         title='Sparse inverse covariance')
plotting.show()
# 可视化二：在web上进行可视化
view = plotting.view_connectome(-estimator.precision_, coords)
# In a Jupyter notebook, if ``view`` is the output of a cell, it will
# be displayed below the cell
view
1
2
3
4
5
6
7
8

偏相关分析是指当两个变量同时与第三个变量相关时，将第三个变量的影响剔除，只分析另外两个变量之间相关程度的过程，判定指标是相关系数的R值。

多个被试的ROI偏相关分析

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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

首先为了方便，定义一个可视化的方法

# 获取数据
from nilearn import datasets
msdl_atlas_dataset = datasets.fetch_atlas_msdl()
rest_dataset = datasets.fetch_development_fmri(n_subjects=n_subjects)
# 提取ROI信号
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,
    high_variance_confounds=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)
    region_ts = masker.transform(func_filename, confounds=confound_filename)
    subject_time_series.append(region_ts)
# 计算group的稀疏逆协方差矩阵
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))
# 显示结果
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()
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56

解读：
1、首先定义了一个可视化的方法，避免重复工作，不过可惜的是，nilearn只能处理最多40个被试。
2、导入数据
3、用masker.fit提取ROI时间序列（多个被试，所以就是list）
4、用GroupSparseCovarianceCV( )代替GraphicalLassoCV( )来计算group的矩阵
5、把稀疏的逆协方差矩阵进行负值化，得到偏相关，然后可视化。