Tensor regression for incomplete observations with application to longitudinal studies¶
作者: Tianchen Xu, Kun Chen, Gen Li
来源: Annals of Applied Statistics
主题: 高维统计 / 随机矩阵
相关性: 5/10
机构绿灯: University of Michigan(US News 前 50,免分进入精读)
链接: https://doi.org/10.1214/23-aoas1830
一、领域脉络与小综述¶
这个方向是什么¶
这个子方向处理的是高维纵向数据在块缺失下的回归与变量选择问题。核心统计挑战在于:当数据具有"样本×特征×时间"的三阶张量结构,且并非所有样本都在所有时间点被观测(块缺失),如何在不插补的前提下,利用所有观测信息进行参数估计与变量选择。该方向目前处于方法成熟期与理论深化期,已有针对张量回归的完备理论与算法,但对缺失机制(尤其是块缺失)的处理仍在不断精细化。
发展脉络¶
根据 introduction 的引用梳理,该领域的发展线索如下:
-
奠基工作(张量回归与降维):
- Zhou et al. (2013) 与 Li et al. (2013):建立了标量-张量回归的基本框架,引入了低秩分解与正则化方法,为高维张量系数估计提供了理论基础。但早期工作多假设数据完整。
- Hung & Wang (2013):提出了针对纵向数据的边际模型,但主要处理单变量或低维情形,未触及高维特征与张量结构。
-
主要进展(缺失数据与纵向结构):
- Xiao et al. (2018):专门处理纵向数据中的块缺失问题,提出了基于边际矩的估计方法,但主要聚焦于均值结构,未充分利用张量协方差结构。
- Luo et al. (2016) 与 Cai & Zhou (2013):在高维协方差矩阵估计与缺失数据补补方面做出了重要贡献,提供了理论上界与极小极大最优性结果,但未直接针对张量回归设定。
-
当前 Frontier(张量缺失与低秩结构):
- Cai et al. (2019) 与 Chen et al. (2019):开始处理张量数据的缺失值问题,利用低秩 Tucker 分解进行填补或估计,但多集中于图像填补或无监督学习,较少涉及回归中的变量选择。
-
本文的位置:
- 本文位于"张量回归"与"块缺失纵向数据"的交叉点。作者指出,现有方法要么要求完整数据,要么无法同时处理高维特征选择与时变效应平滑性。本文试图在不插补的前提下,通过协方差结构直接利用块缺失数据。
子线索聚类¶
被引文献大致落在三条子线索上: - 线索一:张量回归与降维(Zhou et al. 2013, Li et al. 2013, Guhaniyogi et al. 2017):关注如何通过低秩结构降低参数维度,主要处理完整数据。 - 线索二:纵向数据与缺失机制(Xiao et al. 2018, Luo et al. 2016, Ibrahim et al. 2001):关注纵向数据的相关性与缺失值的似然推断,但多限于低维或单变量响应。 - 线索三:高维变量选择与协方差估计(Bickel & Levina 2008, Rothman et al. 2008, Cai & Zhou 2013):提供高维正则化与协方差估计的技术工具,是本文方法的理论基石。
这个方向在追问的核心问题¶
- 如何定义与处理块缺失:传统的逐点缺失假设过于理想,块缺失(某些样本在某些时间点完全缺失)更符合实际,但也带来信息损失与计算困难。
- 如何在缺失下保持估计效率:不插补而直接利用观测矩,能否达到与完整数据相近的估计效率?正则化参数如何选择?
- 如何同时实现变量选择与平滑估计:在特征维度 \(p \gg n\) 且时间效应非参数时,如何保证系数的稀疏性与时间效应的平滑性?
⚠️ 作者的 framing¶
- 作者的说法:作者将缺口 frame 为"现有方法无法同时处理块缺失、高维特征选择与时变效应平滑"。他们强调自己的方法"无需插补"、"利用所有观测"、"同时实现选择与平滑"。
- 被淡化的竞争路线:作者未深入讨论多重插补或似然方法在块缺失下的表现,也未与基于深度学习的填补方法对比。此外,对于缺失机制是否必须 MCAR(完全随机缺失),作者在理论部分假设了 MCAR,但在应用部分未做严格检验,这一点被淡化。
- 缺失的引用:Introduction 中未引用因果推断中处理缺失数据的最新进展(如 Proximal Causal Inference 中的 missing data 处理),也未提及矩阵补全领域的最新低秩填补理论。这可能是作者刻意聚焦于回归设定,但也留下了"因果视角下的块缺失"这一缺口。
张力¶
被引文献之间未见明显对立结论。大多数工作都在各自假设下推进,如 Xiao et al. (2018) 处理块缺失但未做高维选择,Zhou et al. (2013) 做高维张量回归但未处理缺失。本文试图填补这一空白,未见直接冲突。
二、最核心、最简单的例子 / 数学问题¶
第一步:符号、模型、可观测数据¶
在展开论文细节前,先立清楚记号与设定:
-
符号:
- \(n\):样本量。
- \(p\):特征维度(高维,\(p \gg n\))。
- \(T\):时间点数。
- \(Y_i \in \mathbb{R}\):第 \(i\) 个样本的标量响应变量(如神经发育得分)。
- \(\mathcal{X}_i \in \mathbb{R}^{p \times T}\):第 \(i\) 个样本的特征张量(特征×时间),即纵向数据矩阵。
- \(\mathcal{B} \in \mathbb{R}^{p \times T}\):待估系数矩阵(张量),刻画特征在不同时间的效应。
- \(\mathcal{M}_i \in \{0,1\}^{p \times T}\):缺失指示矩阵,\(1\) 表示观测到,\(0\) 表示缺失。本文主要考虑块缺失,即 \(\mathcal{M}_i\) 的某些时间列全为 \(0\)。
-
模型:
- 标量-张量回归模型:
\[Y_i = \langle \mathcal{X}_i, \mathcal{B} \rangle + \epsilon_i = \text{vec}(\mathcal{X}_i)^\top \text{vec}(\mathcal{B}) + \epsilon_i\]其中 \(\langle \cdot, \cdot \rangle\) 为张量内积,\(\epsilon_i\) 为噪声。目标是估计 \(\mathcal{B}\)。
- 结构假设:
- 稀疏性:\(\mathcal{B}\) 只有少部分行(特征)非零,即行稀疏。
- 平滑性:\(\mathcal{B}\) 的每一行(时间效应)随时间平滑变化。
- 低秩性(隐含):通过正则化隐式利用。
- 标量-张量回归模型:
-
可观测数据:
- 研究者能观测到的是 \(\{ (Y_i, \mathcal{X}_i^{\text{obs}}) \}_{i=1}^n\),其中 \(\mathcal{X}_i^{\text{obs}}\) 是 \(\mathcal{X}_i\) 在 \(\mathcal{M}_i=1\) 处的观测值。
- 观测不到的是 \(\mathcal{X}_i\) 在 \(\mathcal{M}_i=0\) 处的缺失值。
- 关键点:本文假设缺失机制为 MCAR(Missing Completely At Random),即 \(\mathcal{M}_i\) 与 \((Y_i, \mathcal{X}_i)\) 独立。这是后续所有无偏估计的基础。
第二步:最小内核¶
为了讲清核心思路,我们考虑一个最简特例:\(p=1\)(单特征)、\(T=2\)(两时间点)、无噪声。
-
问题退化: 此时 \(\mathcal{B} = [\beta_1, \beta_2]^\top\) 是一个二维向量。模型变为 \(Y_i = X_{i1}\beta_1 + X_{i2}\beta_2\)。 假设样本分为两组:
- 组 1(\(n_1\) 个样本):两个时间点都观测到,有 \((X_{i1}, X_{i2}, Y_i)\)。
- 组 2(\(n_2\) 个样本):只有时间点 1 观测到,时间点 2 缺失(块缺失),有 \((X_{i1}, Y_i)\)。
-
传统方法的困境: 如果删除组 2,损失信息;如果插补组 2 的 \(X_{i2}\),引入偏差。 若直接对组 1 做 OLS,得到 \(\hat{\beta}^{(1)}\),方差大(样本少)。 若只用组 2 的 \(X_{i1}\) 与 \(Y_i\),只能估 \(\beta_1\),无法估 \(\beta_2\)。
-
本文的核心想法(最小内核): 利用协方差结构把所有样本"拼"起来。
- 矩条件:
- 对组 1(完整样本):\(E[Y_i X_{i1}] = \beta_1 E[X_{i1}^2] + \beta_2 E[X_{i1}X_{i2}]\)。
- 对组 1 + 组 2(所有样本):\(E[Y_i X_{i1}]\) 可以用全部 \(n_1+n_2\) 个样本精确估计。
- 对组 1(完整样本):\(E[X_{i1}^2]\) 和 \(E[X_{i1}X_{i2}]\) 可以用 \(n_1\) 个样本估计。
- 方程组: 我们有 \(n_1+n_2\) 个样本估计的 \(E[Y_i X_{i1}]\)(左边),和 \(n_1\) 个样本估计的协方差项(右边系数)。这构成了一个关于 \((\beta_1, \beta_2)\) 的线性方程组。
- 推广到高维: 当 \(p \gg n\) 时,这个方程组是欠定的。本文引入正则化(Lasso 型惩罚 + 平滑样条惩罚),在约束方程下寻找稀疏且平滑的解。
- 矩条件:
-
为什么这个内核重要: 它揭示了本文方法的本质:不插补数据,而是插补矩方程。通过构造基于协方差的估计方程,利用不同缺失模式下可用的样本片段,拼凑出完整的估计量。这避免了插补模型设定错误的风险,同时最大化了样本利用率。
三、这篇论文做了什么¶
三句话¶
- 研究了高维纵向数据在块缺失模式下的标量-张量回归问题。
- 核心方法是提出了基于协方差的正则化估计方程,无需数据插补即可利用所有观测样本。
- 主要结论是在 MCAR 假设下,证明了估计量的变量选择一致性,并通过微生物组数据验证了方法的有效性。
关键设定与假设¶
在第二节最小记号的基础上,补全完整设定:
- 块缺失模式: 定义 \(\mathcal{O}_k \subseteq \{1, \dots, T\}\) 为第 \(k\) 类缺失模式(观测到的时间点集合)。假设共有 \(K\) 种缺失模式,样本按缺失模式分为 \(K\) 组。
- 协方差估计方程:
定义 \(\Sigma_{t,t'} = E[X_{it} X_{it'}]\)(协方差矩阵),\(g_t = E[Y_i X_{it}]\)(响应与特征的协方差)。
模型 \(Y_i = \sum_t X_{it} \beta_t + \epsilon_i\) 两边同乘 \(X_{it}\) 并取期望,得:
\[g_t = \sum_{t'} \Sigma_{t,t'} \beta_{t'}\]这是一个线性系统。关键在于,\(g_t\) 和 \(\Sigma_{t,t'}\) 可以利用不同缺失模式下的样本进行估计:
- 若样本在 \(t\) 时刻有观测,则可用于估计 \(g_t\)。
- 若样本在 \(t, t'\) 时刻都有观测,则可用于估计 \(\Sigma_{t,t'}\)。
- 正则化目标函数:
\[\min_{\mathcal{B}} \sum_{k=1}^K w_k \| \hat{g}^{(k)} - \hat{\Sigma}^{(k)} \text{vec}(\mathcal{B}) \|_2^2 + \lambda_1 P_{\text{sparse}}(\mathcal{B}) + \lambda_2 P_{\text{smooth}}(\mathcal{B})\]其中:
- \(\hat{g}^{(k)}, \hat{\Sigma}^{(k)}\) 是基于第 \(k\) 类缺失模式样本的矩估计。
- \(P_{\text{sparse}}\):组 Lasso 惩罚,实现行稀疏(变量选择)。
- \(P_{\text{smooth}}\):融合 Lasso 或样条惩罚,实现时间效应平滑。
- 假设条件:
- MCAR:缺失指示变量独立于潜在数据。
- 子高斯性:特征与误差项服从次高斯分布,保证矩估计的尾概率界。
- 限制特征值条件:设计矩阵需满足稀疏特征值条件,保证 Lasso 的变量选择一致性。
主要结果¶
- 定理 1(估计一致性):
在适当选择正则化参数 \(\lambda_1, \lambda_2\) 下,估计量 \(\hat{\mathcal{B}}\) 与真值 \(\mathcal{B}^*\) 的误差满足:
\[\|\hat{\mathcal{B}} - \mathcal{B}^*\|_F \leq C \sqrt{\frac{s \log(pT)}{n}}\]其中 \(s\) 为非零特征数。这达到了与完整数据下相近的收敛率(仅多出 \(\log T\) 因子)。
- 定理 2(变量选择一致性): 在 irrepresentable condition 或更弱的限制等距条件下,以概率趋于 1,\(\hat{\mathcal{B}}\) 能正确识别非零特征行。
- 技术难点: 难点在于块缺失导致矩估计的方差结构复杂。不同缺失模式下的样本量不同,协方差估计 \(\hat{\Sigma}^{(k)}\) 的精度不同,如何加权组合这些矩方程是关键。作者通过权重 \(w_k\)(通常取各组样本量比例)来平衡方差。
证明路线与技术技巧¶
-
整体路线:
- 矩估计的集中不等式:首先证明 \(\hat{g}^{(k)}\) 和 \(\hat{\Sigma}^{(k)}\) 在次高斯假设下以高概率接近真值。
- 基本不等式:利用 Lasso 的基本不等式,将估计误差转化为经验过程的界。
- 限制特征值条件验证:证明在块缺失下,组合后的设计矩阵仍满足稀疏特征值条件。
- 变量选择:利用 KKT 条件与阈值化技术,证明非零系数的符号一致性。
-
关键跳跃点:
- 引理 1(协方差估计的界):在块缺失下,协方差矩阵估计的误差界不仅依赖于总样本量 \(n\),还依赖于最小缺失模式样本量 \(n_{\min}\)。这是理论中最吃劲的地方——如果某类缺失模式样本太少,对应的协方差估计方差会爆炸,进而拖累整体估计。
- 处理相关性:时间点之间的相关性 \(\Sigma_{t,t'}\) 需要被估计,且估计误差会传导至回归系数。作者通过分块矩阵分析,将时间维度的相关性"解耦",转化为一系列并行的估计问题。
-
技术技巧点名:
- 次高斯尾概率:用于控制高维矩估计的偏差。
- 组 Lasso(Group Lasso):用于实现行稀疏,这是高维变量选择的标准工具。
- 平滑样条惩罚:通过差分惩罚 \(P_{\text{smooth}}(\mathcal{B}) = \sum_{j,t} (\beta_{j,t+1} - \beta_{j,t})^2\) 实现时间平滑。
- 交替方向乘子法(ADMM):用于求解带非光滑惩罚的优化问题。
真实例子与应用¶
- 数据 / 场景:
早产生儿肠道微生物组与神经发育的纵向研究。样本量 \(n \approx 100\),特征 \(p \approx 150\)(微生物属),时间点 \(T=4\)(第 1, 2, 3, 4 周)。响应变量为 2 岁时的神经发育评分。
- 缺失情况:并非所有婴儿在所有时间点都提供了粪便样本,存在块缺失。
- 怎么用: 将微生物丰度数据构造为 \(p \times T\) 的张量,应用本文方法。
- 结果: 识别出若干与神经发育显著相关的微生物属,且其效应随时间变化呈现平滑趋势。相比删除缺失样本的方法,本文利用了更多信息,发现了更多潜在关联。
- 想说明什么: 验证方法在真实数据上的可行性,展示其发现时变效应的能力,并暗示块缺失处理的重要性。
🔎 结论是否比证明窄¶
- MCAR 假设的依赖:理论结果严格依赖于 MCAR 假设。若缺失机制为 MAR(Missing At Random,依赖观测数据)或 MNAR(非随机缺失),矩方程的无偏性不再成立,结论可能失效。作者在文中提及了这一点,但未给出 MAR 下的修正方案。
- 变量选择条件的严格性:变量选择一致性需要 irrepresentable condition,这在相关性强时很难满足。作者在模拟中展示了方法在较弱条件下的表现,但理论证明仍需较强假设。
四、开放问题¶
- MAR 下的推广:本文理论严格依赖 MCAR。若缺失机制为 MAR(缺失依赖观测到的协变量),矩方程 \(E[Y X]\) 不再无偏。能否引入逆概率加权(IPW)或增强逆概率加权(AIPW)来修正矩方程,在 MAR 下实现无偏估计?这扎根于 Introduction 中对"缺失机制"的讨论,以及因果推断中处理缺失数据的经典工具。
- 协方差结构的鲁棒性:方法依赖于对协方差矩阵 \(\Sigma\) 的估计。在高维 \(p \gg n\) 下,协方差矩阵估计本身就很困难。能否引入更鲁棒的协方差估计子(如 thresholding 或 low-rank approximation),并证明其对最终回归估计的影响?扎根于定理证明中对协方差估计误差界的依赖。
- 因果效应的解释:在微生物组例子中,作者将关联解释为"link"。但若存在时变混杂(如饮食、用药),回归系数 \(\mathcal{B}\) 并不具因果解释。能否将此框架扩展至边际结构模型或G-估计,在纵向缺失数据下进行因果推断?这扎根于您对因果推断的兴趣,以及文中对"association analysis"的定位。
- 计算效率的优化:当 \(p\) 极大(如微生物种水平,\(p > 1000\))且 \(T\) 较大时,ADMM 算法可能收敛较慢。能否利用张量结构设计更高效的算法(如随机梯度下降或随机化 SVD)?扎根于您对统计计算的兴趣。
Maintained by 陈星宇 · Homepage · Source on GitHub