文章来源于"脑机接口社区"CCA典型关联分析原理与Python案例​mp.weixin.qq.com

Rose今天分享一下CCA的相关原理以及Python应用，CCA在EEG等脑电数据的特征提取中使用很多，很有必要熟悉其原理。

CCA典型相关分析

CCA(canonical correlation analysis)利用综合变量对之间的相关关系来反映两组指标之间的整体相关性的多元统计分析方法。它的基本原理是：为了从总体上把握两组指标之间的相关关系，分别在两组变量中提取有代表性的两个综合变量U1和V1(分别为两个变量组中各变量的线性组合)，利用这两个综合变量之间的相关关系来反映两组指标之间的整体相关性。

1936年，Hotelling提出典型相关分析。考虑两组变量的线性组合, 并研究它们之间的相关系数p(u,v).在所有的线性组合中, 找一对相关系数最大的线性组合, 用这个组合的单相关系数来表示两组变量的相关性, 叫做两组变量的典型相关系数, 而这两个线性组合叫做一对典型变量。在两组多变量的情形下, 需要用若干对典型变量才能完全反映出它们之间的相关性。下一步, 再在两组变量的与u1,v1不相关的线性组合中, 找一对相关系数最大的线性组合, 它就是第二对典型变量, 而且p(u2,v2)就是第二个典型相关系数。这样下去, 可以得到若干对典型变量, 从而提取出两组变量间的全部信息。

典型相关分析的实质就是在两组随机变量中选取若干个有代表性的综合指标(变量的线性组合), 用这些指标的相关关系来表示原来的两组变量的相关关系。这在两组变量的相关性分析中, 可以起到合理的简化变量的作用; 当典型相关系数足够大时, 可以像回归分析那样, 由- 组变量的数值预测另一组变量的线性组合的数值。

原理描述

案例实现

# 导入工具包

import h5py

import rcca

import sys

import numpy as np

import cortex

zscore = lambda d: (d-d.mean(0))/d.std(0)

第一步：加载数据

请从CRCNS下载数据：http://crcns.org/data-sets/vc/vim-2以下分析假定该数据位于当前目录中名为“ data”的目录中。

data = []

vdata = []

numSubjects = 3

# subjects 是3个受试者列表.

subjects = ['S1', 'S2', 'S3']

# xfms 是Pycortex中变换名称的列表，该名称用于对齐每个受试者的功能和解剖数据。

xfms = ['S1_xfm', 'S2_xfm', 'S3_xfm']

dataPath ="./data/VoxelResponses_subject%d.mat"

for subj in range(numSubjects):

# 打开数据

f = h5py.File(dataPath % (subj+1),'r')

# 获取数据大小

datasize = (int(f["ei"]["datasize"].value[2]),int(f["ei"]["datasize"].value[1]),int(f["ei"]["datasize"].value[0]))

# 从Pycortex获取皮质面罩

mask = cortex.db.get_mask(subjects[subj], xfms[subj], type = 'thick')

# 获取该受试者的训练数据

data_subj = np.nan_to_num(zscore(np.nan_to_num(f["rt"].value.T)))

data.append(data_subj.reshape((data_subj.shape[0],)+datasize)[:, mask])

# 获取受试者的验证数据

vdata_subj = np.nan_to_num(zscore(np.nan_to_num(f["rv"].value.T)))

vdata.append(vdata_subj.reshape((vdata_subj.shape[0],)+datasize)[:, mask])

第二步：定义CCA参数

# 这里设置1e-4和1e2之间的一系列正则化值

regs = np.array(np.logspace(-4, 2, 10))

# 这里考虑3到10之间的成分数量

numCCs = np.arange(3, 11)

# 初始化cca模型

cca = rcca.CCACrossValidate(numCCs=numCCs, regs=regs)

第三步：对数据训练，分析并保存分析结果

"""

说明：

由于数据量大，此分析的计算量很大。

在笔记本中运行它会花费大量时间，因此建议对其进行并行化和/或在计算机群集上运行它，然后加载结果以进行可视化。

"""

# 利用cca训练数据

cca.train(data)

# 利用cca对验证数据进行验证

cca.validate(vdata)

# 计算方差，解释每个体素中的验证响应

cca.compute_ev(vdata)

# 保存分析结果

cca.save("./data/CCA_results.hdf5")

第四步：可视化分析结果

# 导入可视化工具包

%matplotlib inline

import matplotlib.pyplot as plt

# 导入Brewer色彩图以进行可视化

from brewer2mpl import qualitative

nSubj = len(cca.corrs)

nBins = 30

bmap = qualitative.Set1[nSubj]

f = plt.figure(figsize = (8, 6))

ax = f.add_subplot(111)

for s in range(nSubj):

# 绘制所有三个对象的所有体素之间的相关性直方图

ax.hist(cca.corrs[s], bins = nBins, color = bmap.mpl_colors[s], histtype="stepfilled", alpha = 0.6)

plt.legend(['Subject 1', 'Subject 2', 'Subject 3'], fontsize = 16)

ax.set_xlabel('Prediction correlation', fontsize = 20)

ax.set_ylabel('Number of voxels', fontsize = 20)

ax.set_title("Prediction performance across voxels", fontsize = 20)

# p <0.05时的显着性阈值(针对多次比较进行了校正)

# 重要性是使用渐近方法计算的(有关详细信息，请参见论文文本)

thresh = 0.0893

ax.axvline(x = thresh, ymin = 0, ymax = 7000, linewidth = 2, color = 'k')

ax.text(thresh+0.05, 5000, 'p<0.05', fontsize = 16)

ax.set_xticklabels(0.1*np.arange(-8, 11, 2), fontsize = 16)

ax.set_yticklabels(np.arange(0, 10000, 1000), fontsize = 16)

该图显示了皮质图上一个对象的跨学科预测结果，即预测的验证响应和实际的验证响应之间的相关性。不重要的相关性(p <0.05，已针对多个比较进行校正)设置为0。

import cortex

from matplotlib import cm

from copy import deepcopy

subj = 0

subjName = "S1"

subjTransform = "S1_xfm"

corrs = deepcopy(cca.corrs[subj])

# 将所有低于显着性阈值的体素设置为0

corrs[corrs

_ = cortex.quickflat.make_figure(cortex.Volume(corrs, subjName,

subjTransform,

cmap = cm.PuBuGn_r,

vmin = 0., vmax = 1.), with_curvature = True)

参考点击上面动图扫描二维码，关注我们