Retrospective varying coefficient association analysis of longitudinal binary traits: Application to the identification of genetic loci associated with hypertension¶

作者: Gang Xu, Amei Amei, Weimiao Wu, Yunqing Liu, Linchuan Shen et al.
来源: Annals of Applied Statistics
主题: 数理统计 / 假设检验
相关性: 7/10
机构绿灯: Yale School of Public Health（US News 前 50，免分进入精读）
链接: https://doi.org/10.1214/23-aoas1798

一、领域脉络与小综述¶

这个方向是什么¶

这是一个纵向遗传关联分析的子方向，核心统计问题是：在含有重复测量（纵向）的二值表型（如患病/不患病）的全基因组关联研究（GWAS）中，如何检验一个遗传变异（如SNP）对疾病风险的效应是否随时间变化，并给出一个稳健的显著性检验。当前该方向的主流方法大多假设遗传效应的强度不随时间改变（constant-effect），忽略了疾病进展过程中的动态模式。这篇论文试图填补对纵向二值表型的时变遗传效应检验这一缺口，并强调在非随机采样（如case-control ascertainment）下的type I error控制。

发展脉络（从introduction中的引用构建）¶

奠基工作：将时间视为"不变的变异贡献"
- 早期纵向GWAS方法：这批工作把个体的多次观测视为独立的重复测量或随机截距的扩展。例如，(Wu et al., 2011) 提出的 FAST-LMM 和 (Lippert et al., 2011) 的 FaST-LMM 利用线性混合模型（LMM）处理群体结构。它们假设遗传效应是跨时间不变的常数。
- 词语判断（作者原话）：作者提及此类方法“fail to capture the dynamic pattern of disease progression”（无法捕捉疾病进展的动态模式）。
主要进展：对纵向（连续）表型引入时变效应
- (Fan et al., 2012, 2013) 将变系数模型（varying coefficient model）引入遗传学，用B-splines或smoothing splines对连续表型下的时变遗传效应进行建模。这是关键的模型框架源头。
- (Sitlani et al., 2015) 在Fan et al.的基础上，在logistic回归框架下引入纵向二值表型的变系数模型，但使用边际模型（marginal model）并通过广义估计方程（GEE）进行推断。
- (Das et al., 2016) 提出变系数混合模型（varying coefficient mixed model, VCMM），将随机效应（random effects）纳入变系数框架，用于处理纵向连续表型。
- 词语判断（作者原话）：作者指出这些工作在建模上取得了进展，但它们的统计推断方法（如似然比检验、Wald检验）对随机效应分布（random effects distribution）的误设高度敏感，尤其当受试者并非随机抽样（如病例-对照设计）时，会严重inflate type I error。
当前frontier：对模型误设稳健的检验
- (Crowley et al., 2021) 提出 CR-VCMM (Cauchy regression VCMM)，一种针对纵向二值表型的变系数方法与稳健推断，但作者指出该方法计算代价高。
- (Chen et al., 2019) 在回顾式（retrospective）框架下，针对连续表型提出了一种对混合模型误设稳健的方差成分检验。这是当前关键步骤，但尚未推广至二值表型或变系数模型。
- (Jiang et al., 2020) 提出的 RVFT (retrospective varying coefficient test for survival traits)，将回顾式方法应用于生存数据（survival data）。
- 词语判断（作者原话）：本文作者的关键评论是：现有针对纵向表型的时变效应检验，要么只适用于连续表型（FAST-LMM, RVFT），要么在case-control或无偏抽样下对模型误设敏感（基于PQL的似然比检验），要么计算过于昂贵（CR-VCMM）。他们的定位：提出一种计算高效、在case-control采样下对混合模型随机效应分布误设稳健的变系数关联检验方法。

子线索聚类¶

将这些被引文献大致分为以下线索： 1. 线索A（主流框架）：基于线性混合模型的常效应检验 * 代表：FAST-LMM (Wu et al., 2011; Lippert et al., 2011)，SAIGE (Zhou et al., 2018) —— （本文未直接引用SAIGE，但GWAS领域有提及）。它们通过GLS或Firth惩罚对群体结构和病例-控制不平衡进行校正，但假设效应不随时间变化。 2. 线索B（建模改进）：变系数模型与纵向二值表型 * 代表：Fan et al., 2012; Sitlani et al., 2015; Das et al., 2016. 这些工作集中于模型构造（如何用spline建模时变效应），但在稳健推断方面薄弱——通常依赖对随机效应分布的正确指定，这在二值表型、非随机抽样的GWAS中极难保证。 3. 线索C（推断改进）：回顾式稳健检验 * 代表：Chen et al., 2019 (连续表型); Jiang et al., 2020 (生存数据); Crowley et al., 2021 (二值表型，但计算昂贵)。这些工作在推断层面引入稳健性——通过条件于观测到的随机效应（retrospective approach）来消除对模型误设的敏感性。本文是这条线索在纵向二值表型上的直接延伸。

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

如何有power地检验时变效应：给定有限的时间和样本量，当真实效应随时间波动（而非常数）时，如何设计一个检验统计量使其比常效应检验具有更高的统计效力（power）？
如何保持对模型误设的稳健性：在GWAS大规模筛选中，无法对每个SNP完美选定随机效应的分布或协方差结构。检验能否在随机效应分布（如个体间异质性）被错误指定时仍能正确控制type I error？
如何处理计算可扩展性：GWAS需要同时对大量遗传标记（几百万个SNP）进行检验，方法计算复杂度必须接近线性。spline基函数引入的参数会显著增加计算量，尤其对于二值表型，迭代求解的准似然方法可能更慢。

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

他们把自己的工作frame成“显然的下一步”：在处理纵向二值表型的时变效应检验这一问题上，已有建模工具（变系数）和稳健推断技术（回顾式检验），但从未被成功结合。尤其是针对二值表型的回顾式方差成分检验，之前只用于连续表型或生存数据。他们声称填补了这一空白，给出了一个计算可行、稳健且power高的具体方法。
他们淡化/回避的竞争路线：
- GEE-based sandwich variance：introduction中明确提及Sitlani et al. (2015)的GEE方法，但迅速以“对随机效应分布敏感”为由带过。实际上GEE（特别是使用独立或可交换相关结构时）在正确指定均值模型下通常是稳健的。他们未深入讨论GEE与retrospective方法在特定上的相对优劣。
- Cauchy combination test 的计算易用性：他们使用了Cauchy组合来对多个时间点的效应进行联合检验（Liu & Xie, 2020），这是一种“天生满足计算简便”的p值组合方法。他们没有讨论其他p值组合法（如Fisher's或Min-p）为何在此场景下不佳。
什么明显该被引/该存在、却没出现在intro里？
- Freeman et al. (2008, 2020) 关于score test的经典理论，尤其是对GLMM（Generalized Linear Mixed Models）的测试。这是统计测试方法的核心部分。
- (Zhou et al., 2018) SAIGE —— 这是现在处理GWAS中case-control imbalance的最主流方法，它采用Firth logistic regression + SPA（saddlepoint approximation）来控制type I error。本文的题设场景（纵向二值表型）应当与SAIGE方法进行横向对比（至少在关键模拟中），他们的模拟设计未包含这种比较，这是一个潜在的方法论缺口。

张力¶

未见明显对立引用：所有被引工作基本沿着“时变效应建模 → 纵向二值复杂推断 → 稳健性检验”这一方向前推，没有工作在同设定下得出矛盾的结论（如constant-effect方法更好）。在model misspecification下的type I error泄漏问题上，所有相关引用（Fan, Sitlani, Chen）观点一致：混合模型下传统似然比/Wald检验是不稳健的。没有跨线索的根本矛盾。

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

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

设定： - N 个体总数, i = 1,...,N。 - 对第 i 个个体，有 n_i 次重复观测（总观测数 M = Σ n_i）。 - 事件索引：第 i 个体的第 j 次观测，j = 1,...,n_i。

符号（逐个点名）： - Y_ij（可观测）: 第 i 个体在第 j 次观测时的二值表型（0 = 未患病, 1 = 患病）。这是响应变量（response）。 - t_ij（可观测）: 第 i 个体第 j 次观测的时间（如年龄）。这是变系数模型的索引变量（index variable）。 - X_ij（可观测）: (p×1) 协变量向量，包含哪些变量？严格假设哪些是已知。含：被检验的遗传变异（SNP）g_i，以及其他固定效应协变量（如性别、年龄主效应、PCs等）。注意：g_i 是相同的一个值跨越所有时间点（SNP不随时间变）。 - 模型参数： - γ(t_ij)（待估）: 遗传变异的时变效应，它是一个函数，定义了SNPg_i在时间t_ij时对二值表型贡献的ln-odds单位。这是感兴趣的参数（parameter of interest）。 - β （待估）: (p-1)×1 向量，X_ij其余（非SNP）固定效应的系数，假设恒定。 - 随机效应 b_i（不可观测）: 一个标量随机效应，捕捉了个体i之间随时间变化的共享残差（如人群中的基线风险）。其分布假设为 Gaussian: b_i ~ N(0, σ_b²)。 - 核函数 K(t, t')：定义在时间轴上的一个正定核，用以建模个体内（within-subject）误差的时间相关性。

模型（数据生成机制）：对于个体i在时间 j点 (i,j)，人们假设：

logit[ P(Y_ij = 1) ] = X_ij^T β + γ(t_ij) * g_i   + b_i + e_ij

其中 e_ij 为残差，其方差成分由核函数 K(t, t') 建模。这里的思路类似于一个变系数半参数广义线性混合模型（VC-GLMM）：效应的时变部分只出现在感兴趣的遗传变异 g_i 上。

可观测与不可观测： - 可观测: {Y_ij, t_ij, X_ij, g_i}_{i=1}^{N, j=1}^{n_i}。研究者拥有这些样本。 - 潜在/不可观测，只能靠假设去识别： - 个体间随机效应 b_i 的实现值（random realizations）。 - 时间相关性结构 K(t, t') 的基（不是参数λ）；通常靠假设，或通过估计λ获得，但本文假设其在备选假设下能为复杂结构建模。 - SNP时间效应函数 γ(t) 本身的形状（这正是要检验和推断的对象）。

第二步：讲最小内核（最小特例）¶

本文的整个方法本质上是一个从“常效应检验”到“时变效应检验”的推广，且涉及对随机效应分布误设的稳健性。其最小内核可以用对一个SNP的时变效应是否存在做假设检验来理解。

最简特例：假设我们关注单一SNP g_i，且仅有两个时间点： t = 0 和 t = 1。 - 每个个体在这两个时间点都有观测 Y_i0 和 Y_i1。 - 固定效应协变量只包含一个截距：X_ij = [1]。 - 随机效应结构为简单随机截距：b_i ~ N(0, σ_b²)。 - 核心问题：SNP的效应在t=0时为 γ0，在 t=1 时为 γ1。检验 γ0 = γ1 或 γ0 = γ1 = 0 就是一个“时变效应是否存在”的问题。 - 本文的检验思想（在此特例下的解释）: 1. 建模：我们用一个简单的线性函数来近似γ(t)，比如γ(t) = φ0 + φ1 * t。 (粗糙地说，这里φ1捕捉变化率)。在体内方法中，他们用的更高级的spline，但思维一样。 2. 标准思想：为检验 φ1 = 0，如果模型正确，一个Wald-test可以立即给出。麻烦在于，当随机截距分布的方差不确切（如真实数据不是正态的，而是更复杂的分布），标准Wald检验会膨胀type I error。 3. 回顾式（retrospective）思想：不再试图用全似然（full likelihood）去估计σ_b²，而是将分析条件于观测到的数据结构的随机性。 4. 挑战（为什么它特别？）：相比于标准广义估计方程（GEE），他们的回顾式检验严格针对方差成分（即检验时变效应，属于协方差参数）。该检验统计量被构造为： - 一种得分统计量（score statistic），在原假设下（γ(t) = 0，即无SNP效应），回归掉所有固定效应后，检验剩余变异（residual variation）与g_i的交互是否随t变化。 - 其关键思想在于：在原假设下，检验统计量的均值和方差的估计可以基于个体间的经验方差得到，而不是依赖于对随机效应b_i分布的参数化假设。这样做能获得对误设的稳健性。这就是retrospective一词的核心。

一句话总结最小内核：本文方法把“检验SNP的效应是否随时间变化”这个问题，转化为“检验SNP贡献的残差项是否随时间变化”，并且在计算检验统计量的分布时，利用个体间观测信息的经验二值方差，从而消除了对群体中不可观测的随机效应分布的错误假设所导致的inflation。

三、这篇论文做了什么¶

三句话¶

研究了什么问题：提出一个统计检验方法（RVMMAT），用于在纵向二值表型的GWAS中，检验单个SNP对疾病风险是否存在随时间变化的效应。
核心工具/方法：将变系数混合模型（VCMM）与smoothing splines技术结合，以双惩罚准似然（double penalized quasi-likelihood, DPQL）估计参数；采用Cauchy combination对多个时间聚合的检验进行联合；使用回顾式方差成分检验（retrospective variance component test），通过经验方差（empirical variance）替代完全基于模型假定的残差方差，以获得对随机效应分布误设的稳健性。
主要结论：通过与不同竞争方法（常效应检验、基于似然比的变系数检验等）的模拟比较，RVMMAT在case-control采样和无偏采样两种情况下均能良好地控制type I error，而基于似然比的同类方法在case-control采样下显著膨胀；且当真实效应为时变时，RVMMAT相比常效应方法具有更高的统计power。在MESA高血压的真实数据分析中，发现了新的（或已确认但此前未在纵向设定下检测出的）GWAS位点和通路。

关键设定与假设（在第二节基础上补全）¶

模型： logit(P(Y_ij = 1)) = X_ij^T β + g_i * γ(t_ij) + b_i + ε_ij 其中： - γ(t) 用自然三次smoothing spline（NCS）建模：γ(t) = Σ_{k=1}^{K} c_k * B_k(t) + b (t)*，其中 B_k(t) 为B-spline基函数，c_k 为系数，而 b(t)* 是惩罚项。实际上s(b)等价于一个特定的随机效应结构，用τ参数来惩罚曲率。 - b_i ~ N(0, σ_b²)：个体间随机截距。 - ε 的协方差矩阵（个体内）由核 K(t,t') 决定，等价于在GLMM中额外的随机效应（如随机斜率）。

关键假设（与已有文献对比）： A1：Smoothing spline modeling：假设时间效应 γ(t) 是光滑的（二阶导连续的）。这比Fan等直接用B-spline或固定阶多项式更柔性（少指定基点/阶数），但额外引入一个平滑参数τ需要估计。 A2：Mixed model assumptions：假设 b_i 和随机曲线效应是相互独立的，且各自服从均值为零的正态分布。本文的稳健检验正是针对这个假设的放松而设计——retrospective test不依赖于b_i的分布。 A3：Retrospective assumption：检验统计量的零分布是基于个体的表型 Y_i 在条件于个体协变量和群体结构时的随机性，而非基于个体间随机截距的实现。这是与标准似然比检验的核心分歧点：后者需要精确估计σ_b²。 A4：Cauchy combination assumption：使用 Cauchy 分布的和近似作为p值组合器，这种组合对p值间的相关性不敏感，但对于相关性非常极端时仍有轻微保守性。

主要结果¶

1. 稳健性（Type I error）（模拟重点）： - 设计：模拟数据在 无偏采样 和 case-control采样（病例比例50%）下产生。 - 竞争方法：RVMMAT vs. VC-test（标准变系数似然比检验，p-value从参方分布中获得） vs. Constant effect test（如 LMM类）。 - 结果： - 在无偏采样下，所有方法type I error相近，在α=0.05约束。 - 在case-control采样下（临床环境下实质），标准VC-test的type I error从0.05膨胀到约0.068-0.075（即失去控制）。RVMMAT的type I error稳定在0.048-0.053附近（稳健，且接近名义水平）。 - 结论：VaR严格优于替代方法在非随机采样下的稳健性。这是其核心selling point。

2. 统计效力（Power）（模拟重点）： - 条件：当真实效应是随时间线性增加或随时间先增后减（倒U形）时，RVMMAT的相对power如何？ - 结果： - 当效应为常数时，三者的power近似相等（RVMMAT略低，<2个百分点）。 - 当效应为时变（正值平均=0.2）时，RVMMAT的power平均比常效应检验高出10%-15%（因常效应检验的均值检验在效应有正负抵消时失效）。 - 即使在混合效应误设的情形下，RVMMAT依然能保持这一优势。 - 结论：当研究者怀疑存在时变效应时，RVMMAT大幅提升power，而不会在常效应时损失太多（≤2%）。

3. 计算性能： - 在一台标准服务器上模拟：对N=2000代、M=16时间点数据集，RVMMAT关联检验整个基因组（~800k SNP）约需3小时（CPU时间）。对比 CR-VCMM 可能需要数天。无直接模拟对比时间，但本文声称为“计算便宜”（computationally inexpensive）。

证明路线与技术技巧（理论型）¶

本文并非一篇纯理论论文，没有大定理和长脚手架式证明。它的技术贡献在于构造统计量并论证其渐近分布。因此，这里的“证明”是在构造源头上：

1. 整体路线： Step 1: DPQL Estimation (建立零假设下的参数) - 原假设 H0: γ(t) = 0 （SNP无任何时变效应）。 - 将模型简化为一个不包含g_i * γ(t)项的混合模型。 - 用双惩罚准似然（DPQL）估计固定效应 β 和方差成分 (σ_b², τ, λ)。这个过程循环迭代，直到收敛，得到的参数是进行score test的基础。

Step 2: 构造回顾式得分统计量 - 原假设下，γ(t)=0, 残差r = Y - Xβ̂ (标准化)。关键在于构建一个检验统计量 S，它与g_i及时间t的交互相关。 - 使用方差成分（variance component）类检验策略：检验统计量等价于一个二次型 Y^T * V_g * Y / 尺度。这里的矩阵V_g蕴含了“具有时变效应的fitted值对应g_i和K(.)的乘法”。 - 关键：作者构造了一个条件于协变量下的稳健协方差估计量，类似于retrospective框架：Var(S) = 2* tr[ (V_g * Sigma_obs)² ]，其中 Sigma_obs 是基于观测个体水平残差的经验协方差矩阵（而非出自族循环）。tr是迹。

Step 3: Cauchy Combination - 为联合检验不同时间尺度的效应，他们不是拟合一个复杂的二维表面，而是选取一个间隔上数个时间点（例如t1, t2,..., tK）计算该检验的p值，然后通过 Cauchy组合方法 T_com = Σ_T_n * tan( (0.5 - p_i) * π ) 合并成一个最终的联合p值。

2. 关键跳跃点： - 估计稳健方差（retrospective trick）：在混合模型中，Var(S) 通常包含一个来自于对随机效应分布求期望的项（2tr(Ω_V²)），这一项在模型忠实于随机效应时容易计算。但在误设下会出错。他们的关键跳跃是：不通过模型假定计算方差，而是通过 Σ_obs 直接估计观测残差的协方差。这要求解决一个关键计算问题：V_g 是M×M,而Σ_obs也是M×M，直接计算V_g * Σ_obs并求迹是O(M^3)，无法扩展。他们的跳跃（技术技巧1）就是利用Kroneycker结构和谱分解大幅降低了这个复杂度到 O(MlogM) 或 O(N)。

3. 技术技巧点名： * DPQL (Double Penalized Quasi-Likelihood)：这是对标准PQL的推广以嵌入空间平滑。通过对随机效应估计和光滑Penalized回归交替进行，解决了GLMM和非参数平滑的联合估计。 * Retrospective Score Test：核心是构建一个无需正确指定随机效应分布的检验统计量。其原理是将统计量在原假设下的方差估计转向经验个体间观测间的矩。这一技巧在GLM领域被用于稳健性检验。 * Cauchy Combination (Cauchy组合方法)：采用了(Liu & Xie, 2020)的Cauchy组合方法，它保证在(p_i)是独立的或弱相关下，组合统计量是τ分布的混合，但作者利用其尾部近似为Cauchy分布的特性，直接用一个tan变换将所有p值组合，并只需求解一个均值的计算。该方法计算极简单且对相关不敏感——这在需要联合多达K个垂直检验的统计组合中至关重要。

真实例子与应用¶

数据：MESA (Multi-Ethnic Study of Atherosclerosis) 研究中的高血压纵向测量。样本量为N=6678（包含白、非裔、西班牙、华裔美国人），有最多5次重复测量（跨度约10年）。表型为二值：收缩压/舒张压是否超过临界值为高血压。

怎么应用: 1. 准备数据：每个个体有一串重复的高血压状态（1/0），对应测量年龄作为时间轴。协变量包括年龄、性别、种族、前5个遗传PC。SNP是约800万个常染色体变异。 2. 拟合模型：对每个SNP，执行RVMMAT流程（Step 1-3），得到SNP-时间交互项的一个p值。 3. 后处理：对显著SNP（如p<5e-8）做通路富集分析（Pathway Analysis），使用的是DAVID工具。

结果： - 显著位点：发现若干经Bonferroni校正后仍然显著的位点（位点名称未在文章中单独列出，但一般指在曼哈顿图中显示），主要位于染色体4（如rs7698793附近）。 - 通路发现：通路分析识别出与G蛋白信号传导（G-protein signaling）和DNA损伤修复（DNA damage）相关的两个重要通路。作者声称，“这些通路可能提供血压变异随年龄变化的遗传基础”。 - 说明的意义：该结果想说明RVMMAT不仅是一个理论上稳健的模拟方法，而且在实际的、有噪声的、大规模的真实流行病学数据中能定位到生物学合理的位点和通路，体现出其“工具价值”。虽然发现的位点可能在以往的横断面研究中被提及，但本研究的纵向视角（检测时变效应）使结果指向了与“年龄相关疾病进展”相关的遗传标记。

🔎 结论是否比证明窄¶

部分结论的泛化受限于模拟设计。本文严格证明了该方法在case-control采样下的type I error能有效控制。但注意：他们使用的case-control采样是根据最后一次观测事期情况抽取的。如果采样是在整个病程期间（例如病例定义为在随访时间内任何时候都患病），稳健性是否需要额外验证？（文中没有提到这种情景）。
Cauchy组合的严格性。文章里提到了 “a joint test using the Cauchy combination method”。由于(Liu & Xie, 2020)已证明其p值的有效性（即使p值间弱相关），但本文并未在特定纵向相关模型下证明该组合无power损失。作者也有警示：在极端相关（ρ→1）下，Cauchy组合可能是保守的（即loss of power）。这个说明又是一个依赖于具体设置的近似值。
计算的可扩展性：虽然号称便宜（3h per 800k SNP），但这个时间限于特殊硬件（具体没说gpu加速与否）。如果是纯CPU的通用服务器，对不同协变量数目下的扩展还需要额外benchmark。声称的计算复杂度改进（O(M^3) → O(N log M)）对现实混合数据集的适用性也有边界（当个体时间点不规则时）。

四、开放问题（点到为止，扎根具体语句）¶

本文留下以下明确的开放问题（扎根于原文的limitations或未竟之处）：

推广到多种混杂效应：原文的方法假设变系数效应只出现在关注的g_i上。一个显而易见的推广是允许多个协变量（如环境暴露、生物标记）具有时变效应。这对应一个“功能型线性模型”的高维推广。这会引入更高维的参数空间，对DPQL估计和Retrospective test的均方误差必有影响。（扎根：引言末段“an important direction to extend our method to... multiple covariates with time-varying effects”）。
缺失数据下的稳健性：MESA数据有大致固定的观测次数，但标准纵向数据常含非随机截尾（informative dropout）。如果个体的流失（loss to follow-up）与未来的血压状况或基因型相关，RVMMAT中的回顾式检验在“条件”框架下能否依然保持其type I error？这是稳健性推断中一个更严重的问题。（扎根：讨论部分提一句“our results assume missing at random”）。
对罕见变异的power：本模拟和MESA分析中使用的是常见变异（MAF>0.05）。对于低频/罕见变异，污染的二值数据的信息量低，Smoothing spline方法是否需要采用不同的惩罚（如Laplacian平滑）或者对SKAT类方差成分检验做扩展？（扎根：结论“our method can be generalized to rare-variant testing through collapsing...” 但未做具体模拟或理论推导）。
计算开销对于百万SNP级的全量表型大数据集：虽然作者已声称计算高效（针对800k SNP~3h），但如果需要同时进行全基因组-时变互作分析（gene-by-time interaction search for all pairs of time points），其计算代价可能仍难以接受。这是纵向GWAS普遍的挑战。Cauchy组合虽然解决了联合检验，但并未降低扫描次数。（扎根：讨论“computational scalability to even larger biobank-scale data”）。

Maintained by 陈星宇 · Homepage · Source on GitHub