Numerically Perturbed Structural Connectomes from 100 individuals in the NKI Rockland Dataset

This dataset contains the derived connectomes, discriminability scores, and classification performance for structural connectomes estimated from a subset of the Nathan Kline Institute Rockland Sample dataset, and is associated with an upcoming manuscript entitled: <em>Numerical Instabilities in Analytical Pipelines Compromise the Reliability of Network Neuroscience</em>. The associated code for this project is publicly available at: https://github.com/gkpapers/2020ImpactOfInstability. For any questions, please contact Gregory Kiar ([e-mail]) or Tristan Glatard ([e-mail]). Below is a table of contents describing the contents of this dataset, which is followed by an excerpt from the manuscript pertaining to the contained data. impactofinstability_connect_dset25x2x2x20_inputs.h5 : Connectomes derived from 25 subjects, 2 sessions, 2 subsamples, and 20 MCA simulations with input perturbations. impactofinstability_connect_dset25x2x2x20_pipeline.h5 : Connectomes derived from 25 subjects, 2 sessions, 2 subsamples, and 20 MCA simulations with pipeline perturbations. impactofinstability_discrim_dset25x2x2x20_both.csv : Discriminability scores for each grouping of the 25x2x2x20 dataset. impactofinstability_connect+feature_dset100x1x1x20_both.h5 : Connectomes and features derived from 100 subjects, 1 sessions, 1 subsamples, and 20 MCA simulations with both perturbation types. impactofinstability_classif_dset100x1x1x20_both.h5 : Classification performance results for the BMI classification task on the 100x1x1x20 dataset. <strong>Dataset</strong><br> The Nathan Kline Institute Rockland Sample (NKI-RS) dataset [1] contains high-fidelity imaging and phenotypic data from over 1,000 individuals spread across the lifespan. A subset of this dataset was chosen for each experiment to both match sample sizes presented in the original analyses and to minimize the computational burden of performing MCA. The selected subset comprises 100 individuals ranging in age from 6 – 79 with a mean of 36.8 (original: 6 – 81, mean 37.8), 60% female (original: 60%), with 52% having a BMI over 25 (original: 54%). Each selected individual had at least a single session of both structural T1-weighted (MPRAGE) and diffusion-weighted (DWI) MR imaging data. DWI data was acquired with 137 diffusion directions; more information regarding the acquisition of this dataset can be found in the NKI-RS data release [1]. In addition to the 100 sessions mentioned above, 25 individuals had a second session to be used in a test-retest analysis. Two additional copies of the data for these individuals were generated, including only the odd or even diffusion directions (64 + 9 B0 volumes = 73 in either case). This allows an extra level of stability evaluation to be performed between the levels of MCA and session-level variation. In total, the dataset is composed of 100 diffusion-downsampled sessions of data originating from 50 acquisitions and 25 individuals for in depth stability analysis, and an additional 100 sessions of full-resolution data from 100 individuals for subsequent analyses. <strong>Processing</strong><br> The dataset was preprocessed using a standard FSL [2] workflow consisting of eddy-current correction and alignment. The MNI152 atlas was aligned to each session of data, and the resulting transformation was applied to the DKT parcellation [3]. Downsampling the diffusion data took place after preprocessing was performed on full-resolution sessions, ensuring that an additional confound was not introduced in this process when comparing between downsampled sessions. The preprocessing described here was performed once without MCA, and thus is not being evaluated. Structural connectomes were generated from preprocessed data using two canonical pipelines from Dipy [4]: deterministic and probabilistic. In the deterministic pipeline, a constant solid angle model was used to estimate tensors at each voxel and streamlines were then generated using the EuDX algorithm [5]. In the probabilistic pipeline, a constrained spherical deconvolution model was fit at each voxel and streamlines were generated by iteratively sampling the resulting fiber orientation distributions. In both cases tracking occurred with 8 seeds per 3D voxel and edges were added to the graph based on the location of terminal nodes with weight determined by fiber count. <strong>Perturbations</strong><br> All connectomes were generated with one reference execution where no perturbation was introduced in the processing. For all other executions, all floating point operations were instrumented with Monte Carlo Arithmetic (MCA) [6] through Verificarlo [7]. MCA simulates the distribution of errors implicit to all instrumented floating point operations (flop). MCA can be introduced in two places for each flop: before or after evaluation. Performing MCA on the inputs of an operation limits its precision, while performing MCA on the output of an operation highlights round-off errors that may be introduced. The former is referred to as Precision Bounding (PB) and the latter is called Random Rounding (RR). Using MCA, the execution of a pipeline may be performed many times to produce a distribution of results. Studying the distribution of these results can then lead to insights on the stability of the instrumented tools or functions. To this end, a complete software stack was instrumented with MCA and is made available on GitHub through https://github.com/gkiar/fuzzy. Both the RR and PB variants of MCA were used independently for all experiments. As was presented in [8], both the degree of instrumentation (i.e. number of affected libraries) and the perturbation mode have an effect on the distribution of observed results. For this work, the RR-MCA was applied across the bulk of the relevant libraries and is referred to as Pipeline Perturbation. In this case the bulk of numerical operations were affected by MCA. Conversely, the case in which PB-MCA was applied across the operations in a small subset of libraries is here referred to as Input Perturbation. In this case, the inputs to operations within the instrumented libraries (namely, Python and Cython) were perturbed, resulting in less frequent, data-centric perturbations. Alongside the stated theoretical differences, Input Perturbation is considerably less computationally expensive than Pipeline Perturbation. All perturbations were targeted the least-significant-bit for all data (t=24and t=53in float32 and float64, respectively [7]). Simulations were performed between 10 and 20 times for each pipeline execution, depending on the experiment. A detailed motivation for the number of simulations can be found in [9]. <strong>Evaluation</strong><br> The magnitude and importance of instabilities in pipelines can be considered at a number of analytical levels, namely: the induced variability of derivatives directly, the resulting downstream impact on summary statistics or features, or the ultimate change in analyses or findings. We explore the nature and severity of instabilities through each of these lenses. Unless otherwise stated, all p-values were computed using Wilcoxon signed-rank tests. <strong>Direct Evaluation of the Graphs</strong><br> The differences between simulated graphs was measured directly through both a direct variance quantification and a comparison to other sources of variance such as individual- and session-level differences. Quantification of Variability – Graphs, in the form of adjacency matrices, were compared to one another using three metrics: normalized percent deviation, Pearson correlation, and edgewise significant digits. The normalized percent deviation measure, defined in [8], scales the norm of the difference between a simulated graph and the reference execution (that without intentional perturbation) with respect to the norm of the reference graph. The purpose of this comparison is to provide insight on the scale of differences in observed graphs relative to the original signal intensity. A Pearson correlation coefficient was computed in complement to normalized percent deviation to identify the consistency of structure and not just intensity between observed graphs. Finally, the estimated number of significant digits for each edge in the graph was computed. The upper bound on significant digits is 15.7 for 64-bit floating point data. The percent deviation, correlation, and number of significant digits were each calculated within a single session of data, thereby removing any subject- and session-effects and providing a direct measure of the tool-introduced variability across perturbations. A distribution was formed by aggregating these individual results. Class-based Variability Evaluation – To gain a concrete understanding of the significance of observed variations we explore the separability of our results with respect to understood sources of variability, such as subject-, session-, and pipeline-level effects. This can be probed through Discriminability [10], a technique similar to ICC which relies on the mean of a ranked distribution of distances between observations belonging to a defined set of classes. Discriminability can then be interpreted as the probability that an observation belonging to a given class will be more similar to other observations within that class than observations of a different class. It is a measure of reproducibility, and is discussed in detail in [10]. This definition allows for the exploration of deviations across arbitrarily defined classes which in practice can be any of those listed above. We combine this statistic with permutation testing to test hypotheses on whether differences between classes are statistically significant in each of these settings. With this in mind, three hypotheses were defined. For each setting, we state the alternate hypotheses, the variable(s) which will be used to determine class membership, and the remaining variables which may be sampled when obtaining multiple observations. Each hypothesis was tested independently for each pipeline and perturbation mode, and in every case where it is possible the hypotheses were tested using the reference executions alongside using MCA. Individual Variation<br> HA: Individuals are distinct from one another.<br> Class definition: Subject ID.<br> Experiments: Session (1 subsample), Direction (1 subsample), MCA (1 subsample, 1 session). Session Variation<br> HA: Sessions within an individual are distinct.<br> Class definition: Session ID | Subject ID.<br> Experiments: Subsample, MCA (1 subsample). Subsample Variation<br> HA: Direction subsamples within an acquisition are distinct.<br> Class definition: Subsample | Subject ID, Session ID.<br> Experiments: MCA. As a result, we tested 3 hypotheses across 6 MCA experiments and 3 reference experiments on 2 pipelines and 2 perturbation modes, resulting in a total of 30 distinct tests. <strong> Evaluating Graph-Theoretical Metrics</strong><br> While connectomes may be used directly for some analyses, it is common practice to summarize them with structural measures, which can then be used as lower-dimensional proxies of connectivity in so-called graph-theoretical studies [11]. We explored the stability of several commonly-used univariate (graphwise) and multivariate (nodewise or edgewise) features. The features computed and subsequent methods for comparison in this section were selected to closely match those computed in [12]. Univariate Differences – For each univariate statistic (edge count, mean clustering coefficient, global efficiency, modularity, assortativity, and mean path length) a distribution of values across all perturbations within subjects was observed. A Z-score was computed for each sample with respect to the distribution of feature values within an individual, and the proportion of "classically significant" Z-scores, i.e. corresponding to p &lt; 0.05, was reported and aggregated across all subjects. The number of significant digits contained within an estimate derived from a single subject were calculated and aggregated. Multivariate Differences – In the case of both nodewise (degree distribution, clustering coefficient, betweenness centrality) and edgewise (weight distribution, connection length) features, the cumulative density functions of their distributions were evaluated over a fixed range and subsequently aggregated across individuals. The number of significant digits for each moment of these distributions (sum, mean, variance, skew, and kurtosis) were calculated across observations within a sample and aggregated. <strong> Evaluating A Complete Analysis</strong><br> Though each of the above approaches explores the instability of derived connectomes and their features, many modern studies employ modeling or machine-learning approaches, for instance to learn brain-behavior relationships or identify differences across groups. We carried out one such study and explored the instability of its results with respect to the upstream variability of connectomes characterized in the previous sections. We performed the modeling task with a single sampled connectome per individual and repeated this sampling and modelling 20 times. We report the model performance for each sampling of the dataset and summarize its variance. BMI Classification – Structural changes have been linked to obesity in adolescents and adults [13]. We classified normal-weight and overweight individuals from their structural networks (using for overweight a cutoff of BMI &gt; 25 [14]). We reduced the dimensionality of the connectomes through principal component analysis (PCA), and provided the first N-components to a logistic regression classifier for predicting BMI class membership, similar to methods shown in [14], [15]. The number of components was selected as the minimum set which explained &gt; 90% of the variance when averaged across the training set for each fold within the cross validation of the original graphs; this resulted in a feature of 20 components. We trained the model using k-fold cross validation, with k = 2, 5, 10, and N (equivalent to leave-one-out; LOO).本数据集包含源自内森·克莱因研究所洛克兰样本（Nathan Kline Institute Rockland Sample, NKI-RS）子集的结构连接组（structural connectome）的衍生连接组、可判别性分数与分类性能，关联一篇即将发表的论文，标题为《分析流水线中的数值不稳定性损害网络神经科学研究的可靠性》（Numerical Instabilities in Analytical Pipelines Compromise the Reliability of Network Neuroscience）。本项目的关联代码已公开部署于：https://github.com/gkpapers/2020ImpactOfInstability。如有任何疑问，请联系Gregory Kiar（邮箱：[e-mail]）或Tristan Glatard（邮箱：[e-mail]）。 以下为本数据集的内容目录，其后附有所包含数据相关的论文节选。 impactofinstability_connect_dset25x2x2x20_inputs.h5 : 源自25名受试者、2次扫描会话、2次子采样及20次输入扰动蒙特卡罗算术（Monte Carlo Arithmetic, MCA）模拟的连接组。 impactofinstability_connect_dset25x2x2x20_pipeline.h5 : 源自25名受试者、2次扫描会话、2次子采样及20次流水线扰动MCA模拟的连接组。 impactofinstability_discrim_dset25x2x2x20_both.csv : 25x2x2x20数据集各分组的可判别性分数。 impactofinstability_connect+feature_dset100x1x1x20_both.h5 : 源自100名受试者、1次扫描会话、1次子采样及20次双扰动类型MCA模拟的连接组与特征。 impactofinstability_classif_dset100x1x1x20_both.h5 : 针对100x1x1x20数据集的BMI分类任务的分类性能结果。 <strong>数据集概况</strong><br> 内森·克莱因研究所洛克兰样本（NKI-RS）数据集[1]包含覆盖全生命周期的1000余名个体的高保真影像与表型数据。本研究为每个实验选取该数据集的一个子集，既匹配原始分析中的样本量，又尽可能降低执行蒙特卡罗算术（Monte Carlo Arithmetic, MCA）的计算负担。所选子集包含100名年龄介于6至79岁的个体，平均年龄为36.8岁（原始数据集年龄范围为6至81岁，平均年龄37.8岁），其中女性占比60%（原始数据集占比60%），52%的个体体质指数（BMI）高于25（原始数据集占比54%）。每名入选个体至少拥有1次结构T1加权成像（MPRAGE）与弥散加权磁共振成像（DWI）数据。DWI数据采用137个弥散方向采集；有关该数据集采集的更多信息可参阅NKI-RS数据发布文档[1]。 除上述100次扫描会话外，另有25名个体提供了第二次扫描会话，用于重测分析。为这些个体额外生成了两份数据副本，仅分别包含奇数或偶数序号的弥散方向（两种情况均为64个弥散方向+9个B0图像帧，总计73帧）。这一设置可在MCA扰动与扫描会话变异层面之间额外增加一层稳定性评估。总体而言，本数据集包含两类数据：一类是源自50次采集、25名个体的100次弥散下采样扫描会话数据，用于深度稳定性分析；另一类是来自100名个体的100次全分辨率扫描会话数据，用于后续分析。 <strong>数据预处理流程</strong><br> 本数据集采用标准FSL[2]工作流进行预处理，包含涡流校正与配准步骤。将MNI152标准脑图谱配准至每一次扫描会话的数据，并将得到的变换矩阵应用于DKT脑区分割模板[3]。弥散数据的下采样操作在全分辨率会话完成预处理后执行，以确保在比较下采样会话时不会引入额外的混杂因素。本文所述的预处理流程仅在未引入MCA的情况下执行一次，因此不纳入本次评估。 我们采用Dipy[4]工具库中的两种标准流水线从预处理数据中生成结构连接组：确定性流水线与概率性流水线。在确定性流水线中，我们使用恒定立体角模型估计每个体素内的张量，随后通过EuDX算法[5]生成纤维束流线。在概率性流水线中，我们对每个体素拟合约束球面反卷积模型，并通过迭代采样得到的纤维取向分布生成纤维束流线。两种流水线均采用每个3D体素8个采样种子的方式进行追踪，并根据终端节点的位置为连接图添加边，边的权重由纤维数量决定。 <strong>扰动设置</strong><br> 所有连接组均通过一次无任何处理扰动的参考执行流程生成。其余所有执行流程中，所有浮点运算（flop）均通过Verificarlo[7]工具接入蒙特卡罗算术（MCA）[6]。MCA用于模拟所有接入的浮点运算中隐含的误差分布。对于每个浮点运算，MCA可在运算求值前或求值后两个位置引入扰动：在运算输入侧执行MCA会限制其计算精度，而在运算输出侧执行MCA则会突出可能引入的舍入误差。前者被称为精度边界（Precision Bounding, PB），后者被称为随机舍入（Random Rounding, RR）。 通过MCA，我们可多次执行流水线以生成结果分布。对这些结果分布的研究可帮助我们理解接入工具或函数的稳定性。为此，我们将完整软件栈接入MCA，并通过https://github.com/gkiar/fuzzy公开部署于GitHub。本研究所有实验均独立使用MCA的RR与PB两种变体。如文献[8]所述，接入程度（即受影响的库的数量）与扰动模式均会对观测结果的分布产生影响。 在本研究中，我们在绝大多数相关库中应用RR-MCA，称之为流水线扰动（Pipeline Perturbation），此时绝大多数数值运算均会受到MCA的影响。反之，我们在少量库的运算中应用PB-MCA，称之为输入扰动（Input Perturbation），此时仅接入库（即Python与Cython）内的运算输入会被扰动，从而产生频率更低、以数据为中心的扰动。除上述理论差异外，输入扰动的计算成本远低于流水线扰动。 所有扰动均针对所有数据的最低有效位（float32与float64格式分别对应t=24与t=53[7]）。根据实验需求，每个流水线执行流程的模拟次数为10至20次。有关模拟次数的详细设计动机可参阅文献[9]。 <strong>评估方法</strong><br> 流水线不稳定性的程度与重要性可从多个分析层面考量：一是直接观测衍生结果的诱导变异，二是下游汇总统计量或特征受到的影响，三是最终分析或研究发现的变化。我们通过这三个层面探究不稳定性的本质与严重程度。除非另有说明，所有p值均通过威尔科克森符号秩检验计算。 <strong>直接评估连接组图</strong><br> 我们通过直接方差量化以及与其他方差来源（如个体间与会话间差异）的比较，直接测量模拟连接组图之间的差异。 变异量化——以邻接矩阵形式存在的连接组图之间通过三项指标进行比较：归一化百分比偏差、皮尔逊相关系数与逐边有效数字。归一化百分比偏差指标[8]将模拟连接组与参考执行流程（无故意扰动）之间的差的范数，相对于参考连接组的范数进行缩放。该比较旨在帮助我们理解观测到的连接组差异相对于原始信号强度的尺度。我们同时计算皮尔逊相关系数，以补充归一化百分比偏差，从而识别观测连接组之间不仅是强度，更是结构层面的一致性。最后，我们计算图中每条边的估计有效数字数量。对于64位浮点数据，有效数字的上限为15.7。 我们在单次扫描会话数据内计算百分比偏差、相关系数与有效数字数量，从而消除个体与会话效应，直接测量不同扰动下工具引入的变异。我们通过聚合这些个体结果得到整体分布。 基于类别的变异评估——为具体理解观测变异的显著性，我们针对已知的变异来源（如个体、会话与流水线层面的效应）探究结果的可分离性。这可通过可判别性（Discriminability）[10]进行分析，该方法与组内相关系数（Intraclass Correlation Coefficient, ICC）类似，依赖于属于指定类别的观测值之间距离的排序分布均值。可判别性可被解释为：属于某一类别的观测值与该类别内其他观测值的相似度，高于与其他类别观测值的相似度的概率。该指标是可重复性的量化方式，详细讨论可参阅文献[10]。这一定义允许我们探索任意定义的类别之间的偏差，实际应用中可包括上述所有类别。我们将该统计量与置换检验结合，以检验各类别之间的差异在各设置下是否具有统计学显著性。 基于此，我们定义了三项假设。对于每个设置，我们将给出备择假设、用于确定类别的变量，以及获取多次观测值时可采样的其余变量。我们针对每种流水线与扰动模式独立检验每项假设，并在可行的情况下同时使用参考执行流程与MCA执行流程进行假设检验。 个体变化<br> HA: 个体之间存在显著差异。<br> 类别定义：受试者ID。<br> 实验设置：会话（1次子采样）、弥散方向（1次子采样）、MCA（1次子采样、1次会话）。 会话变化<br> HA: 同一受试者的不同扫描会话之间存在显著差异。<br> 类别定义：会话ID | 受试者ID。<br> 实验设置：子采样、MCA（1次子采样）。 子采样变化<br> HA: 同一次采集中的不同弥散方向子采样之间存在显著差异。<br> 类别定义：子采样 | 受试者ID、会话ID。<br> 实验设置：MCA。 综上，我们针对2种流水线、2种扰动模式下的6次MCA实验与3次参考实验，共检验了3项假设，总计完成30次独立检验。 <strong>图论指标评估</strong><br> 尽管连接组可直接用于部分分析，但常见做法是通过结构指标对其进行汇总，随后在图论研究中作为低维连接性代理[11]。我们探究了数种常用单变量（全图级）与多变量（节点级或边级）特征的稳定性。本节计算的特征与后续比较方法均与文献[12]保持高度一致。 单变量差异——对于每项单变量统计量（边数量、平均聚类系数、全局效率、模块化度、同配性与平均路径长度），我们观测了受试者内所有扰动下的取值分布。针对每个样本，我们根据个体内特征值分布计算Z分数，并报告并聚合所有受试者中“经典显著”的Z分数占比（即对应p < 0.05）。我们同时计算并聚合单个受试者的估计值中包含的有效数字数量。 多变量差异——对于节点级（度分布、聚类系数、介数中心性）与边级（权重分布、连接长度）特征，我们在固定范围内评估其分布的累积密度函数，并随后在所有受试者间聚合。我们计算样本内观测值的各分布矩（和、均值、方差、偏度与峰度）的有效数字数量，并进行聚合。 <strong>完整分析流程评估</strong><br> 尽管上述每种方法均探究了衍生连接组及其特征的不稳定性，但许多现代研究采用建模或机器学习方法，例如探索脑-行为关联或识别组间差异。我们开展了一项此类研究，并针对前文所述的连接组上游变异，探究其结果的不稳定性。我们为每名个体采样一条连接组以执行建模任务，并将该采样与建模流程重复20次。我们报告数据集每次采样的模型性能，并汇总其变异。 BMI分类——结构改变与青少年及成人的肥胖相关[13]。我们根据个体的结构网络对正常体重与超重个体进行分类（超重的BMI cutoff值为>25[14]）。我们通过主成分分析（PCA）降低连接组的维度，并将前N个主成分输入逻辑回归分类器以预测BMI类别，该方法与文献[14][15]中的方法类似。主成分的数量选择为：在原始连接组的交叉验证中，每个折的训练集上平均解释>90%方差的最小集合，最终得到20个主成分的特征集。我们采用k折交叉验证训练模型，k分别取2、5、10与N（等价于留一法，LOO）。