On Ignorability of Preferential Sampling in Geostatistics

讲者: Ganggang Xu
会场: Time Series and Longitudinal Data and Bayesian Methods
报告题目: On Ignorability of Preferential Sampling in Geostatistics
链接: arXiv
来源: JCSDS 2026 · 返回会议总览

---

一、领域脉络与小综述

• 这个方向是什么：地理统计学中的“优先采样”（preferential sampling）问题。根本问题是：当观测位置（采样点）的分布与所观测的空间过程（如温度、树木直径）存在依赖时，如何无偏地估计空间回归系数、空间协方差函数和预测？标准地统计模型假设采样位置是确定性的或独立于空间过程，但这一假设在许多实际场景（如动物携带传感器、生态调查中倾向于去物种丰富区域）中不成立。该子方向当前成熟度较高，已有大量似然方法，但计算昂贵且依赖对采样机制的参数化假设。

• 发展脉络（history）：

  • 奠基工作：Diggle et al. (2010) 首次将优先采样形式化，提出用对数高斯Cox过程（LGCP）建模采样位置，用高斯过程建模标记（mark），通过两个高斯随机场之间的参数化关系（如 X(s) = γY(s)）捕捉依赖。这为后续所有工作提供了标准框架。
  • 主要进展（似然方法）：
    • Pati et al. (2011) 发展了贝叶斯方法估计标记过程均值。
    • Zidek et al. (2014) 将框架扩展到时空背景。
    • Ferreira & Gamerman (2015) 将其用于指导新采样位置的选择（最优设计）。
    • Dinsdale & Salibian-Barrera (2019a) 指出早期蒙特卡洛似然近似可能不准确，并展示了TMB（Template Model Builder）等数值方法在中等复杂度下优于蒙特卡洛方法。
    • Dinsdale & Salibian-Barrera (2019b) 将方法扩展到移动传感器（如海豹携带的探头）产生的数据。
    • Pennino et al. (2019) 在物种分布模型中应用LGCP框架。
    • Amaral et al. (2024) 考虑了空间变化的采样强度。
  • 当前frontier（非似然/加权方法）：
    • Schliep et al. (2023) 和 Hsiao & Waller (2025) 采用复合似然方法，引入与强度相关的权重来改进标记过程参数估计。Hsiao & Waller (2025) 特别提出了逆采样强度加权（ISIW）方法，并发现其在模型误设下预测性能优于标准优先采样方法，且指出“准确的参数估计与预测性能相关性很小”。
  • 本文的位置：本文挑战了该领域的一个隐含共识——即必须显式建模采样机制才能纠正偏差。作者证明，在Diggle et al. (2010) 的框架及其扩展下，一些忽略优先采样的非似然方法（最小二乘、核平滑）仍然可以产生无偏且一致的估计量。这为“何时可以忽略优先采样”提供了理论条件，并提出了一个无需参数化采样机制的计算高效方法。

• 子线索聚类：

  1. 全似然方法（Full Likelihood）：以Diggle et al. (2010) 为起点，包括Pati et al. (2011), Zidek et al. (2014), Dinsdale & Salibian-Barrera (2019a,b), Pennino et al. (2019)。核心是联合建模点过程和标记过程，通过MCMC或TMB进行推断。优点是理论上最有效，缺点是计算昂贵、对模型误设敏感。
  2. 加权/复合似然方法（Weighted/Composite Likelihood）：以Schliep et al. (2023) 和 Hsiao & Waller (2025) 为代表。通过引入基于采样强度的权重来修正估计方程，避免联合建模。计算上更可行，但理论保证（如一致性）尚未完全建立（本文指出“no theoretical guarantees have yet been established”）。
  3. 标记点过程框架下的扩展：针对更复杂的数据结构，如重复标记点过程（Xu et al. 2024, Yin et al. 2021）、函数型数据（Gervini & Baur 2020）。这些工作关注如何利用复制或函数型结构来建模共享潜场。

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

  1. 识别：在优先采样下，空间回归系数β和协方差参数θ是否可识别？需要什么假设？
  2. 估计：如何构造计算可行且对模型误设稳健的估计量？似然方法 vs. 矩方法 vs. 加权方法的优劣。
  3. 推断：如何为估计量构造有效的置信区间？渐近理论是否成立？
  4. 预测：在优先采样下，克里金（kriging）预测是否仍然最优？Hsiao & Waller (2025) 的结果暗示参数估计精度与预测性能可能脱钩。

• ⚠️ 作者的 framing：

  • 作者把缺口 frame 成：现有似然方法需要指定采样机制的参数形式，易受模型误设影响，且计算昂贵。他们发现一个“令人惊讶的发现”（surprising finding）：在Diggle et al. (2010) 的框架下，一些忽略优先采样的方法仍然一致。因此，他们的工作成为“显然的下一步”：系统地研究可忽略的条件，并开发一个无需参数化采样机制、计算高效且稳健的方法。
  • 被淡化或回避的竞争路线：作者淡化了Hsiao & Waller (2025) 的ISIW方法。虽然引用了它，但强调其“尚未建立理论保证”，而本文提供了理论。作者也回避了与全似然方法在效率上的直接比较，承认自己的方法方差更大，但用“计算速度快数千倍”和“对模型误设稳健”来平衡。
  • 什么明显该被引/该存在、却没出现在intro里？：本文的intro没有引用任何关于“因果推断中处理选择偏差”的文献，例如工具变量、倾向性评分加权、双重稳健估计等。虽然领域不同，但“当选择机制依赖于结果时如何做推断”是因果推断的核心问题，其方法论（如逆概率加权）与本文的“逆采样强度加权”有深刻联系。这是一个值得研究者去查的潜在连接点。

• 张力：未见明显对立引用。所有被引工作都认同优先采样是一个需要处理的问题，分歧在于如何处理（似然 vs. 加权 vs. 矩方法）。本文的“可忽略性”发现是第一个挑战“必须处理”这一共识的工作，但作者并未声称自己的方法在所有情况下都优于似然方法，而是指出了在特定条件下可以忽略。

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

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

• 符号：

  • s ∈ S ⊂ R²: 空间位置。
  • Z(s): 在位置 s 处的观测值（标记，mark），是一个标量随机变量。
  • w(s) ∈ R^p: 在位置 s 处的 p 维空间协变量。
  • β ∈ R^p: 回归系数，是要估的参数。
  • Y(s): 一个零均值、平稳、各向同性的高斯随机场，代表空间相关性。方差为 σ²_Y，协方差函数为 C_Y(||s-t||)。
  • e(s): 独立同分布的高斯噪声（块金效应，nugget effect），方差为 σ²_e。
  • N: 一个定义在 S 上的对数高斯Cox过程（LGCP），其实现就是观测到的采样位置集合 {s_1, ..., s_n}。
  • λ(s): LGCP的潜强度函数，λ(s) = λ₀(s) exp[X(s)]。λ₀(s) 是基线强度，X(s) 是另一个零均值高斯随机场。
  • X(s): 控制采样强度的潜高斯随机场。其协方差函数为 C_X(||s-t||)。
  • C_XY(||s-t||): X(s) 和 Y(t) 之间的交叉协方差函数。这是捕捉“优先采样”的关键。如果 C_XY(r) ≠ 0，则采样位置和观测值相关。
  • ρ(s) = E[λ(s)]: 点过程 N 的一阶强度函数。
  • ρ₂(s,t) = E[λ(s)λ(t)]: 点过程 N 的二阶强度函数。
  • |S_n|: 观测区域 S_n 的面积。
  • |N|: 观测到的采样点个数。

• 模型：

  1. 数据生成模型（标记过程）：Z(s) = w(s)^T β + Y(s) + e(s)。这是一个经典的地统计模型，其中 Y(s) 和 e(s) 独立。
  2. 采样位置生成模型（点过程）：采样位置 {s_i} 是LGCP N 的一个实现。这意味着给定潜高斯场 X(s)，采样位置是一个强度为 λ(s) = λ₀(s) exp[X(s)] 的非齐次泊松过程。
  3. 依赖结构：X(s) 和 Y(s) 是联合高斯的，其依赖由交叉协方差函数 C_XY(||s-t||) 完全刻画。Diggle et al. (2010) 假设了特殊形式 X(s) = γY(s)，本文将其推广到任意各向同性交叉协方差。

• 可观测数据：

  • 可观测：位置 s_i 和该位置上的观测值 Z(s_i)，以及位置上的协变量 w(s_i)。即 { (s_i, Z(s_i), w(s_i)) }。
  • 不可观测/潜在：潜高斯场 Y(s) 和 X(s) 的完整实现，以及噪声 e(s_i)。我们只能通过观测到的 Z(s_i) 和 s_i 的分布来推断 Y(s) 和 X(s) 的协方差结构。

第二步：讲最小内核

本文的核心发现可以浓缩为一个最简特例：当采样位置由一个LGCP生成，且标记过程是高斯过程时，即使存在优先采样（即 X(s) 和 Y(s) 相关），最小二乘估计量 ˆβ 对非截距项**的回归系数仍然是无偏的。

最简特例：假设模型为 Z(s) = β₀ + β₁ w(s) + Y(s) + e(s)，其中 w(s) 是标量协变量。采样位置来自LGCP，其潜强度为 λ(s) = λ₀ exp[X(s)]，且 X(s) 与 Y(s) 相关（C_XY(0) ≠ 0）。

核心思路：考虑最小二乘估计量 ˆβ = (∑ w(s)w(s)^T)^{-1} ∑ w(s)Z(s)。其期望可以通过Campbell定理和Stein引理计算。

1. Campbell定理：对于点过程 N 上的函数 f(s)，有 E[∑_{s∈N} f(s)] = ∫_S E[λ(s) f(s)] ds。这里 λ(s) 是随机强度。
2. Stein引理：对于联合高斯的 (X, Y)，有 E[Y exp(X)] = Cov(X, Y) E[exp(X)]。

现在，计算 E[∑ w(s)Z(s)]： E[∑ w(s)Z(s)] = ∫ E[λ(s) w(s) Z(s)] ds = ∫ w(s) E[ exp(X(s)) * (w(s)^T β + Y(s) + e(s)) ] λ₀(s) ds = ∫ w(s) w(s)^T β E[exp(X(s))] λ₀(s) ds + ∫ w(s) E[Y(s) exp(X(s))] λ₀(s) ds 应用Stein引理到第二项：E[Y(s) exp(X(s))] = Cov(X(s), Y(s)) E[exp(X(s))] = C_XY(0) E[exp(X(s))]。 所以： E[∑ w(s)Z(s)] = [∫ ρ(s) w(s) w(s)^T ds] β + C_XY(0) ∫ ρ(s) w(s) ds。

类似地，E[∑ w(s)w(s)^T] = ∫ ρ(s) w(s) w(s)^T ds。

因此，E[ˆβ] = β + C_XY(0) * [∫ ρ(s) w(s) w(s)^T ds]^{-1} ∫ ρ(s) w(s) ds。

关键观察：如果 w(s) 的第一个元素是1（截距项），那么 ∫ ρ(s) w(s) ds 的第一个元素是 ∫ ρ(s) ds，而其他元素是 ∫ ρ(s) w_j(s) ds。矩阵 [∫ ρ(s) w(s) w(s)^T ds]^{-1} 乘以这个向量，结果是一个向量，其第一个元素非零，但其他元素为零。这是因为对于非截距项，w(s) 是中心化的或与截距正交的，但这里更直接的论证是：偏差项 C_XY(0) * [∫ ρ(s) w(s) w(s)^T ds]^{-1} ∫ ρ(s) w(s) ds 是一个常数向量，它只影响截距项 β₀ 的估计，而不影响斜率 β₁。

结论：在这个最简特例下，ˆβ₁ 是 β₁ 的无偏估计量，而 ˆβ₀ 有偏差 C_XY(0) * α，其中 α 是一个可计算的常数。这个偏差可以通过估计 C_XY(0) 来校正。整个论文的一般化工作就是将这个核心观察推广到更一般的协方差结构、更复杂的参数估计（如协方差参数）和渐近理论。

三、这篇论文做了什么

• 三句话：

  1. 研究了什么问题：在地统计学的优先采样框架下，研究了何时可以“忽略”优先采样，即使用忽略采样机制的经典方法（最小二乘、核平滑）仍能得到一致估计。
  2. 核心工具/方法：利用Campbell定理和Stein引理，在Diggle et al. (2010) 的LGCP框架下，推导了最小二乘估计量 ˆβ 和核平滑估计量 ˆV_Y(r), ˆC_XY(r) 的期望和渐近性质。提出了基于最小对比度（MC）和复合似然（CL）的协方差参数估计方法。
  3. 主要结论：在一般各向同性交叉协方差假设下，最小二乘估计对非截距回归系数无偏；核平滑估计对半变异函数和交叉协方差函数一致。基于此，提出了无需参数化采样机制的估计和推断方法，并在模拟和真实数据中展示了其计算效率和稳健性。

• 关键设定与假设：

  • 模型设定：Z(s) = w(s)^T β + Y(s) + e(s)，其中 Y(s) 是平稳各向同性高斯过程，e(s) 是独立高斯噪声。
  • 采样机制：采样位置来自LGCP，其潜强度 λ(s) = λ₀(s) exp[X(s)]，X(s) 是平稳各向同性高斯过程。
  • 依赖结构：X(s) 和 Y(s) 是联合高斯的，其依赖由各向同性交叉协方差函数 C_XY(||s-t||) 刻画。这是对Diggle et al. (2010) 中 X(s) = γY(s) 这一线性假设的显著放宽。
  • 渐近框架：采用“递增域”（increasing-domain）渐近，即观测区域 S_n 随 n 扩大至 R²。条件 (C1) 保证区域在各方向均匀增长。
  • 技术条件：
    • (C2): 强混合条件，保证点过程的空间依赖随距离衰减足够快（多项式速率），LGCP满足此条件。
    • (C3)-(C6): 协变量、基线强度、协方差函数和交叉协方差函数的有界性和可积性条件。这些是保证矩存在和积分收敛的标准正则性条件。
    • (C7): 设计矩阵的期望可逆，是回归分析的标准条件。
    • (C8)-(C10): 对参数化半变异函数和权重函数的正则性条件，保证参数估计的唯一性和可识别性。

• 主要结果：

  • Theorem 1（回归系数）：在条件(C1)-(C7)下，最小二乘估计量 ˆβ_n 渐近正态，收敛到 β*₀ = β₀ + C_XY(0)α，其中 α 只影响截距项。这意味着非截距回归系数被一致估计，截距偏差可通过估计 C_XY(0) 校正。该定理还给出了渐近协方差矩阵 Σ_n 的显式表达式，为推断提供了基础。
  • Theorem 2（协方差函数）：在条件(C1)-(C6)下，基于矩的块金估计 ˆω、核平滑的半变异函数估计 ˆV_Y(r) 和交叉协方差函数估计 ˆC_XY(r) 都是一致的。收敛速度分别为 O_p(|S_n|^{-1/2}) 和 O_p(h_n + (|S_n|h_n)^{-1/2})。这是一个“令人惊讶的发现”：即使存在优先采样，这些经典的非参数估计量仍然有效。
  • Theorem 3（参数协方差）：在条件(C1)-(C6)和(C8)-(C10)下，基于最小对比度（MC）和复合似然（CL）的协方差参数估计量 ˆθ_n,MC 和 ˆθ_n,CL 都是一致的，收敛速度为 O_p(|S_n|^{-1/2})。这证明了即使不建模采样机制，也能正确估计空间协方差结构的参数。

• 证明路线与技术技巧：

  • 整体路线：

    1. 证明 ˆβ 的渐近性质（Theorem 1）：
       • Step 1: 偏差分析：利用Campbell定理和Stein引理（Lemma 1）计算 E[∑ w(s)Z(s)] 和 E[∑ w(s)w(s)^T]，得到 E[ˆβ] 的表达式，从而识别出偏差仅存在于截距项。
       • Step 2: 渐近正态性：将 ˆβ 视为一个估计方程的解。证明估计方程的期望为零（在 β*₀ 处），并计算其方差。利用递增域下的强混合条件（C2）和分块技巧（Lyapunov CLT），证明估计方程的中心极限定理。最后通过Slutsky定理得到 ˆβ 的渐近分布。
    2. 证明非参数估计的一致性（Theorem 2）：
       • Step 1: 分解：将核平滑估计量 ˆV_Y(r) 分解为 A_n(r), B_n(r), C_n(r) 三部分，分别对应真实残差项、ˆβ 的估计误差项和交叉项。
       • Step 2: 主项分析：利用Campbell定理和Lemma 1计算 A_n(r) 的期望和方差，证明其依概率收敛到目标积分。方差的主要部分来自二阶矩，其阶为 O(1/(|S_n|h_n))，在 |S_n|h_n → ∞ 时趋于0。
       • Step 3: 余项分析：证明 B_n(r) 和 C_n(r) 依概率收敛到0，且收敛速度比主项更快（因为 ˆβ 是 √n 一致的）。
       • Step 4: 偏差分析：通过泰勒展开，量化核平滑带来的渐近偏差，其阶为 O(h_n)。
       • Step 5: 合并：结合主项收敛、余项消失和偏差阶，得到 ˆV_Y(r) 的收敛速度和一致性。ˆC_XY(r) 的证明类似。
    3. 证明参数协方差估计的一致性（Theorem 3）：
       • Step 1: 无偏估计方程：证明在真实参数 θ₀ 下，估计方程 U*_n(θ₀) 的期望为零。这再次依赖于Lemma 1，证明 E[ ( [Z*(s)-Z*(t)]² - 2ζ(||s-t||; θ₀) ) exp(X(s)+X(t)) ] = 0。
       • Step 2: 方差衰减：计算 U*_n(θ₀) 的方差，证明其阶为 O(|S_n|³)，因此 Var[U*_n(θ₀)/|S_n|²] = O(1/|S_n|) → 0。这需要计算复杂的四阶矩，并利用Isserlis定理和Lemma 1进行化简。
       • Step 3: 泰勒展开：对估计方程进行泰勒展开，结合Step 1和Step 2，得到 ˆθ_n 的 √n 一致性。

  • 关键跳跃点：

    • Lemma 1（Stein引理的推广）：这是整个证明的基石。它给出了 E[Y exp(X)] 和 E[Y² exp(X)] 的显式表达式，其中 (X, Y) 是联合高斯的。这个引理将复杂的期望计算简化为协方差和矩的计算，是连接点过程（exp(X) 项）和标记过程（Y 项）的桥梁。
    • 证明 E[U*_n(θ₀)] = 0：这是证明参数协方差估计一致性的核心。它依赖于一个关键事实：[Z*(s)-Z*(t)]² 的条件期望（给定 X 和 Y 的联合分布）恰好等于 2ζ(||s-t||; θ₀)，而这个条件期望与 X 无关。因此，当对 X 取期望时，exp(X(s)+X(t)) 项被抵消，使得整体期望为零。这解释了为什么即使存在优先采样，基于残差平方和的估计方程仍然无偏。

  • 技术技巧点名：

    • Campbell定理：用于将点过程上的求和期望转化为空间积分，是处理随机采样位置的标准工具。
    • Stein引理（Lemma 1）：用于计算高斯随机变量指数矩的期望，是连接点过程和标记过程依赖的关键。
    • Isserlis定理（Wick定理）：用于计算高维高斯随机变量的高阶矩，在证明Theorem 3的方差项时用于展开复杂的四阶矩。
    • 强混合条件（C2）与分块技巧：用于在递增域渐近下证明中心极限定理，处理空间依赖数据。
    • Lyapunov中心极限定理：用于证明分块后的估计方程之和的渐近正态性。
    • 泰勒展开与Delta方法：用于从估计方程的一致性推导参数估计量的一致性。

• 真实例子与应用：

  • 数据：巴罗科罗拉多岛（Barro Colorado Island）50公顷森林动态样地中的“Trichilia tuberculata”树种数据。观测的是树木位置和胸径（DBH）。
  • 方法应用：
    1. 将 log(DBH - 9) 作为标记 Z(s)。
    2. 将地形坡度的平方根作为空间协变量 w(s)。
    3. 假设标记过程具有指数协方差函数 C_Y(r) = σ²_Y exp(-r/φ_Y)。
    4. 使用本文提出的MC和CL方法估计回归系数 β 和协方差参数 (σ²_Y, φ_Y, σ²_e)。
    5. 同时，使用非参数估计量（4）-（6）估计协方差和交叉协方差函数。
  • 结果：
    • MC和CL方法给出的参数估计非常接近，而忽略优先采样的MLE给出的估计值差异较大（特别是 β₁ 和 φ_Y）。
    • TMB方法在此数据集上遇到计算问题，作者推测是因为其假设 X(s) = γY(s) 被违反。
    • 估计的交叉协方差函数 ˆC_XY(r) 为负值，且其依赖范围与协方差函数 ˆC_Y(r) 明显不同。这支持了本文放宽线性依赖假设的必要性，并定量证实了“树木密度高的区域胸径较小”的生态学发现。
  • 这个例子想说明：本文方法在实际复杂数据中可行，且能处理TMB等参数化方法无法处理的更灵活的依赖结构。它展示了方法对模型误设的稳健性，并提供了比MLE更合理的生态学解释。

• 🔎 结论是否比证明窄：

  • Theorem 1 的结论是 ˆβ_n 收敛到 β*₀，其中截距有偏。作者在正文中声称“remains unbiased for all regression coefficients in β, except for the intercept term”。这是准确的，与定理陈述一致。
  • Theorem 2 的结论是 ˆV_Y(r) 和 ˆC_XY(r) 一致。作者在正文中声称“can still be consistently estimated”。这是准确的。
  • Theorem 3 的结论是 ||ˆθ_n - θ₀|| = O_p(|S_n|^{-1/2})。作者在正文中声称“establishes consistency”。这是准确的，但需要注意，这个一致性是在假设参数化半变异函数 ζ(r; θ) 正确指定（条件C8）的前提下成立的。如果参数模型误设，估计量会收敛到某个伪真值（pseudo-true value），而非 θ₀。作者在模拟中展示了在模型误设下（Scenario 2）方法仍优于TMB，但理论上的鲁棒性并未在定理中覆盖。这是一个值得注意的窄化。

四、开放问题

1. 扩展到非LGCP采样机制：本文的理论严格依赖于采样位置来自LGCP（对数高斯Cox过程）。作者在结论中提到了“repulsive point processes”（如Matérn硬核过程）作为未来方向。扎根于：Section 7, “it would be interesting to study the scenarios where the underlying point process moves beyond an LGCP, e.g. repulsive point processes”。这是一个明确的开放问题：当采样位置具有排斥性时，本文的“可忽略性”结论是否仍然成立？证明需要新的技术工具来处理非泊松的矩结构。

2. 空间变系数模型：本文假设回归系数 β 是常数。作者在结论中提到了“spatially varying regression coefficients”。扎根于：Section 7, “extend our proposed method to specialized geostatistical models under more complex preferential sampling mechanism, such that those with spatially varying regression coefficients”。这是一个自然的推广：当 β(s) 是空间变化的函数时，最小二乘估计的偏差结构会如何变化？如何用非参数方法估计 β(s) 并同时校正优先采样偏差？

3. 效率提升与最优权重：本文的方法（特别是CL和MC）比MLE/TMB方差更大。作者在模拟中承认了这一点。一个开放问题是：能否在不建模采样机制的前提下，通过选择最优权重函数 w(s,t) 来提高估计效率？扎根于：Section 4, “A good choice of w(s,t) can improve the efficiency of the resulting estimators” 以及 Section 5.1, “Our method exhibits slightly higher variance compared to MLE and TMB”。这指向一个半参数效率问题：在给定 C_XY(r) 非参数的情况下，估计 θ 的有效界是什么？本文的CL和MC估计量是否达到了这个界？

4. 与因果推断中“可忽略性”的深层联系：本文的标题和核心概念“Ignorability”直接借用了因果推断的术语。在因果推断中，“可忽略性”（ignorability）或“无混杂性”（unconfoundedness）是指给定协变量后，处理分配独立于潜在结果。本文的“可忽略性”是指给定潜高斯场后，采样位置独立于标记过程？还是指在特定估计方法下，采样机制可以被忽略？扎根于：本文的整个理论框架。一个值得探索的问题是：能否将本文的“可忽略性”条件（即 C_XY(r) 的某种结构）与因果推断中的标准可忽略性假设（Y(1), Y(0) ⟂ T | X）建立形式化的联系？这可能会为两个领域之间的方法迁移打开大门。

---

Maintained by 陈星宇 · Homepage · Source on GitHub

评论