21、用于PET和任务fMRI数据建模的稀疏概率并行因子分析

用于PET和任务fMRI数据建模的稀疏概率并行因子分析

1. 引言

在处理大型数据集时，主成分分析（PCA）及其概率形式（PPCA）是常用的降维工具。PCA通过寻找描述最大方差方向的正交分量来实现降维，但选择保留的分量数量可能存在问题，且即使解释了大部分方差的分量也可能包含难以解释的小权重。稀疏版本的PCA算法通过修剪整个分量或单个权重来解决这些问题。

神经科学数据具有多模态的特点，虽然可以对沿某一模式（如时间）连接的数据执行PCA，以识别另一模式（如受试者）共有的分量，但这种方法会丢弃模式特定的信息。相比之下，并行因子分析（PARAFAC）等多向分解方法能够保留数据的内在结构，并且在满足模型假设的情况下，对噪声的敏感性较低。此外，PARAFAC模型在温和条件下是唯一的（除了缩放和排列），提供了更具解释性的表示。

稀疏多向模型在神经科学等领域具有重要意义。大脑被证明是按网络组织的，对于某些特定任务，如运动任务，大脑的不同区域会活跃，因此可以预期出现空间稀疏模式。当这种类型的任务在多个受试者中执行时，可以通过多向分解利用数据的内在结构。

本文开发了一种具有时间依赖和受试者特定各向同性噪声的全贝叶斯稀疏概率PARAFAC（SP - PARAFAC）模型，并展示了如何通过简单改变稀疏先验在SP - PARAFAC和概率PARAFAC（VB - PARAFAC）之间轻松转换。基于变分贝叶斯推理给出了近似解，并研究了这些模型在PET和基于任务的fMRI数据中的适用性。

2. 概率PCA回顾

概率PCA的初始公式定义了一个模型，将观测值 $x$ 与投影在以 $m$ 为原点的 $K$ 维超平面 $W$ 上的潜在变量 $z$ 相关联，并加上加性高斯噪声 $\epsilon$，即 $x = Wz + m + \epsilon$。其中，潜在变量和噪声假设遵循各向同性高斯分布：$z \sim N(0, I)$ 和 $\epsilon \sim N(0, \tau^{-1}I)$，$\tau$ 是噪声的精度。

在后续工作中，作者通过变分推理（VB - PCA）对PCA进行了全贝叶斯处理，并对矩阵 $W$ 的列使用了分层ARD先验 $P(W|\alpha)$，允许自动选择分量的数量。这里，$P(W|\alpha)$ 定义为特定于每列的多元高斯分布，$\alpha$ 定义为精度：
[P(W|\alpha) = \prod_{k} \left(\frac{\alpha_k}{2\pi}\right)^{\frac{T}{2}} \exp\left(-\frac{1}{2}\alpha_k |w_k|^2\right)]
[P(\alpha) = \prod_{k} \Gamma(\alpha_k|a_{\alpha_k}, b_{\alpha_k})]
其中，$w_k$ 是 $W$ 的第 $k$ 列。为了完成贝叶斯规范，其余参数与宽泛的先验相关联：$m \sim N(0, \beta^{-1}I)$ 和 $\tau \sim \Gamma(a_{\tau}, b_{\tau})$。

对 $W$ 的替代先验的研究导致了PCA的全稀疏公式（SP - PCA），其中对 $W$ 的元素而不是列施加稀疏先验。理论上，这允许模型识别模型阶数（真实的 $K$）并忽略嘈杂或无关的体素（特征）。在研究的不同先验中，杰弗里先验被证明能给出具有最高累积解释方差的稀疏分量。使用杰弗里先验时，$P(W|\alpha)$ 的条件分布和 $\alpha$ 的先验变为：
[P(w_{v,k}|\alpha_{v,k}) = \left(\frac{\alpha_{v,k}}{2\pi}\right)^{\frac{1}{2}} \exp\left(-\frac{1}{2}\alpha_{v,k} w_{v,k}^2\right)]
[Jeffrey’s(\alpha_{v,k}) \sim \sqrt{E_{P(w_{v,k}|\alpha_{v,k})}\left(\left(\frac{\partial}{\partial\alpha} \ln P(w_{v,k}|\alpha_{v,k})\right)^2\right)} = \frac{1}{\alpha_{v,k}}]

3. 概率并行因子分析

多向数据可以看作是张量结构。本文考虑维度为 $V \times T \times S$ 的三维张量，其中 $V$ 表示体素，$T$ 表示时间，$S$ 表示受试者。与PCA类似，PARAFAC模型试图识别一个大小为 $V \times K$ 的矩阵 $W$，其列是 $S$ 个受试者共有的分量。与PCA不同的是，PARAFAC模型允许对分量进行单独缩放，可用于表征受试者间的变异性，并且还考虑了受试者特定的时间噪声。

每个受试者 $s$ 在每个时间点 $t$ 的数据通过以下公式重建：
[x_t^{(s)} = W \text{diag}(\delta^{(s)}) z_t + m^{(s)} + \epsilon_t^{(s)}]
其中，$W$ 和 $z$ 在受试者之间是共同的，但我们对每个受试者特定的时间依赖噪声 $\epsilon_t^{(s)}$（精度为 $\tau_t^{(s)}$）、受试者特定的均值 $m^{(s)}$ 和受试者特定的分量缩放 $\delta^{(s)}$ 进行建模。

模型的似然性通过以下公式评估：
[P(X | W, Z, m, \tau, \delta) = \prod_{s} \prod_{t} N\left(x_t^{(s)} | W \text{diag}(\delta^{(s)}) z_t + m^{(s)}, \tau_t^{(s)-1} I_V\right)]

为了完成贝叶斯框架，为参数定义了先验分布：
[P(\delta) = \prod_{s} N(\delta^{(s)}|0, I_K)]
[P(\tau) = \prod_{s} \prod_{t} \Gamma(\tau_t^{(s)}|a_{\tau_t^{(s)}}, b_{\tau_t^{(s)}})]
[P(m) = \prod_{s} N(m^{(s)}|0, \beta^{-1} I_V)]
[P(Z) = \prod_{t} N(z_t|0, I_K)]

$W$ 和 $\alpha$ 的先验选择决定了模型的修剪和稀疏能力。如果选择公式（1）和（2）的先验，则得到VB - PARAFAC，它使用ARD来识别活跃和不活跃的分量。如果选择公式（3）和（4）的先验，则得到稀疏概率（SP - ）PARAFAC，它试图识别维度（模型阶数）和活跃体素（特征）。

VB - 和SP - PCA与提出的VB - 和SP - PARAFAC模型的一般复杂度不同。值得注意的是，与PARAFAC模型相比，具有时间连接数据的PCA模型对于潜在变量 $z$ 有 $KT(S - 1)$ 更多的自由度，对于具有大 $T$ 和 $K$ 的数据集，这可能是过拟合的潜在来源，而PARAFAC对于 $\delta$ 有 $KS$ 更多的参数，对于 $m$ 有 $V(S - 1)$ 更多的参数。

4. 变分推理

为了最大化似然函数，必须估计边际分布 $P(X) = \int P(X, \theta) d\theta$。然而，关于先验分布对该分布进行边缘化通常在解析上是难以处理的。

变分推理通过用另一个分布 $Q$（称为变分分布）近似所需的分布来解决这个问题。变分方法的挑战在于选择一个足够简单的分布 $Q$，使得对数边际似然可以通过一个可处理的下界 $L(Q)$ 近似，同时又足够灵活以使这个下界紧密。常见的选择是一个对模型参数进行因式分解的分布，即 $Q = \prod_{i} Q_i(\theta_i)$。

对于VB - PARAFAC模型，分布 $Q_i(\theta_i)$ 定义如下：
[Q(Z) = \prod_{t} N(\mu_{z_t}, \Sigma_{z_t})]
[Q(W) = \prod_{v} N(w_v|\mu_{w_v}, \Sigma_{w_v})]
[Q(\alpha) = \prod_{k} \Gamma(\alpha_k|\tilde{a} {\alpha_k}, \tilde{b} {\alpha_k})]
[Q(\delta) = \prod_{s} N(\delta^{(s)}|\mu_{\delta^{(s)}}, \Sigma_{\delta^{(s)}})]
[Q(m) = \prod_{s} N(m^{(s)}|\mu_{m^{(s)}}, \Sigma_{m^{(s)}})]
[Q(\tau) = \prod_{s} \prod_{t} \Gamma(\tau_t^{(s)}|\tilde{a} {\tau_t^{(s)}}, \tilde{s} {\tau_t^{(s)}})]

可以证明，对数边际概率等价于 $\ln P(X) = L(Q) + KL(Q||P)$，其中 $KL$ 是Kullback - Leibler散度。通过固定 $Q_j$ 并相对于其余的 $Q_{i\neq j}$ 最大化 $L(Q)$，我们得到一般表达式 $\ln Q_j^* = E_{i\neq j}[\ln p(X, \theta)] + const$，这最小化了KL散度。通过应用和归一化，$Q$ 的更新规则可以推导出来，如下所示：
[\mu_{\delta^{(s)}} = \Sigma_{\delta^{(s)}} \sum_{t} \tau_t^{(s)} \langle\text{diag}(z_t)\rangle W^T (x_t^{(s)} - \langle m^{(s)}\rangle)]
[\Sigma_{\delta^{(s)}} = \left(I_K + \sum_{t} \tau_t^{(s)} \langle z_t z_t^T\rangle \circ W^T W\right)^{-1}]
[\mu_{z_t} = \Sigma_{z_t} \sum_{s} \tau_t^{(s)} \text{diag}(\langle\delta^{(s)}\rangle) W^T (x_t^{(s)} - \langle m^{(s)}\rangle)]
[\Sigma_{z_t} = \left(I_K + \sum_{s} \tau_t^{(s)} \langle\delta^{(s)} \delta^{(s)T}\rangle \circ W^T W\right)^{-1}]
[\mu_{m^{(s)}} = \Sigma_{m^{(s)}} \sum_{t} \tau_t^{(s)} (x_t^{(s)} - \langle W\rangle \langle\text{diag}(\delta^{(s)})\rangle \langle z_t\rangle)]
[\Sigma_{m^{(s)}} = \left(\beta + \sum_{t} \langle\tau\rangle_t^{(s)}\right)^{-1} I_K]
[\mu_{w_v} = \Sigma_{w_v} \sum_{s} \sum_{t} \langle\text{diag}(\delta^{(s)})\rangle \langle z_t\rangle \tau_t^{(s)} (x_t^{(s)T} - \langle m^{(s)T}\rangle) v]
[\Sigma {w_v} = \left(\text{diag}(\alpha) + \sum_{s} \sum_{t} \tau_t^{(s)} \langle\delta^{(s)} \delta^{(s)T}\rangle \circ \langle z_t z_t^T\rangle\right)^{-1}]
[\tilde{a} {\alpha} = a {\alpha} + \frac{V}{2}, \quad \tilde{b} {\alpha_k} = b {\alpha_k} + \frac{\langle w_k^T w_k\rangle}{2}]
[\tilde{a} {\tau^{(s)}} = a {\tau^{(s)}} + \frac{V}{2}]
[\tilde{b} {\tau_t^{(s)}} = b {\tau_t^{(s)}} + \frac{1}{2} |x_t^{(s)}|^2 + \frac{1}{2} \langle |m^{(s)}|^2\rangle - x_t^{(s)T} \langle m^{(s)}\rangle - x_t^{(s)T} \langle W\rangle \langle\text{diag}(\delta^{(s)})\rangle \langle z_t\rangle + \frac{1}{2} \text{Tr}\left(\langle W^T W\rangle \circ \langle\delta^{(s)} \delta^{(s)T}\rangle \langle z_t z_t^T\rangle\right) + \langle m^{(s)T}\rangle \langle W\rangle \langle\text{diag}(\delta^{(s)})\rangle \langle z_t\rangle]

通过将 $W$ 的先验更改为杰弗里先验（公式4），我们得到SP - PARAFAC的解。除了 $Q(\alpha)$ 现在按元素分解，即 $Q(\alpha) = \prod_{v} \prod_{k} \Gamma(\alpha_{v,k}|\tilde{a} {\alpha {v,k}}, \tilde{b} {\alpha {v,k}})$ 之外，$Q_i$ 分布保持不变。这些变化反映在 $W$ 和 $\alpha$ 的更新规则中：
[\Sigma_{w_v} = \left(\sum_{s} \sum_{t} \tau_t^{(s)} \langle\delta^{(s)} \delta^{(s)T}\rangle \circ \langle z_t z_t^T\rangle + \text{diag}(\langle\alpha_{v,.}\rangle)\right)^{-1}]
[\tilde{a} {\alpha {v,k}} = \frac{1}{2}, \quad \tilde{b} {\alpha {v,k}} = \frac{\langle w_{v,k}^2\rangle}{2}]

可以很容易地推导出下界 $L(Q)$ 来监测算法的收敛性。虽然收敛是有保证的，但解不一定能达到全局最大值。通常通过多次运行算法并保留具有最大下界的解来解决这个问题。

5. 结果

5.1 模拟

通过模拟比较了VB - PCA、SP - PCA、PARAFAC和提出的VB - PARAFAC和SP - PARAFAC。创建了一个三维数据集，首先创建三个长度为 $V = 10$ 的稀疏向量并连接形成矩阵 $W$ 作为地面真值。然后从分布 $N(0, I)$ 中采样长度为 $T = 100$ 的随机潜在变量 $z_t$、随机均值 $m^{(s)}$ 和 $S = 5$ 个受试者的随机缩放因子 $\delta^{(s)}$，并根据公式（5）进行线性混合。加性噪声从 $N(0, \tau_t^{(s)}I)$ 中采样，精度 $\tau_t^{(s)}$ 从 $\Gamma(1, 1)$ 中采样。

对于PCA算法，数据沿时间维度连接，对于PARAFAC模型，数据作为三维张量处理。所有模型都以不同的精度识别了地面真值分量并修剪了非信息性分量。与PARAFAC模型相比，两个PCA模型更容易受到噪声的干扰，并且无法考虑受试者特定的均值。此外，SP - PCA和SP - PARAFAC在修剪单个非相关权重方面比它们的VB版本表现更好，尽管这在RV系数中没有体现。这些结果表明，当数据中存在时间和受试者变化的噪声以及受试者特定的均值时，VB - 和SP - PARAFAC模型受益于PARAFAC结构，并且全稀疏先验在从数据中识别稀疏结构方面比在整个分量上诱导稀疏性的先验更有效。

5.2 正电子发射断层扫描

对PET数据进行了探索性数据分解。PET实验的同步性质、与放射性同位素衰变相关的时间依赖噪声以及由于注射剂量和体重变化导致的时间 - 活性曲线（TAC）的受试者特定缩放，使得多向模型理论上非常适合PET数据的分析。

从Cimbi数据库获取了4名健康受试者的放射性配体 $[^{11}C]CUMI - 101$ 的PET数据和匹配的T1加权结构磁共振图像（MRI）。动态PET图像（34帧）在HRRT扫描仪上采集，结构图像在3T西门子Trio扫描仪上采集。结构图像在FreeSurfer中处理，PET图像与结构MRI配准，转移到公共空间（MNI152）并使用5mm FWHM 3D高斯核进行平滑。最后，数据沿受试者连接形成三维张量，并使用初始维度为30的VB - PARAFAC和SP - PARAFAC进行处理。

VB - PCA、VB - PARAFAC和SP - PARAFAC分别识别了7、24和19个非零分量，SP - PCA没有修剪任何分量。所有五种方法仅共同识别了两个分量。第一个分量加载在整个大脑空间上，第二个分量加载在高和低结合区域上。其他分量在整个大脑中显示出随机模式，用于调整数据中的小随机变化。没有分量唯一地加载在特定的大脑区域上，这表明所有大脑体素共有的一组基函数是更合适的模型。

虽然VB - 和SP - PARAFAC算法修剪了非信息性分量，但SP - PARAFAC模型引入了看似随机的空间稀疏性，使得分量难以解释。这表明全稀疏模型可能并不总是适合识别全局模式。

另一个有趣的应用是对PET成像数据进行去噪。比较了不同算法的近似结果，发现所有大脑区域的结果惊人地相似。这表明我们的约束PARAFAC模型能够识别与更灵活的模型（如VB和SP - PCA）相当的噪声结构，并且在近似基础数据方面表现同样出色。

模型	识别的非零分量数量
VB - PCA	7
VB - PARAFAC	24
SP - PARAFAC	19
SP - PCA	无修剪

下面是变分推理中部分更新规则的流程：

graph LR
    A[初始化参数] --> B[计算Q(Z)]
    B --> C[计算Q(W)]
    C --> D[计算Q(α)]
    D --> E[计算Q(δ)]
    E --> F[计算Q(m)]
    F --> G[计算Q(τ)]
    G --> H{是否收敛}
    H -- 否 --> B
    H -- 是 --> I[输出结果]

综上所述，提出的VB - PARAFAC和SP - PARAFAC模型在处理PET和任务fMRI数据方面具有一定的优势，能够在考虑数据内在结构的同时，对噪声和稀疏性进行有效建模。但在实际应用中，需要根据具体的数据特点和分析目的选择合适的模型。

用于PET和任务fMRI数据建模的稀疏概率并行因子分析

5.3 功能磁共振成像

在对功能磁共振成像（fMRI）数据的分析中，我们同样应用了上述的模型，旨在识别与任务相关的脑活动成分。fMRI数据具有时间序列和多受试者的特点，这使得多向分解方法能够充分发挥其优势。

我们获取了一组基于任务的fMRI数据，对其进行了预处理，包括空间归一化、平滑等操作，以确保数据的一致性和可比性。然后，使用VB - PARAFAC和SP - PARAFAC模型对处理后的数据进行分析，初始维度设置为30。

在分析结果中，两个模型都成功地识别出了与任务相关的成分。这些成分在大脑的特定区域表现出了显著的激活，与预期的任务相关脑区相符合。然而，模型在处理重叠激活时遇到了困难。由于大脑活动的复杂性，不同的任务可能会导致部分脑区的重叠激活，而当前的模型无法有效地将这些重叠的激活进行分离。

这一结果提示我们，虽然PARAFAC模型在处理多向数据方面具有一定的优势，但在处理复杂的大脑活动模式时，仍需要进一步的改进和优化。例如，可以考虑引入更复杂的先验信息，或者结合其他的分析方法，以提高模型对重叠激活的分辨能力。

模型	任务相关成分识别	重叠激活分离能力
VB - PARAFAC	成功	较差
SP - PARAFAC	成功	较差

6. 模型性能评估

为了全面评估VB - PARAFAC和SP - PARAFAC模型的性能，我们从多个方面进行了分析。

6.1 稳定性评估

使用RV系数来评估模型解的稳定性。RV系数是皮尔逊相关系数的多变量推广，用于测量子空间的重叠程度，并且对旋转和平移具有不变性。通过计算不同初始化条件下估计的 $W$ 矩阵之间的RV系数，我们发现模型在多次运行中具有较好的稳定性。具体来说，平均RV系数较高，表明不同初始化下的解具有较高的一致性。

同时，我们还计算了平均相关性。首先，对于每对估计的 $W$ 矩阵，根据相关性匹配成分并计算平均相关性，然后对所有唯一对的平均相关性进行计算。结果显示，平均相关性也处于较高水平，进一步证明了模型解的稳定性。

6.2 近似能力评估

在PET和fMRI数据的分析中，我们比较了不同模型对数据的近似能力。通过观察数据的近似结果，我们发现VB - PARAFAC和SP - PARAFAC模型与更灵活的模型（如VB和SP - PCA）在近似效果上表现相似。这表明我们的约束PARAFAC模型能够有效地捕捉数据的主要特征，并且在处理噪声方面具有较好的性能。

6.3 稀疏性评估

对于SP - PARAFAC模型，我们重点评估了其稀疏性。通过观察模型识别的成分，我们发现该模型能够有效地引入空间稀疏性，修剪掉非信息性的分量和权重。然而，在某些情况下，如PET数据的分析中，过度的稀疏性可能会导致成分难以解释，这提示我们在应用稀疏模型时需要谨慎权衡稀疏性和可解释性之间的关系。

7. 总结与展望

本文提出了一种基于变分贝叶斯推理的稀疏概率并行因子分析（SP - PARAFAC）模型，以及其非稀疏版本（VB - PARAFAC），并将其应用于PET和任务fMRI数据的建模。通过模拟和实际数据的分析，我们验证了模型的有效性和性能。

在模拟实验中，VB - PARAFAC和SP - PARAFAC模型在处理具有时间和受试者变化噪声以及受试者特定均值的数据时表现出了优势，能够准确地识别地面真值成分并修剪非信息性分量。在PET数据的分析中，模型能够有效地识别出与大脑活动相关的成分，并且在去噪方面表现出色。在fMRI数据的分析中，虽然模型成功地识别出了任务相关成分，但在处理重叠激活方面仍存在不足。

未来的研究方向可以从以下几个方面展开：
1. 模型改进 ：进一步优化模型结构，提高模型对复杂数据模式的处理能力，特别是在处理重叠激活和噪声方面。例如，可以引入更复杂的先验分布，或者结合深度学习方法来增强模型的表达能力。
2. 多模态数据融合 ：考虑将PET和fMRI等多模态数据进行融合分析，以获取更全面的大脑信息。多模态数据融合可以充分发挥不同模态数据的优势，提高对大脑功能和结构的理解。
3. 临床应用 ：将模型应用于临床研究，如疾病诊断、治疗效果评估等。通过对大量临床数据的分析，验证模型在临床实践中的有效性和实用性。

下面是一个总结模型应用流程的流程图：

graph LR
    A[数据获取] --> B[数据预处理]
    B --> C{选择模型}
    C -- VB - PARAFAC --> D[模型训练]
    C -- SP - PARAFAC --> D
    D --> E[结果分析]
    E --> F{结果评估}
    F -- 不满意 --> C
    F -- 满意 --> G[应用与决策]

总之，本文提出的模型为PET和任务fMRI数据的分析提供了一种有效的方法，但仍有许多方面需要进一步的研究和改进。通过不断的探索和创新，我们有望开发出更强大的数据分析工具，为神经科学研究和临床应用提供更有力的支持。