Accounting for overdispersion and clustering in binomial data from N-of-1 trials

技术摘要：N-of-1 试验中二项分布数据的过度离散与聚类处理

问题陈述
N-of-1 试验是针对特定慢性病患者设计的以患者为中心的、随机、双盲、交叉实验。虽然其主要目标是评估特定患者的疗效，但通过汇总一组 N-of-1 试验的结果来获得群体层面的疗效估计并评估治疗个体化是否成功，具有显著价值。然而，用于汇总二项分布数据的标准统计方法往往无法充分考虑 N-of-1 设计中固有的两个关键特征：

1. 过度离散（Overdispersion）： 二元结果的方差往往超过了标准二项分布所规定的均值-方差关系。
2. 层级聚类（Hierarchical Clustering）： 由于患者内部的重复治疗阶段，测量值在患者内部是相关的（患者内相关性）。

虽然针对此类情况（特别是对于二元结果）已提出了贝叶斯方法，但目前缺乏能够同时严格处理过度离散和聚类的频率学方法。本文通过提出并比较两种用于汇总 N-of-1 试验二项分布数据的频率学分析策略来填补这一空白。

方法论
作者提出了两种不同的频率学方法来估计治疗偏好或成功率的汇总比例 ($\theta$)。

• 方法 1：拟似然（Quasi-Likelihood, QL）法（仅处理过度离散）
  该方法使用 Beta-Binomial 分布对数据进行建模，以解释过度离散 ($\rho$)，但在处理聚类结构时将患者视为独立的单位。

  • 模型： $X_j \sim \text{BetaBinom}(n_j, \theta, \rho)$，其中 $X_j$ 是患者 $j$ 在 $n_j$ 个阶段中的成功次数。
  • 估计： 使用拟似然估计量计算汇总比例 $\theta$：$\hat{\theta}_{QL} = \frac{\sum c_j X_j}{\sum c_j n_j}$，其中权重 $c_j$ 取决于过度离散参数 $\rho$。
  • 方差估计： 标准误通过 Jackknife 方法导出，以确保一致性。
  • 估计 $\rho$： 本文评估了若干种用于估计组内相关系数（ICC）参数 $\rho$ 的估计量，包括方差分析（AOV）、Fleiss-Cuzick（FC）、Pearson（PE）、扩展拟似然（EQL）、伪似然（PL）和矩估计法（MM）。作者指出，某些估计量（如 AOV）的闭式方差公式可能会产生负值或偏差，因此他们依赖于基于 Jackknife 的方差估计进行比较。

• 方法 2：随机效应法（处理过度离散与聚类）
  该方法改编自 Molenberghs 等人的研究，利用层级随机效应模型同时处理过度离散和患者内部重复测量的相关性。

  • 模型： $X_j \sim \text{Binom}(n_j, \theta_j)$，其中 $\theta_j = \pi_j \kappa_j$。这里，$\pi_j$（通过 Beta 分布建模）解释了剩余的过度离散，而 $\kappa_j$（源自正态分布的随机效应）则解释了层级聚类。
  • 估计： 参数通过数值积分进行极大似然（ML）估计。最终的汇总估计量由随机成分乘积导出：$\hat{\theta} = \hat{E}(\pi_j) \times \hat{E}(\kappa_j)$。
  • 近似： 正态随机变量的 Sigmoid 函数期望值使用带有缩放因子 ($\lambda \approx 0.588$) 的标准正态累积分布函数进行近似。

模拟研究与结果
作者进行了一项包含 1,000 次迭代的模拟研究，旨在比较这些方法在不同真实比例 ($\theta$) 和过度离散 ($\rho$) 情景下，以及不同样本量（$J$ 名患者和 $n_j$ 个阶段）下的表现。

• $\rho$ 估计量的表现： 在方法 1 中，扩展拟似然（EQL）估计量在 $\theta = 0.50$ 时表现出较高的均方误差（MSE）和偏差。然而，尽管 EQL 对 $\rho$ 的估计表现不佳，其得到的汇总比例估计值 ($\hat{\theta}$) 的表现却与其他估计量相当，甚至有时更好。作者将其归因于一种抵消效应，即 $\hat{\rho}$ 的偏差降低了 Delta 方法展开式中 $\text{Var}(\hat{\theta})$ 中 $\rho$ 项的方差贡献。
• 方法 1 与方法 2 的对比：
  • MSE： 方法 1 的估计量通常表现出稳定的 MSE。方法 2 在某些情景下具有竞争力的 MSE，但在 $\theta = 0.75$ 时表现较差。
  • 覆盖概率： 方法 1 在所有情景下都能保持接近名义 95% 水平的覆盖概率。方法 2（随机效应法）在 $\theta = 0.75$ 时显示出显著的覆盖率下降（低至约 78%）。作者将其归因于广义线性混合模型（GLMMs）中方差估计值的偏态抽样分布，这使得对称置信区间在这些特定情景下并不适用。
• 样本量效应： 增加患者数量 ($J$) 比增加每位患者的阶段数 ($n_j$) 能更有效地降低 MSE。

实际数据应用
这些方法被应用于三个现实世界的 N-of-1 试验数据集：

1. 阿米替林治疗纤维肌痛： 所有方法得出的偏好汇总估计值相似（$\hat{\theta} \approx 0.67$）。未检测到显著的过度离散（$\rho$ 的置信区间包含零）。
2. 非甾体抗炎药对比对乙酰氨酚治疗骨关节炎： 结果在各方法间保持一致（总评分 $\hat{\theta} \approx 0.42$）。未发现显著的过度离散。
3. 茶碱治疗不可逆性 CAL： 汇总估计值保持一致（$\hat{\theta} \approx 0.51$）。同样，没有显著证据表明存在过度离散。

在所有实际数据案例中，尽管过度离散参数的估计值有所不同，但汇总比例的点估计在不同方法之间表现得非常一致。

意义与主张
本文声称提供了首个专门处理二项分布数据中过度离散和层级聚类的频率学方法。

• 方法论贡献： 作者引入了一个用于处理过度离散的拟似然框架，并改编了一个用于同时处理过度离散和聚类的随机效应框架。
• 解释性区别： 本文强调了两种方法在概念上的区别。方法 1 得出的是“群体平均”（Population-averaged, PA）估计，而方法 2 得出的是“受试者特异性”（Subject-specific, SS）估计。鉴于 N-of-1 试验本质上关注个体差异，作者认为方法 2（SS）的解释可能更符合 N-of-1 的语境，尽管其计算复杂度较高且在偏态情景下偶尔会出现覆盖率问题。
• 实际效用： 研究表明，虽然过度离散参数 ($\rho$) 的估计可能不稳定（尤其是使用 EQL 时），但由此产生的汇总比例估计值 ($\hat{\theta}$) 在不同的估计策略下都是稳健的。作者得出结论，在实际应用中，方法的选择取决于推断目标是群体层面还是个体层面；当担心方法 2 的覆盖特性时，方法 1 提供了一个更简单且稳健的替代方案。