Burt Variogram Matching

基本思想

Fig. 1
如图所示,Burt 的方法主要分为两部分,首先生成一个随机的map,再通过平滑等操作增强其空间自相关性(Spatial Autocorrelation,SA),得到保留SA的map。其中的核心思想是让permutation得到的map和target map之间的variogram趋势一致。

Variogram

  • 对于特定的map,variogram是衡量不同距离下点对值的变异大小的曲线,反映了map中的SA整体情况。
    • 距离:点和点之间的距离,在surface上可以使用geodesic distances(根据surface走形等得到),在volume上可以使用欧氏距离。
  • γ(h)=i1Nji+1NWijViji1Nji+1NWij\gamma(h)=\frac{\sum_{i-1}^N\sum_{j-i+1}^{N}W_{ij}V_{ij}}{\sum_{i-1}^N\sum_{j-i+1}^{N}W_{ij}}
    • WijW_{ij}为随距离衰减的指数函数,在两点间距离等于hh时达到最大值。Wij=e(2.68(hdij))22b2W_{ij}=e^{-\frac{(2.68(h-d_{ij}))^2}{2b^2}}
    • VijV_{ij}为两点间值的差异。Vij=12(xixj)2V_{ij}=\frac{1}{2}(x_i-x_j)^2

步骤

Step 1. 打乱target map点序,生成random map (x0x_0)

Step 2. 平滑random map得到smooth map (xkx_k),计算variogram

  • xk,i=j=1kK(dij)x0,jj=1kK(dij)x_{k,i}=\frac{\sum_{j=1}^{k}K(d_{ij})x_{0,j}}{\sum_{j=1}^{k}K(d_{ij})}K(d)K(d)是一个有关距离的平滑核,让random map的每一点的值与周围点的值有关。

Step 3. 计算target map的variogram

Step 4. 对smooth map进行线性变换,作为surrogate map

  • 将smooth map的variogram线性变化到target map的variogram,γ[x]=βγ[xk]+α\gamma[x]=\beta \gamma[x_k]+\alpha
  • 选出变换后SSE最小的map xkx_{k*},将其变换得到的αk\alpha_{k*}βk\beta_{k*}用于变换xkx_{k*}的值,作为surrogate map。xs=βk12xk+αk12N(0,1){x_{s}}=|\beta_{k*}|^{\frac{1}{2}}x_{k*}+|\alpha_{k*}|^{\frac{1}{2}}N(0,1)

参数

  • knn:在距离矩阵中,对每一个点只选择前knn个邻近点的距离,其他点的距离置零。
  • kernel:第2步中使用的平滑核类型。
    • 一般选择指数核或者高斯核,如果选择uniform则不受到距离影响。
  • delta:第2步中使用的平滑核大小(包含百分之多少的点)。
  • ns:在进行第2、3步时,只选择ns个点建立variogram。
  • nh:第2、3步中建立variogram时的距离取值数量(variogram精度)。
  • pv:在fit variogram时只对齐前(距离短的)pv的点。
    • 如果值不是大块大块或者遍布全脑,可以把这个参数调小一点

Moran spectral randomization

基本思想

对距离矩阵进行特征根分解。将target map投射到每一个特征向量,并根据特征值进行缩放,得到SA。并让得到的randomization map的SA与target相同。

SA指标

Step 1. 距离矩阵的构建和分解

  • 对于N个voxel或vertex的map,可以计算得到NNN*N的距离矩阵。对于距离矩阵取倒数并将对角线置零得到矩阵W。
  • 对于W进行特征根分解,根据(Griffith 2000)的结论,一定会有一个特征值非常非常接近0 (具体为啥我也不知道...)。所以可以得到N-1个特征值和对应的特征向量。

Step 2. 向量投射与SA指标

  • 计算taget map x与每个特征向量的相关系数的平方rxvi2r^2_{xv_{i}}。这个系数相当于x在这个特征向量上的能量谱。
    • 这个概念和傅里叶变换中的能量谱非常接近,只不过傅里叶变换中得到的是每个频段的能量密度,这里得到的是每个特征向量上的能量密度。
    • rxv2=1r^2_{xv}=1,所有特征向量上的能量谱和等于1。
  • SA指标:Ix=rxvT2mI_x=r^2_{xv^T}m
    • m是缩放后的特征值,m=hλkm=h\lambda_{k},其中h是nijwij\frac{n}{\sum_{ij}w_{ij}}

三种标准

  • 如何评估自相关系数是否一致

    • 根据上面的公式,自相关系数其实就是target map在每个特征向量的power,并根据特征值进行缩放得到的结果。因此,在permutation的时候,我们希望得到的map在这些特征向量上的power也尽量和target的一样。因此,需要permutation map的得到的rxvr_{x'v}满足以下条件:
    • 标准1:rxv2=1r_{x'v}^2=1
    • 标准2:rxvT2m=Ixr^2_{x'v^T}m=I_{x}
    • 标准3: 对于每个特征向量,rxviT2=rxviT2r^2_{x'v_i^T}=r^2_{xv_i^T}
  • 在算法上,提供了三种评估target map和permutation map间SA一致性的标准,singleton最严,pair,triplet依次越来越松。

singleton

  • 满足标准1,2,3。

pair

  • 满足标准1,2。
  • 对于标准3,放松至rxviT2+rxvjT2=rxviT2+rxvjT2r^2_{x'v_i^T}+r^2_{x'v_j^T}=r^2_{xv_i^T}+r^2_{xv_j^T}

triplet

  • 满足标准1,2。
  • 对于标准3,放松至rxviT2+rxvjT2+rxvkT2=rxviT2+rxvjT2+rxvkT2r^2_{x'v_i^T}+r^2_{x'v_j^T}+r^2_{x'v_k^T}=r^2_{xv_i^T}+r^2_{xv_j^T}+r^2_{xv_k^T}

优势和劣势

  • 优势:这个方法运算速度非常快,因为只要找到一组rxvr_{x'v},就能找到很多permutation map。而且不需要对得到的每一个map做fit,相比burt的方法非常节约时间。
  • 劣势:如果数据跟一个特征向量非常相关,得到的map就会有点问题。详细可以参考Burt的文章(Burt et al. 2020 https://www.sciencedirect.com/science/article/pii/S1053811920305243; supplementary figure 6)

参考文献

Wagner, H. H., & Dray, S. (2015). Generating spatially constrained null models for irregularly spaced data using Moran spectral randomization methods. Methods in Ecology and Evolution, 6(10), 1169–1178. https://doi.org/10.1111/2041-210X.12407

Vos de Wael, R., Benkarim, O., Paquola, C., Lariviere, S., Royer, J., Tavakol, S., Xu, T., Hong, S. J., Langs, G., Valk, S., Misic, B., Milham, M., Margulies, D., Smallwood, J., & Bernhardt, B. C. (2020). BrainSpace: A toolbox for the analysis of macroscale gradients in neuroimaging and connectomics datasets. Commun Biol, 3(1), 103. https://doi.org/10.1038/s42003-020-0794-7

Burt, J. B., Helmer, M., Shinn, M., Anticevic, A., & Murray, J. D. (2020). Generative modeling of brain maps with spatial autocorrelation. NeuroImage, 220, 117038. https://doi.org/10.1016/j.neuroimage.2020.117038

工具包

https://brainsmash.readthedocs.io/en/latest/

https://brainspace.readthedocs.io/en/latest/