2
2
2
·2771·
113
129.
[3 ] Hartigan JA. Clustering Algorithm [ M] . New York : Wiley Press ,1975.
[4 ] Talairach J , Tournoux P. Coplanar Stereotaxic Atlas of the Human Brain
[M] . Stuttgart : Thieme Medical ,1988.
[5 ] Hu ZH , Wu YG, Chen FY, et al. AD RAT model applied in fMRI[J ] .
中国医学影像技术 2004 年第 20 卷第 11 期 Chin J Med Imaging Technol ,2004 ,Vol 20 ,No 11
Acta Psychologica Sinica ,2002 ,34 (Suppl) :94
胡正珲 ,吴义根 ,陈飞燕 ,等。AD RAT 模型的 fMRI 研究 [J ] . 心理学
报 ,2002 ,34 (增刊) :94
95.
95.
[6 ] DeGroot MH. Probability and Statistics [ M ] . Chicago : Addison
Wesley
Inc ,1975.
Basic principle of SPM : an introduction
- Part Ⅱ: applications to PET and fMRI
WU Yi
gen , LI Ke
( Key Laboratory of Nuclear Analysis Techniques , Institute of High Energy Physics ,
Chinese Academy of Sciences , Beijing 100039 , China)
[ Abstract] The basic mathematic principle of the SPM was introduced in the first part of this essay. Then some useful application models
in the PET and fMRI data set were emphasized in the second part of the paper.
[ Key words] SPM; Tomography , emission
computed ; Functional magnetic resonance imaging
SPM 软件包数据处理原理简介
———第二部分 :应用于 PET 及 fMRI
(中国科学院高能物理研究所核分析技术重点实验室 ,北京 100039)
吴义根 ,李 可
[摘 要 ] 在本文第一部分介绍了 SPM 软件包对数据进行预处理的基本数学原理后 ,此第二部分将重点介绍它在 PET 和
fMRI 数据处理中的一些具体应用模型 。
[关键词 ] 统计参数图 ; 体层摄影术 ,发射型计算机 ; 功能磁共振成像
[中图分类号 ] TP391 ; R445. 2 [文献标识码 ] A [文章编号 ] 1003
3289 (2004) 11
1772
04
用 SPM 软件做数据分析中 ,最重要的工作之一就是如何
建立 (或选择) 线性模型 ,下面分正电子发射断层扫描 ( PET)
实验和功能磁共振成像 (fMRI) 实验 ,介绍脑功能研究中常用
的几种模型 。
1 正电子发射断层扫描实验模型[1 ]
1. 1 整体效应 正电子发射断层扫描实验中 ,由于需要对多
个被试以及一个被试在不同时间点的多次测量 ,整体脑血流
(gCBF , PET/ H2O15 实验) 或整体脑葡萄糖代谢率 ( gCMRGlc ,
PET/ FDG实验) 将不相同 ;同时 ,一般由于不做定量测量 ,数
据值可用局部放射性活度 (rA) 来代替局部脑血流 (rCBF) 或者
局部脑葡萄糖代谢率 (rCMRGlc) ,此时的整体活度 (gA) 还取决
于注射的放射性标记药物剂量以及它们在脑部积累的比例 ,
这种整体效应会影响到局部放射性活度值 ,数据处理中需要
消除它 。这里 ,第 i 次扫描的整体效应定义为局部脑血流
rCBF、rCMRGlc 或 rA 对全脑中所有的像素取平均 : gi = Yi
·=
(
k = Yi
k) / K ,其中 K 是像素总数 。
K
k = 1 Yi
消除整体效应的最简单方法就是将不同扫描的数据乘以
一个比例因子 ( scaling model ,正比模型) ,使得所有扫描的 gi
都统一到某个常数 ,SPM 中选择了 rCBF 的典型值 50 ml/ (min·
dl) ,注意到这个数值可以任意的 ,对于 PET/ FDG 实验等依然
适用 。这样 Yi
,由于有比例因子 ,
如果假设 Y’i
k 的误差就不一样 ,需要
用带权重的回归方法去拟合 , F 检验方法等就不能适用 ,从这
点讲这种方法不再是简单的了 。更精确的方法是协方差分析
模型 (ANCOVA model) ,就是将整体效应作为一个影响因素而
放在线性模型中 :
Yi
k 将变为 Y’i
k 的误差对 i 相同 ,则 Yi
(12)
k 是其他因素的总和 , ζk 是反映整体效应的参数 , g·
式中 μi
是对所有扫描整体效应的平均 ,εk 是误差项 。
k + ζk ( gi - g·) +εk
k/ ( gi/ 50)
k = μi
[ 基金项目 ] 国家重点基础研究发展规划资助项目 ( G1999054006) 。
[作者简介 ] 吴义根 (1964 - ) ,男 ,江苏东台人 ,硕士 ,副研究员。研究方
向 :脑功能成像原理及应用。E
[ 收稿日期 ] 2004
mail : wuyg @ihep . ac. cn
07
24
选择何种方法目前还没有一个标准 ,它取决于具体的实
验 。在定量测量实验中 ,对于正常状态下的正常被试 ,gCBF
或 gCMRGlc 通常变化很小 ,并且其数值也远离零 ,这时协方差
分析模型比正比模型能更好地反映局部与整体的关系 。在非
2
中国医学影像技术 2004 年第 20 卷第 11 期 Chin J Med Imaging Technol ,2004 ,Vol 20 ,No 11
·3771·
定量测量实验中 ,即使被试的 gCBF 或 gCMRGlc 非常稳定 ,但
由于注射药物剂量 、它们在脑部积累的比例以及机器的灵敏
度等的不同 ,gA 的变化将是较大的 。同时 ,对于 gCBF 或 gCM
RGlc 为常数时 ,局部与整体活度之间有着正比关系 。而更重
要的 ,由于这类实验是计数 (count) 测量 ,数据服从泊松 ( Pois
k 其误差也大 ,乘上比例因子后 ,不同扫描
son) 分布 ,大的 Yi
k 的误差的差异性反而变得小了 (但没有完全消除) ,此时正
Yi
比模型中的缺点自然变为了优点 ,这种情况下 ,正比模型明显
优于协方差分析模型 。当然 ,如果注射药物剂量等因素大致
相同 ,协方差分析模型依然是很好的方法 。
上述整体效应假设了与研究任务无关 ,因而在线性模型
中将被看成是不感兴趣的部分而不参与后面统计分析 。对于
与任务有关的情况 ,需要考虑用其他的实验来判断它们之间
的关系 。
1. 2 单个被试激活模型 (activation design) 这是最简单的一
种实验模型 ,它是用来研究单个被试在不同外部刺激条件下
大脑活动差异 ,每种条件下需要多个扫描 (如果可以假定不同
条件下实验数据的误差相同 ,那么可以只在其中一种条件下
做多次扫描而其他条件下只做单次扫描 。但只要实验条件许
可 ,每种条件下的多次扫描既可以提高实验精度 ,有可以减小
模型的偏差) 。如果采用 ANCOVA 模型来消除整体效应 ,并
假设可以多于两种实验条件 (如 Q 种条件) ,类似于方程 (12) ,
有下列模型 :
qi
qi = ξk
Yk
q + μk + ζk ( gqi - g··) +εk
(13)
qi 代表第 k 像素在条件 q 下第 i 次扫描的数据 , ξk
q 是第 k 像
Yk
素条件 q 对 Yk
qi 的贡献 , ( gqi - g··) 是在条件 q 下第 i 次扫描
的整体效应与其所有扫描的平均值之差 , ξk 代表整体效应差
qi 的贡献 , μk 代表第 k 像素所有扫描的平均值 (不是直
异对 Yk
接求平均 ,而是作为未知参数待求) 。后一项的引入一方面是
为了突出ξq
k 的影响 ,而没有任何实际意义 ;另一方面 ,它的引
入导致对ξq
k 加了约束条件 :
Q
q =1
ξq
k = 0 ,如果选取新的矩阵
元 ,可以减少模型中的一个感兴趣参数 ,如对于两种条件的实
验 ,如果对应条件一的矩阵元取 - 1 ,对应条件二的取 1 , ξq
k
可以只为一个ξk
,从而减少参与后面统计检验的参数 。由于
这项的引入 ,设计矩阵将不是满秩的 ,参数拟合时可以用降秩
方法去处理 。和ζk 一样 , μk 也属于不感兴趣的部分而不参
于其后的统计检验 ,后面的模型里也会引入这样的不感兴趣
项 。
图 1a 是该模型设计矩阵图 (5 列 ,12 行) 的一个举例 ,它
是三种条件 ,每种条件各做 4 次扫描的实验 。12 行表示 12 个
扫描 ,5 列中 ,第一 、二和三列对应于条件一 、二和三的矩阵元
(加某种条件时相应位置上的矩阵元为 1 ,否则为 0) ,第四列
对应μk 的矩阵元 (均为 1) ,第五列对应于 ( gqi - g··) (每个扫
描的整体效应差异值) 。需要注意的是 ,这个图上的扫描是按
照前 4 个 、中间 4 个和后 4 个分别对应于条件 1 、2 和 3 的顺
序 ,实际上也可以是随机的 ,只要放对每个扫描的相应条件就
可以了 。在 SPM 软件中 ,该模型对应于“single subject : condi
tions and covariates”中的 covariates 个数为 0 的情况 。
1. 3 单个被试参数模型 (parametric design) 和上个模型一
样 ,它也是用来研究单个被试在不同外部刺激条件的大脑活
动差异 ,但不同的是 ,上面的模型里 ,条件被简单分成了几类 ,
但没有量化它们 ,如果实验中也能对它们做定量测量 (形象化
为‘打分’,这里用 s 代表) ,比如手指运动研究中的运动频率 ,
心理研究中任务的难易程度等 ,那么可以用此参数模型 :
k
i = ρk ( si - s·) + μk + ζk ( gi - g·) +εi
Yk
(14)
式中ρk 是分数变化单位量对 Yk
qi 的贡献 , ( si - s·) 是第 i 次扫
描的分数与所有扫描平均分数之差 。与公式 (13) 相比 ,与任
务有关的贡献不再简单地分为有或无 (对应 1 和 0) ,而可以
根据任务的量值 (分数) 来连续取值 ,公式 (14) 可以定量地求
出大脑活动对外部刺激条件的依赖关系 ,而不仅仅是不同条
件间的活动差异 。图 1b 是该模型设计矩阵图 (3 列 ,12 行) 的
一个举例 ,它是一个参数 ,共做 12 次扫描的实验 。12 行表示
12 个扫描 ,3 列中 ,第一列对应于不同扫描时任务的分数 (距
平差) ,第二列对应 μk 的矩阵元 (均为 1) ,第三列对应于 ( si -
s·) (每个扫描的整体效应差异值) 。在 SPM 软件中 ,该模型
对应于“single subject : covariates only”的情况 。
公式 (14) 只对应于单因素的简单情况 ,对于多个因素 ,可
以将其简单外推 ,类似于 ρk
,增加相应因素的贡献 。同时 ,如
qi 与分数 s 之间不是线性关系 ,比如指数关系 (神经生理
果 Yk
和心理活动中常见的规律) 等 ,如果可以通过变换函数将其转
换成参数的线性关系 ,如指数关系可以通过求对数 : 对应于
si
,用 ln ( si) 代替 ,那么公式 (14) 也可以适用 。
另外 ,还有上述二种模型的混合 ,即外部刺激条件中既含
有可定量的部分 (简称为参数部分) 又含有只简单区分有无的
部分 (简称为条件部分) ,如果二者之间没有关联 ,公式可以参
照 (13) 和 (14) 来建立 。在 SPM 软件中 ,该模型对应于“single
subject : conditions and covariates”中的 covariates 个数不为 0 且
covariates 没有相互作用的情况 。
1. 4 单个被试协变量关联模型设计矩阵图 如果上面提到
的混合模型里 ,参数部分和条件部分之间有关联 ,并假设有 Q
种条件且每种条件下参数取值个数相同 ( R 个) ,那么在 q ( q
= 1 ,2 , …, Q) 条件下 ,当参数取第 r( r = 1 ,2 , …, R) 个值时 ,
第 i 次扫描的数据为 :
qri = ξk
Yk
q +ρk
q ( sqr - s··) + μk + ζk ( gqri - g···) +εk
qri
q 为 q 条件下的平均贡献 ,ρk
(15)
式中ξk
q 为 q 条件下参数的平均微
分效应 。图 1c 是该模型设计矩阵图 (6 列 ,12 行) 的一个举
例 ,它是两种条件 ,一个参数 ,每种条件各做 6 次扫描 (12 次
扫描中参数取值不全相同) 的实验 。12 行表示 12 个扫描 ,5
列中 ,第一 、二列对应于第一 、二种条件下平均贡献的矩阵元 ,
第三 、四列对应于第一 、二种条件下的参数的分数 (距平差) ,
第五列对应 μk 的矩阵元 (均为 1) ,第三列对应于 ( gi - g·)
(每个扫描的整体效应差异值) 。在 SPM 软件中 ,该模型对应
于“single subject : conditions and covariates”中的 covariates 个数不
为 0 且 covariates 有相互作用的情况 。
2
2
·4771·
中国医学影像技术 2004 年第 20 卷第 11 期 Chin J Med Imaging Technol ,2004 ,Vol 20 ,No 11
图 1 单个被试激活模型设计矩阵图 a :简单激活模型 ; b :参数模型 ; c :协变量相互作用模型
covariates”中的 covariates 个数为 0 情况 。
(2) 多个被试协变量关联模型 (图 2b)
jqri = ξk
Yk
q ( sjqr - s···) + γk
j + ζk
q +ρk
j ( gjqri - g····) +εk
jqri
(17)
在 SPM 软件中 ,该模型对应于“multi
subject : conditions and co
variates”中的 covariates 个数不为 0 且 covariates 有相互作用的
情况 。
(3) 多个被试配对 t 检验模型 (图 2c)
j +ε′k
Y′k
jqi = ξ′k
q + γ′k
jqi
(18)
它与方程 (16) 类似 ,只是整体效应用的是正比模型 (所以将
jqi 等换成了 Y′k
jqi ) 。在 SPM 软件中 ,该模型对应于“population
Yk
main effect 2 cond’s , 1 scan/ cond. (paired t
test) ”的情况 ,这是
一种比较常用的模型 。
(4) 多个被试多组多种条件模型 (图 2d)
图 2 多个被试激活模型设计矩阵图 a :简单激活模型 ;b :协变量相
互作用模型 ;c :配对 t 检验模型 ;d :两组两种条件模型 (其中 a ,b 中整
体效应是一个被试、一个被试去处理的 ,c 中整体效应采用了正比的
方法)
1. 5 多个被试模型 对于 PET 实验 ,特别是用 FDG做示踪剂
的实验 ,单个被试的多次测量一般不容易实现 ,同时 ,有的实
验需要反映某个群体的整体特性 ,这些情况下需要多个被试
的研究 。此时 ,前面提到的模型中需要添加与个体差异效应
有关的项 ,为了节省篇幅 ,下面只给出几个相应的方程 、设计
矩阵图及如何从 SPM 软件中选用 ,不再给出具体的解释 。
(1) 多个被试激活模型 (图 2a)
jqi = ξk
Yk
(16)
在 SPM 软件中 ,该模型对应于“multi - subject : conditions and
j ( gjqi - g··) +εk
q + γk
j + ζk
jqi
mjqri
j + ζk
mq + μk
mjqri = ξk
Yk
j ( gmjqri - g·····) +εk
(19)
它也与方程 (16) 类似 ,只是将被试根据需要而分组 (用 m 代
表组) 。在 SPM 软件中 ,该模型对应于“multi
group : conditions
and covariates”中的 covariates 个数为 0 情况 。
2 功能磁共振成像实验
2. 1 高通滤波器 (high pass filter) 由于有较高的时间分辨 ,
周期性的生理脉动 (如呼吸和心跳) 对功能核磁信号的影响需
要考虑 ,SPM 软件中采用高通滤波器来消除它的影响 。高通
滤波器是将周期性的生理脉动归为时间序列功能核磁信号中
,由于它具有周期性 ,可将它用系列余弦
不感兴趣项 YP ( ts)
βrf r ( ts) , s = 1 ,2 , …, N ,其中 N 为
函数展开 : YP ( ts) =
实验总的扫描次数 , f r ( ts) = cos( rπ ts -
tN -
号不被过滤 ,通常高通滤波器中的最大频率分量的周期 (临界
值 ,cut
off) 取实验周期 (trial ,对 block design 为 on 和 off 的时间
和 ,对 event
related design 为两个静息后首次刺激的间隔) 的两
倍 (或略多一点) ,这样 K 就取实验周期数的两倍 (或少于此
) 。为了使有用信
t1
t1
r =1
K
2
2
中国医学影像技术 2004 年第 20 卷第 11 期 Chin J Med Imaging Technol ,2004 ,Vol 20 ,No 11
·5771·
3. 2 F 检验[2 ] 假设模型中的参数β可以分成两部分 :βT =
[β1 ⁝β2 ] T ,同时我们希望检验 H :β1 = 0 ,相应地 ,设计矩阵
也分成两部分 : X = [ X1 ⁝X2 ] , 那么完整模型就是 : Y =
+ε ,当假设 H 正确时 ,此模型就简化为一个
值) ,特殊情况 ,如果实验只有一个周期 ,此时 K 为 1 。低频生
理噪声被归结为 YP ( ts) 后 ,此项在统计推断中没被考虑 ,这
样低频生理噪声将被过滤 ,信号完好被保留 ,这也是高通滤波
器的由来 。特别指出的是 ,SPM 软件中 ,用户只需根据实验特
点给出高通的 cut
2. 2 低通滤波器 (low pass filter) 功能核磁信号时间序列反
映的是脑血流等的特性 ,因而具有血流响应函数 (hrf) 形式 。
数据处理中 ,为了优化信号以及减少白噪声 ,需要将此时间序
列用高斯函数或 hrf 去卷积 ,也就是做时间上的平滑 ,这样高
频噪声将被过滤掉 。
off , K 值就会自动选取 。
对时间序列 做时间平滑后 ,它们将不再是独立的 ,此时
方程 (8) 变为 :
KYk = KXβk +εk
K 是平滑矩阵 。定义 X
= KX ,方程 (20) 的解为 :
βk = ( X
TX
) - 1 X
TKYk
而其误差矩阵为 :
(20)
(21)
TKKT ( X
TX
) - 1
TX
) - 1 X
Var{βk} = σ2
k ( X
(22)
2. 3 固定响应模型 它是将功能核磁信号中感兴趣的部分
用一个固定的响应函数来表示 。对于 block design ,这个响应
函数可以是半正弦函数 、方波形式 ( box
car ,1 表示激活 ,0 或
- 1 表示基线) 以及它们与脑血流响应函数的卷积 ; 对于
related design ,这个响应函数可以就是脑血流响应函数或
event
者加起始时间偏移 (也是待定参数) 。这类模型适用于只研究
激活区而不研究响应曲线的实验 ,好处是待定参数少 ,回归计
算快速 ,缺点是响应模式固定可能代表不了真实情况 。
2. 4 无固定响应模型 它是将功能核磁信号中感兴趣的部
分用系列基准函数展开 (也就是它们的线性组合 ,系数为待定
参数) ,这些基准函数可以是傅立叶级数或γ函数 。它适用于
研究响应曲线和对相应特性不是很清楚的只研究激活区实
验 ,特点是响应特性可以是任意的 (仅在理论上可以这样说 ,
事实上由于计算能力以及数据量的有限性 ,基准函数不可能
取太多 ,也就是这种任意性是有限制的) ,因而具有普适型 ,缺
点是参数可能太多影响回归速度 。
3 统计推断和结果
3. 1 t 检验 在建立模型并估计出其参数后 ,需要对参数做
统计推断 。如果方程 (8) 的设计矩阵 X 是满秩 ,那么参数的
估计值将符合多维正态分布 : ^β~ N [β,σ2 ( XTX) - 1 ] ,它们的
线性组合也服从一维正态分布 : cT^β~ N [ cTβ,σ2 cT ( XTX) - 1 c ]
。由于方差σ2 也需要通过同样的数据拟合得到 (即 ^σ2 ,见公
式 (10) ) ,用 ^σ2 代替σ2 后 , cT^β将不服从正态分布 。根据 fisher
定律 , ^β与^σ2 是独立的 ,参数组合量
cT^β - d
^σ2 cT ( XTX) - 1 c
将服从 t
分布 tn- p ,其中 n 是总的扫描次数 , p 是 X 的秩 , d 是 cTβ的真
值 (待检验量) 。检验 cT^β可以用
cT^β - d
^σ2 cT ( XTX) - 1 c
代替 ,这是
一个检验 H : cTβ = d 的 t 检验问题 。对于 X 非满秩情况 ,
( XTX) - 1 无法直接求出 ,要用求伪逆的方法 (pseudo
inverse) 计
算 ,但上述方法仍然适用 。
[ X1 ⁝X2 ]
β1
…
β2
约化模型 (reduced model) : Y = X2β2 +ε。设完整模型和约化
模型的残差分别 S (β) 和 S (β2)
,那么由于缺少 β1 使得约化
模型比完整模型的残差多了额外的一项 S (β1 | β2) = S (β2)
- S (β) 。在假设 H 成立的情况下 , S (β1 | β2) ~σ2χ2
q ,且与
S (β) 独立 ,其中自由度 q 是 X 与 X2 的秩之差 ,如果 H 不成
立 , S (β1 | β2) 是一个偏心的 χ2 分布 ,但仍然与 S (β) 独立 。
由此有 :
F =
( S (β2) - S (β) ) / ( p - p2)
S (β) / ( n - p)
~ Fp- p2 , n- p
式中 p、p2 分别是 X 与 X2 的秩 , n 是总扫描次数 。 F 值越大
表示 H :β1 = 0 越不成立 。
3. 3 统计参数图 完成了上述统计推断后 ,根据 t 或 F 值以
及相应的阈值 ,将小于阈值的像素的 t 或 F 值置为 0 ,大于阈
值的像素的数值 t 或 F 值保留 ,将得到于此阈值相对应置信
度下的脑功能激活图 ,由于它是对模型的参数做统计分析而
得到的 ,并且是加了阈效应的关于 t 或 F 值的图像 ,所以被称
为统计参数图 。可根据选用的 t 或 F 值记为 SPM{ t} , SPM
{ F} 。同时 ,由于在给出 t 或 F 值后仍需要相应的自由度才
能给出置信度 ,SPM 还提供了与 t 或 F 值及它们相应的自由
度所对应置信度的{0 ,1}正态分布的参数值 Z ,这样也得到关
于 Z 图像 ,记为 SPM{ Z} 。脑功能激活图可以叠加到脑结构
图(T1 加权 MRI 图或 CT 图) 上 ,这个融合图像能清楚地看出
哪些脑区参与了所要研究的功能 。此外 ,根据这些激活图 ,
SPM 还可以输出脑血流等生理信号的响应曲线以及不同状态
的响应特性等 。
4 结语
脑功能成像是神经科学中的一个重要的研究方向 , 其数
据处理在此研究中与实验设计具有同等重要的作用 。作为脑
功能成像研究的有力分析工具之一 ,功能强大的 SPM 软件包
不但能处理 fMRI 信号 ,也可以适用于 PET、SPECT 数据 ,而且
能在这几种不同模态之间的实现相互融合处理 ,对比结果 。
本文在介绍了 SPM 的基本原理之后 ,讲述了其在 PET、fMRI
处理中的典型应用 ,目的是为了更加深入理解此分析工具 ,便
于在以后的研究工作中能够根据自己的实验设计更加灵活地
来分析数据 。
[ 参考文献 ]
[1 ] Members & collaborators of the Wellcome Department of Imaging Neuro
science. Statistical Parametric Mapping Introduction [ EB/ OL ] . http :/ /
www. fil. ion. ucl. ac. uk/ spm/
[2 ] Draper NR , Smith H. Applied Regression Analysis [ M] . 2nd ed. New
York : Wiley Press ,1981.