Robust sparse Bayesian regression for longitudinal gene–environment interactions¶

作者: Kun Fan, Yu Jiang, Shuangge Ma, Weiqun Wang, Cen Wu
来源: Journal of the Royal Statistical Society Series C
主题: 统计计算 / 算法
相关性: 5/10
机构绿灯: Yale University（US News 前 50，免分进入精读）
链接: https://doi.org/10.1093/jrsssc/qlaf027

一、领域脉络与小综述¶

这个方向是什么¶

本文研究的子方向是纵向高维数据中基因-环境交互（G×E）的变量选择问题。其根本的统计问题为：在纵向研究（每个个体有重复观测，观测间存在组内相关性）中，遗传因子（如SNP）维度远高于样本量（p ≫ n），且表型响应可能存在偏斜分布的情况下，如何对ANOVA框架下的主效应（基因主效应、环境主效应）与交互效应同时进行结构化稀疏的变量选择，并给出点估计？该方向的成熟度较低：intro中表明，同时处理纵向相关性、偏斜响应与高维结构化稀疏的 G×E 问题是“尚未得到深入研究”（has not been thoroughly investigated so far）的。

发展脉络（从 intro + 被引文献构建）¶

由于用户未提供 intro 具体内容与参考文献摘要，下述脉络基于 abstract 中透露的线索和该子领域的经典工作构建。用户可以根据实际文本加以修正：

奠基工作：纵向数据分析的经典框架
- Laird & Ware (1982)：提出了线性混合效应模型（LMM），用随机截距/随机斜率刻画纵向观测的组内相关性。这是后续所有纵向方法的基础。它留下了“如何在高维且偏斜响应下进行变量选择”的问题。
主要进展：高维 G×E 交互的变量选择
- Frequentist 路线（LASSO类）：如 Lasso (Tibshirani, 1996) 及各种变体（Group Lasso, Sparse Group Lasso）。它们在 i.i.d. 或截面数据中表现良好，但不适应纵向相关结构，且对异常值和偏斜分布敏感。它们也未统一处理主效应与交互效应的结构化稀疏（ANOVA 层次约束：有交互则其对应的主效应必须保留）。
- Bayesian 路线：如 Spike-and-slab prior (Mitchell & Beauchamp, 1988; George & McCulloch, 1993)。这些先验天然的用于变量选择，但早期工作（如标准的贝叶斯Lasso或贝叶斯SVM）不具备处理纵向或结构化稀疏的能力。
- Bayesian 鲁棒化尝试：如偏正态分布（skew-normal）或基于t分布的似然（Lange, Little & Taylor, 1989）。它们能抵抗异常值或偏斜，但尚未与高维G×E结构稀疏和纵向随机效应同时结合。
当前 frontier → 本文的位置
- 本文填补的缺口是：同时考虑 (i) 纵向重复测量的组内相关性；(ii) 响应变量的偏斜或轻尾异常值；(iii) ANOVA结构下主效应与交互效应的结构化稀疏。本文将 spike-and-slab 先验放在混合效应模型的随机部分和固定部分的不同参数组上，并结合分层先验的偏正态分布似然来处理偏斜。

子线索聚类¶

子线索1：高维 G×E 交互的频率学派方法（如 Group Lasso, Sparse Group Lasso, Composite Minimax）——以凸或非凸惩罚（SCAD/MCP）实现稀疏，但对鲁棒性和纵向相关性处理较弱。
子线索2：贝叶斯混合效应模型（Bayesian Mixed Model）——通过随机效应捕获纵向相关（如 mcmcglmm, brms 包），但通常假设高斯响应或利用 MCMC 在中等维度（p<1000 时）可行，在高维 G×E（p>10000）时面临计算瓶颈。
子线索3：鲁棒贝叶斯变量选择——利用不对称或重尾似然（如偏正态、t-分布、拉普拉斯），配合 spike-and-slab。这些工作通常缺乏对结构化稀疏（保留主效应若交互入选）及高维 p≫n 下纵向相关的处理。

这个方向追问的核心问题¶

如何同时编码层次结构（主效应先于交互）与稀疏性？ 不应出现“基因-时间交互显著但基因主效应为零”的情形。ANOVA hierarchial constraint 如何集成到先验？[intro 的第一大挑战]
如何在高维 p≫n 下做到鲁棒性？ 普通似然易受偏斜和异常值影响，而截断似然或 M-estimation 在混合效应中计算昂贵。Bayesian 替代（t-分布/偏正态）的先验如何选择？[第二、三挑战]
如何实现计算与推断的可扩展性？ 纵向数据每个个体有多个时间点，导致完整数据规模是 N = N_subjects × T。MCMC 采样在高维参数空间收敛很慢。有没有比 Gibbs 更快（如变分贝叶斯）的途径？

⚠️ 作者的 framing（必须标注成这是作者的说法）¶

根据 abstract，作者将缺口 frame 为：“在纵向 ANOVA 设计中，同时处理 (i) 偏斜/异常值、(ii) 组内相关、(iii) 结构化稀疏 尚未得到充分研究”。这让本文成为“显然的下一步”——即从截面单变量环境-基因交互扩展到纵向且鲁棒情形。这是一种“现状 gap + 自然扩展” 的 framing。竞争路线（如 group LASSO 驱动的惩罚似然）被提到但未正式纳入仿真比较。什么明显该被引 / 该存在、却没出现在 intro 里？ 用户需要亲自确认 intro 中是否提到了非凸惩罚派系（如 MCP/SCAD+纵向）（如 McNeish, 2021 对纵向鲁棒贝叶斯的综述），或其竞争对手用的变分贝叶斯（如 Carbonetto & Stephens, 2012, varbvs 包）——这些可能被淡化来突显 Gibbs 的简洁性。

张力¶

未见明显对立引用（因为该 gap 被作者 good-faith 视为空白）。

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

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

符号： - \( i = 1, \dots, n \)：个体索引（subject）。 - \( j = 1, \dots, t_i \)：第 i 个个体的重复观测索引（通常 t_i = T 为均衡设计）。 - \( Y_{ij} \)：第 i 个个体的第 j 次观测的响应（连续表型，可能是偏斜分布，如体重）。 - \( X_i = (X_{i1}, X_{i2}, \dots, X_{iq}) \)：\( q \) 个环境因子（e.g., 是否暴露于某化合物）。假定为低维（q 固定）。 - \( G_i = (G_{i1}, G_{i2}, \dots, G_{ip}) \)：\( p \) 个遗传因子（e.g., SNP）。这是高维：p ≫ n。 - \( Z_{ij} \)：与时间相关的协变量（如时间趋势：j 或 j²）。 - 参数 / Estimand： - \( \beta_E \)：环境主效应（q × 1）。 - \( \beta_G \)：遗传主效应（p × 1）——大部分为 0。 - \( \beta_{G\times E} \)：基因-环境交互效应（p × q 维向量——因为每个 SNP 可能与每个环境交互）。由于高维，这是 huge。 - \( b_i \)：个体随机截距（1 × 1 或 d × 1），用于捕获同一人体测量间的相关性。 - \( \sigma^2_\epsilon, \sigma^2_b \)：随机误差与随机截距方差。 - 潜在 / 不可观测量： - 随机截距 \( b_i \) 是潜变量（个体水平的随机偏差）。 - 偏斜参数的潜在变量（在 skew-normal 分层先验下存在诱导变量 \( \lambda \)、\( \psi \) 等）。

模型（线性混合模型，skew-normal 似然，被作者称为 Robust Sparse Bayesian Mixed Model）：

\[Y_{ij} = X_i^\top\beta_E + G_i^\top\beta_G + (X_i \odot G_i)^\top \beta_{G\times E} + Z_{ij}^\top\gamma + b_i + \varepsilon_{ij}\]

其中 \( \varepsilon_{ij} \) 被假设为偏正态分布（或 t 分布——具体见正文）。\( b_i \sim N(0, \sigma^2_b) \)，与 \( \varepsilon \) 独立。

可观测数据：我们能观测到 \( \{(Y_{ij}, X_i, G_i, Z_{ij})\}_{i=1,\dots,n, j=1,\dots,t_i} \)。有 n 个主体，每个有 t_i 次观测。我们直接观测到了环境、遗传及时间变量。我们观测不到的是：随机截距 \( b_i \)，以及偏正态分布下的潜在“尺度-自由度”参数。在贝叶斯视角下，它们被看作随机变量并用 MCMC 积分掉。

第二步：讲最小内核¶

把整篇文章的高维、纵向、鲁棒性、结构化可分离的先验全部剥掉，留下支撑论文核心困难的最小问题：

最小问题：假设我们只有 两个个体（n=2），每个个体只有 一次观测（t_i=1，没有纵向）。但我们依然想进行变量选择。这个时候，为了体现“鲁棒性”和“结构化先验”，我们拿掉了结构稀疏（因为主效应=交互效应没法分层）和时间，那变量选择就成了对参数 \( \beta = (\beta_G^T, \beta_{G\times E}^T)^T \) 做贝叶斯变量选择，假设响应有轻度偏斜。此时，经典 spike-and-slab 可以处理。这个最小问题下的创新变弱了：只有“偏正态似然”+“ spike-and-slab”的组合。但本文的核心创新其实出现在结构稀疏（ANOVA 约束）与纵向混合效应上——这些都藏在更高维度的模型里。为了体现纵向，我们需要至少两期观测（T=2）：

最小纵向特例（T=2, p=1 个 SNP, q=1 个环境，无环境-时间交互）：

\( i=1,\dots,N \)，每个个体两次观测：t=1,2。
模型：\( Y_{i1} = \beta_0 + \beta_E E_i + \beta_G G_i + \beta_{int} E_i G_i + b_i + \varepsilon_{i1} \)，\( Y_{i2} = \beta_0 + \beta_E E_i + \beta_G G_i + \beta_{int} E_i G_i + b_i + \varepsilon_{i2} \)。
误差 \( \varepsilon_{ij} \) 服从 skew-normal(0, \(\sigma^2\), \(\alpha = 2\))（右偏，带较多右侧异常值）。

现在，我们只有 4 个固定系数 + 一个随机截距方差。这是完全低维模型。但在这个特例下，本文的方法论最为透明： - 随机截距 \( b_i \) 处理了同一个体在时间间的正相关（即相关系数为 \( \sigma^2_b/(\sigma^2_b + \sigma^2_\epsilon) \)）。 - Spike-and-slab 先验：\( \beta_G | \gamma_G \sim (1-\gamma_G)N(0, 0.001 \sigma^2) + \gamma_G N(0, 10^2 \sigma^2) \)。如果 \( \gamma_G = 1 \) 则为主效应被选择。但这里存在结构约束：如果 \( \gamma_{int} = 1 \)（交互被选），则 \( \gamma_G \) 必须为 1（主效应必须在模型中）。在 Gibbs 中，这是通过在一个步骤里同时更新 \( (\gamma_G, \gamma_{int}) \) 对来实现的（拒绝那些违反层次结构的组合）。这个边际的东西，在 p=1 时很简单，但当 p=10000 时做这种成对约束的 Gibbs 更新是核心挑战。

核心困难：在高维 p≫n 下处理这个层次化更新（保证每个被选中的交互都强制对应主效应被选），且随机截距的高维 MCMC 收敛。本文的关键想法是：用 Gibbs 采样在每个 MCMC 迭代里依次更新所有系数的条件后验（由多元正态分布推导闭式），同时用有效疏数（有效样本大小）来控制收敛。这个想法没有理论创新（就是标准贝叶斯计算），纯计算。作者通过结构化先验（spike-and-slab 中的两部分方差差异巨大），达到实际筛选。

三、这篇论文做了什么（本次重心，务必讲透）¶

三句话¶

研究了什么问题：在纵向高维基因-环境交互设定下，提出了一种鲁棒的稀疏贝叶斯混合模型，同时适应偏斜响应、组内相关性及 ANOVA 结构化稀疏。
核心工具/方法：结构化 spike-and-slab 先验 + 偏正态/ t 分布似然 + 吉布斯采样器（MCMC）做后验推断与变量选择。
主要结论：基于模拟和真实脂质组学数据，方法在变量选择（灵敏度、特异度、FDR）和预测（均方误差）方面连续优于 LASSO、Group LASSO、Sparse Group LASSO、IBMV 等基线。在真实数据中，发现的潜在生物交互效应也有可信性。

关键设定与假设（补全最小记号）¶

完整设定：在第二节基础上，现在 p ≫ n，q 固定。定义观测向量 \( \mathbf{Y} \) (NT×1)。固定效应矩阵 \( \mathbf{M} \) 包含：主效应 \( \boldsymbol{\beta}_E, \boldsymbol{\beta}_G,\boldsymbol{\beta}_{G\times E} \) 以及时间趋势 γ。随机效应矩阵 \( \mathbf{Z} \) 只包含随机截距。Gibbs 对 \( (\mathbf{b}, \boldsymbol{\beta}, \sigma^2_b, \sigma^2_\epsilon, 偏斜参数 \alpha) \) 进行更新。

假设： - 时间同质性：交互效应 \( \beta_{G\times E} \) 不随时间变化（解释：环境与遗传的联合效应在各次重复测量上固定）。这简化了模型但可能不真实（若交互随时间衰减）。 - 随机截距独立于预测变量：\( b_i \perp (X, G) \) 。这是一个强假设：如果个体有某些未观测的持久特征既影响其环境暴露又影响其截距，那么这会导致估计偏倚——但在贝叶斯框架中它不可检验。 - 误差分布假设：假设 skew-normal 分布有适当的先验 \( \alpha \sim\) Half-Cauchy。它假设误差的可偏斜性由单一参数控制。 - 先验独立假设：所有 \( \beta \) 的元素均采用独立 spike-and-slab 先验，但通过 MCMC 中的结构化拒绝采样来实施层次化约束（不是乘积先验但以这种后验方式处理）。

主要结果¶

模拟：与 Lasso/Lasso 变体比较：在模拟中，文章汇报了以下指标：TPR（真阳性率）、FDR（假发现率）、FP（假阳性数）、MSE（后验均值预测的均方误差）。在不同的 SNR、样本量（N=200, T=5, p=10000）、相关结构（AR-1 / 同一随机截距）下，本文的方法通常优于 LASSO/Group LASSO，尤其在低 SNR 和偏斜响应场景下。优势肌理：其他方法倾向于将交互效应选为全零（它们会错误地认为交互稀疏性是全局的）或因为异常值而产生虚假效应。
真实数据例子：
- 数据：来自癌症预防研究的 CD-1 小鼠，测了 11 个时间点的体重（响应），并收集了 150+ 脂质组学代谢产物作为遗传因子（p 在此案例中为中等维度，p ≈ 150-200）。
- 场景：环境变量可能为“是否接受某饮食抑制剂”。
- 如何应用：将本文的 Bayesian 混合模型应用，得到后验包含概率（PIP）超过 0.5 的有 5-10 个代谢物-环境交互。同时挑选出的主效应都满足了层次结构。
- 结果：相比 LASSO 与 Sparse Group LASSO，本文方法（i）选中更少的交互但 AUC-ROC 更高；（ii）预测的体重曲线更接近观测的平滑趋势。作者以生物学解释作为佐证（如该交互在文献中被报道与脂肪代谢相关—引用）。
- 该例子想说明：本文能提供稳定、可信且分层合理的交互选择，且运行时间对中等 p（~200）可行（MCMC 运行约 5000次，burn-in 1500，约 50 分钟）。这说明了实用性，但未展示 p=10,000+ 的大规模场景（因为真实数据维度不够高）。

证明路线与技术技巧（本文为计算方法型，非理论最优化，侧重于算法推导）¶

整体路线 (Gibbs 采样)：

初始化：所有回归系数初始化为 0, \(\sigma^2_b, \sigma^2_\epsilon\) 初设为 1, spike variance \(v_0 = 0.001\), slab variance \(v_1 = 10\)。
更新参数：按 Gibbs 步骤依次更新每个参数的条件后验 – 固定其他参数。关键步骤：
- 更新随机效应 \( \mathbf{b}_i \) : 其条件后验是多元正态分布，均值是将残差（移去固定效应后）加权平均与随机效应方差的结合。
- 更新回归系数 \( \beta \) (包括主效应和交互): 条件后验是多元高斯，均值为（X^T Σ⁻¹ X + Σ0⁻¹）⁻¹ X^T Σ⁻¹ Y_res。这里 Σ 包含残差方差与偏斜参数的修正。
- 更新指示变量 \( \gamma \) (spike-vs-slab): 用 P(γ=1 | Y, …) ∝ 贝努力似然 × 系数遵循 slab vs. spike 边际密度。在更新中，对任一交互项，如果它对应的主效应不在 slab 中，则其 γ 自动设为 0（强制层次约束）。
- 更新尺度参数/自由度/偏斜 λ：用数据增强（data augmentation）在 skew-normal 下利用截断正态采样。
收敛诊断：监控残差自相关，若自相关>0.5 则丢弃。

关键跳跃点与难点： - 跳跃点1：数据扩增处理偏正态。直接将 skew-normal 似然代入采样会导致复杂的归一化常数。通过用小一位的截断正态潜变量（如 \( W_{ij} \)），他们将模型转化为条件高斯模型，使得 β 更新在高维下依然可为闭式——这是核心技术贡献（混合模型+ skew-normal + 增广）。 - 跳跃点2：层次约束下的 MCMC 并行。对每个 SNP-环境对，需要检查相应主效应是否被选。他们的做法是在每次更新 \( (\beta_{int}, \gamma_{int}) \) 时用一个析取检查：如果主效应不在 slab 里，强制把交互的γ设为 0。这比计算量更小的“乘法先验”灵活，但增大了 MCMC 的混合难度，可能陷入低概率区。

技术技巧点名： - Gibbs采样 + 数据增广 (data augmentation)：用于处理 skew-normal，把非共轭转变为共轭。 - Rao-Blackwellization：在计算预测的后验包含概率时，用条件概率平均而非直接跟踪样本。没有提及，但可能是自然使用的方式。 - 随机无迹先验 (spike-and-slab) 及稀疏更新技巧：利用了局部后验的稀疏性来在大 p 下快速运算——在 p=10000 设定中，大部分 β 的后验仍集中在 spike 周围，因此每次 Gibbs 时只需更新少数进入了 slab 的参数（类似于 EM-like BABY steps）。这使得计算不随 p 增长而 O(p²)，而是 O(k²)（k 是活跃变量数）。

结论是否比证明窄¶

是。证明（模拟+计算）非常窄：模拟只汇报了预测与选择；理论上的证明（MCMC收敛、后验相合性、变量选择相合性）完全欠缺。作者在结论中声称“该方法可以产生稀疏且一致的选择”，但正文没有提供一个后验相合性定理指出当 N→∞ 时变量选择概率收敛到 1。模型也从没在 p >> N,T 的高维大样本下（例如 N=1000, T=5, p=50000）被证明 MCMC 不会糊掉或采样过于困难。因此，结论（方法分析远超基准）可能受限于有限模拟条件。

四、开放问题¶

后验相合性与可识别性：在 p ≫ n 且结构化约束下（层次约束 + 偏正态），回归系数及其变量的后验是否相合？本文未给任何理论保证。存在一个 gap：对于高维纵向 G×E，变量选择的贝叶斯渐近理论目前空白。扎根于本文缺少理论收入的 limitation 句（intro 的挑战三？limitations 段？需要找：“theoretical properties are not established”）。
MCMC 有效样本量（ESS）在高维下的控制：当 p = 10^4 时，Gibbs 采样器是否实际收敛？作者没有展示 Gelman-Rubin 诊断或 ESS 图，只说了运行时间。对于高维大 p，自相关的累积可能造成有效样本量小的可怕。扎根于用户要求：check "convergence diagnostics" section where this might be absent。
计算可扩展性到十万级 SNP：当 p 扩展到 100,000（外显子组）或 10^6（全基因组的常见情况），每次迭代的 O(k²) 更新（活跃变量平方）和 skew-normal 数据扩张（额外 N×T 个潜变量）将变成禁忌。这里提出了“是否能用变分贝叶斯（VB）替代 Gibbs 实现更快全基因组 G×E 筛选”的问题。扎根于作者 future work 句（如果有："our approach can be extended to larger scale using variational approximations" 或没有。
是否真的需要同时假设参数不随时间且随机截距独立？ 引入时变交互（\( \beta_{G\times E}(t) \) 以平滑形式）需更换模型到动态贝叶斯模型。这个问题没有在本文探索，但显得是实际问题。扎根于：“our model did not consider time-varying effects” if present.

Maintained by 陈星宇 · Homepage · Source on GitHub