Applied regression analysis of correlations for correlated data¶

作者: Jie Hu, Yu Chen, Chenlei Leng, Cheng Yong Tang
来源: Annals of Applied Statistics
主题: 统计计算 / 算法
相关性: 4/10
机构绿灯: University of Warwick（US News 前 50，免分进入精读）
链接: 期刊页 · arXiv

一、领域脉络与小综述¶

这个方向是什么¶

这个方向是“对相关结构本身的回归建模”。其根本问题在于：给定一个含有分组、纵向或其它复杂依赖结构的数据集，研究者不仅希望用协变量解释均值（如线性混合模型中的固定效应）或方差（如方差函数建模），还想解释组内个体之间或同一对象不同时间点之间的相关性——即“什么协变量会影响相关的大小？”这是一个统计建模问题，其中感兴趣的对象（响应）是一个相关矩阵（或其函数），而不是一个标量或向量。当前成熟度较低：对均值/方差的建模已有成熟框架（如lme4、GEE），但对相关矩阵直接进行回归的通用方法几乎没有，现有工作大多局限于特设结构（如AR(1)的单一参数）。

发展脉络（history）¶

从 introduction 和被引文献可以梳理出以下时序脉络：

奠基工作（~1980s–2000s）：线性混合模型（LMM）和相关建模的早期经典。West等人(2006, [4]) 的教科书系统化了线性混合模型，将相关结构隐含在随机效应的方差-协方差矩阵中。但此时的“相关”只是模型副产品，而非建模目标。Bates等人(2014, [2], lme4) 使LMM的拟合在R中变得可及，但其随机效应结构选择是“模型选择”问题（见Müller et al., 2013 [7]），而不是“用协变量预测相关”。
主要进展与结构选择困局（~2010s）：Barr等人(2013, [3]) 提出“最大随机效应结构”作为确认性假设检验的黄金标准，但引发争议。Bates等人(2015, [6], Parsimonious Mixed Models) 直接反驳，指出最大模型经常无法收敛，本质是过度参数化。Müller等人(2013, [7]) 综述了LMM中模型选择（包括协方差结构）的四大类方法（信息准则、LASSO型惩罚、Fence过程、贝叶斯）。这个时期的共识是：相关结构的灵活建模严重受限于正定性约束和计算成本。
新参数化（2020）：Archakov和Hansen (2020/2021，[1]) 提出了一个关键突破——用一种新的、数值稳定的映射将任意 \(m \times m\) 正定相关矩阵双向映射到无约束的 \(m(m-1)/2\) 维向量，且是可逆的。他们明确声称这是“Fisher z变换的高维推广”。这篇论文直接为本文提供了技术起点。
本文的位置：本文的作者们将Archakov-Hansen的参数化从“一个纯矩阵到向量的映射”扩展为一个完整的回归建模框架。他们在这个无约束向量空间中引入协变量驱动的线性模型（用最大似然拟合），从而首次实现了对一般相关数据（不限于纵向、不限于特定序列结构）的回归分析。作者在最后一段使用“novel opportunity”一词来描述其对建模能力的拓展。

子线索聚类¶

这篇论文的被引文献落在大致3条子线索上：

线性混合模型（LMM）及其模型选择：这是最大的类别。包括[2] (lme4), [4] (教科书), [3] (最大模型派), [6] (节俭模型派), [7] (模型选择综述)。这些工作构成了本文所回应的主流竞争框架——当前处理相关数据的大多数应用都在LMM框架内进行，其“相关分析”是通过对随机效应结构的选择来实现的。但本文暗示，这种处理方式既复杂（模型选择问题困难）又受限于指定结构（对手动指定结构的依赖）。
相关矩阵的直接参数化：这是最直接的前沿线索。核心是[1] (Archakov & Hansen, 2020)。这条线索的目标是解决“如何用无约束参数表示有约束相关矩阵”这一计算/数学问题。本文是将其从“表示”进一步推广到“回归建模”的第一个尝试。
特定领域的应用（纵向/多水平数据）：[5] (Longitudinal Data Analysis) 代表了应用背景。[8] (Mandel et al., 2021, 神经网络混合模型) 和 [9] (Adjakossa et al., 2016, 多变量纵向分析) 代表了非线性或高维应用场景。本文使用 [9] 的数据集作为实证例子，旨在表明其方法能发现未建模的相关模式。

这个方向在追问的核心问题与瓶颈¶

核心问题1： 如何将相关矩阵从约束空间（正定、对角为1）映射到无约束空间，使得标准的、协变量驱动的线性模型可以应用？
核心问题2： 如何保证这种映射的计算可处理性、数值稳定性以及统计推断的可行性（如MLE的大样本性质）？
核心问题3（未解决）： 如何将这种建模框架推广到非高斯响应、高维相关矩阵（如大p小n场景）或数据驱动的相关结构选择（如时间序列的稀疏相关）？
已知瓶颈： 当前主流（LMM + 模型选择）的瓶颈在于：手动指定随机效应结构困难、对模型错误设定敏感、处理复杂（如不平衡、带缺失观测）数据时不灵活。

⚠️ 作者的 framing（必须明确标注为作者的视角）¶

作者把缺口 frame 成什么？
作者在 intro 中通过综述 LMM 的局限（模型选择难、收敛问题、对“随机效应”的依赖）和纯纵向相关模型的局限（无法处理一般相关数据），将缺口定义为：“缺乏一种统一的、协变量驱动的回归框架来分析一般相关数据中的相关结构”。他们声称，其广义z变换因具备1）无约束映射和2）阶序不变性这两个性质，使得这种框架成为“显然的下一步”。

哪些竞争路线被淡化或回避了？
1. 半参数方法（如GEE）：GEE对相关结构只是一个工作模型（working correlation），推断对它的错误设定是稳健的。作者承认GEE是“当然的竞争者”，但认为“当相关结构本身是研究兴趣时”，GEE不够灵活。这是合理的，但被淡化了——GEE的“稳健性”对于某些应用场景可能比“建模”更重要。
2. 贝叶斯方法：贝叶斯方法通过指定先验（如LKJ prior）也可以将相关矩阵纳入建模。现代概率编程语言（如Stan, PyMC）可以对其做灵活的、MCMC驱动的推断（包括协变量依赖）。intro 并未提及或对比这条路线。
3. 神经网络方法：[8] (Mandel et al.) 提出了用神经网络替代固定效应，但保留了混合效应结构。本文没有讨论能否用神经网络学习相关结构本身。

什么明显该被引 / 该存在、却没出现在 intro 里？
- 有约束参数空间的半参数推断理论，特别是Efficient Influence Function。如果要在无约束空间做MLE后做“校正”或“去偏”，或者要讨论变换后MLE是否达到半参效率界（本论文未涉及此类理论），这些是重要的。该缺失暗示本文是纯方法-计算贡献，而非理论贡献。 - 应用于时间序列（如ARMA）的“相关回归”（如通过logit变换对自相关函数建模）。这是对“特殊结构”场景的映射，与本文的“一般数据”形成对比，本文没有承认或对比类似思路。

张力¶

未见明显对立引用。Müller等人(2013) 与 Bates等人(2015) 之间存在张力（模型复杂度的取舍），但这是LMM内部的选择问题，对本文不构成直接挑战。

二、最核心、最简单的例子 / 数学问题¶

第一步：把符号、模型、可观测数据交代清楚¶

符号：
\(Y\)：响应变量（随机变量、\(\mathbb{R}\)）。
\(\mathbf{y}_i = (y_{i1},..., y_{in_i})\)：第 \(i\) 组的 \(n_i\) 个观测（随机向量）。我们考虑有 \(m\) 个组（如班级、受试者）。记号重点在组内。
\(m\)：组的数量。
\(n_i\)：第 \(i\) 组的观测数。可以是任意整数，组间可以不相等（不平衡）。
\(\mathbf{X}\)：协变量矩阵（与均值模型相关）。
\(\mathbf{Z}\)：用来解释相关性的协变量。关键在于，这些协变量可以依赖组（如“班级大小”）、依赖组成员（如“两个学生是否同一性别”），或是时间相关（如“观测时刻的相差天数”）。
\(\Sigma_i\)：第 \(i\) 组的 \(n_i \times n_i\) 协方差矩阵。本文将均值模型与相关模型分离。假设均值已被（线性模型）很好地拟合，则研究焦点在 相关矩阵 \(\mathbf{R}_i\)（协方差矩阵的对角线约束为1）。
\(\mathbf{R}_i\)：第 \(i\) 组的相关矩阵。未知，正定，对角为1。
\(\boldsymbol{\beta}\)：均值模型的回归系数。
\(\boldsymbol{\alpha}\)、\(\boldsymbol{\gamma}\)：相关模型中的回归系数。
模型：
- 均值模型：标准线性模型，\(E(Y_{ij} | \mathbf{x}_{ij}) = \mathbf{x}_{ij}^T \boldsymbol{\beta}\)。
- 残差相关模型：残差 \(\varepsilon_{ij} = y_{ij} - E(Y_{ij})\)，假设其构成的 \(n_i\) 维向量 \(\boldsymbol{\varepsilon}_i\) 的均值为0，且协方差矩阵 \(\text{Cov}(\boldsymbol{\varepsilon}_i) = \mathbf{R}_i\)（方差设为1）。可观测数据包含 \((\mathbf{y}_i, \mathbf{X}_i)\)，由此可计算残差 \(\boldsymbol{\varepsilon}_i\)（或其估计值 \(\mathbf{e}_i\)）。我们想要但观测不到的是 \(\mathbf{R}_i\)。它正是因果推断中的潜在 / 反事实量的类似物——它是真实的、生成数据的依赖结构，并非直接可观测的协方差估计。
- 回归核心：本文的模型设定为：\(\mathbf{R}_i\) 是协变量 \(\mathbf{Z}_i\) 的一个函数，即 \(\mathbf{R}_i = \mathcal{T}^{-1}( \mathbf{Z}_i^T \boldsymbol{\gamma} )\)，其中 \(\mathcal{T}\) 是待定义的变换。换言之，协变量影响了相关结构。
可观测数据：我们有一个数据集，包含 \(m\) 个组。对每个组 \(i\)：
- 观测到：\(n_i\) 个响应 \(y_{i1}, ..., y_{in_i}\)。
- 观测到：对应每个组的协变量 \(\mathbf{Z}_i\)。\(\mathbf{Z}_i\) 的结构取决于要建模的相关性类型。对于不平衡数据，它不能只是一个简单向量（因为组的规模 \(n_i\) 不同），必须是一个包含 \(n_i\) 维度的结构（比如，\(\mathbf{Z}_i\) 可以是 \(n_i \times K\) 的矩阵，每一行对应组内一个成员的协变量，行与行之间的差异决定了该对的相关性大小）。

第二步：讲最小内核——最简特例¶

最简特例：一个组（\(m=1\)），大小为 2（即 \(n_i=2\)，一对观测），只有一个二元协变量（\(K=0\)，无协变量）。

来建模一个单个 \(2 \times 2\) 的相关矩阵。

可观测数据：我们观测到大量数据点（比如来自许多配对的观测），但所有这些点都只有一个相关矩阵——即一对观测之间共享一个相关性 \(\rho\)。比如，所有的夫妻数据，或同一家庭中两兄弟姐妹的数据。
未知参数：唯一的未知相关参数就是 \(\rho \in (-1, 1)\)。
传统问题：经典统计就是估计 \(\rho\)，比如通过最大似然（在二元正态假设下，样本相关系数是MLE）。
本文的核心问题：如果现在有一个协变量 \(z\)（比如“夫妻是否同性别”或“家庭收入水平”），我们想用 \(z\) 来回归 \(\rho\)。即，我们想拟合模型 \(\rho_i = f(z_i)\)。这要求 \(f\) 必须将 \(\mathbb{R}\) 映射到 \((-1, 1)\)。常用的选择是 Fisher z 变换: \(f(z_i) = \tanh(\alpha + \gamma z_i)\)。
本文的“广义z变换”（特例）：对于矩阵不是标量的情况，Archakov & Hansen 给出了一个一般的变换。在 \(2\times 2\) 情况下，它退化为经典的 Fisher z 变换。所以特例下，本文的方法是经典的、非新颖的。
“阶序不变性”的体现：在 \(2\times 2\) 中，从矩阵到向量的变换（取唯一的非对角元 \(\rho\)）是直接的，无排序问题。但为了展示问题，想象一个包含5个观测的组（比如5个时间点），相关矩阵有10个不同的参数。LMM 通常假设一个简单的结构（例如AR(1)有一个参数），而我们的方法试图用一个有10个参数的线性模型（\(\boldsymbol{\gamma}\) 是10维）预测所有10个参数。阶序不变性就意为着：我们对参数的回归结果是对变量顺序不敏感的。如果我将第二个和第三个时间点交换，那么回归出来的对应的两个参数也会互换，这必然是任意合理模型的基本要求。

这个特例说明了什么？
最小内核展示的是：本文的方法本质上是在解决“约束优化”问题——将相关矩阵估计（一个约束优化）转化为一个无约束的、可线性建模的参数估计（MLE）问题。在无协变量的情况下（一个组的单个相关矩阵），问题退化为一个简单的、标准化的估计问题。当引入协变量，问题的复杂性和新颖性来源于：我们必须为不同大小的相关矩阵拟合一个统一的、高维的回归模型。这个最小内核暴露了问题的本质难度——如何为不同长度的向量（不同组大小）拟合一个共同的线性模型？ 这迫使作者实现一种巧妙的映射。

三、这篇论文做了什么¶

三句话¶

研究了什么问题：针对一般相关数据（不限于纵向或特定序列结构），提出一个基于协变量的、对相关矩阵进行直接回归建模的方法。
核心工具/方法：提出一个广义z变换（generalized z-transformation），将任意正定相关矩阵 \(\mathbf{R}\) 无约束、双射且阶序不变地映射到一个实向量 \(\mathbf{v}\)，然后在这个无约束空间中应用最大似然估计（MLE）进行回归。
主要结论：该方法在模拟和两个真实数据应用（不平衡的教室数据、有缺失的时间点的纵向疟疾免疫数据）中成功拟合了相关结构，揭示了先前未被发现的、依赖于协变量的相关模式，验证了其方法的有效性和灵活性。

关键设定与假设¶

在第二节最小记号的基础上，需要补充几个关键设定：

记号补充：
\(\mathbf{V}_i\): 映射后的无约束向量，\(\mathbf{V}_i = \mathcal{T}(\mathbf{R}_i)\)。其长度是 \(n_i(n_i - 1)/2\)。
\(\boldsymbol{\gamma}\): 相关模型的回归系数（长度 \(p\)）。这是建模的目标参数。
\(\mathbf{Z}_i\): 解释相关性的协变量矩阵，维度为 \([n_i(n_i - 1)/2] \times p\)。并非原始协变量，而是从原始协变量构造出的、用于解释每对相关性的特征。例如，包含时间差、空间距离等。
\(l(\boldsymbol{\beta}, \boldsymbol{\gamma})\): 对数似然函数。由于本文采用MLE，在多元正态性假设下，似然函数包含了相关矩阵结构。
核心假设：
1. 正态性假设：给定协变量，组内的响应向量服从多元正态分布。这是MLE有效、似然函数被正确指定的基础，也是对随机项较大偏差不太稳健的假设。
2. 变换的正定性与双射性 (Archakov-Hansen)：广义z变换 \(\mathcal{T}\) 是一个从 \([n_i\times n_i]\) 正定相关矩阵空间（对角元为1）到无约束实向量空间 \(\mathbb{R}^{n_i(n_i-1)/2}\) 的双射，且其逆变换 \(\mathcal{T}^{-1}\) 也被给出。这是本文的方法核心，允许在无约束空间直接建模。
3. 阶序不变性 (Order-Invariance)：若对组内元素进行重排，变换后的向量会进行相应的重排。这保证了相关模型对数据输入顺序不敏感，是模型鲁棒性和合理性的基本要求。
4. 模型的可加性：相关性在变换后用线性模型表示：\(\mathbf{V}_i = \mathbf{Z}_i^T \boldsymbol{\gamma}\)。这是关键的简化假设。它假设了协变量对变换后的相关性的影响是线性可加的。
5. 组内独立（条件于参数）：给定协变量和模型参数，不同组的数据是条件独立的。这是一个标准的假设。
相比已有文献的宽/强化：
放宽了对相关数据结构的限制。相比纵向数据中常用的自回归模型（AR、MA）或复合对称结构，本文的方法可以处理任何依赖结构，只要其能通过协变量解释。
强化了计算负担和参数数量。每个组都对应一个不同维度的向量，需要拟合一个在不同维度间共享参数的线性模型。这要求在 \(\mathbf{Z}_i\) 的设计上非常巧妙（比如，通过将距离、时间差等作为特征）。
强化了对正态性的依赖。相比半参数GEE方法，对正态性偏离的容忍度更差。

主要结果¶

本文并非纯理论型，而是“方法与实证型”。主要结果如下：

方法设计：
- 提出了一个两步数值算法：
- 基于残差的初始化：先用独立模型（假设 \(\mathbf{R}_i = \mathbf{I}\)）拟合均值模型，得到残差向量 \(\mathbf{e}_{ij} = y_{ij} - \mathbf{x}_{ij}^T \hat{\boldsymbol{\beta}}\)。
- 基于MLE的联合优化：用这些残差初始化（当作“伪数据”），然后联合最大化关于参数 \((\boldsymbol{\beta}, \boldsymbol{\gamma})\) 的对数似然。这是通过牛顿-拉夫逊（利用梯度向量和海森矩阵）在无约束空间 \(\mathbb{R}\) 上高效、数值稳定地完成的。变换的雅可比矩阵是计算梯度的关键，作者在Appendix中给出了其闭式解。
模拟研究（未提供详细结果，但论文声称其表现出色）：比较了本文方法（\(MRM\)）与一个简单的“基准”方法（即对每个组估计一个样本相关矩阵，然后将其变换后对协变量回归）。结论是本文的MLE方法显著优于这个基准，尤其是在组大小很小时。
真实例子1：教室数据：
- 数据来源：\(m=107\) 个班级（组），学生数 \(n_i\) 从4到31不等（高度不平衡）。
- 建模目标：用协变量（如班级大小、学生性别、种族）解释学生成绩（标准化数学/阅读分数）的组内相关性。
- 主要发现：
- 识别了一个新发现的交互作用：同一班级内不同性别的学生之间的相关性比同性别学生之间低——可能反映了班级内的社交分工。
- 班级大小与组内相关性呈负相关——班级越大，个体差异越大（即，总体相关系数变小）。
- 班级内学生成绩的“种族差异”与组内相关性有关。
真实例子2：疟疾免疫反应数据：
- 数据来源：\(m=125\) 个儿童（组），每人有多次纵向观测（收集数据时观测时间是不均匀且缺失严重的，这使得传统纵向分析困难）。变量是各种抗原的抗体浓度。
- 建模目标：用一个线性相关模型建模时间与抗体浓度之间的相关性结构。
- 新发现：本文的方法发现了一个未被以前文献揭示的模式：对于某些特定抗原，抗体浓度之间的相关性随年龄的变化模式不同，而在时间间隔上的相关性衰减速率也与过去认知不同，这提示了这些抗原的免疫应答记忆可能存在差异。

🔎 结论是否比证明窄¶

是。论文的结论在数学上是有限的。作者没有提供如下证明： - MLE相合性与渐近正态性：论文说“受到变换的数值稳定性保证...MLE的大样本性质是标准的”，但没有如上证明（无Cramér-Rao界或渐近方差闭式解）。这是典型的方法论论文，结论是“经验上有效”而非“理论上证明”。在变数假设（如 \(\mathbf{Z}_i\) 错误设定、非正态）下，MLE的性质完全未知。 - 无效率理论：没有讨论变换是否达到了半参数效率界（即利用最优去偏信息）。 - 无高维保障：无任何Lasso或收缩的手，所以参数个数 \(p\) 必须显著小于总样本量。

四、开放问题（点到为止）¶

【可解释性与模型选择】 本文的线性模型 \(\mathbf{V}_i = \mathbf{Z}_i^T \boldsymbol{\gamma}\) 虽然灵活，但如何选择 \(\mathbf{Z}_i\)？例如，如果一组有10个成员，协变量矩阵 \(\mathbf{Z}_i\) 需要为成员对的每种组合都包含特征。这会引发一个严重的高维模型选择问题。目前论文没有提供任何正则化或模型选择方法。扎根于：文中所有例子都是通过明确的、人为设计好的特征（比如性别差异指示符）来定义 \(\mathbf{Z}_i\)，没有处理未知特征的情形。
【效率理论】 变换后MLE是否达到了对应假设下的半参数效率界？若否，能否用“去偏化”或“有效影响函数”来校正？这个问题对于将本文方法纳入“半参数推断”的理论框架至关重要，特别是在应用半参数理论（moderately_familiar 领域）时。扎根于：论文只提到了估计，没有讨论任何效率性质或渐近方差。
【高维设置】 当相关矩阵很大（比如组内有 \(n_i = 1000\)）时，广义z变换涉及 \(n_i \times n_i\) 的谱分解，计算复杂度是 \(O(n_i^3)\)，很快变得不可行。如何处理高维情况？能否用稀疏相关模式（如图模型诱导稀疏性）或低秩近似来逼近广义z变换，从而将计算复杂度降低到与组大小线性或多项式相关？扎根于：作者的例子中组大小是最多31名学生的班级。
【非参数推广】 能否将变换后的线性模型 \(\mathbf{V}_i = \mathbf{Z}_i^T \boldsymbol{\gamma}\) 替换为非参数回归（如高斯过程、样条）或者深度学习模型？这能极大地扩展模型的表达能力，尤其是在协变量-相关关系高度非线性时。扎根于：线性可加性是论文的核心假设之一，限制了灵活性。

Maintained by 陈星宇 · Homepage · Source on GitHub