之前用神经网络炼丹的时候,受制于硬件条件,训练起来极为痛苦,尤其是一开始我头铁一波直接用全脑的4D图像作为输入数据,网络复杂度根本上不去,只要结构稍微复杂一点就必定炸显存。毕竟单独一个时间点上的3D数据大小就能达到几十万(如果体素大小3毫米的话normalise之后就是61×73×61)。那就要考虑一些数据降维的方法了,先行提取一些特征出来(答辩的时候还有个老师问我为啥不直接用全脑数据来着,我说穷啊,要是硬件条件允许我也想试试)。
功能连接(Functional Connectivity)矩阵
先说说
功能连接矩阵吧,这应该是用得最多而且也比较简单的一种方法了。思路也比较简单:
fMRI是四维的时间序列,也就是在一段时间内不同时间点上的三维图像。那么,只要提取出某一个体素在全部时间上的数据也就得到了
这个体素的时间序列。这个
时间序列也就代表着这个体素之内
血氧水平随时间的变化,也就反映了功能的表现。 所以,从
时间序列之间的相关性,就可以得到功能的相关性。(当然这里存在问题,这也是功能连接以及下面有效连接这两种方法的缺点,因为这两种方法两两计算脑区连接度,那么就考虑不到是否存在相互作用的第三点的问题。)
计算功能连接一般使用
皮尔森相关系数。也有研究使用其他方法,以期得到更稀疏的功能连接(但我忘了具体咋算的了)。两组时间序列的皮尔森相关系数由下式得,也就是这两段时间序列的
协方差与除以两段时间序列的标准差之积:
如果用Python写的话可以直接使用numpy提供的函数:
corr = np.corrcoef(X, Y)
不过很多工具箱也提供了直接计算功能连接矩阵的功能(如
DynamicBC和
GRETNA,这些同样也可以用来计算有效连接等其他脑网络),就不需要自己手动写脚本提取了。
皮尔森相关系数的结果
取值范围为
。当
时,表示两组时间序列正相关,反之则为负相关。当等于零时,代表两组时间序列相互独立,无相关关系。由此可以得出某两组脑区在某一段时间内的功能是相互协同的还是相互拮抗的。
通过
计算各个脑区的平均时间序列,
两两计算相关系数,就可以得到这段时间内全脑的相关性矩阵,即功能连接矩阵。
功能连接矩阵,使用AAL90模板,共90个脑区,因此连接矩阵为90x90
另外,还可以使用动态窗口法在不同的时间窗内计算功能连接矩阵,得到随时间变化的功能连接矩阵。
有效连接(Effective Connectivity)矩阵
有效连接与功能连接相似,但是计算的是
脑区之间的因果关系,也就是脑区之间的互相影响。与功能连接不同,
功能连接是
无向有权图,得到的矩阵是
对称矩阵;而
有效连接使用因果关系,是有方向的,也就是说会得到
有向有权图,所以有效连接矩阵是
不对称的。
常用的方法有
动态因果模型(Dynamic Causal Modelling, DCM)和
格兰杰因果关系(Granger Causality)。我常用的方法是格兰杰因果。其定义为:对于两个平稳的时间序列X和Y,如果通过X和Y过去的值来预测X的当前值要比仅使用X自身的过去值预测X的当前值要更准确,则认为Y是X的格兰杰原因。
打个比方就是,
半个小时前天气很好很暖和,但是突然开始下雪,大雪下了半个小时,现在天气很冷。假设
气温是X,
天气状况是Y,那么只使用
X的过去值(暖和)来预测
现在(冷)明显是不准确的。那么
同时考虑X和Y,因为
下雪所以
使暖和的温度变冷了,这就要更为准确,所以我们可以认为
天气状况Y是气温X的原因。
那么构建因果关系首先要构建X和Y的自回归模型:
其中a和b是模型的自回归系数,是常数;p为模型阶数,即使用当前时刻前1到p个时间点来预测当前值;ε是预测误差。因为现在的自回归模型是
单变量的自回归模型,即只使用X来预测X(或只使用Y来预测Y),因此
预测误差只取决于X(或Y)的过去值。我们可以使用预测误差的方差来评估模型精度。
接下来咱们考虑
X和Y互相影响的情况,这时候在上式中X加入一项Y,在Y中加入一项X,构造双变量自回归模型:
c和d同样是自回归系数。这时
预测误差ε会同时受到X和Y的影响。同样使用预测误差的方差来评估模型精度。那么,
当双变量自回归模型的精度更高时,例如:
。
这个关系的度量通过下式计算:
反之,考虑X对Y的因果关系也同理。
计算的方法与功能连接相似,提取出脑区的平均时间序列两两计算。得到的矩阵也类似。
独立成分分析(Independent Component Analysis,ICA)
这个方法就不太一样了,
独立成分分析可以把全脑信号分解成
若干个相互独立的信号源产生的
信号的叠加。例如,某个脑区同时受到了多个脑区产生的信号的影响,此时这个区域的信号是这些原始信号的叠加,那么我们可以通过这种方法分离出原始信号源。
全脑的 fMRI 数据可以认为是由若干组相互独立的信号的线性叠加。因此利 用 ICA 方法可以提取出空间上独立的信号源。(也可以进行时间独立的分解。)
独立成分分析的问题可以描述为如下:
简单来说,
就是找到一个矩阵W,能够将原始信号分解为若干独立的信号。
那么,用什么方法可以找到W呢?在这里介绍一种
Infomax方法,按照信息极大化原则进行。它的思想是,当
系统的输出熵极大的时候(熵极大也就代表
输出中包含的信息量极大),这个时候输出的
各分量之间的互信息最小(也就是各分量中包含的
重复信息得到了去除),这个时候可以认为各个输出成分是互相独立的。
Infomax算法原理图
Infomax的实现和神经网络非常相似,也就是
用W将初始信号分解为若干分量后,把这若干个输出经过一个
非线性函数得到结果,
计算结果的总熵。然后在不断的迭代中
调节W的值,使得总熵极大。
当然,我们还是可以使用MATLAB工具箱(例如Group ICA of fMRI Toolbox,GIFT)来帮助我们进行分析。
可以使用工具箱先对数据进行成分数量的估计,然后按照估计的结果设置要分离出多少个独立成分。那么最后我们可以得到一个尺寸为
独立成分数x时间序列长度的二维矩阵(也就是分离出的
各个独立信号源脑区的平均时间序列)。还可以通过GIFT工具箱查看各个独立成分的详细信息,例如位置、成分内各个体素位置的信号等等。
总结
以上大概介绍了三种特征提取的方式(也基本就是我毕设里用的玩意了),其实适用于fMRI数据的降维手段还有很多,例如张量的方法,像Tucker分解或者CP分解(张量方法我还不熟悉,可能以后会再了解了解),当然也有一些研究选择了直接用全脑的数据去炼丹……这样会带来巨大的时空开销,我也没实际仔细调试过,不好评价效果如何。
总之就是根据实际情况,我个人觉得ICA的方法更灵活,但也要更加麻烦一些。三种方法应用都很广泛,根据实际情况选择吧~