非线性混合效应模型在神经影像学中分析纵向变形

care because it involves performing calculations while re-specting the constraints inherent to the space where oursamples come from [1]. This is challenging when the ge-ometry of the space is not very well understood.Fortunately, in certain cases, the canonical representa-tions of the structured objects of interest lead to spaces thatare relatively better characterized. For instance, when thespace is a Riemannian manifold, the differential geometryprovides a suitable Riemannian metric which serves as abuilding block for inference models. Within the last tenor so years, numerous authors have reported how respect-ing the geometry/structure of the data in tasks as simple assummarizing the observed variability in the empirical dis-tribution [2] yields good empirical performance. There isa rapidly increasing body of work showing how manifold-based formulations yield powerful algorithms for problemsincluding affine motion tracking [3, 4], human detection [5],modeling spatial distribution of earthquakes [6], classifica-tion tasks [7], robot localization [8], registration [9], andnumerous settings in medical imaging [10, 11, 12].Symmetric positive definite matrices (SPD) are a goodrepresentative example of structured data that lie on a Rie-mannian manifold. In computer vision, SPD matrices areoften encountered due to covariance matrices which serveas feature descriptors in tracking [4] as well as for textureanalysis and classification [5]. Diffusion tensor magneticresonance imaging, a type of neuroimaging modality, yieldsan SPD matrix at each image voxel as the local diffusiontensor. But even aside from the SPD manifold, probabilitydensities, shape spaces, and certain warping functions areother examples of Riemannian manifolds which commonlymanifest themselves in vision applications [12, 13, 14, 15].The need for metrics, estimators, hypothesis tests and othertools for inference tasks has led to much interest in Rie-mannian analogs of classical vision and machine learningmethods. Consequently, we now have manifold versions ofnearest neighbor queries, dimensionality reduction, dictio-nary learning, sparse coding, kernel methods, boosting, av-eraging and regression. These developments were recently125400Riemannian Nonlinear Mixed Effects Models: AnalyzingLongitudinal Deformations in Neuroimaging0Hyunwoo J. Kim † Nagesh Adluru † Heemanshu Suri †0Baba C. Vemuri § Sterling C. Johnson † Vikas Singh †0† University of Wisconsin-Madison § University of Florida0http://pages.cs.wisc.edu/˜hwkim/projects/riem-mem0摘要0在视觉中，对流形值数据进行操作的统计机器学习模型正在广泛研究，其动机是应用于活动识别、特征跟踪和医学成像等领域。尽管非参数方法在文献中已经得到相对充分的研究，但对于参数模型（在小样本情况下可能具有优势）的高效表达方式直到最近才出现。到目前为止，流形回归模型（如测地线回归）仅限于对横截面数据的分析，即统计学中所谓的“固定效应”。但在大多数“纵向分析”（例如，当参与者提供多次测量时）中，固定效应模型的应用存在问题。为了满足这一需求，本文将非线性混合效应模型推广到响应变量为流形值的情况，即 f: Rd →M。我们推导出底层模型和估计方案，并展示了这种模型在纵向脑成像数据上可以提供的直接优势，无论是在群体水平还是个体水平上的分析。我们的结果的直接后果是，可以以可计算的方式进行流形值测量的纵向分析（特别是对称正定流形）。01. 引言0计算机视觉和统计机器学习中的许多应用涉及对“结构化”对象的操作，例如概率密度、特殊类别的矩阵/张量、树和图。文献中普遍认为，在现成的机器学习模型中将这些对象以向量化形式处理，往往会导致不令人满意的结果，因为欧几里得几何的基本条件，包括距离和角度，可能无法满足。在这种数据上定义机器学习模型需要25410编制在手册[16]中，该手册提供了迄今为止取得的成就以及尚未解决的技术问题的概述。基于流形的计算机视觉算法中的大部分技术贡献都集中在基本问题的非参数形式上。最近，参数形式的版本开始受到关注，这是由于在生命科学中样本量较小的应用。例如，考虑一个简单的参数模型，该模型旨在识别预测变量（特征）x与流形值响应变量y之间的关联，f：x→y。最近几年，针对这个问题开发了高效的参数算法，该模型被用于识别“基于形状”的流形变量（例如分割脑区域的形状）与受试者特定变量x（例如年龄）之间的关系。仅在几个月前，该模型的推广用于多维预测变量——通过多个预测变量（或特征）x∈Rd的线性组合来拟合响应变量y，甚至在预测变量是流形值的情况下。这样的构造还用于识别一组协变量之间的关联，同时控制一组“干扰”变量[18,19]。流形值回归。在疾病进展的纵向神经影像学研究中，统计模型需要捕捉随时间变化的形态学变化，同时控制来自同一受试者的重复测量的依赖性。考虑以下分析：我们有两组受试者，“对照组”和“疾病组”，分别对应健康对照组和高风险疾病个体。我们想要识别关于时间（或疾病进展）的组间差异。已知体素X的解剖学变化可以通过变形图的雅可比矩阵J(X)来捕捉，例如第一个时间点和第二个时间点之间相隔两年的变化。最常用的“变形”特征是雅可比矩阵的对数行列式log(det(J(X)))——一个标量体素值，用于捕捉体积/解剖学变化。所谓的柯西变形张量（CDT）[20]表示为...0J(X)TJ(X)是J(X)的更丰富表示，是SPD（3）流形上的对象，参见图1。为了理解图像中的每个特定体素与预测变量（例如年龄或疾病状态）的关联，可以在每个体素上进行预测变量和CDT之间的回归，以识别受年龄或疾病影响最大的脑区域。这可以通过最近提出的测地线回归算法[17,19,18]或流形上的核回归[21]来实现。纵向分析中的随机效应。现在，让我们考虑一个稍微复杂一些的设置，每个受试者在多个时间点提供数据，相隔几年。在这样的纵向设置中，我们获得一个CDT...0在每个连续时间点（即对）之间的图像。标准线性回归（或其流形值模拟）对时间样本的依赖性是不可知的。由于受试者在研究中被多次检查，来自同一受试者的重复测量作为受试者特定的“随机效应”。这种依赖性违反了固定效应模型（例如，广义线性模型）的独立同分布假设，包括流形版本[17，18]。固定效应模型假设所有数据都是来自具有随机噪声的响应变量Y的相同底层生成函数的独立同分布样本。如图2所示，每个受试者可能具有不同的趋势。例如，受试者A具有早期疾病发生（截距）。受试者B显示更快的疾病进展（斜率）。此外，根据参与者的年龄范围，受试者之间的变异性可能比受试者内的变异性更大。因此，对于图2中的数据，固定效应线性模型更容易适应人群水平的变异性（黑色）而不是每个受试者的轨迹（红色）。事实上，这种受试者特定的随机效应可以通过更一般的混合效应模型来建模。本研究的总体目标和贡献是为纵向响应变量集合Y是流形值变量的情况推导出公式/算法，并拟合线性（或非线性）混合效应模型。0相关工作。如上所述，视觉领域中有越来越多的工作将数据的特定几何/结构直接纳入估计问题中。各种统计构造已经推广到了黎曼流形：包括回归[22，23]，基于边界和提升的分类器[7]，分类[24]，核方法[25]，滤波[10]和字典学习[26，27]。文献中还包括与射影维数降低方法相关的思想。例如，通过所谓的PGA对主成分分析（PCA）进行推广[28]，测地线PCA[29]，精确PGA[30]，流形上的CCA[31]，带有框架丛的水平维数降低[32]，以及PGA的推广到黎曼流形的乘积空间，即张量场[24]。我们应该注意到，早期关于流形上的单变量线性模型（与测地线回归相关）的重要工作是使用了不严格的黎曼度量来研究微分同胚群的度量[33，34]。尽管有这些发展，将参数回归模型成功推广到黎曼流形的第一个结果相对较新：测地线回归[17，35]，多项式回归[36]和具有流形值响应的多变量线性回归[18]。最后，我们指出，最近在[37]中关于混合效应模型的一个独立结果处理了一维流形[0，1]，它是实线R上的单位区间，具有特定设计的度量来捕捉S形函数样式。这项工作不涉及实际的流形值变量（例如SPD）；此外，对于超过数百个体素的情况，它在计算上是不切实际的。相比之下，我们将分析的3DCDT图像超过1M+个体素。405060708090�ij��ij��ij�25420图1. 形态学研究中生成的一个数据面板示例。（a，d）移动和固定的脑图像。 （b）将（a）移动到（d）的扭曲空间网格。 （c）局部变形的矢量场。（e，f）变形场的行列式的映射。 （g，h）柯西变形张量场（CDTs）（√0J T J）。在可以分析的脑形态学的不同特征中，CDTs是本文的重点。0年龄0脑萎缩0混合效应0早发性0更快的进展0晚发性0图2.这个图展示了我们感兴趣的关键效应。每个受试者的脑萎缩进展速率（加速效应）和变化的起始时间（时间偏移）不同。常规的广义线性模型（GLM）具有固定效应无法捕捉这些效应，而在混合效应模型中包括随机效应（个体特定的斜率和截距）可以捕捉这些效应。与实际的流形值变量（例如SPD）相比，这项工作不涉及实际的流形值变量；此外，对于超过数百个体素的情况，它在计算上是不切实际的。相比之下，我们将分析的3D CDT图像超过1M+个体素。02. 初步概念和符号0我们首先简要回顾线性混合效应模型及其估计方法。然后，我们总结一些我们将使用的基本微分几何符号。本文的扩展版本和[38，39]描述了更多细节。此外，我们引入CDT和雅可比矩阵来捕捉纵向形态学脑变化。02.1. 线性混合效应模型0一般来说，回归模型（如线性/多项式）的估计假设数据来自一个无...0基本模型是具有独立同分布噪声的模型；因此，协变量/特征的效应与整个样本相关。这些模型被称为固定效应模型。例如，线性回归模型也是一个固定效应模型，如下所示：0y = β0 + β1x1 + ... + βpxp + ε，(1)0在这里，y ∈ R，x ∈ Rp，β = [β0,...,βp]T ∈Rp+1。我们可以看到，系数是'固定的'，在整个人群中都是相同的。然而，在纵向研究中（见图2），来自同一受试者的重复测量不再是独立的。我们需要更灵活的规范-通常协变量/特征对个体受试者（或群体）有不同的影响，这被称为随机效应。例如，脑萎缩和疾病进展的速率可以因个体而异，如下所示：0y_i = u_1iz_1 + ... + u_qiz_q + �_i, (2)0其中z是一个已知向量，指定样本属于哪个主体（或组），u_qi是第i个主体（或组）的第q个随机效应，记作u_i。这种固定效应和随机效应的组合形成了混合效应模型[40]。当模型是线性的时候，我们得到线性混合效应模型，下面我们将介绍它的非线性模拟。非线性混合效应模型是推导我们最终的流形数据模型的一个中间（但必要）步骤，该模型在第3.2节中介绍。模型的规范化。令y_i =[y_ij]_nij=1为主体i的n_i个重复观测的响应变量集合。这里y_i是一个n_i维的向量，垂直堆叠了主体i的y_ij响应。记号[i,j]简单地表示主体i的第j个观测。类似地，令大小为n_i×p的主体特定矩阵X_i为：0j=1，我们收集为Zi =�z1�ij� z2�ij� . . . zq�ij��25430主体i的所有p个测量值作为行组成的矩阵X_i。矩阵Z_i提供了每个主体（设计矩阵）的纵向测量次数信息。与X_i类似，我们通过指定行来定义Z_i，如下所示：0j=1。这些对应于集合0对于第i个主体，有p和q个变量（特征），我们希望估计X_i的固定效应和Z_i的随机效应对y_i的影响。在经典的设置中，线性混合效应模型（[40]）给出如下形式：0y_ij = β_0 + β_1x_1ij + ... + β_px_pij +0u_1iz_1�ij� + ... + u_qiz_q�ij� + ��ij�,0其中β_1,...,β_p是整个人群共享的固定效应，u_1i,...,u_qi是第i个主体的（特定于主体的）随机效应。随机效应u_i =[u_1i u_2i ∙∙∙u_qi]T被假设服从多变量正态分布。未解释的随机误差�_ij来自正态分布N(0,Σ^2�)。我们可以使用矩阵表示来简洁地写出模型：0y_i = X_iβ + Z_iu_i + �_i. (3)0让'vstack(∙)'表示参数的垂直堆叠。通过表示y = vstack(y_1, y_2, ...,y_N)，以及类似地表示X、Z、u，所有N个主体的最终模型可以表示为：0y = Xβ + Zu + σ^2�I,0其中u�N(0,˜Σ)且˜Σ=diag(Σ_1,Σ_2,...,Σ_N)=Σ�I（当Σ_i=Σ�i），且Z=diag(Z_1,Z_2,...,Z_N)。一般来说，线性混合效应模型的估计没有解析解，除非˜Σ和σ^2�已知。02.2. 基本微分几何符号0设M是一个可微（光滑）流形，维数任意。可微流形M是一个局部类似于欧几里得空间并具有全局定义的微分结构的拓扑空间。Riemann流形(M,g)是一个带有平滑变化的内积g的可微流形M。所有切空间上的内积族称为Riemann度量，它在曲面上定义了各种几何概念，如曲线的长度等。测地线是一条局部最短路径，类似于R^d中的直线，这样的测地线将定义我们在SPD矩阵空间中的协方差矩阵轨迹。与欧几里得空间不同，注意在曲面上两点之间可能存在多条测地线。因此，M上两点之间的测地距离定义为连接两点的最短测地线的长度（即SPD矩阵）。测地距离有助于测量我们轨迹估计的误差（类似于Frobenius或0在欧几里德设置中，从 y i 到 y j 的测地线由以 y i为锚点的切空间上的切向量参数化，其中 Exp ( y i , ∙ ) : T yi M → M 是指数映射，其逆映射是对数映射，Log ( y i , ∙ ): M → T y iM。这两个操作在流形和切空间之间来回移动。与上述符号分开，矩阵指数（和对数）简单地是 exp( ∙ )（和 log( ∙)）。最后，平行传输是流形上的广义平行平移。给定一个可微曲线 γ : I → M，其中 I 是一个开区间，沿着曲线 γ平行传输的 v 0 ∈ T γ ( t 0 ) M可以解释为在流形上平行平移 v 0，保持 v ( t ) 和 γ之间的角度。从 y 到 y ′ 的平行传输是 Γ y → y ′v（更多细节请参见扩展版本）。03. CDT图像的纵向分析0设 I i,j 表示从主体 i 在时间点 j 获取的图像。给定连续访问(j, j + 1) 的图像 I i,j 和 I i,j +1，我们可以计算将两个图像对齐的变形[41, 42]。设 I i, 1给出主体特定的坐标系，表示为 Ω i。这将为主体 i在时间上经历的变形（即CDT图像）提供（中间的）共同坐标系，j = 1, 2, ∙ ∙ ∙ , n i。表示每个主体 i 的所有 (n i - 1)个时间变形（即CDT图像）的全局模板表示为Ω。然后，对于每个图像（更准确地说，对于每个 (I i,j + 1,I i,j) 对），给出了体素（空间位置）vox ∈ Ω的非线性变形 Φ( vox )0Φ : I i,j + 1 → I i,j0Φ( vox + d vox ) = Φ( vox ) + J ( vox ) d vox + O ( d vox 2 )0其中 J ( vox ) 表示位置 vox处变形的雅可比矩阵。CDT的一个好的特性是它保持雅可比矩阵的行列式，因为 det( J ( vox )) >0。因此，CDT表示在不影响体积变化信息的情况下，很好0J T J )（证明在扩展版本中）。由体素位置 vox组成的CDT“图像”与 I 1,1 大小相同，并由一个 3 × 3SPD矩阵√派生自黑盒微分同胚求解器0J T J是每个体素的CDT图像。它提供了主体的两个纵向图像之间的变形场。各种结果已经描述了CDT图像在分析中的好处；我们实验的一个例子在第6节中。请注意，对于一对实数 a和 b，线性模型“y = ax +b”的预测始终是一个实数。但是，如果 a 和 b是SPD矩阵（即CDT），则不再成立：“y = ax +b”可能不是一个有效的SPD矩阵（例如，考虑 x < 0的情况）。除了预测的有效性之外，在经典线性回归中，拟合的残差 (y - (Ax + b)) 可以直接估计̸25440在流形值数据上，SPD流形上的误差（残差）必须通过流形上的测地距离来测量预测值 ˆ y 和真实响应 y之间的距离。因此，任何需要扩展到流形值数据的机器学习模型都必须特别关注这些问题，并且通常需要专门的参数估计优化方法[17, 18,21]。考虑到这些问题，我们首先描述了具有特定主体截距的线性混合效应模型，然后过渡到流形值数据的非线性混合效应模型。03.1. 具有特定主体截距的模型0我们知道，在任何纵向数据集中，重复测量的误差/噪声是相关的。为了考虑这一方面，一种常见的方法是将随机效应表示为干扰参数。如果集合 { i = 1, i = 2, ∙ ∙ ∙ , i = N }索引列，我们可以将设计矩阵 Z 写成 diag ( � n 1 , � n 2 , ∙ ∙ ∙, � n N )，其中 � n i = [1 ∙ ∙ ∙ 1] T ∈ R n i 。然后，模型 (3)变为0y = β0 + βTx + ui, (5)0其中y, β0, ui ∈ R，β, x ∈ Rp。注意，z ∈ RN从矩阵Z中恢复出与受试者i的访问j对应的特定行，将其与u的点积给出受试者特定的随机效应，ui =ziju。该模型存在两个问题。1）对于整个人群，它具有相同的斜率β，而研究中的受试者可能具有不同的疾病进展速度；2）另一个问题是对ui的解释，它被视为y空间或x空间中受试者特定的偏移，即取决于我们将其移动到(5)式中的等号的左侧还是右侧。在医学应用中，模型的可读性对于理解疾病非常重要。我们的解决方案涉及明确地为x和y添加受试者特定的偏移。03.2. 具有ψi(x)的非线性混合效应模型0基于上述动机，我们可以通过一个受试者特定的随机函数ψi(∙)来扩展线性混合效应模型，如下所示0y = β0 + zijui + β.ψi(xij) (6)0根据ψi(∙)的形式，(6)可以是非线性混合效应模型（NLMM）。当ψi是Identity时，我们得到一个线性混合效应模型。在我们的分析中，我们使用ψi(x):=αi(x−τi−t0)+t0，灵感来自于[43]，其中每个受试者可以有自己的疾病进展速度（αi）和不同的起病时间τi，但(β0，β，t0)对于整个人群是共同的。然后我们有0y = β0 + zijui + β(αi(xij−τi−t0) + t0). (7)0注意，这种扩展与广义线性混合效应模型[44]不同，例如，y = h−1(xβ +0zijui), ui � N(0, Σi).接下来，我们将(5)和(7)中的混合效应模型扩展到流形值数据。04. 混合效应模型在流形上0线性混合效应模型（LMM）可以根据加法的顺序和解释方式在流形设置中以多种方式进行扩展。例如，回想一下欧几里得空间中加法的结合性(a + b) + c = a + (b +c)在流形上并不直接转化，即Exp(Exp(a, b), c) ≠Exp(Exp(a, c′),b′)，其中b′和c′分别是b和c的平行传输切向量，使得它们位于正确的切空间中。LMM在（3）中的自然扩展可以写成0y = Exp(Exp(Exp(B, Vx), Uiz), ε), (8)0其中y, B, Bi ∈ M，V ∈ TBMp，Ui ∈ ThijMq，hij =Exp(B, Vxij)，xij ∈ Rp和zij ∈Rq。请注意，流形M上的基点B类似于（5）中的截距β0，而V（和Ui）对应于斜率β（和随机效应ui）。不幸的是，上述模型与Ui相关的问题。请注意，在hij处，Ui在不同的切空间中使用。此外，特别是在具有GL不变度量的SPD（n）流形上，切向量的范数随着相应切空间的基点B的函数而变化，即∥U∥2B = �U, U�B =tr(UB−1UB−1)。因此，相应的尺度可能不同。为了解决这个问题，我们改变指数映射的顺序，并提出了一种在流形上具有受试者特定截距（y偏移）的混合效应模型。此外，与欧几里得空间不同，通常情况下，x的偏移与y的偏移之间没有等价关系。因此，我们可以明确地添加x的偏移，表示为τi。然后，我们在流形上的表述如下0y � ij � = Exp(Exp(Bi, ΓB→Bi(V)(x � ij � − τi)), εij), (9)0Bi = Exp(B, Ui z � ij �), (10)0其中τi ∈ Rp，Bi ∈ M，V ∈ TBMp，Ui ∈TBMq，其余变量与之前相同。与标准混合效应模型一样，假设Ui服从正态分布。有了这个构造，我们现在可以在发病时间中包括一个受试者特定的时间偏移（类似于（7）），并假设疾病的进展具有相同的整体模式，但只有其速度/速率和发病时间在受试者之间有所变化。该模型由0y � ij � = Exp(Exp(Bi, ΓB→Bi(V)αi(x � ij � − τi − t0), εij)))0Bi = Exp(B, Ui), (11)0其中αi是受试者特定的加速度，αi < 1（和αi >1相应）表示比人群慢（和快相应）。此外，τi是受试者特定的发病时间偏移25450时间：τi >0表示晚发病时间，而t0是发病时间的全局偏移。最后，Ui（或Bi）是表征响应变量空间中受试者特定偏移的切向量（或基准点）（参见图2中的欧几里得情况）。与经典设置一样，我们可以假设随机效应服从正态分布，即ΓB→IUi �NSYM(0，σ2U)，αi � N(1，σ2α)，τi � N(0，σ2τ)。05. 参数估计过程0在流形上估计混合效应模型是具有挑战性的；可证明准确的估计方法几乎不可能在整个脑图像上运行。即使对于欧几里得响应变量，非线性混合效应模型的高效估计方法仍在文献中积极研究中，例如交替算法[44]、拉普拉斯和自适应高斯积分算法[46]，以及具有MCMC的广义EM算法[47]。不幸的是，这个问题在流形设置中只会变得更糟。即使对于实线上的单变量流形，拟合非线性混合效应模型也需要大约一天的时间[37]。在我们的数据集中，体素的数量超过100万，对于整个大脑进行精确分析是不现实的。因此，我们提出了基于模型的某种几何解释的近似算法。05.1. RNLMM的估计0我们观察到我们模型的主要构建块，Riemannian非线性混合效应模型（RNLMMs），是一个流形值多元通用线性模型（MMGLM）。该模块具有称为Log-Euclidean框架的高效参数估计。众所周知[19]，在实践中，可以在响应变量Y的Frechet均值的切空间中以一个居中的X近似估计，即B≈¯Y，τ≈¯X。与全局流形值线性模型[19]一样，参数V将对应于完整的数据集；然而，我们允许基准点B和τ的受试者特异性变异性通过Bi(r)和τi(r)来实现，其中r ∈R可以被视为共享全局V的局部模型之间的混合率。这由0y � ij � = Exp(Exp(Bi(r), ΓB→Bi(r)(V)(x � ij � − τi(r))), ε). (12)0换句话说，r ∈ R 是用于全局平均人口特定基准点 Bi(r)和时间偏移量τi(r)的权重——所有受试者共享固定效应V，但每个受试者对应于自己的偏移量τi(r)和Bi(r)在x和y空间中。当r =0时，模型简化为[18]中仅具有全局截距的模型（请参见扩展版本）。我们对（12）的估计总结在Alg中。1中，其中y �� ij �是通过以下方式获得的切向量：将响应y � ij �映射到Bi(r)处的切空间，并将该映射平行传输到TBM。我们现在简要地解释一下0算法1 流形混合效应模型1：计算每个受试者的均值¯yi，0¯yi=argminy∈M0j=1 d(y, y�ij)2。（13）0类似地计算整个种群的¯y。2：给定r，通过以下方式解决Bi(r)（¯yi和¯y的插值）：0Bi(r)=¯y(¯y−1¯yi)r=¯yi(¯y−1i¯y)1−r，0≤r≤1。03：y�ij=Γ¯Bi(r)→BLog(¯Bi(r), y�ij)。04：通过群作用将y�ij转运到I。5：通过x�ij(r)=(1−r)(x�ij−¯x)+r(x�ij−¯xi)对x进行居中。06：使用MMGLM[19]在传输的y�ij和x�ij(r)上计算V�。7：预测由以下给出：0ˆyij=Exp(Bi(r), ΓB→¯Bi(r)(V�)xij(r))。（14）0在描述如何高效进行估计时，我们首先在步骤2中解决两个SPD矩阵的线性插值，这是使用SPD流形上的测地距离的解析解的形式[48]（请注意，当样本数较大时，存在递归方案[49]）。在步骤4中，我们使用群作用从B传输切向量到I，反之亦然，这被认为比平行传输更高效，但等效[50]。05.2. 用ψi(x)估计RNLMM0使用特定于受试者的随机函数ψi(∙)估计模型（11）涉及一些额外的技术挑战。为了降低问题的复杂性，我们首先找到主要的纵向变化方向η，控制受试者特定的随机效应¯Yi和¯Xi（因为Ui和τi是随机效应）。该方案在算法2中描述。0算法2 计算纵向变化方向01：计算响应的种群Frechet平均值¯y。2：计算每个受试者的Frechet平均值¯yi。3：解决y�ij=Γ¯yi→ILog(¯yi, y�ij)。04：解决x�ij=xij−¯xi，其中¯xi=Ej[xij]。05：收集X�=[x�1;...x�N]和Y�=[y�1,...y�N]。6：通过最小二乘估计计算纵向变化方向η，η=((X�)TX�)−1((X�)TY)∈TIM。0一旦估计出纵向变化方向η（整个种群的固定效应），我们一次解决一组参数。此过程在算法3中描述，我们在给定η的估计的情况下解决所有参数。请注意，对参数的先验偏差可能会降低统计功率。因此，我们对所有参数使用了非信息性先验。虽然算法3使用了非信息性先验，但稍作修改，我们可以轻松地将正态分布先验纳入其中，请参阅扩展版本。�,�,RNLMMs on longitudinal CDTs. We now present re-sults using our Riemannian nonlinear mixed effects models(RNLMM) using subject specific transformation functionsψi(x�ij�) (6). Here, x�ij� ∈ R is used to represent the ageof each subject at the previous visit and y�ij� ∈ M, (theCDT image calculated from scans at two points). For theseresults, we used data from subjects who had at least threevisits. We estimated our model at each voxel in the brain(1.3M+) using a total of N = 228 participants that had atleast two CDT images. The maps for acceleration (αi), spa-tial shift (Ui) and time shift (τi) for each of the subjects offerunique advantages. For instance, these maps are not offeredby standard linear mixed effects models where only a sub-ject specific slope or intercept is used as the random-effects(independently noted in [37]). Fig. 4 shows four representa-tive subject-specific acceleration maps. The regions wherethis specific individual has a faster (slower) aging (or dis-ease progression) compared to the population average rateare colored in yellow (and blue) color-scales respectively.These RNLMM maps can be used to perform additional“downstream” statistical tests using parametric tests. Here,we cover two specific examples. In Fig. 4, we show thekind of results our model can offer at the individual level.Fig. 4(a)-(d) shows four results, each pertaining to a dif-ferent participant in the study. Fig. 4(a)-(b) show maps for25460算法3 具有ψi(x)的流形混合效应模型01：计算种群的Frechet平均值¯y∈M。2：计算每个受试者的Frechet平均值¯yi∈M。3：通过算法（2）计算主要纵向变化方向η。4：计算受试者特定的基准点（随机效应）Bi=Exp(B, Ui)，其中Ui=argminUi d(¯yi, Exp(B,Ui))2+λUi∥Ui∥2B。5：y�ij=ΓBi→ILog(Bi,yij)。6：当直到收敛时执行以下步骤：7：通过固定所有其他变量计算共同的变化速度V=cη和共同的时间截距t0=b/c。0空0ij p T ij q i 空格0ij p T ij q i空格0ij p T ij p ij0空格b c0空格=0ijqij0ijqTiTij0ijpTijy�ij0其中b := t0c，qi := η(1−αi)，pij := η(αixij−αiτi)。08:给定V，t0，通过广义最小二乘估计计算受试者特定的加速度αi和时间偏移τi，其中αi和τi的先验为di/αiΣjWTijWij−ΣjWTijVΣjWTijV−ΣjVT V0αidi0=0ΥTijWijΥTijV0其中Υij := y�ij−Vt0，Wij := V(Xij−t0)和di = αiτi。9: end while0Alg. 3与Alg.1有许多共同的步骤。在第7步中，我们通过固定所有其他变量（c是一个虚拟变量）来估计固定效应V和t0。在第8步中，我们通过固定V和t0来估计受试者特定的随机效应αi和τi（di是虚拟变量）。更多细节，包括推导，可以在扩展版本中找到。06. 实验0目标。我们实验的总体目标是评估所提出的公式是否可以作为驱动神经影像学图像数据纵向分析的核心模块。为此，在进行特定疾病的纵向数据分析时，该过程应该为群体分析提供有意义的结果，例如，当将人群分为具有分层变量（例如性别或疾病风险因素）的子群时，对于受试者/体素特定的“随机”效应（特别是加速度和空间偏移）的统计显著差异的“图”应该是科学可解释的，同时与基线一致。我们下面的实验展示了模型满足这一要求的程度。数据。CDT图像（表示受试者特定的变形）来自于一项早期阿尔茨海默病（AD）的纵向神经影像学研究。纵向变形（或变换）是通过在两次连续访问之间进行受试者内部的T1加权图像配准获得的，即Φi,j：Ii,j→Ii,j+1。体素级的CDTs是从空间导数�Φi,j(vox)中得到的。0变形场的细节。由于空间限制，CDT与Jacobian行列式的比较。我们首先进行了一个激励实验，以证明使用CDTs而不是Jacobian行列式的合理性，即CDTs是否实际上携带更多信息？我们使用CDTs测试了中年人群与老年人群之间大脑纵向变化的群体差异，并将这些结果与通过行列式获得的结果进行比较。为了避免这种比较中的混淆因素，我们使用Cramér的检验，这是一种非参数检验，适用于单变量和流形值数据，因为它不需要对零分布进行任何规定[51