跳转至

Uncertainty intervals for multilevel models with missing not at random data

作者: Minna Genbäck
主题: 因果推断
相关性: 7/10
链接: https://arxiv.org/abs/2606.24228

一、领域脉络与小综述

这个方向是什么

本文处理的根本问题是：在纵向数据（longitudinal data）中，当数据缺失机制不可忽略（Missing Not At Random, MNAR）时，如何对线性混合效应模型（multilevel / mixed-effects model）中的固定效应参数进行推断。标准的多层模型软件（如 lme4）默认假设数据为随机缺失（Missing At Random, MAR），即缺失概率仅依赖于已观测到的变量。然而，在老龄化、健康等纵向研究中，这一假设常被违反（例如，健康状况更差的个体更易退出研究），导致估计偏倚。本文的核心贡献是提出一种敏感性分析方法，通过引入一个量化偏离MAR程度的敏感性参数，将点识别放松为部分识别（partial identification），从而在比MAR更弱的假设下给出参数的估计区间（uncertainty intervals）。

发展脉络（history）

  • 奠基工作：缺失数据理论的基石是 Little & Rubin (2002) 的专著，系统定义了MAR、MNAR等机制。Manski (2003) 提出了部分识别（partial identification）的框架，为放弃点识别、转而寻求识别集（identification set）提供了理论基础。
  • 主要进展（MNAR建模）：处理MNAR的主流框架分为两类：
    • 选择模型（Selection Model）：将缺失/脱落作为因变量建模。Carpenter & Kenward (2013) 的专著是这一方向的经典参考。Enders (2011, 2023) 综述了包括全信息最大似然（FIML）、多重插补（MI）在内的现代缺失数据处理方法，并讨论了MNAR模型。
    • 模式混合模型（Pattern Mixture Model）：将缺失模式作为预测变量。Josefsson et al. (2016) 将倾向得分匹配与模式混合模型结合，处理纵向因果推断中的非随机脱落。Hammon & Zinn (2020) 提出了针对二元多层MNAR数据的多重插补方法，使用删失双变量probit模型。
    • 这些方法都依赖于额外的、不可检验的分布假设（如联合正态性、特定的参数结构）来实现点识别，一旦假设错误，结论可能具有误导性。
  • 当前 frontier（部分识别与敏感性分析）：为规避强假设，部分识别路线被提出。Vansteelandt et al. (2006) 系统阐述了“不确定性区间”（uncertainty intervals）的推断框架，通过指定敏感性参数的范围来得到覆盖真实参数的区间。Genbäck et al. (2015) 和 Genbäck & de Luna (2019) 将这一思路应用于非多层回归模型和因果推断中的未观测混杂问题，使用相关性作为敏感性参数。
  • 本文的位置：本文是上述部分识别路线在多层模型这一重要但尚未被充分覆盖的设定下的直接推广。作者明确指出，此前类似的分析（Genton, 2004; González-Farías et al., 2004）仅推导了“每个簇要么全部观测、要么全部缺失”这一不现实的特殊情形下的偏倚表达式，而本文给出了一般性的偏倚解析表达式（Theorem 1），并首次在模拟和真实数据中进行了系统评估。

子线索聚类

  1. 缺失数据理论与方法：Little & Rubin (2002), Molenberghs et al. (2015), Enders (2023)。提供基础框架和综述。
  2. 多层模型与软件：Bates et al. (2015) (lme4), Demidenko (2004)。提供模型拟合工具。
  3. MNAR建模（点识别）：Carpenter & Kenward (2013), Enders (2011), Josefsson et al. (2016), Hammon & Zinn (2020)。依赖强分布假设实现点识别。
  4. 部分识别与敏感性分析：Manski (2003), Vansteelandt et al. (2006), Genbäck et al. (2015), Genbäck & de Luna (2019), Imai et al. (2010)。放弃点识别，通过敏感性参数获得识别集。本文属于此线索。

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

  • 核心问题1：如何在放弃MAR假设后，仍能对参数进行有意义的推断？
  • 核心问题2：如何量化偏离MAR的程度，并将其与估计偏倚联系起来？
  • 核心问题3：如何将部分识别框架从简单的回归模型推广到更复杂的结构（如多层模型、纵向数据）？
  • 已知瓶颈：在多层模型中，偏倚的解析表达式复杂，计算上需要处理高维截断正态分布的期望，计算负担大。

⚠️ 作者的 framing

  • 作者的缺口frame：作者将缺口frame为“多层模型下的MNAR敏感性分析缺乏一般性的偏倚解析表达式和系统评估”。他们声称此前的工作（Genton, 2004; González-Farías et al., 2004）只处理了不现实的特殊情形，而Grilli & Rampichini (2010) 只讨论了特例的偏倚，没有给出解析表达式。因此，本文的Theorem 1被呈现为“显然的下一步”。
  • 被淡化/回避的竞争路线：作者明确将本文定位为与“点识别”路线（如Josefsson et al., 2016; Hammon & Zinn, 2020）的替代方案。他们淡化了这些方法在特定假设下的有效性，强调其“不可检验”的风险。作者也回避了与多重插补（MI） 框架下处理MNAR的复杂方法的直接比较（如Wang et al., 2025 被提及但未深入对比）。
  • 什么明显该被引/该存在、却没出现在intro里？：本文的敏感性参数是相关性，这与Genbäck et al. (2015) 和 Genbäck & de Luna (2019) 一脉相承。然而，在因果推断的敏感性分析中，另一种常见参数化是E值（E-value）（VanderWeele & Ding, 2017）或混杂强度。作者没有讨论为何选择相关性而非其他参数化方式，也未与E值方法进行任何对比。这是一个值得研究者去查的问题：在多层模型设定下，E值方法是否适用？与相关性参数化相比有何优劣？

张力

未见明显对立引用。所有被引工作基本认同“MAR假设在健康老龄化研究中可能不成立”这一前提，分歧在于如何应对：是采用更强的分布假设实现点识别，还是放弃点识别进行敏感性分析。本文属于后者，并与前者形成互补而非对立。

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

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

  • 符号：

    • i = 1, ..., N: 簇（cluster）的索引，例如家庭（household）。
    • k = 1, ..., n_i: 簇 i 内观测的索引，例如个体 j 在时间点 k 的观测。
    • y_i: n_i × 1 向量，簇 i 的可观测结果变量（例如记忆得分）。
    • X_i: n_i × m 矩阵，簇 i 的可观测协变量（包括时间、年龄、教育等）。
    • β: m × 1 向量，感兴趣的参数（固定效应系数）。
    • η_i: n_i × 1 向量，结果模型的随机成分（包含随机效应和误差）。
    • z*_i: n_i × 1 潜在变量向量，决定缺失机制。z_ik = I(z*_ik > 0) = 1 表示 y_ik 被观测到，否则缺失。
    • δ: m × 1 向量，缺失模型中的固定效应系数。
    • u_i: n_i × 1 向量，缺失模型的随机成分。
    • Σ_i: n_i × n_i 矩阵，η_i 的协方差矩阵。
    • Ω_i: n_i × n_i 矩阵，u_i 的协方差矩阵。
    • ∆_i: n_i × n_i 矩阵，敏感性参数矩阵，即 cov(η_i, u_i)。当 ∆_i = 0 时，缺失机制为MAR。
    • y_i^obs: y_i 中实际被观测到的子向量。
    • X_i^obs: X_i 中对应于 y_i^obs 的行组成的子矩阵。
    • Σ_i^obs: Σ_i 中对应于 y_i^obs 的行和列组成的子矩阵。
    • v_i: 一个随机向量，定义为 u_i 在给定观测模式下的条件分布（见Theorem 1）。

  • 模型：

    • 结果模型（线性混合效应模型）：y_i = X_i β + η_i，其中 η_i ~ MVN(0, Σ_i)。
    • 缺失模型（广义线性混合效应模型，probit链接）：z*_i = X_i δ + u_i，其中 u_i ~ MVN(0, Ω_i)。观测指示变量 z_ik = I(z*_ik > 0)。
    • 联合模型：(η_i, u_i) 服从联合多元正态分布，协方差矩阵为 [Σ_i, ∆_i; ∆_i^T, Ω_i]。∆_i 是未知的，是敏感性分析的对象。

  • 可观测数据：

    • 可观测：对于每个簇 i，我们能观测到 X_i（所有协变量，无论结果是否缺失），以及 y_i^obs（结果变量的观测子集）和对应的缺失指示变量 z_i。
    • 不可观测/潜在：y_i 中缺失的部分（y_i^mis），以及潜在变量 z*_i。η_i 和 u_i 也是不可观测的随机效应。核心困难：我们无法从数据中检验 ∆_i = 0（即MAR假设）是否成立。

第二步：讲最小内核

本文的核心思路可以用一个最简特例来理解：单层线性回归，且缺失只发生在结果变量上，协变量完全观测。

  • 最简设定：

    • 没有簇结构（N 个独立个体）。
    • 结果模型：y_i = x_i^T β + η_i，η_i ~ N(0, σ^2)。
    • 缺失模型：z*_i = x_i^T δ + u_i，u_i ~ N(0, 1)（probit模型，方差固定为1以识别）。
    • 联合分布：(η_i, u_i) ~ MVN(0, [σ^2, ρσ; ρσ, 1])。这里，敏感性参数 ρ 是一个标量，即 cor(η_i, u_i)。当 ρ=0 时，缺失机制为MAR。

  • 核心问题：我们想估计 β，但 y_i 只在 z*_i > 0 时被观测到。如果我们忽略缺失机制，直接用观测到的 y_i^obs 做OLS回归，会得到有偏的估计 \tilde{β}。

  • 核心思路（Theorem 1 的特例）：

    1. 写出偏倚表达式：在给定 x_i 和 z_i=1（即被观测到）的条件下，y_i 的条件期望为： E[y_i | x_i, z_i=1] = x_i^T β + E[η_i | u_i > -x_i^T δ]。 由于 (η_i, u_i) 联合正态，有 E[η_i | u_i > -x_i^T δ] = ρσ * E[u_i | u_i > -x_i^T δ]。 这里 E[u_i | u_i > -x_i^T δ] 是一个截断正态分布的期望，可以解析计算（即逆米尔斯比率，Inverse Mills Ratio）。
    2. 偏倚形式：因此，用观测数据做OLS得到的 \tilde{β} 的偏倚为： bias(\tilde{β}) = (∑_{i: z_i=1} x_i x_i^T)^{-1} ∑_{i: z_i=1} x_i * [ρσ * E(u_i | u_i > -x_i^T δ)]。 这个表达式与Theorem 1完全对应：(X^obs)^T (Σ^obs)^{-1} X^obs 退化为 ∑ x_i x_i^T，∆_i Ω_i^{-1} E(v_i) 退化为 ρσ * E(u_i | ...)。
    3. 如何做敏感性分析：
      • 首先，在MAR假设（ρ=0）下，用 lme4 或 glm 拟合结果模型和缺失模型，得到 \tilde{β}、\hat{δ}、\hat{σ} 等估计。
      • 然后，指定一个敏感性参数 ρ 的值（例如 ρ=0.3）。这个值代表“未观测到的结果与缺失倾向的相关性”，是先验知识或专家判断。
      • 接着，将 \hat{δ} 和 ρ 代入偏倚表达式，计算 \hat{bias}。
      • 最后，得到偏倚校正后的估计：\hat{β} = \tilde{β} - \hat{bias}。
      • 通过让 ρ 在一个合理区间内（如 [0, 0.5]）变化，可以得到一系列 \hat{β} 和对应的置信区间，它们的并集就是不确定性区间。

  • 为什么这个特例抓住了本质：这个特例清晰地展示了本文的核心数学操作：将MNAR导致的偏倚表达为敏感性参数（ρ）与一个可估计的截断正态期望的乘积。多层模型的一般化只是将这个标量 ρ 扩展为矩阵 ∆_i，并将截断正态期望从一元扩展到多元（E(v_i)），但核心的“偏倚 = 敏感性参数 × 可估计项”这一结构保持不变。

三、这篇论文做了什么

三句话

  1. 研究了什么问题：在线性多层（混合效应）模型中，当数据缺失机制为MNAR时，如何对固定效应参数 β 进行敏感性分析，以评估结论对MAR假设偏离的稳健性。
  2. 核心工具/方法：通过联合建模结果和缺失过程（均为多层模型），引入协方差矩阵 ∆_i 作为敏感性参数，推导出广义最小二乘（GLS）估计量 \tilde{β}_{GLS} 的偏倚解析表达式（Theorem 1）。通过指定 ∆_i 的值（或由少数标量敏感性参数构造），可估计并校正该偏倚，进而构建不确定性区间。
  3. 主要结论：模拟研究表明，当敏感性参数被正确指定时，偏倚校正方法在有限样本下表现良好，能显著降低估计偏倚并恢复接近名义水平的覆盖率。应用于SHARE数据的实例分析表明，即使考虑MNAR，关于孤独感和体力活动对记忆轨迹影响的主要结论与MAR分析基本一致。

关键设定与假设

  • 设定：N 个簇，每个簇 i 有 n_i 个观测。结果模型为线性混合效应模型（公式1），缺失模型为带probit链接的广义线性混合效应模型（公式2）。两者共享相同的协变量矩阵 X_i。
  • 假设：
    1. 联合正态性：(η_i, u_i) 服从联合多元正态分布（公式3）。这是推导偏倚解析表达式的核心分布假设。
    2. 模型正确指定：结果模型和缺失模型的固定效应部分（X_i β, X_i δ）和随机效应结构（Σ_i, Ω_i）被正确指定。这是所有参数估计（\tilde{β}, \hat{δ}, \hat{Σ}_i, \hat{Ω}_i）的基础。
    3. 敏感性参数 ∆_i 已知或可指定：这是敏感性分析的核心。作者建议将复杂的 ∆_i 参数化为少数几个有意义的标量相关性（如 ρ_hh, ρ_ind, ρ_err），并假设其位于一个先验合理的区间内。
    4. 与已有文献的对比：相比Genton (2004) 和 González-Farías et al. (2004) 的“全有或全无”缺失假设，本文的假设更一般，允许任意缺失模式。相比Grilli & Rampichini (2010) 的数值讨论，本文给出了解析表达式。相比Josefsson et al. (2016) 和 Hammon & Zinn (2020) 的点识别方法，本文放弃了额外的分布约束，转而进行部分识别。

主要结果

  • Theorem 1（核心理论结果）：给出了GLS估计量 \tilde{β}_{GLS} 在MNAR下的偏倚解析表达式。该表达式将偏倚分解为两部分：一个由观测数据决定的权重矩阵 (∑ (X_i^obs)^T (Σ_i^obs)^{-1} X_i^obs)^{-1} 和一个包含敏感性参数 ∆_i 的项 ∑ (X_i^obs)^T (Σ_i^obs)^{-1} {∆_i Ω_i^{-1} E(v_i)}^obs。关键：E(v_i) 是一个双重截断多元正态分布的期望，其截断点由 X_i δ 决定，这构成了计算的主要难点。
  • 模拟研究（核心量化结论）：
    • 偏倚校正效果：在4种设计下（嵌套/非嵌套，不同相关性强度），未校正的MAR估计（\tilde{β}）在MNAR下存在显著偏倚（例如，Design 2中 β_{time} 的偏倚为7.1%）。应用偏倚校正后（假设已知真实 ρ 值），偏倚大幅降低（β_{time} 偏倚降至0.8%），MSE也显著减小（见表5）。
    • 覆盖率：未校正估计的95%置信区间覆盖率极低（例如，Design 2中 β_{time} 覆盖率仅0.6%）。校正后，覆盖率恢复到接近名义水平（87%-95%），尽管对于 time 变量在嵌套设计下仍偏低（约40%），但作者指出这是由于标准误估计不准所致。
    • 模型选择：模拟发现，非嵌套模型（unnested） 在数值稳定性和偏倚校正效果上均优于嵌套模型（nested）。嵌套模型在拟合时频繁出现收敛错误（见表B1），且对 time 变量的标准误估计严重偏低（仅70-80%），导致校正后覆盖率仍不理想。非嵌套模型则表现稳健。
  • 真实数据例子（SHARE）：
    • 数据：欧洲健康、老龄化和退休调查（SHARE）第4-9波数据，27,313名65岁以上个体。
    • 方法应用：拟合多层模型，分析基线孤独感和体力活动与记忆轨迹的关联。敏感性参数设定为 ρ_hh, ρ_ind, ρ_err ∈ [0, 0.5]，假设记忆衰退与脱落风险正相关。
    • 结果：MAR分析显示“经常感到被冷落”与记忆得分平均下降0.59分相关，“每周至少一次剧烈运动”与平均上升0.41分相关。MNAR敏感性分析给出的不确定性区间（见表2）虽然略宽，但方向一致，结论稳健。例如，“经常感到被冷落”的估计范围是 (-0.77, -0.46)，仍显著为负。
    • 例子想说明什么：验证了所提方法在实际纵向研究中的可用性，并展示了如何通过不确定性区间来评估MAR结论对MNAR假设的敏感程度。

证明路线与技术技巧

  • 整体路线（Theorem 1 证明）：
    1. 写出GLS估计量的期望：E[\tilde{β}_{GLS}] = E[ (∑ (X_i^obs)^T (Σ_i^obs)^{-1} X_i^obs)^{-1} ∑ (X_i^obs)^T (Σ_i^obs)^{-1} y_i^obs ]。
    2. 条件期望：对 X 取条件，将期望移到内部：E[\tilde{β}_{GLS} | X] = (∑ (X_i^obs)^T (Σ_i^obs)^{-1} X_i^obs)^{-1} ∑ (X_i^obs)^T (Σ_i^obs)^{-1} E[y_i^obs | X]。
    3. 计算 E[y_i^obs | X]：这是关键跳跃点。由于 y_i 的观测与否取决于 z_i，我们需要计算 E[y_i | X, z_i]。利用Lemma 1，将 η_i 分解为与 u_i 独立的部分和与 u_i 相关的部分：η_i = w_i + ∆_i Ω_i^{-1} u_i。由于 w_i 与 u_i 独立且均值为0，E[η_i | X, z_i] = ∆_i Ω_i^{-1} E[u_i | X, z_i]。
    4. 代入并完成证明：将 E[y_i^obs | X] = X_i^obs β + {∆_i Ω_i^{-1} E[u_i | X, z_i]}^obs 代入步骤2，得到 E[\tilde{β}_{GLS} | X] = β + bias，其中 bias 即为Theorem 1中的表达式。
  • 关键跳跃点：Lemma 1 是整个证明的核心。它巧妙地利用了多元正态分布的性质，通过构造一个与 u_i 独立的 w_i，将复杂的条件期望 E[η_i | u_i 的截断] 转化为一个可处理的项 ∆_i Ω_i^{-1} E[u_i | 截断]。这个转化将偏倚与敏感性参数 ∆_i 直接联系起来。
  • 技术技巧点名：
    • 多元正态分布的条件分布与独立性：利用 w_i = η_i - ∆_i Ω_i^{-1} u_i 构造独立成分，这是处理相关正态随机变量的标准技巧。
    • 双重截断多元正态分布：E(v_i) 的计算需要数值积分，作者使用了 MomTrunc R包。这是计算瓶颈所在。
    • GLS估计：使用GLS而非OLS，以处理簇内相关性（Σ_i 非对角）。

🔎 结论是否比证明窄

  • 窄结论1：Theorem 1 的证明严格依赖于联合正态性假设（公式3）。然而，作者在结论和讨论中并未强调这一假设的敏感性。如果真实数据生成过程偏离正态性，偏倚表达式可能不再准确。作者在模拟中使用了正态生成的数据，未测试对非正态的稳健性。
  • 窄结论2：模拟研究中的偏倚校正假设敏感性参数 ρ 是已知的（即使用了数据生成时的真实值）。在实际应用中，ρ 是未知的，只能指定一个范围。模拟并未评估当 ρ 被错误指定时（例如，真实 ρ=0.3 但假设范围是 [0, 0.1]），不确定性区间的覆盖率会如何变化。作者在论文中声称“如果敏感性参数的真值位于预先指定的范围内，则不确定性区间将以至少 (1-α)% 的概率覆盖 β”，但模拟并未直接验证这一声称在 ρ 被错误指定时的表现。
  • 窄结论3：模拟中使用的标准误是 lmer 返回的默认标准误，未考虑偏倚校正带来的额外不确定性。作者在4.1节承认了这一点，并认为“任何修正都不足以证明额外的计算成本是合理的”。这意味着最终的不确定性区间可能偏窄，其实际覆盖率可能低于名义水平。模拟结果（表5）中，校正后估计的覆盖率在部分设计下（如嵌套设计的 time 变量）仍远低于95%，部分原因就在于此。

四、开放问题

  1. 放松联合正态性假设：Theorem 1 的推导依赖于 (η_i, u_i) 的联合正态性。能否在更弱的分布假设（如椭圆分布、或仅假设前两阶矩）下推导出类似的偏倚表达式或界？这扎根于本文的公式(3) 和Lemma 1 对正态性的依赖。
  2. 标准误的校正：本文使用的标准误是未校正的 lmer 输出，忽略了偏倚估计 \hat{bias} 本身的抽样变异性。如何推导出考虑了这一变异性的正确标准误，并构建渐近有效的置信区间？这扎根于本文4.1节的讨论：“we choose to ignore this potential issue”。
  3. 计算可扩展性：当簇内观测数 n_i 很大时，计算双重截断多元正态分布的期望 E(v_i) 在计算上变得不可行。是否存在更高效的近似方法（如基于MCMC的采样、或利用稀疏结构）来克服这一瓶颈？这扎根于本文第5节的讨论：“the computational burden becomes substantial when clusters contain many observations”。
  4. 敏感性参数的选择与校准：本文使用相关性作为敏感性参数。是否存在其他更直观或更易从先验知识中获得的参数化方式（如E值、或直接对偏倚大小进行参数化）？如何系统性地为这些参数指定合理的范围？这扎根于本文2.2节的讨论：“it may be difficult to specify and interpret in practice”。

Maintained by 陈星宇 · Homepage · Source on GitHub

评论