PyMVPA_manual 4.8 event related 设计的数据分析 笔记

4.8时间相关的实验设计

与前几节不同的是，事件相关设计的数据分析时 将 time 时间这一因素考虑进去了

这里介绍两个方法：

1）实验条件的建模和模型参数估计的分析

2）时间空间数据样本的提取。

4.8.1 事件相关设计数据的预处理不是事件相关的

对于与事件相关的分析，大多数处理都是对从一组事件中派生出来的数据样本进行的。其余的数据可以被认为是无关的。然而，某些预处理只有在对完整的时间序列而不是对分段的事件样本执行时才有意义。一个例子是，通常需要对原始的连续时间序列进行去趋势处理。

在当前的形状中，我们的数据集由表示连续fMRI体积的样本组成。在这个阶段，我们可以很容易地进行线性去趋势。

>>> from mvpa2.tutorial_suite import *

>>> ds = get_raw_haxby2001_data(roi=(36,38,39,40))

>>> poly_detrend(ds, polyord=1, chunks_attr='chunks')

>>> orig_ds = ds.copy()

4.8.2设计说明

对于任何与事件相关的分析，我们都需要一些关于实验设计的信息：什么时候用什么刺激了多长时间(也许用什么强度)。在pymvpa中，这是通过编译事件定义列表来完成的。在许多情况下，事件由onset, duration和其他一些额外的属性组成的，例如刺激条件或记录的session number.

find_events() 这个函数的作用是  将事件相关的事件一个个提取出来，onset,duration 这两个属性是默认获取的，其他属性需要传递给这个函数，并且这个函数会根据所提供的属性值将属性值一样的几个事件合成一个事件，这就要求输入给这个函数的属性需要真的将我们想要的事件区分开，又不至于将不同的事件合并。所以传入的参数很重要。这个例子中传入的属性只有 targets和chunks  实际上如果两个target一样的事件就会合在一起，所以记住传入真正能分开的参数吧

>>> events = find_events(targets=ds.sa.targets, chunks=ds.sa.chunks)

>>> print len(events)

204

>>> for e in events[:4]:

... print e

{'chunks': 0, 'duration': 6, 'onset': 0, 'targets': 'rest'}

{'chunks': 0, 'duration': 9, 'onset': 6, 'targets': 'scissors'}

{'chunks': 0, 'duration': 6, 'onset': 15, 'targets': 'rest'}

{'chunks': 0, 'duration': 9, 'onset': 21, 'targets': 'face'}

我们不仅向函数提供targets，而且还向提供chunks信息，因为我们不希望有跨越多个记录sessions的事件。 find_events() 依次解析所有提供的属性，并在任何属性中的值发生变化时记录事件。生成的事件定义是包含以下内容的字典：

1.Onset 事件的开始作为序列中的索引(在本例中，这是 volume 的id)

2.Duration “序列元素数目”中事件的持续时间(即volume数)

3.该事件的属性组合，即特定位置上所有给定属性的实际值。

现在让我们把自己限制在face 和house的刺激block上。我们可以很容易地过滤掉所有其他事件。

>>> events = [ev for ev in events if ev['targets'] in ['house', 'face']]

>>> print len(events)

24

>>> for e in events[:4]:

... print e

{'chunks': 0, 'duration': 9, 'onset': 21, 'targets': 'face'}

{'chunks': 0, 'duration': 9, 'onset': 63, 'targets': 'house'}

{'chunks': 1, 'duration': 9, 'onset': 127, 'targets': 'face'}

{'chunks': 1, 'duration': 9, 'onset': 213, 'targets': 'house'}

4.8.3响应建模

当我们必须处理数据时，多个并发信号在时间上是重叠的，例如在快速事件相关的fMRI研究中，通常有必要对数据进行适当的模型拟合，并对模型参数估计进行分析，而不是对原始数据进行分析。

pymvpa可以利用nipy的GLM建模功能。它期望刺激事件的信息作为实际的时间戳，而不是数据样本索引，因此我们必须转换我们的事件列表。

>>> # temporal distance between samples/volume is the volume repetition time

>>> TR = np.median(np.diff(ds.sa.time_coords))

>>> # convert onsets and durations into timestamps

>>> for ev in events:

... ev['onset'] = (ev['onset'] * TR)

... ev['duration'] = ev['duration'] * TR

现在，我们可以对所有相关刺激条件下的血流动力学反应模型进行拟合。函数eventreated_dataset()为我们执行所有操作。对于给定的输入数据集，我们需要提供事件列表、带有每个示例时间戳的属性的名称以及我们希望建模的条件的信息。后者被指定为condition_attr参数。这可以是一个属性名，在这种情况下，所有唯一的值都将用作条件。它也可以是多个属性名称的序列，属性的所有唯一值的组合都将用作条件。在下面的示例('targets','chunks')中，我们希望为每个示例数据集(chunks)的每个刺激条件(targets)建立一个单独的模型。

>>> evds = fit_event_hrf_model(ds,

... events,

... time_attr='time_coords',

... condition_attr=('targets', 'chunks'))

>>> print len(evds)

24

这将为每个块chunks的每个目标值target value生成一个参数估计样本。

在进行分类分析之前，我们仍然需要对每个特征进行规范化(此时每个体素的GLM参数估计)。

>>> zscore(evds, chunks_attr=None)

其余的是直截了当的：我们用选择的分类器建立一个交叉验证分析，并运行它：

>>> clf = kNN(k=1, dfx=one_minus_correlation, voting='majority')

>>> cv = CrossValidation(clf, NFoldPartitioner(attr='chunks'))

>>> cv_glm = cv(evds)

>>> print '%.2f' % np.mean(cv_glm)

0.04

还不错!让我们将其与一种更简单的方法进行比较，这种方法也适用于像这样的分块设计实验。

>>> zscore(ds, param_est=('targets', ['rest']))

>>> avgds = ds.get_mapped(mean_group_sample(['targets', 'chunks']))

>>> avgds = avgds[np.array([t in ['face', 'house'] for t in avgds.sa.targets])]

我们对关于其余条件的所有体素进行规范化。这为所有刺激条件产生了一些粗略的“激活”分数。随后，我们在每次运行中平均每个条件的所有样本。这将产生与GLM建模相同大小的数据集。我们可以重新使用交叉验证设置。

>>> cv_avg = cv(avgds)

>>> print '%.2f' % np.mean(cv_avg)

0.04

也不错。然而，值得重复的是，这种简单的平均采样方法限于具有所有感兴趣信号的清晰的时间间隔的块设计，而HRF建模更适合于具有快速刺激交替的实验。

4.8.4从时间序列到时空样本

现在我们想尝试一些不同的东西。与其将所有时间信息压缩成一个单一的模型参数估计，我们还可以考虑整个时空信号在我们感兴趣的区域和整个刺激块的持续时间。换句话说，我们可以进行灵敏度分析(见分类模型参数-灵敏度分析)，揭示分类相关信息的时空分布。

在开始进行事件提取之前，我们希望对每个特征进行规范化(即此时是一个体素)。在这种情况下，我们再次使用实验“静止条件”的平均值和标准差来对它们进行Z评分，得到的值可能被解释为“激活分数”。

>>> zscore(ds, chunks_attr='chunks', param_est=('targets', 'rest'))

对于这种分析，我们不需要将事件开始信息转换为时间戳，而是可以对样本索引进行操作，因此我们再次从原始事件列表开始。

>>> events = find_events(targets=ds.sa.targets, chunks=ds.sa.chunks)

>>> events = [ev for ev in events if ev['targets'] in ['house', 'face']]

我们所有的事件都是相同的长度，连续9个fMRI卷。稍后，我们想查看从前到刺激阻滞后的时间敏感性曲线，因此我们应该延长事件的持续时间。

>>> event_duration = 13

>>> for ev in events:

... ev['onset'] -= 2

... ev['duration'] = event_duration

下一步也是最重要的一步是将原始时间序列数据集实际分割为与事件相关的样本。pymvpa提供s eventrelated_dataset()作为函数来执行此转换。让我们这样做，它只需要原始数据集和我们的事件列表。

>>> evds = eventrelated_dataset(ds, events=events)

>>> len(evds) == len(events)

True

>>> evds.nfeatures == ds.nfeatures * event_duration

True

此时，需要查看数据集的映射器--尤其是在转换为事件期间添加的链映射器中的最后两个项目。

>>> print evds.a.mapper[-2:]

我们分析的其余部分都照常进行，很快就完成了。 我们希望对svm分类器进行交叉验证分析.我们主要不是对它的性能感兴趣，而是对它分配给特性的权重感兴趣。请记住，现在每个特性都是时点体素，因此我们有机会查看数据中与分类相关的信息的时空轮廓。然而，我们将允许计算一个混淆矩阵，这样我们就可以保证分类器的性能相当好，因为只有一个泛化模型值得检查，否则它就会变得不合适，分配的权重也就没有意义了。

>>> sclf = SplitClassifier(LinearCSVMC(),

... enable_ca=['stats'])

>>> sensana = sclf.get_sensitivity_analyzer()

>>> sens = sensana(evds)

4.8.5绘图示例

我们已经使用一些MRI查看器的应用检查了敏感性的时空剖面，但我们也可以在这里组装一个信息丰富的图形。让我们组成一个图表，显示原始的周期刺激时间序列，归一化的效果，以及经过训练的svm分类器的相应的灵敏度剖面。我们将对两个例子做这件事，我们可以通过检查完整的地图得到它们的坐标。

>>> example_voxels = [(28,25,25), (28,23,25)]

绘图将由流行的matplotlib包完成。

首先，在初始去趋势后绘制原始信号。为此，我们将同样的时间序列分割应用于原始的去趋势数据集，并为我们的示例体素绘制所有脸和房子事件的平均信号。下面的代码将使用matplotlib的pylab接口(导入为pl)创建绘图。如果您熟悉MATLAB的绘图工具，这应该不难理解。