9 发表偏倚

来源页面：https://doing-meta.guide/pub-bias

回顾前几章，我们已经介绍了大量 Meta 分析技术。我们不仅学习了如何合并效应量，也已经知道如何评估结果的稳健性、考察异质性模式，以及检验效应为何会彼此不同。

所有这些方法都能帮助我们从 Meta 分析中得出有效结论。然而，这些结论实际上依赖于一个此前尚未被我们认真质疑的隐含前提：也就是我们收集到的数据是全面的，至少也是对目标研究领域具有代表性的。

在前文关于研究检索的部分，我们提到 Meta 分析通常试图纳入所有可获得证据，以便得到一个能够恰当描述整个研究领域的综合效应量。从统计角度看，我们或许可以容忍分析中遗漏少量研究，但前提只能是这些研究是出于偶然原因而“缺失”的。

遗憾的是，Meta 分析往往无法纳入所有现有证据。更糟的是，也有充分理由认为，我们收集到的数据中有些研究并不是“完全随机地”缺失。现实世界并不完美，支配科学实践的激励与“规则”也同样如此。这意味着，某些系统性偏倚会决定一项研究最终是否进入我们的 Meta 分析。

一个很好的例子来自药物治疗研究中一则不算遥远的轶事。早在 20 世纪 90 年代，人们普遍认为抗抑郁药物，例如选择性 5-羟色胺再摄取抑制剂（selective serotonin re-uptake inhibitors, SSRIs），对抑郁患者是有效的。这一结论的大部分证据来自已发表药物治疗试验的 Meta 分析，在这些试验中，抗抑郁药被拿来与安慰剂药片比较。考虑到抗抑郁药市场价值数十亿美元且持续增长，这一问题的重要性不言而喻。

这也就有助于理解 Irving Kirsch 及其同事发表题为 The Emperor’s New Drugs 的文章时，为何会引起巨大震动。该文认为，事情也许远没有人们想象得那样乐观。

Kirsch 及其同事依据《信息自由法》获得了此前未发表的抗抑郁药试验数据，这些数据原先由制药公司提交给美国食品药品监督管理局（FDA）。他们发现，一旦把这些未发表数据也纳入考虑，抗抑郁药相较安慰剂的获益充其量也只是极小，而且在临床上几乎可以忽略。Kirsch 及其同事认为，其原因在于公司只发表了结果有利的研究，而把那些“令人失望”的研究搁置不发。

围绕这一问题的争论随后愈演愈烈，直到今天 Kirsch 的主张仍然存在争议。我们之所以选择这个例子，并不是为了站队，而是为了说明：缺失研究可能对 Meta 分析推断的有效性构成严重威胁。在 Meta 分析文献中，这类问题通常被概括为发表偏倚（publication bias）。

发表偏倚问题提醒我们：Meta 分析中的任何发现，其质量都不会高于其所依据的数据本身。Meta 分析技术只能处理手头已有的数据。因此，如果收集到的数据本身就是扭曲的，那么再好的统计模型也只会复制这种内在偏倚。你也许还记得，本书一开始在讨论“文件抽屉问题”时，我们已经提到过这一根本局限。事实上，“文件抽屉问题”和“发表偏倚”常常被当作同义词使用。

发表偏倚及相关问题对 Meta 分析结果的影响可能非常巨大。它可能导致我们高估治疗效果、忽视不良副作用，或者强化那些实际上并不成立的理论。

因此，在本章中，我们将讨论发表偏倚扭曲研究结论的多种形式；也会介绍若干 Meta 分析者可用于评估数据中发表偏倚风险的方法，以及如何尽量减轻这类偏倚。

9.1 什么是发表偏倚？

发表偏倚是指：一项研究被发表的概率会受到其结果影响。大量证据表明，如果研究结果达到统计学显著，或者支持原始假设，那么它更可能进入公共领域。

在检索符合条件的研究时，我们通常只能依赖某种形式上已经公开的证据，例如同行评议论文、预印本、书籍，或其他可获取报告。一旦存在发表偏倚，这不仅意味着我们的数据集中缺失了一些研究，更意味着这些缺失的研究很可能恰恰是不利结果的研究。

Meta 分析技术能够帮助我们得到总体平均效应量的无偏估计。但如果我们的样本本身就是扭曲的，那么即便这个效应估计在统计意义上是“正确”的，它也无法代表现实。这就像试图估计一座冰山的大小，却只测量其露出水面的尖端：即使我们对水面以上的高度测得极其精确，最终结论仍然不可避免地是错误的。

严格说来，发表偏倚其实只是众多未充分报告偏倚中的一种。还有若干其他因素同样会扭曲我们在 Meta 分析中得到的证据，包括：

• 引文偏倚（citation bias）：即便已经发表，结果为阴性或不确定的研究也较少被相关文献引用，因此更难通过参考文献追溯等方式被发现。
• 时间滞后偏倚（time-lag bias）：结果为阳性的研究往往比结果不利的研究更早发表。这意味着，近期完成的研究中，阳性结果常常已经可见，而不显著结果尚未公开。
• 多重发表偏倚（multiple publication bias）：结果“成功”的研究更可能在多篇期刊文章中重复报告，因此至少有一篇被检索到的概率更高。这种把同一研究结果拆分到多篇论文中的做法，也被称为 “salami slicing”。
• 语言偏倚（language bias）：在大多数学科中，英语是发表研究证据的主要语言。其他语言的出版物更不容易被识别，尤其当研究者自己无法在不借助翻译的情况下理解其内容时更是如此。如果英文研究系统性地不同于其他语言研究，这也会带来偏倚。
• 结局报告偏倚（outcome reporting bias）：许多研究，尤其是临床试验，会测量不止一个结局。有些研究者会只报告那些得到正向结果的结局，而把不支持假设的结局删去。这同样会导致偏倚：从技术上看，研究确实发表了，但由于不利结果未被报告，它仍会在我们的 Meta 分析中“缺席”。

这些未充分报告偏倚可以被看作使我们更难发现既有证据的系统性因素。然而，即便我们真的能纳入所有相关研究，结果也未必就没有问题。偏倚也可能来自研究者在分析和报告研究时所采用的可疑研究实践（questionable research practices, QRPs）。

前文讨论“研究者自由度”时，我们已经提到过这一概念。所谓 QRP，可以理解为研究者滥用这些自由度，使结果朝自己想要的方向“弯曲”。遗憾的是，目前并没有关于什么构成 QRP 的完全一致定义，但有几个例子经常被提到。

最著名的一种 QRP 是 p-hacking，也就是不断调整分析，直到结果达到传统显著性阈值 \(p < 0.05\)。这可能涉及离群值剔除方式、亚组分析方式，或者缺失数据处理方式等。

另一种 QRP 是 HARKing，即 hypothesizing after the results are known，意思是“在知道结果之后才提出假设”。其中一种形式是：把探索性分析中偶然发现的结果，事后包装成研究一开始就设定好的先验假设。比如，研究者可能对同一数据集做了许多不同检验，然后再围绕那些恰好显著的结果“发明”出理论假设。这种做法存在严重缺陷，会抬高研究的假阳性发现率，从而增加虚假发现的风险。另一种 HARKing 则是删除所有未被数据支持的原始假设，最终也可能导致结局报告偏倚。

9.2 在 Meta 分析中处理发表偏倚

很明显，发表偏倚、其他报告偏倚以及 QRPs，都会对 Meta 分析的有效性造成强烈而有害的影响。之所以构成重大挑战，是因为我们通常几乎不可能准确知道偏倚到底有多大，甚至连偏倚是否存在都难以断定。

在 Meta 分析中，我们可以采用一些技术，在一定程度上降低发表与报告偏倚以及 QRPs 带来的扭曲风险。其中有些方法属于研究检索层面，另一些则是统计方法。

• 研究检索。前文关于研究检索的部分已经讨论了如何搜索符合条件的研究。如果存在发表偏倚，这一步尤为关键，因为仅检索已发表文献可能得到一个并不能充分代表全部证据的数据集。我们可以通过进一步检索灰色文献来缓解这个问题，灰色文献包括学位论文、预印本、政府报告或会议论文等。值得庆幸的是，预注册在许多学科中也越来越普遍。这使得我们可以在 ICTRP 或 OSF Registries 等研究注册平台中寻找带有未发表数据的研究，并联系作者索取尚未公开的数据。（例如 Mahmood 及其同事 2014 年的论文就详细介绍了如何开展全面的灰色文献检索以及其中可能遇到的挑战，而且该文可在线开放获取。）灰色文献检索可能既繁琐又令人沮丧，但它是值得的。一项大型研究发现，纳入灰色文献和未发表文献有助于避免高估真实效应。
• 统计方法。我们也可以借助统计程序来考察发表偏倚的存在。这些方法没有一种能够直接“识别”发表偏倚本身，但它们可以检验数据中某些可能暗示偏倚存在的性质。有些方法还可以在校正发表偏倚后，估计真实的总体效应。

本章接下来将展示若干常见的统计方法，用于评估和控制发表偏倚。我们先从聚焦于小样本研究效应（small-study effects）的方法开始。它们的共同点在于：通过考察研究精确度与观察到的效应量之间的关系，来寻找发表偏倚的线索。

9.2.1 小样本研究效应方法

Meta 分析中有多种小样本研究效应方法可用于评估并校正发表偏倚。这些技术中不少已经沿用了很多年。顾名思义，这些方法尤其关注小样本研究。从统计角度看，这通常意味着标准误较高的研究。小样本研究效应方法假定，小研究更容易受到发表偏倚影响。

这种假设基于三个核心想法：

1. 大型研究往往需要大量资源和时间投入，因此无论结果是否显著，都更有可能被发表。
2. 中等规模研究面临更高的未发表风险。不过，即便统计功效只是中等水平，也往往足以产生显著结果。因此，只有部分这类研究会因为产生了“不理想”的非显著结果而不被发表。
3. 小样本研究最容易产生非显著结果，因此也最可能停留在“文件抽屉”里。对小研究而言，只有非常大的效应才可能达到显著，因此最终发表出来的往往只有那些效应量特别高的小研究。

可见，这一假设背后的机制其实相当简单：它认为发表偏倚之所以存在，是因为只有显著结果才被发表。既然获得显著结果的概率会随着样本量增大而上升，那么发表偏倚就会不成比例地影响小研究。

9.2.1.1 漏斗图

前文讨论效应量时我们已经知道，研究样本量与标准误密切相关。标准误越大，效应量的置信区间就越宽，也越不容易达到统计学显著。因此，很自然可以推断：小样本研究效应主要会影响标准误较大的研究。

假设我们收集到的数据受到发表偏倚影响。在这种情况下，可以预期：标准误较大的研究，其效应量通常高于标准误较小的研究。原因在于，那些效应较低的小研究往往没有达到显著水平，因此从未进入发表流程，自然也不会被纳入我们的 Meta 分析。

通常，人们通过漏斗图（funnel plot）来检视小样本研究效应。漏斗图本质上是一个散点图：横轴是研究的观察效应量，纵轴是其标准误的某种度量。通常纵轴会倒置，也就是说，图中纵轴位置越高，代表标准误越小。

如果不存在发表偏倚，数据点大致应构成一个对称的倒漏斗形状，这也是其名称的由来。位于图上部的研究（标准误较低）应该彼此更接近，也更接近合并效应量；而在下部，随着标准误增大，漏斗会逐渐“张开”，效应量在合并效应的左右两侧分布得更分散。

如果回想前文关于效应量行为的讨论，以及固定效应模型的介绍，就更容易理解为什么研究会形成漏斗。标准误反映的是研究的精确度：标准误越小，我们越期待观察效应量能够更好地估计真实效应量；而标准误较大时，效应量精度较低，因此更可能偏离总体中的真实效应。

下面我们用实际代码生成一个漏斗图。在 {meta} 包中，meta::funnel 函数可用于为一个 Meta 分析对象绘制漏斗图。（在本章中我们始终写成 meta::funnel。这并非绝对必要，但在实践中有时能减少报错和混淆，因为 {metafor} 包中也有另一个 funnel 函数。）这里我们基于 m.gen 这一 Meta 分析对象来绘图。我们还使用了两个参数：xlim 用于控制横轴范围，studlab 用于显示研究标签。在运行 meta::funnel 后再调用 title，就可以为图形添加标题。

代码如下：

# 加载 'meta' 包
library(meta)

# 绘制漏斗图
meta::funnel(m.gen,
             xlim = c(-0.5, 2),
             studlab = TRUE)

# 添加标题
title("Funnel Plot (Third Wave Psychotherapies)")

如前所述，这个漏斗图横轴显示每项研究的效应量（这里是标准化均值差，SMD），纵轴显示标准误（从大到小）。为了便于解释，图中还叠加了一个理想化的漏斗轮廓，表示我们希望研究遵循的分布模式。漏斗中央的垂直线表示平均效应量。由于 m.gen 是在随机效应模型下得到的，因此漏斗图也采用随机效应估计值。

如果不存在小样本研究效应，研究点大致应围绕该漏斗形状分布。在这个例子中是否如此呢？坦率地说，并不完全是。虽然标准误较低的研究确实更集中地围绕估计的真实效应，但整体模式看起来并不对称。原因在于，图右下角有三项效应量非常高的小研究（Shapiro、Kang 和 Danitz-Orsillo）。

然而，图左下角并没有与之对应的研究来“平衡”这种模式。也就是说，我们没有看到那些效应量非常低甚至为负的小研究。另一个令人担忧的细节是：样本中精确度最高的研究——de Vibe——也似乎并不符合理想漏斗模式，它的效应量明显比预期更小。

总体来看，这个数据集在漏斗图中呈现出一种不对称模式，可能提示存在发表偏倚。也许这三项小研究只是因为效应足够大、恰好达到显著，所以才得以发表；而另一些标准误类似但效应更小、因此不显著的小研究，则可能仍然躺在文件抽屉中，因而没有被纳入我们的 Meta 分析。

为了考察不对称模式与统计显著性之间的关系，一个很好的做法是绘制轮廓增强漏斗图（contour-enhanced funnel plot）。这类图有助于区分发表偏倚和其他形式的不对称。轮廓增强漏斗图用不同颜色标出图中各研究所处的显著性区域。在 meta::funnel 中，只需向 contour 参数提供所需的显著性阈值即可。通常取 0.9、0.95 和 0.99，分别对应 \(p < 0.1\)、\(p < 0.05\) 和 \(p < 0.01\)。通过 col.contour 可以指定这些轮廓区域的颜色；随后再用 legend 添加图例，说明不同颜色的含义。

代码如下：

# 定义轮廓区域填充色
col.contour = c("gray75", "gray85", "gray95")

# 生成漏斗图（这里不显示研究标签）
meta::funnel(m.gen, xlim = c(-0.5, 2),
             contour = c(0.9, 0.95, 0.99),
             col.contour = col.contour)

# 添加图例
legend(x = 1.6, y = 0.01,
       legend = c("p < 0.1", "p < 0.05", "p < 0.01"),
       fill = col.contour)

# 添加标题
title("Contour-Enhanced Funnel Plot (Third Wave Psychotherapies)")

现在，我们看到漏斗图中出现了三个着色区域。我们尤其关注 \(p < 0.05\) 和 \(p < 0.01\) 区域，因为落在这些区域内的效应量传统上被视为显著。

添加轮廓区域后，图形更具启发性：它表明那三项小研究虽然标准误很大，但其效应都达到了统计学显著；而标准误相近但不显著的研究只有一项。如果我们想通过“填补”左下角缺失研究来增加对称性，那么这些假想研究大多会落在图中的非显著区域，甚至可能具有显著的负向效应。

对于较大研究，模式略有不同。我们确实看到若干研究的 \(p > 0.05\)，而且效应分布也没有那么偏向一侧。不过，潜在问题在于：虽然它们并未严格达到显著，但除一项外，其余研究都非常接近显著性阈值，也就是落在 \(0.1 > p > 0.05\) 区域。可能的解释是，这些研究在原始论文中使用了不同的效应量计算方法，从而得到了显著结果；也可能只是“趋势显著”就已经足以让研究发表。

总之，对轮廓增强漏斗图的检查支持了我们最初的判断：漏斗图存在不对称，而这可能由发表偏倚造成。不过，切忌草率下结论，必须谨慎解释漏斗图，因为发表偏倚只是漏斗图不对称的众多可能原因之一。

漏斗图不对称的其他解释

虽然发表偏倚会导致漏斗图不对称，但也存在一些较为“无害”的原因，可能产生类似模式：

• 不对称也可能由研究间异质性造成。漏斗图假定效应量分散主要来自抽样误差，但并未控制各研究可能在估计不同真实效应这一事实。
• 也可能是小研究在实施过程上与大研究不同，从而产生更高效应。例如在临床研究中，样本量较小时更容易保证每位参与者都严格接受既定治疗；而在大型研究中，这种治疗忠实度（treatment fidelity）可能下降，从而导致效应变小。此时，检查纳入研究的具体特征就很有意义。
• 低质量研究往往更容易报告较大效应，因为其偏倚风险更高。大型研究投入更大，方法学通常也更严谨。这同样可能造成漏斗图不对称，即便并不存在发表偏倚。
• 最后，漏斗图不对称也完全可能只是偶然现象。

可见，单靠目测（尤其是轮廓增强）漏斗图，确实已经可以发现一些“红旗”，提示结果可能受到发表偏倚影响。

然而，仅靠肉眼解释漏斗图显然也有局限。并不存在一个明确规则告诉我们什么程度的不对称算“过度不对称”，因此漏斗图推断总带有一定主观性。于是，采用定量方法来评估漏斗图不对称就很有帮助，这通常就是通过 Egger 回归检验完成的。

9.2.1.2 Egger 回归检验

Egger 回归检验是常用的定量方法，用于检验漏斗图中的不对称。与目测漏斗图一样，它只能识别小样本研究效应，不能直接证明发表偏倚存在。该检验基于一个简单线性回归模型：

这个模型中的响应变量 \(y\) 是 Meta 分析中观测效应量 \(\hat\theta_k\) 除以其标准误之后的结果。所得值等价于 \(z\) 分数，可直接告诉我们效应是否显著：当 \(z \ge 1.96\) 或 \(z \le -1.96\) 时，就意味着 \(p < 0.05\)。该响应变量回归于研究标准误的倒数，也就是研究的精确度。

不过，在 Egger 检验中，我们关注的并不是回归系数 \(\beta_1\) 的大小和显著性，而是截距 \(\beta_0\)。为了评估漏斗图不对称，我们要考察 \(\hat\beta_0\) 的大小，以及它是否显著偏离零。若显著偏离，则 Egger 检验提示存在漏斗图不对称。

为什么回归截距能告诉我们漏斗图是否不对称？因为在线性回归中，截距表示当其他预测变量都为零时 \(y\) 的取值。在这里，预测变量是研究的精确度，因此截距表示当精确度为零时（也就是标准误无限大时）预期的 \(z\) 分数。

下面几幅图展示了漏斗图不对称如何影响 Egger 检验背后的回归斜率与截距。

接下来看看，如果把这种回归模型拟合到 m.gen 的数据上，会得到怎样的结果。在 R 中，我们可以从 m.gen 中提取原始数据，构造响应变量 y 与预测变量 x。下面的代码用到了管道和 {tidyverse} 中的 mutate，最后通过线性模型函数 lm 让 \(z\) 分数 y 回归于精确度 x，再输出 summary。

# 加载所需包
library(tidyverse)

m.gen$data %>%
  mutate(y = TE/seTE, x = 1/seTE) %>%
  lm(y ~ x, data = .) %>%
  summary()

## [...]
##
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)   4.1111     0.8790   4.677 0.000252 ***
## x            -0.3407     0.1837  -1.855 0.082140 .
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
##
## [...]

从结果可见，回归模型的截距为 \(\hat\beta_0 = 4.11\)。它显著大于零（\(t = 4.677\), \(p < 0.001\)），提示漏斗图中的数据确实存在不对称。总体而言，这支持了我们关于存在小样本研究效应的初步判断。不过仍需再次强调：这并不能确定该模式就是由发表偏倚造成的。

更方便的做法，是直接使用 {meta} 包中的 metabias 函数来执行 Egger 截距检验。该函数只需要一个 Meta 分析对象，并把 method.bias 设为 "linreg" 即可。如果应用到 m.gen 上，将得到与前面相同的结果。

metabias(m.gen, method.bias = "linreg")

如何报告 Egger 检验结果

通常只需报告截距值、其 95% 置信区间，以及 \(t\) 值和 \(p\) 值即可。{dmetar} 包中提供了一个便捷函数 eggers.test，它是 metabias 的封装，能以更适合报告的格式返回结果。若你没有安装 {dmetar}，其源码可在线查看：https://raw.githubusercontent.com/MathiasHarrer/dmetar/master/R/eggers.test.R。

eggers.test(m.gen)

项目	截距	置信区间	t	p
Egger’s test	4.111	2.347-5.875	4.677	0.00025

m.gen 使用的效应量指标是经小样本偏倚校正的 SMD（Hedges’ \(g\)）。已有研究指出，直接对 SMD 运行 Egger 检验可能会提高假阳性率，因为标准化均值差与其标准误并不独立。

这一点可以从组间 SMD 标准误的公式直接看出。该公式中包含 SMD 本身，这意味着研究的标准误会随观察效应量大小而变化，也就是说，SMD 与其标准误之间存在一种人为相关。

Pustejovsky 和 Rodgers提议，在检验标准化均值差的漏斗图不对称时，使用一个修正标准误。在这个版本中，只保留原标准误公式的第一部分，因此观测效应量本身会从公式中消失：

其中，\(SE^*_{\text{SMD}_{\text{between}}}\) 表示修正后的标准误。检查采用这一改进后 Egger 检验是否仍给出相同结论，是很有意义的。下面的代码中，我们把每项研究实验组与对照组的样本量加入数据集，计算修正标准误，并据此重新运行分析。

# 添加实验组 (n1) 和对照组 (n2) 样本量
n1 <- c(62, 72, 44, 135, 103, 71, 69, 68, 95,
        43, 79, 61, 62, 60, 43, 42, 64, 63)

n2 <- c(51, 78, 41, 115, 100, 79, 62, 72, 80,
        44, 72, 67, 59, 54, 41, 51, 66, 55)

# 计算修正后的标准误
ThirdWave$seTE_c <- sqrt((n1+n2)/(n1*n2))

# 用修正后的标准误重新运行 metagen，得到新的 Meta 分析对象
m.gen.c <- metagen(TE = TE, seTE = seTE_c,
                   studlab = Author, data = ThirdWave, sm = "SMD",
                   fixed = FALSE, random = TRUE,
                   method.tau = "REML", hakn = TRUE,
                   title = "Third Wave Psychotherapies")

# Egger 检验
metabias(m.gen.c, method = "linreg")

可以看到，虽然具体数值有所不同，但结果解释并没有改变，这说明我们先前的发现具有一定稳健性。

直接在 metabias 中使用 Pustejovsky-Rodgers 方法

在 {meta} 的较新版本中，metabias 也可以通过把 method.bias 设为 "Pustejovsky"，直接使用 Pustejovsky 与 Rodgers提出的修正标准误公式来执行 Egger 检验。

不过，这只有在 {meta} 的 Meta 分析对象中已经包含实验组与对照组样本量（分别保存在 n.e 和 n.c 元素中）时才可行。像本例中的 metagen 对象，通常并不会自带这些值，因此需要手动补上。例如：

m.gen$n.e = n1; m.gen$n.c = n2
metabias(m.gen, method.bias = "Pustejovsky")

需要注意的是，在这一设置下，metabias 所使用的模型写法与前述公式等价，但并不完全相同。metabias 把修正标准误作为模型预测变量，并用逆方差作为权重；而我们上面示范的做法，则是在等式两侧都使用了修正标准误。因此，两种方式得到的结果不会完全一致。

9.2.1.3 Peters 回归检验

效应量与标准误之间的依赖关系并不仅存在于标准化均值差中。基于二分类结局数据的效应量，例如（对数）比值比、风险比或率，也存在类似的数学关联。

为了避免在二分类结局数据上使用 Egger 检验时出现假阳性膨胀，我们可以改用 Peters 及其同事提出的另一种回归检验。Peters 检验中，对数转换后的效应量回归于样本量的倒数：

这里，\(\log\psi_k\) 可以代表任何基于二分类结局数据的对数效应量（如 log odds ratio），而 \(n_k\) 是第 \(k\) 项研究的总样本量。

在拟合这一回归模型时，每项研究 \(k\) 还会被赋予一个不同的权重 \(w_k\)，该权重取决于样本量和事件计数。因此，这是一种加权线性回归，与 Meta 回归类似，但并不完全相同。权重公式为：

其中，\(a_k\) 为治疗组事件数，\(c_k\) 为对照组事件数；\(b_k\) 与 \(d_k\) 分别为治疗组和对照组的非事件数。与 Egger 检验不同，Peters 检验使用的是 \(\beta_1\) 而不是截距来判断漏斗图不对称。如果统计检验表明 \(\beta_1 \neq 0\)，我们就可以认为数据中存在不对称。

如果我们的 Meta 分析是使用 metabin 或 metaprop 得到的，那么可以通过 metabias 执行 Peters 检验，只需提供合适的 Meta 分析对象，并将 method.bias 设为 "peters"。下面，我们对前文中基于风险比构建的 m.bin 进行检验：

metabias(m.bin, method.bias = "peters")

输出结构与 Egger 检验非常相似。结果会告诉我们这是一个“基于样本量的回归检验”，也就是使用了 Peters 方法。该检验不显著（\(t = -0.08\), \(p = 0.94\)），因此没有发现漏斗图不对称。

漏斗图不对称检验的统计功效

一般只建议在 Meta 分析包含足够多研究时，才检验漏斗图不对称。当研究数量太少时，Egger 或 Peters 检验的统计功效可能不足以识别真实存在的不对称。通常建议至少在 \(K \ge 10\) 时才进行这类检验。

默认情况下，当研究数量少于这一阈值时，metabias 会直接报错。不过，虽然并不推荐，仍可通过降低函数中的 k.min 参数来强制执行。

9.2.1.4 Duval 与 Tweedie 的 trim-and-fill 方法

到这里，我们已经了解了几种用于检查 Meta 分析中小样本研究效应的方法。知道数据中可能存在发表偏倚固然重要，但我们真正更关心的是这种偏倚到底有多大。我们想知道，它只是轻微扭曲了综合效应估计，还是已经大到足以改变结论解释。

也就是说，我们需要一种能计算偏倚校正后真实效应量的方法。不过我们已经知道，发表偏倚本身并不能被直接测量。我们只能借助小样本研究效应这一“替代指标”，它至多只能指向发表偏倚的可能性。

因此，我们真正能校正的是小样本研究效应，而不是发表偏倚“本身”。如果效应量不对称的确是由发表偏倚造成的，那么校正这种不平衡，就能得到一个在纳入全部证据时更接近真实情况的估计。

最常见的校正漏斗图不对称的方法之一，就是 Duval 与 Tweedie 的 trim-and-fill 方法。它的核心思想很简单：不断“填补”假想缺失研究，直到漏斗图恢复对称。由此扩展后的数据集所得到的合并效应量，就被看作校正后的小样本研究效应估计。其算法包含两个步骤：

• Trimming（修剪）：首先，算法识别漏斗图中的离群研究。在我们前面的例子中，这些就是散落在图右侧的小研究。一旦识别出这些研究，就将其从分析中移除，并重新计算合并效应。这一步通常使用固定效应模型。
• Filling（填补）：接下来，把重新计算得到的合并效应视为所有效应量的中心。对于每一个被修剪掉的研究，都在漏斗另一侧添加一个镜像研究。例如，如果重新估计的平均效应是 0.5，而某被修剪研究效应量为 0.8，那么其镜像研究的效应量就设为 0.2。对所有被修剪研究完成该操作后，漏斗图将大致恢复对称。随后在包括原始研究、修剪研究与填补研究在内的全部数据上再次计算平均效应（通常使用随机效应模型），其结果就被看作校正后的综合效应估计。

需要特别注意的是，当研究间异质性很大时，trim-and-fill 方法并不能给出可靠结果。因为如果各研究并不共享同一个真实效应，那么即使是大样本研究也可能明显偏离平均效应。这样一来，这些研究也可能被算法当作需要“修剪”和“填补”的对象，但它们实际上未必真的受到发表偏倚影响，这显然会导致无效结果。

在 {meta} 中，我们可以用 trimfill 函数把这一算法应用到数据上。该函数默认设置已经相当合理，所以通常只需把 Meta 分析对象传给它即可。下面仍以 m.gen 为例。不过在开始之前，先看一下这个 Meta 分析中的 \(I^2\) 异质性大小：

m.gen$I2

我们会发现，\(I^2 = 63\%\)，异质性已经相当可观。结合 trim-and-fill 在高异质性数据中的局限，这一点值得警惕。

因此，我们会做两次 trim-and-fill 分析：一次使用全部研究；另一次作为敏感性分析，去除前文识别出的两个离群值（第 3 和第 16 项研究）。结果分别保存为 tf 和 tf.no.out。

# 使用全部研究
tf <- trimfill(m.gen)

# 去除离群值后再分析
tf.no.out <- trimfill(update(m.gen,
                             subset = -c(3, 16)))

先看看包含全部研究的分析结果：

summary(tf)

## Review:     Third Wave Psychotherapies
##                              SMD             95%-CI %W(random)
## [...]
## Filled: Warnecke et al.   0.0520 [-0.4360;  0.5401]        3.8
## Filled: Song & Lindquist  0.0395 [-0.4048;  0.4837]        4.0
## Filled: Frogeli et al.    0.0220 [-0.3621;  0.4062]        4.2
## Filled: Call et al.      -0.0571 [-0.5683;  0.4541]        3.8
## Filled: Gallego et al.   -0.0729 [-0.5132;  0.3675]        4.0
## Filled: Kang et al.      -0.6230 [-1.2839;  0.0379]        3.3
## Filled: Shapiro et al.   -0.8277 [-1.4456; -0.2098]        3.4
## Filled: DanitzOrsillo    -1.1391 [-1.8164; -0.4618]        3.3
##
## Number of studies combined: k = 26 (with 8 added studies)
##
##                         SMD           95%-CI    t p-value
## Random effects model 0.3428 [0.1015; 0.5841] 2.93  0.0072
##
## Quantifying heterogeneity:
##  tau^2 = 0.2557 [0.1456; 0.6642]; tau = 0.5056 [0.3816; 0.8150];
##  I^2 = 76.2% [65.4%; 83.7%]; H = 2.05 [1.70; 2.47]
##
## [...]
## Details on meta-analytical method:
## - Inverse variance method
## - Restricted maximum-likelihood estimator for tau^2
## - Q-profile method for confidence interval of tau^2 and tau
## - Hartung-Knapp adjustment for random effects model
## - Trim-and-fill method to adjust for funnel plot asymmetry

可以看到，trim-and-fill 一共补入了 8 项研究。被修剪并填补的研究不仅包括我们识别出的离群值，还包括几项效应相对较高的小研究。填补出来的效应量都很低，其中一些甚至是高度负向的。输出还给出了校正后的综合效应估计：\(g = 0.34\)。它仍然显著，但比最初 m.gen 得到的 \(g = 0.58\) 小了很多。

再把它与去除离群值后的结果进行比较：

summary(tf.no.out)

## Review:     Third Wave Psychotherapies
## [...]
##
## Number of studies combined: k = 22 (with 6 added studies)
##
##                         SMD           95%-CI    t p-value
## Random effects model 0.3391 [0.1904; 0.4878] 4.74  0.0001
##
## Quantifying heterogeneity:
##  tau^2 = 0.0421 [0.0116; 0.2181]; tau = 0.2053 [0.1079; 0.4671];
##  I^2 = 50.5% [19.1%; 69.7%]; H = 1.42 [1.11; 1.82]
## [...]

在这里，结果几乎完全一样，仍然是 \(g = 0.34\)。总体而言，trim-and-fill 方法提示：我们 Meta 分析中最初得到的 \(g = 0.58\) 很可能因小样本研究效应而被高估，真实效应可能明显更小。这种高估很可能与发表偏倚有关，但也无法排除其他解释，因此 trim-and-fill 的校正结果本身也未必一定有效。

此外，我们还可以绘制包含补入研究的漏斗图。只需把 trimfill 的输出交给 meta::funnel 即可。下面代码为两次 trim-and-fill 分析（含离群值与去除离群值）分别绘制轮廓增强漏斗图，并用 par 将两幅图并排显示：

# 定义轮廓和颜色
contour <- c(0.9, 0.95, 0.99)
col.contour <- c("gray75", "gray85", "gray95")
ld <- c("p < 0.1", "p < 0.05", "p < 0.01")

# 用 par 设置一行两图
par(mfrow=c(1,2))

# 轮廓增强漏斗图（完整数据）
meta::funnel(tf,
             xlim = c(-1.5, 2), contour = contour,
             col.contour = col.contour)
legend(x = 1.1, y = 0.01,
       legend = ld, fill = col.contour)
title("Funnel Plot (Trim & Fill Method)")

# 轮廓增强漏斗图（去除离群值）
meta::funnel(tf.no.out,
             xlim = c(-1.5, 2), contour = contour,
             col.contour = col.contour)
legend(x = 1.1, y = 0.01,
       legend = ld, fill = col.contour)
title("Funnel Plot (Trim & Fill Method) - Outliers Removed")

在这些漏斗图中，被补入的研究以空心圆表示。

9.2.1.5 PET-PEESE

Duval 与 Tweedie 的 trim-and-fill 方法历史较长，而且可以说是最常见的小样本研究效应校正方法之一。不过正如前面所提到的，它远非完美，也不是估计偏倚校正综合效应的唯一办法。近年来，一个名为 PET-PEESE 的方法越来越受欢迎，特别是在经常使用 SMD 作为结局指标的研究领域，例如心理学或教育学。和前面的方法一样，PET-PEESE 针对的仍是小样本研究效应，而后者通常被看作发表偏倚的潜在线索。

PET-PEESE 实际上是两个方法的组合：precision-effect test（PET）与 precision-effect estimate with standard error（PEESE）。先看 PET。它基于如下简单回归模型，将研究效应量回归于其标准误：

和 Peters 检验一样，这里使用加权回归。研究权重 \(w_k\) 由方差的倒数给出，这与常规固定效应 Meta 分析中的权重相同：

需要指出的是，PET 使用的回归模型与 Egger 检验在本质上是等价的。主要区别在于：在 PET 中，\(\beta_1\) 量化漏斗图不对称，而在 Egger 检验中，不对称是通过截距体现的。

不过，PET 真正关心的并不是由 \(\beta_1\) 反映的不对称程度，而是截距 \(\beta_0\)。因为在上式中，截距代表所谓的极限效应（limit effect），即标准误为零时的预期效应量。这可以理解为完全没有抽样误差时的观察效应量。在其他条件相同的前提下，不受抽样误差影响的效应量本身就应最接近真实总体效应。

PET 的思想是：通过把标准误作为预测变量纳入模型，来控制小研究对结果的影响。理论上，这样得到的截距 \(\beta_0\) 就代表了校正所有小样本研究效应后，Meta 分析中的真实效应：

PEESE 的公式与之非常相似，唯一的区别是它把标准误平方（即效应量方差 \(s_k^2\)）作为预测变量：

研究权重公式 \(w_k\) 保持不变。之所以要对标准误平方，是因为小研究尤其容易报告被严重高估的效应，而统计功效较高研究中的这种问题则被认为轻得多。

当真实效应 \(\beta_0\) 等于零时，PET 往往表现更好；而当真实效应不等于零时，PEESE 的表现更优。因此，Stanley 和 Doucouliagos提出把两者结合起来，以平衡各自优缺点，这就是 PET-PEESE。它用 PET 或 PEESE 中某一个模型的截距 \(\beta_0\) 作为校正后真实效应的估计。

到底选 PET 还是 PEESE，取决于 PET 截距的大小。如果 PET 的截距 \(\beta_{0_{\text{PET}}}\) 在单侧检验中显著大于零（\(\alpha = 0.05\)），就采用 PEESE 的截距作为真实效应估计；如果 PET 截距并未显著大于零，就保留 PET 估计。

在 R 中，大多数回归模型默认使用双侧检验，因此如果要模拟单侧 \(\alpha = 0.05\) 的规则，在效应方向为正向有利时，只要 PET 截距的估计值大于零且其 \(p < 0.1\)，我们就把它视为“显著”，从而转用 PEESE。（如果负效应代表有利结果，则这一规则的方向应反过来。）

于是，PET-PEESE 的决策规则可以写成：

这条 if-else 规则初看确实有些绕，但一个实际例子会更清楚。下面我们基于 m.gen 看看 PET-PEESE 会给出怎样的真实效应估计。

目前 {meta} 中并没有直接实现 PET-PEESE 的现成函数，因此我们用 lm 自己来写。拟合 PET 和 PEESE 模型之前，先准备所需变量，并把数据框命名为 dat.petpeese。首先当然是标准化均值差。无论最初 Meta 分析是用 metacont 还是 metagen 完成的，每项研究的 SMD 都会保存在 Meta 分析对象的 TE 中。

# 构建数据集，先放入效应量
dat.petpeese <- data.frame(TE = m.gen$TE)

接下来需要效应量的标准误。对于 PET-PEESE，通常也建议使用 Pustejovsky 与 Rodgers提出的修正标准误；James Pustejovsky 还曾在其博客中推荐这一做法，并把这种替代方案称作 “SPET-SPEESE”： https://www.jepusto.com/pet-peese-performance/。

因此，我们使用前面相同的修正公式来计算 seTE_c，以避免标准误与效应量本身相关。同时再加入一个变量 seTE_c2，即标准误平方，因为 PEESE 需要它作为预测变量。

# 实验组 (n1) 与对照组 (n2) 样本量
n1 <- c(62, 72, 44, 135, 103, 71, 69, 68, 95,
        43, 79, 61, 62, 60, 43, 42, 64, 63)

n2 <- c(51, 78, 41, 115, 100, 79, 62, 72, 80,
        44, 72, 67, 59, 54, 41, 51, 66, 55)

# 计算修正标准误
dat.petpeese$seTE_c <- sqrt((n1+n2)/(n1*n2))

# 添加标准误平方（即方差）
dat.petpeese$seTE_c2 <- dat.petpeese$seTE_c^2

最后，还要计算每项研究的逆方差权重 w_k。这里同样用修正标准误平方来近似方差。

dat.petpeese$w_k <- 1/dat.petpeese$seTE_c^2

现在，dat.petpeese 已包含拟合 PET 和 PEESE 所需的全部变量。下面拟合两个加权线性回归模型，并直接用 summary 打印系数：

# PET
pet <- lm(TE ~ seTE_c, weights = w_k, data = dat.petpeese)
summary(pet)$coefficients

# PEESE
peese <- lm(TE ~ seTE_c2, weights = w_k, data = dat.petpeese)
summary(peese)$coefficients

要决定使用 PET 还是 PEESE，首先要看 PET 的结果。这里我们发现，PET 的极限估计为 \(g = -1.35\)。虽然该效应显著（\(p < 0.10\)），但它明显小于零，因此按规则应采用 PET 估计本身。

然而，\(g = -1.35\) 作为偏倚校正后的真实效应，显然并不可信。它意味着这种干预在真实情况下不仅没有好处，反而对结局产生了高度负面影响，似乎相当“有害”。这看起来极不合理。某种“真诚实施”的干预没有效果是可能的，但真正有害的干预并不常见。

事实上，这正是 PET-PEESE 的常见局限之一：它有时会对数据偏倚进行过度校正。我们的例子似乎就是如此：尽管所有观察效应量都是正向的，校正后的效应量却变成了高度负值。再看 PEESE 的结果，也会发现其估计同样为负，虽然绝对值略小（\(g = -0.44\)）。

当出现这种情况时，最好不要把截距解释为真实效应的点估计。更稳妥的做法是说：PET-PEESE 提示，在校正小样本研究效应后，这种干预类型没有效应。也就是说，我们更适合把 \(\hat\theta_{\mathrm{PET-PEESE}}\) 视为 0，而不是直接接受模型给出的负效应估计。

PET-PEESE 的局限

PET-PEESE 不仅可能系统性地对综合效应量过度校正，有时甚至会在根本不存在发表偏倚时高估真实效应。总体来看，当纳入研究数量较少（\(K < 20\)）、样本量偏小且异质性很高（例如 \(I^2 > 80\%\)）时，PET-PEESE 的表现较差。

遗憾的是，这样的 Meta 分析在现实中并不少见。因此，我们不建议把 PET-PEESE 作为唯一的小样本研究效应校正方法。不过，了解这一方法及其实现方式仍然很重要，因为它在某些研究领域正变得越来越常见。

用 rma.uni 代替 lm 实现 PET-PEESE

在上面的示例中，我们使用 lm 和研究权重实现了 PET-PEESE。但这一做法并非毫无争议。

通过 lm 实现的加权回归模型，与 rma.uni 等 Meta 回归函数之间存在一个细微却关键的差异：lm 采用的是乘性误差模型，而 Meta 分析函数通常采用加性误差模型。这里不展开技术细节，感兴趣的话可以阅读 Wolfgang Viechtbauer 关于该问题的优秀说明文档：https://www.metafor-project.org/doku.php/tips:rma_vs_lm_lme_lmer。

核心结论是：由于 lm 假定抽样误差方差存在一个比例常数，因此它并不是处理 Meta 分析数据的完美工具。至少作为敏感性分析，使用 rma.uni 而非 lm 来实现 PET-PEESE 是有道理的。实际操作上，就是在 mods 中把 \(SE_{\theta_k}^{(2)}\) 作为调节变量，例如 PET 可写成：

rma.uni(TE, seTE^2, mods = ~seTE, data = dat, method = "FE")

9.2.1.6 Rücker 的极限 Meta 分析方法

另一种计算校正后效应量估计的方法，是执行 极限 Meta 分析（limit meta-analysis），这是 Rücker 及其同事提出的方法。与 PET-PEESE 相比，这种方法更复杂，计算也更精细。因此，这里我们重点把握它的总体思想，而把繁重计算交给 R 完成。

Rücker 方法的核心思想是：构建一个显式考虑小样本研究效应偏倚的 Meta 分析模型。回顾一下，随机效应 Meta 分析可写为：

其中，\(\hat\theta_k\) 是第 \(k\) 项研究的观察效应量，\(\mu\) 是真实总体效应量，\(\epsilon_k\) 是抽样误差，\(\zeta_k\) 则量化研究间异质性带来的偏离。

在极限 Meta 分析中，我们扩展了这个模型。因为一旦存在小样本研究效应，我们就不能再认为效应量与标准误是彼此独立的。理由在于：如果发表偏倚特别影响小研究，那么小研究的效应量通常就会比大研究更大。Rücker 方法通过引入一个新的项 \(\theta_{\text{Bias}}\) 来把这种偏倚加入模型，并假设它与 \(\epsilon_k\) 以及 \(\zeta_k\) 相互作用，而且会随着 \(\epsilon_k\) 增大而增强：

这里要特别注意：式中的 \(\mu_*\) 已不再代表通常意义上的总体真实效应，而是一个在“标准”随机效应 Meta 分析中没有直接对应物的全局均值（除非 \(\theta_{\text{Bias}} = 0\)）。

接下来的思路与 PET-PEESE 类似。我们设想研究的效应量估计越来越精确，也就是其抽样误差 \(\epsilon_k\) 逐渐趋近于零。此时，\(\epsilon_k\) 会从模型中“消失”：

这里的 \(\mathrm{E}(\hat\theta_k)\) 表示当 \(\epsilon_k\) 趋于零时，\(\hat\theta_k\) 的期望值。这就是“极限 Meta 分析”的基本形式：它给出了在去除高标准误研究的扭曲影响之后，效应量的调整后估计。由于 \(\zeta_k\) 通常可由研究间异质性方差 \(\tau^2\)（或其平方根 \(\tau\)）来表示，我们可以进一步写成：

其中，\(\hat\theta_{*}\) 表示在校正小样本研究效应后，合并效应量的估计。Rücker 方法使用极大似然来估计该式中的参数，包括校正后的真实效应 \(\hat\theta_*\)。此外，它还可以为每项研究 \(k\) 给出一个“收缩后”的研究层面效应估计：

这里的 \(SE^2_k\) 是研究 \(k\) 标准误的平方（即观察方差），\(\hat\theta_k\) 是原始观察效应量。（这一公式可以从式 9.11 的更完整形式推导出来，技术性说明可见 Rücker 等人的原文公式 2.4 到 2.6。）

与 PET-PEESE 相比，Rücker 极限 Meta 分析的一大优点，是它把异质性方差 \(\tau^2\) 显式纳入了模型。另一个更实际的优点是，它可以在 R 中直接实现，所用函数是 limitmeta，该函数属于 {metasens} 包。

由于 {metasens} 与 {meta} 出自同一研究团队，二者通常配合得相当顺畅。例如，要对 m.gen 进行极限 Meta 分析，只需把它作为 limitmeta 的第一个参数：

# 先安装并加载 'metasens'
library(metasens)

# 运行极限 Meta 分析
limitmeta(m.gen)

## Results for individual studies
## (left: original data; right: shrunken estimates)
##
##                           SMD        95%-CI      SMD        95%-CI
## Call et al.              0.70 [ 0.19; 1.22]    -0.05 [-0.56; 0.45]
## Cavanagh et al.          0.35 [-0.03; 0.73]    -0.09 [-0.48; 0.28]
## DanitzOrsillo            1.79 [ 1.11; 2.46]     0.34 [-0.33; 1.01]
## de Vibe et al.           0.18 [-0.04; 0.41]     0.00 [-0.22; 0.23]
## Frazier et al.           0.42 [ 0.13; 0.70]     0.13 [-0.14; 0.42]
## Frogeli et al.           0.63 [ 0.24; 1.01]     0.13 [-0.25; 0.51]
## Gallego et al.           0.72 [ 0.28; 1.16]     0.09 [-0.34; 0.53]
## Hazlett-Stevens & Oren   0.52 [ 0.11; 0.94]    -0.00 [-0.41; 0.40]
## Hintz et al.             0.28 [-0.04; 0.61]    -0.05 [-0.38; 0.26]
## Kang et al.              1.27 [ 0.61; 1.93]     0.04 [-0.61; 0.70]
## Kuhlmann et al.          0.10 [-0.27; 0.48]    -0.29 [-0.67; 0.08]
## Lever Taylor et al.      0.38 [-0.06; 0.84]    -0.18 [-0.64; 0.26]
## Phang et al.             0.54 [ 0.06; 1.01]    -0.11 [-0.59; 0.36]
## Rasanen et al.           0.42 [-0.07; 0.93]    -0.25 [-0.75; 0.25]
## Ratanasiripong           0.51 [-0.17; 1.20]    -0.48 [-1.17; 0.19]
## Shapiro et al.           1.47 [ 0.86; 2.09]     0.26 [-0.34; 0.88]
## Song & Lindquist         0.61 [ 0.16; 1.05]     0.00 [-0.44; 0.44]
## Warnecke et al.          0.60 [ 0.11; 1.08]    -0.09 [-0.57; 0.39]
##
## Result of limit meta-analysis:
##
##  Random effects model     SMD            95%-CI     z     pval
##     Adjusted estimate -0.0345 [-0.3630; 0.2940] -0.21   0.8367
##   Unadjusted estimate  0.5771 [ 0.3782; 0.7760] -0.21 < 0.0001
## [...]

输出首先展示了每项研究的原始效应估计（左）与收缩后的效应估计（右）。可以看到，调整后的效应量明显小于原始观察值，有些甚至变成了负值。输出第二部分给出了调整后的合并效应估计，这里为 \(g = -0.03\)，也就是说，在校正小样本研究效应后，总体效应基本接近于零。

如果这些小样本研究效应确实是由发表偏倚导致的，这个结果就相当令人沮丧了。它意味着我们最初的发现可能完全是虚假的，是选择性发表掩盖了这种治疗其实并无效果这一事实。不过，我们仍然很难证明：小样本研究效应在我们的数据中是否仅仅由发表偏倚造成。

对于极限 Meta 分析，也可以绘制相应的漏斗图：只需把 limitmeta 的结果交给 funnel.limitmeta。这种图看起来和 meta::funnel 生成的漏斗图几乎一样，只是多了一条灰色曲线。这条曲线一方面表示当纵轴上的标准误为零时，校正后的平均效应；另一方面也象征着随着标准误增大，小样本研究效应所带来的偏倚如何逐渐增强。

如果把 shrunken = TRUE，还可以在图中显示各研究收缩后的研究层面效应估计。代码如下：

# 创建 limitmeta 对象
lmeta <- limitmeta(m.gen)

# 带曲线的漏斗图
funnel.limitmeta(lmeta, xlim = c(-0.5, 2))

# 带曲线和收缩后研究效应的漏斗图
funnel.limitmeta(lmeta, xlim = c(-0.5, 2), shrunken = TRUE)

需要注意的是，limitmeta 并不只适用于 SMD；任何 {meta} 包生成的 Meta 分析对象都可以使用。为说明这一点，我们再看看前文基于风险比的 m.bin：

limitmeta(m.bin)

## Result of limit meta-analysis:
##
##  Random effects model     RR           95%-CI    z     pval
##     Adjusted estimate 2.2604 [1.8066; 2.8282] 7.13 < 0.0001
##   Unadjusted estimate 2.0217 [1.5786; 2.5892] 7.13 < 0.0001

在这个例子中，原始估计与调整后估计几乎一致。考虑到前面 Peters 检验已经提示该 Meta 分析中小样本研究效应不明显，这一点并不令人意外。

9.2.2 P-curve

前面我们介绍的各种方法，都是通过考察小样本研究效应来评估发表偏倚风险。尽管具体实现不同，它们都建立在同一个思想之上：选择性报告会让研究效应量依赖于样本量。也就是说，标准误更高、精确度更低的小研究，平均上会比大研究表现出更高的效应量，因为只有那些效应特别大的小研究才会被发表，其他则停留在文件抽屉中。

不过，这种“理论”虽然直观，却也可以说多少有些没抓住重点。小样本研究效应方法假设，发表偏倚是通过效应量驱动的；而更现实的看法可能是，发表偏倚实际上是通过 \(p\) 值 运作的。现实中，研究结果通常只有在 \(p < 0.05\) 时才被认为值得发表。

正如前面提到的，科研是由人类完成的，因此同样受到金钱、声望和激励结构的影响。“significant \(p\), or no PhD” 这句恶名昭彰的话，非常准确地概括了这一点。研究者往往承受巨大外部压力，被要求“产出”小于 0.05 的 \(p\) 值。他们知道，能否过这一阈值，往往决定研究能否发表，也决定它是否被视为“成功”。这或许解释了为什么阴性或不显著结果越来越少出现在已发表文献中。

可以说，小样本研究效应方法只是间接捕捉了发表偏倚的机制。的确，选择性报告会导致小研究更容易显示较大效应；但那只是因为更大的效应提高了获得 \(p < 0.05\) 的检验统计量的概率。对小样本研究效应方法而言，\(p = 0.049\) 与 \(p = 0.051\) 的研究几乎没有差别；但在现实中，这样微小的差别可能对研究者来说是天壤之别。

接下来我们介绍一种聚焦于 \(p\) 值的发表偏倚方法：p-curve。这种方法的特别之处在于，它只关注显著效应及其 \(p\) 值分布。它既可以判断 Meta 分析数据背后是否存在真实效应，也可以估计其大小。更重要的是，它还明确控制了诸如 p-hacking 之类的可疑研究实践，而这一点是小样本研究效应方法通常做不到的。

P-curve 是一种相对较新的方法。它是在近年来社会科学所经历的“可重复性危机”背景下发展起来的。许多看似牢固的研究发现事实上并不稳健，无法被系统重复验证。这促使人们重新关注检测发表偏倚的方法，因为选择性报告很可能正是重复失败的重要解释之一。若不充分控制选择性发表，Meta 分析很可能只是复制了已发表文献中原本就存在的偏倚。

P-curve 的提出也部分是为了回应传统发表偏倚方法，尤其是 Duval 与 Tweedie trim-and-fill 方法的不足。Simonsohn 及其同事发现，trim-and-fill 往往只会带来很小的下调，而且经常无法识别这样一种情况：被分析的数据背后其实根本没有真实效应。

顾名思义，p-curve 是基于一条 \(p\) 值曲线。它有点像直方图，显示 Meta 分析中那些 \(p < 0.05\)、\(p < 0.04\)、\(p < 0.03\) 等研究分别有多少。P-curve 方法认为，这个显著 \(p\) 值直方图的形状取决于研究样本量，而且更重要地，取决于数据背后的真实效应量。

为了说明这一点，作者模拟了 9 组 Meta 分析结果。为了让模式足够清晰，每组想象中的 Meta 分析都包含极其庞大的 \(K = 10^5\) 项研究。在这 9 组模拟中，单项研究样本量从 \(n = 20\) 到 \(n = 100\) 不等，真实效应量则从 \(\theta = 0\) 到 0.5 不等。模拟假定每个 Meta 分析中的研究共享同一个真实效应，也就是服从固定效应模型。随后，只保留这些模拟中达到显著的效应量对应的 \(p\) 值，并绘制其直方图，结果如下图所示。

图 9.1：不同研究样本量和真实效应下的 P-curve。

第一行展示的是在不存在真实效应时，显著 \(p\) 值的分布。可以看到，无论单项研究样本量多大，模式都是相同的：\(p = 0.04\) 这种“勉强显著”的值与 \(p = 0.01\) 这种高度显著的值，看起来几乎一样常见。这种平坦的 p-curve 会在数据中没有真实效应、即零假设 \(\theta = 0\) 成立时出现。

当零假设为真时，\(p\) 值被假定服从均匀分布：每一个 \(p\) 值出现的概率都相同。即使效应量真的为零，偶然得到显著结果仍然可能发生，这就是假阳性或 \(\alpha\) 错误。但这种情况有多大概率，我们是知道的：既然 \(p\) 值在效应为零时服从均匀分布，那么所有 \(p\) 值中大约有 5% 会小于 0.05，这正是我们在假设检验中常用来拒绝零假设的显著性水平 \(\alpha = 0.05\)。

第二行和第三行的 p-curve 看起来则完全不同。此时零假设为假，数据中存在真实效应，因此显著 \(p\) 值会形成右偏分布。如果数据确实包含真实效应，那么高度显著的结果（例如 \(p = 0.01\)）就会比勉强显著的结果（例如 \(p = 0.049\)）更常见。而随着真实效应增大、样本量增大，这种右偏还会变得越来越明显。

不过，即便 Meta 分析中的研究统计功效严重不足，比如每项研究只有 \(n = 20\) 人却试图探测一个仅有 \(\theta = 0.2\) 的小效应，我们仍然可以看到右偏的 p-curve。这说明 p-curve 对真实效应量的变化相当敏感。只要确实存在真实效应，我们通常只需查看显著 \(p\) 值的分布，就有机会识别出来。

现在，想象一下如果研究者对结果进行了 p-hacking，p-curve 会是什么样子。通常，分析者只有在结果不显著但接近显著时才会开始 p-hacking；他们会调整分析细节，直到 \(p < 0.05\) 为止。既然一旦达到阈值通常就足以发表，因此在这之后往往不会继续折腾。可想而知，广泛的 p-hacking 会导致左偏 p-curve：略低于 0.05 的 \(p\) 值被过度代表，而高度显著的结果则相对不足。

总之，p-curve 可以被视为一种诊断工具，用来评估发表偏倚与 p-hacking 的可能性。下面我们将讨论 p-curve analysis，也就是基于经验 p-curve 的一组统计检验。重要的是，这些检验并不是直接针对发表偏倚本身，而是着眼于数据中是否存在证据价值（evidential value）。这其实正是我们在 Meta 分析中最关心的：我们想知道，所估计的效应是不是虚假的，是不是只是选择性报告造成的幻象。P-curve 正是为此而设计的。它能帮助我们判断，研究发现背后是否真的存在现实中的效应，还是说它不过是“一场声嘶力竭却毫无所指的喧哗”。

9.2.2.1 证据价值检验

为了评估证据价值是否存在，p-curve 使用两类检验：右偏检验（test for right-skewness）和33% 功效检验（也可视为对 p-curve 平坦性的检验）。先从右偏检验开始。正如前面所说，p-curve 的右偏程度取决于研究样本量及其背后的真实效应。因此，如果一种检验能够确认 Meta 分析对应的 p-curve 显著右偏，那就非常有帮助。若显著右偏，则说明这些显著 \(p\) 值很可能是由真实效应驱动的。

9.2.2.1.1 右偏检验

为检验右偏，p-curve 首先使用二项检验。二项检验适用于服从二项分布的数据。所谓二项分布，指的是数据可划分为两类（例如成功/失败、正面/反面、是/否），其中 \(p\) 表示其中一类结果的概率，\(q = 1-p\) 表示另一类结果的概率。

在这里，我们首先把 p-curve 分为两个区间：一类是 \(p < 0.025\) 的值；另一类是显著但大于 0.025 的值。由于 p-curve 只包含 0 到 0.05 之间的显著 \(p\) 值，这相当于把横轴一分为二。如果 p-curve 的确右偏，我们就应当看到这两类值的数量不相等，因为获得 \(p < 0.025\) 的概率将明显大于获得 \(p > 0.025\) 的概率。

设想我们的 p-curve 包含 8 个值，其中 7 个小于 0.025。此时可以使用 R 中的 binom.test 来检验：在零假设“高 \(p\) 值与低 \(p\) 值出现概率相同”成立时，得到这种数据有多大概率。由于我们假定较小的 \(p\) 值更常见，因此可以把 alternative 设为 "greater" 来执行单侧检验。

k <- 7   # p<0.025 的研究数
n <- 8   # 显著研究总数
p <- 0.5 # 零假设下的概率

binom.test(k, n, p, alternative = "greater")$p.value

结果会显示该二项检验显著（\(p < 0.05\)）。这意味着，在这个例子中，较小的 \(p\) 值显著多于较大的 \(p\) 值，因而 p-curve 呈右偏，也就提示数据背后存在真实效应。

二项检验的缺点在于：它需要把原本连续的 \(p\) 值人为二分。为避免信息损失，我们还需要一种不必把数据离散化的方法。

P-curve 通过为每个显著 \(p\) 值计算一个新的 \(pp\) 值来实现这一点。\(pp\) 值表示：当 p-curve 是平坦的（即不存在真实效应）时，得到一个至少和当前 \(p\) 值一样极端结果的概率。它给出了在“只考虑显著值”的条件下，该 \(p\) 值出现的概率。由于当 \(\theta = 0\) 时，\(p\) 值服从均匀分布，因此对于连续结局而言，\(pp\) 值实际上就是把显著 \(p\) 值线性映射到 \([0,1]\) 区间，例如：\(p = 0.023 \times 20 = 0.46 \rightarrow pp\)。

拿到每项显著研究 \(k\) 的 \(pp_k\) 之后，就可以用 Fisher 方法检验右偏。Fisher 方法是一种“古老”的 Meta 分析技术，由 R. A. Fisher 在 20 世纪早期提出。它可以合并多个研究中的 \(p\) 值，用于检验其中是否至少有一项测量到了真实效应。其做法是对每个 \(pp\) 值取对数，求和后再乘以 -2，所得统计量服从自由度为 \(2 \times K\) 的 \(\chi^2\) 分布（其中 \(K\) 是 \(pp\) 值总数）：

作为简单示例，假设 p-curve 包含 5 个显著 \(p\) 值：0.001、0.002、0.003、0.004 和 0.03。先把它们转换为 \(pp\) 值：

p <- c(0.001, 0.002, 0.003, 0.004, 0.03)
pp <- p*20

# 查看 pp 值
pp

然后按上式计算 \(\chi^2\)：

chi2 <- -2*sum(log(pp))
chi2

结果为 \(\chi^2 = 25.96\)。由于共有 5 项研究，因此自由度为 \(\text{d.f.} = 2 \times 5 = 10\)。在 R 中，可用 pchisq 检验在零假设下出现该 \(\chi^2\) 值的概率：

pchisq(26.96, df = 10, lower.tail = FALSE)

得到的 \(p\) 值约为 0.0026，这意味着零假设很不可能成立，因此应被拒绝。换言之，这些 \(p\) 值确实呈现右偏，可以被解释为数据背后存在证据价值。

9.2.2.1.2 平坦性检验

右偏检验能够帮助我们判断显著 \(p\) 值分布是否暗示真实总体效应的存在。但问题在于，这一检验本身依赖于数据的统计功效。因此，如果右偏检验不显著，并不能自动说明没有证据价值。原因可能有二：一是确实没有真实效应；二是 p-curve 中值的数量太少，导致即使数据实际上右偏，检验也没有足够功效显著。

因此，当右偏检验不显著时，我们还必须排除“统计功效不足”这一解释。右偏检验的零假设是“不存在证据价值”，换言之，我们试图通过证明经验 p-curve 不是平坦的，来拒绝这一零假设。

平坦性检验则把逻辑反过来。新的零假设并不是“完全没有效应”，而是“数据中存在一个很小的效应，因此 p-curve 应该是轻微右偏的”。在平坦性检验中，我们的目标就是证明：观察到的 p-curve 甚至连这种轻微右偏都没有。如果能做到这一点，就意味着即便是非常小的效应也可以被排除，数据中很可能根本不存在任何证据价值。

P-curve 通过 33% 功效检验来实现这一点。它的思想是：对每一项显著研究 \(k\)，在假定真实效应非常小、恰好只能以 33% 的功效被检测出来时，计算该研究对应的预期 \(pp\) 值。33% 这一阈值在某种程度上带有任意性，是 p-curve 发明者将其作为“几乎可以忽略的极小效应”的经验界限。

我们不展开 33% 功效下 \(pp\) 值的具体统计推导，但要知道：它涉及非中心分布，例如非中心 \(F\) 分布、非中心 \(t\) 分布和非中心 \(\chi^2\) 分布，具体取决于结局类型。（在本章稍后介绍的 R 版 p-curve 实现中，效应量会先通过 \(z = \frac{\hat\theta}{SE_{\hat\theta}}\) 转成 \(z\) 分数；在此基础上，33% 功效检验所需的预期 \(pp\) 值是通过自由度为 1 的非中心 \(\chi^2\) 分布计算的。）

总的来说，平坦性检验首先为每一个显著 \(p\) 值计算基于 33% 功效的 \(pp\) 值。如果这一假设与观察到的 \(p\) 值分布吻合，那么这些 33% 功效下的 \(pp\) 值将近似服从均匀分布。这样，我们就能用和右偏检验相同的方法来处理它们，只不过这一次输入的是 33% 功效下的 \(pp\) 值。

9.2.2.1.3 P-curve 结果的解释

到这里，我们已经介绍了若干用于分析经验 p-curve 的检验。即使其中有些统计概念一时难以完全掌握，也不用担心。理解 p-curve 的方法学确实需要一点时间，而下面的实际示例会进一步帮助你。最重要的是把握这些检验背后的思路，以及它们的结果该如何解读。

在解释 p-curve 时，我们实际上要同时面对 4 类结果：二项式右偏检验与平坦性检验的结果，以及基于 \(pp\) 值的右偏检验和平坦性检验结果。更复杂的是，p-curve 还引入了另外两项我们尚未详细讲解的检验：基于半个 p-curve 的右偏检验和平坦性检验。这两项检验与前述基于 \(pp\) 值的检验原理相同，只不过只使用较小的显著 \(p\) 值（即 \(p < 0.025\)）。

引入 half p-curve 检验，是为了防范一种更“激进”的 p-hacking。虽然相对不那么常见，但确实有可能出现这样的情况：研究者不断调参数，直到结果不只是勉强显著，而是变得“高度显著”。这样做会扭曲 p-curve 的形状，使得即使在没有真实效应时，曲线看起来也不一定左偏，甚至可能轻微右偏。基于 half p-curve 的检验能够在一定程度上控制这一问题，因为即便是很激进的 p-hacker，要想得到非常小的 \(p\) 值（例如 \(p < 0.01\)）也越来越困难，除非数据背后确实存在真实效应。由于 half p-curve 只包含小于 0.025 的值，因此这里不再进行二项检验。

在解释 p-curve 结果时，本质上要回答两个问题。

1. p-curve 是否提示存在证据价值？这主要看右偏检验。
2. 如果我们无法确认存在证据价值，那么证据价值是否缺失或不足？这主要看平坦性检验。

在实践中，可以使用如下判断规则：

• 存在证据价值：half p-curve 的右偏检验显著（\(p < 0.05\)），或者 full 与 half 两条曲线的右偏检验都达到 \(p < 0.1\)。
• 证据价值缺失或不足：full curve 的平坦性检验显著（\(p < 0.05\)），或者 half curve 的平坦性检验与二项检验都达到 \(p < 0.1\)。

如何理解 “no-no” 情形

每一次 p-curve 分析最终都会落入三种结果之一。若右偏检验显著，我们得出“存在证据价值”的结论。若右偏检验不显著、但平坦性检验显著，则表示证据价值缺失，或者效应非常非常小。

第三种也是最棘手的一种情况，可以称为 “no-no” 情形：既无法确认存在证据价值，也无法确认其缺失，也就是右偏检验和平坦性检验都不显著。

在解释上，这意味着：我们既不能确认真实效应存在，也不能排除一个相对较小但真实的效应。通常这类情形发生在 p-curve 中只包含很少研究的时候。它固然令人失望，但其真正传达的信息是：仅凭现有 p-curve，我们还不知道真实效应是否存在，需要更多证据来澄清。

9.2.2.2 在 R 中进行 P-curve 分析

到目前为止，我们已经了解了 p-curve 分析背后的不少理论，现在是时候在真实示例中使用它了。幸运的是，p-curve 的发明者 Simonsohn、Simmons 和 Nelson开发了一个可以自动完成上述检验并返回关键结果的应用，即 p-curve app，在线版本可见：http://p-curve.com/。

在 R 中，我们可以使用 pcurve 函数。这个函数模仿网页版应用的行为，并且专门为 {meta} 包生成的 Meta 分析对象设计。

pcurve 函数

pcurve 函数收录在 {dmetar} 包中。安装并加载 {dmetar} 后即可直接使用。如果你没有安装 {dmetar}，可以按以下方式手动使用：

1. 在线获取该函数源码：https://raw.githubusercontent.com/MathiasHarrer/dmetar/master/R/pcurve2.R
2. 将完整源码复制到 RStudio 控制台中运行，让 R “学会”这个函数。
3. 确保 {stringr} 和 {poibin} 已安装并加载。

使用 pcurve 非常简单，只需提供一个之前生成的 {meta} Meta 分析对象即可。在本例中，我们继续使用前面通过 metagen 得到的 m.gen。

不过，在正式运行分析之前，我们先去除前文已识别的两个离群研究（第 3 和第 16 项）。稍后会解释为什么这么做是合理的。

library(dmetar)
library(meta)

# 更新 m.gen，排除离群值
m.gen_update <- update(m.gen, subset = -c(3, 16))

# 运行 p-curve 分析
pcurve(m.gen_update)

## P-curve analysis
##  -----------------------
## - Total number of provided studies: k = 16
## - Total number of p<0.05 studies included into the
##   analysis: k = 9 (56.25%)
## - Total number of studies with p<0.025: k = 8 (50%)
##
## Results
##  -----------------------
##                     pBinomial  zFull pFull  zHalf pHalf
## Right-skewness test     0.020 -3.797 0.000 -2.743 0.003
## Flatness test           0.952  1.540 0.938  3.422 1.000
## Note: p-values of 0 or 1 correspond to p<0.001 and p>0.999,
## respectively.
##
## Power Estimate: 66% (31.1%-87.8%)
##
## Evidential value
##  -----------------------
## - Evidential value present: yes
## - Evidential value absent/inadequate: no

运行 pcurve 后，会得到两部分输出：一部分是 p-curve 分析结果，另一部分是观察到的 p-curve 图形。

输出告诉我们，这个 Meta 分析中有 \(k = 9\) 个显著效应被纳入 p-curve，其中大多数（\(k = 8\)）都高度显著，即 \(p < 0.025\)。在 Results 部分，可以看到三项右偏检验全部显著：二项检验（pBinomial; \(p = 0.02\)）、完整 p-curve 的 \(pp\) 值检验（pFull; \(p < 0.001\)），以及基于 half p-curve 的检验（pHalf; \(p = 0.003\)）。这里 full 与 half 检验对应的统计量是 \(z\) 分数，这是因为 pcurve 函数使用的是 Stouffer 方法而不是 Fisher 方法来合并结果；两者本质上非常接近。

按照前面的解释规则，这些结果说明我们的数据具有证据价值。下一行的三项平坦性检验均不显著，\(p\) 值从 0.938 到 1 不等，这同样合乎逻辑地说明：证据价值既不缺失，也不算不足。输出中 Evidential value 部分也给出了同样的解释。

函数生成的 p-curve 图包含三条线：实线代表基于数据的经验 p-curve；虚线表示在假定 33% 功效时的预期 \(p\) 值分布；点线表示在不存在效应时应有的均匀分布。可以明显看到，实线呈右偏形态，这正是数据中存在真实效应时我们期待看到的模式。

总体来看，这些结果表明：数据中确实存在证据价值，也就是说，背后存在一个真实的非零效应。我们仍然不能排除发表偏倚对 Meta 分析结果有影响；但至少根据 p-curve 的结果，我们可以说，这个综合效应并非完全虚假，不只是选择性报告制造出来的“海市蜃楼”。

有意思的是，这一结论与前面某些小样本研究效应方法并不完全一致。无论是 PET-PEESE，还是极限 Meta 分析，它们都把校正后的平均效应估计在大约零附近。细心的读者也许已经注意到：这些小样本研究效应方法在前面的示例中并未去除离群值，这可能是结果差异的一部分原因。作者也重新在去除离群值的条件下跑过 PET-PEESE 和极限 Meta 分析，结果与之前大体相同。

9.2.2.3 P-curve 的效应量估计

至此，我们已经看到经验 p-curve 的分析如何帮助判断 Meta 分析中是否存在证据价值。不过，即便这一步得到清楚结论，我们获得的信息仍然有限。知道数据背后存在真实效应当然很有帮助，但如果还能知道这个真实效应有多大，就更理想了。幸运的是，p-curve 也可以回答这个问题。只要知道各研究的样本量，就可以寻找一个最能拟合 p-curve 形状的真实效应量。

要理解这为什么可行，首先必须讨论非中心分布（non-central distribution）这一概念，以及它与效应量之间的关系。为说明非中心分布，我们从最常见的统计检验之一——独立样本 \(t\) 检验——开始。\(t\) 检验通常用于判断两组均值是否不同。

\(t\) 检验的零假设是两组均值 \(\mu_1\) 与 \(\mu_2\) 完全相同，因此其差值为零。当零假设成立时，\(t\) 统计量服从中心 \(t\) 分布。中心 \(t\) 分布与标准正态分布相似，并以零为中心。

如果零假设不成立，中心 \(t\) 分布显然就不再适合描述数据。若两组均值存在真实差异，\(t\) 统计量就不会围绕零分布，而会服从非中心 \(t\) 分布。这种分布通常不对称，且离散程度更大。最关键的是，它的中心会偏离零，而偏离的幅度由非中心参数 \(\delta\) 控制。\(\delta\) 越大，非中心分布的峰值就离零越远。

下图展示了这一点：左边是 \(\delta = 0\) 的中心 \(t\) 分布，右边是 \(\delta = 5\) 的非中心 \(t\) 分布。可以看到，右边曲线明显不对称，且峰值大约位于 5，而左侧的中心分布则对称地围绕零。

左边的中心分布表示零假设为真时预期出现的 \(t\) 值；右边则表示备择假设为真时预期出现的 \(t\) 值。

换一种说法：左边代表没有效应时的 \(t\) 分布，右边代表存在效应时的 \(t\) 分布。中心分布例如可以对应 SMD = 0，而右侧非中心分布则可代表效应量为 SMD = 1.3 时的预期 \(t\) 值分布。效应越大，非中心参数 \(\delta\) 就越大，非中心分布也就离零越远。

在常规统计教材中，非中心分布之所以不常见，是因为在一般假设检验里我们通常不需要它。多数情况下，我们的备择假设是非特异性的：例如做两样本 \(t\) 检验时，我们真正关心的只是零假设“组间没有差异”是否成立。一旦数据与零假设不符，我们就拒绝它，并认为“某种效应存在”。至于效应到底多大，并不是那一步的重点。

而具体的备择假设通常只在进行功效分析时才需要。规划实验时，我们通常希望样本量足够大，使得假阴性的概率不超过 20%。为了做到这一点，我们必须假定一个真实效应量，因为统计功效取决于样本量与真实效应的共同作用。

把这些拼在一起，我们会发现：非中心 \(t\) 分布的形状，也就是统计功效，实际上只取决于两件事：样本量与真实效应。同样地，p-curve 的形状也由样本量与真实效应共同决定。这一点非常重要：它意味着只要知道 p-curve 中各研究的样本量，我们就有可能反推出其真实效应量。

为了具体说明这一点，下面给出一个小例子。对独立样本两组 \(t\) 检验（假定两组方差相等）而言，\(t\) 值等于组间均值差 \(\text{MD}_{\text{between}}\) 除以其标准误 \(SE_{\text{MD}_{\text{between}}}\)：

代入组间均值差及其标准误公式后，可得：

可以看到，这里 \(t\) 的自由度等于两组总样本量 \(n_1+n_2\) 减 2。

假设某研究实验组有 \(n_1 = 30\) 人，对照组有 \(n_2 = 20\) 人。研究报告两组均值分别为 13 和 10，且两组标准差都为 5。我们可以这样计算 \(t\) 值：

# 计算均值差
md <- 13-10

# 计算均值差的标准误
n1 <- 30
n2 <- 20
s1 <- s2 <- 5
s_pooled <- sqrt((((n1-1)*s1^2) + ((n2-1)*s2^2))/
                   ((n1-1)+(n2-1)))

se <- s_pooled*sqrt((n1+n2)/(n1*n2))

# 计算 t 值（等价于方差相等时的两样本 t 检验）
md/se

计算结果为 \(t_{48} = 2.078\)。这个结果是否支持零假设“组间均值无差异、没有效应”呢？

可以用 pt 函数回答。它给出的，是在自由度 48 且零假设为真时，获得大于 2.078 的 \(t\) 值的概率，也就是单侧 \(t\) 检验的 \(p\) 值：

pt(2.078, df = 48, lower.tail = FALSE)

结果为 \(p = 0.02\)，说明检验显著。因此我们拒绝“实验组效应为零或负”的零假设，接受备择假设：实验组具有正向效应（假定更高分数代表更好结局）。

既然零假设所对应的中心 \(t\) 分布并不适合这组经验数据，那么更合适的应是某个非中心 \(t\) 分布。问题只是：到底是哪一个？此时我们还不知道哪个真实效应量、也就是哪个非中心参数 \(\delta\)，最能代表产生这个经验 \(t\) 值的总体。

作为第一次猜测，我们可以假定真实效应量是 \(\theta = 0.6\) 的标准化均值差。这意味着，我们之所以观察到 \(t = 2.078\)，是因为实验中存在一个中等偏大的真实效应。根据这一 \(\theta\) 值，可用下式计算非中心参数：

在本例中：

theta <- 0.6
delta <- theta/sqrt((n1+n2)/(n1*n2))

# 查看 delta
delta

为了看看 \(\delta = 2.078\) 的非中心 \(t\) 分布长什么样，我们做一个小模拟。用 rt 函数抽取一百万个自由度为 48 的随机 \(t\) 值，两次：一次设非中心参数为 0（即零效应），一次设为刚刚得到的 \(\delta\)（即真实效应 SMD = 0.6）。然后绘制直方图：

# “1 后面 6 个 0”在 R 中也可以写成 1e6
# 记得先加载 tidyverse（因为这里用到了管道）
rt(n = 1e6, df = 48, ncp = 0) %>%
  hist(breaks = 100,
       col = "gray50",
       xlim = c(-4,8),
       ylim = c(0, 40000),
       xlab = "t-value",
       main = NULL)

rt(n = 1e6, df = 48, ncp = delta) %>%
  hist(breaks = 100,
       col = "gray95",
       xlim = c(-4,8),
       ylim = c(0, 40000),
       add = TRUE)

图中左边是中心 \(t\) 分布（无效应），右边是非中心分布（\(\theta = 0.6\)）。由于样本量已经不算太小（\(N = 50\)），非中心分布看起来不像前面的图那样明显右偏，但仍可以看到它整体向右平移，峰值大约位于 \(\delta\)。

关键问题当然是：如果这一备择分布确实正确，也就是说真实效应就是 \(\theta = 0.6\)，那么获得大于 \(t_{48} = 2.078\) 的值的概率到底是多少？要回答这个问题，可以再次使用 pt，但这次加入 ncp = delta：

# 别忘了这里的 t=2.078
pt(2.078, df = 48, ncp = delta, lower.tail = FALSE)

结果表明，在这一备择分布下，得到比经验值更大的 \(t\) 值的概率大约为 50%。这意味着，大约一半的值会比我们观察到的更大，另一半更小。换言之，自由度 48 且真实效应为 0.6 的非中心 \(t\) 分布与我们的经验结果相当吻合，看起来总体真实效应确实很可能接近 SMD = 0.6。

这其实就是 p-curve 用于估计真实效应量的基本思想。对于 p-curve 中每一个显著 \(p\) 值，它都会计算：在给定某个候选效应量/非中心参数、给定研究自由度（可由样本量推出），并且只考虑显著结果（\(p < 0.05\)）的条件下，观察到一个至少如此极端的统计量的概率。于是，对研究 \(k\) 来说，其 \(pp\) 值可以写为：

因为研究的自由度通常已知，式中唯一未知的是 \(\delta\)，也就是对应的真实效应量 \(\theta\)。不过，由于当候选真实效应量“猜对了”时，得到的 \(pp\) 分布应当接近均匀分布，因此这个未知量是可以通过搜索找到的。就像当数据符合零假设时，\(p\) 值服从均匀分布一样；当结果符合正确的非中心分布时，\(pp\) 值也会服从均匀分布。

因此，我们只需尝试很多很多个候选效应量，把由其决定的 \(\delta\) 代入上式，然后考察所得 \(pp\) 分布与均匀分布相差多少。那个最接近均匀分布的候选效应量，就被看作真实效应的最佳估计。P-curve 用 Kolmogorov-Smirnov（KS）检验中的 \(D\) 距离度量这一偏离程度。

P-curve 的效应量估计同样已经实现在 pcurve 函数中。只需把 effect.estimation = TRUE，并提供每项研究的样本量 N。还可以用 dmin 和 dmax 控制候选效应量的搜索范围。下面我们让 pcurve 在 Cohen’s \(d = 0\) 到 1 之间搜索。注意，dmin 必须大于等于 0，也就是说 p-curve 能识别的最小值就是“无效应”。

# 添加实验组 (n1) 与对照组 (n2) 样本量
# 这里已删除第 3 和第 16 项研究
n1 <- c(62, 72, 135, 103, 71, 69, 68, 95,
        43, 79, 61, 62, 60, 43, 64, 63)

n2 <- c(51, 78, 115, 100, 79, 62, 72, 80,
        44, 72, 67, 59, 54, 41, 66, 55)

# 带效应量估计的 p-curve 分析
pcurve(m.gen_update,
       effect.estimation = TRUE,
       N = n1+n2,
       dmin = 0,
       dmax = 1)

## P-curve analysis
##  -----------------------
## [...]
##
## P-curve's estimate of the true effect size: d=0.389

现在，输出中多了两项内容：真实效应量估计，以及一张显示搜索结果的图。图中，不同候选效应量形成一条平滑的 V 形曲线，其顶点位于 \(d = 0.389\)。在这个位置，计算得到的 \(pp\) 分布与均匀分布之间的差异（纵轴上的 \(D\) 值）最小，因此它被视为最好的真实效应估计。

需要强调的是，只有当这张图呈现出类似这里这样良好的 V 形 时，p-curve 的效应量估计才更值得信赖。如果图形更为紊乱，则说明 p-curve 可能没有找到真正的最小值；如果图形像一条平滑下降的斜线，则可能表示搜索范围过窄，此时应当增大 dmax 再重新运行。

总体来看，p-curve 给出的 \(d = 0.389\) 稍低于去除离群值后 Meta 分析得到的综合效应（\(g = 0.45\)），但仍然足以支持前面的结论：这些研究确实包含证据价值。

如何报告 p-curve 结果

在报告 p-curve 分析时，至少建议包含三项右偏检验与三项平坦性检验的 \(p\) 值，以及对这些结果的解释。如果同时估计了真实效应量，也应一并报告。一个简洁的表格可以写成这样：

检验	pBinomial	zFull	pFull	zHalf	pHalf	证据价值存在	证据价值缺失/不足	\(\hat d\)
Right-Skewness Test	0.020	-3.80	<0.001	-2.74	0.003	yes	no	0.39
Flatness Test	0.952	1.54	0.938	3.42	>0.999	yes	no	0.39

P-curve 的开发者还强烈建议为每次分析建立一张披露表（disclosure table），说明结果是从原始文章的哪个部分提取的，以及其原始报告形式。相关示例和其他实践建议可见 Simonsohn 等人的论文。

P-curve 与研究间异质性

这里还需要解释一下，为什么我们在 p-curve 分析中排除了离群研究。包含离群值时，Meta 分析的异质性为 \(I^2 = 63\%\)，已经相当高。问题在于，已有研究发现：当研究间异质性很高时，p-curve 并不是一个稳健的真实效应量估计方法。

Van Aert 等人因此建议，最好只在异质性较小到中等时使用 p-curve，他们提出的经验阈值是 \(I^2 = 50\%\)。如果 Meta 分析的异质性高于此值，一种折中办法就是像本例这样，先去掉离群值再做 p-curve；更好的办法则是在合理的亚组内分别进行分析。

9.2.3 选择模型

本章最后介绍的一类发表偏倚方法，是所谓的选择模型（selection models）。尽管选择模型用于评估选择性发表影响的想法早已存在，但近年来，人们对其应用兴趣明显增加。

前面介绍的所有发表偏倚方法，本质上都基于某种“理论”，用来解释选择性发表为什么以及如何影响 Meta 分析结果。比如，小样本研究效应方法假定研究未被发表的风险与其样本量和效应量成比例；p-curve 则假定 0.05 的 \(p\) 值像一道“魔法阈值”，使得 \(p \ge 0.05\) 的结果通常比显著结果更容易从数据中消失。

而选择模型可以看作这些方法的一种广义化版本。它允许我们把任何一种我们认为可能塑造发表偏倚过程的机制，都显式建模出来。这使其非常灵活：无论我们的假设非常简单，还是高度复杂，选择模型都可以据此构建。

所有选择模型的共同思想，是先指定一种分布，用以预测某项研究根据其结果被发表（即被“选中”）的概率。通常这里的“结果”就是研究的 \(p\) 值，因此选择模型可以被理解为一个函数：对不同的 \(p\) 值返回不同的发表概率。定义好这样的选择函数后，就可以据此“剥离”假设中的选择性发表偏倚，并据此估计校正后的真实效应量。

当然，这个校正后的效应只有在我们所设定的选择模型本身大体正确时才有意义。必须始终记住：我们所定义的模型只是众多可能解释中的一种，即便它看起来与数据很匹配，也未必就是真实的发表机制。选择性发表到底如何塑造了数据，终究无法被完全知晓。但即便如此，选择模型仍然非常有价值，因为它们至少能帮助我们粗略评估：发表偏倚是否可能影响了数据，以及可能是怎样影响的。

本章中，我们将介绍两类基于阶梯函数（step function）的简单选择模型。因此，先从阶梯函数本身说起。

9.2.3.1 阶梯函数选择模型

要进行任何一种选择模型分析，都需要两部分：一个效应量模型，以及选择模型本身。可以把它们都理解为函数：输入某个值 \(x\)，返回该值出现的概率。

效应量模型由函数 \(f(x_k)\) 表示，它与随机效应模型本质相同：假定观察效应量 \(\hat\theta_k\) 围绕平均效应 \(\mu\) 呈正态分布，并因抽样误差与研究间异质性方差 \(\tau^2\) 而偏离。只要知道 \(\mu\)、\(\tau^2\)、研究标准误，并假定效应量服从正态分布，函数 \(f(x_k)\) 就能告诉我们：在不存在发表偏倚时，观察到某个效应量 \(x_k\) 的概率有多大。

但是一旦存在发表偏倚，这个效应量分布以及 \(f(x)\) 本身就不再是真实世界的恰当刻画。因为选择性发表导致某些研究被过度代表——通常是那些样本量小、效应量又出乎意料地高的研究。换言之，在不知不觉中，我们实际上给了这些研究更高的权重。因此，我们需要构造一个更“现实”的版本 \(f^*(x)\)，把这种发表机制也纳入进去。

这正是通过权重函数 \(w(p_k)\) 实现的。权重函数告诉我们，第 \(k\) 项研究根据其 \(p\) 值被选中的概率是多少。有了它之后，就能定义一个同时考虑效应量分布与发表机制的修正函数：

式中分母是 \(w(p_k)f(x_k)\) 的积分。这里的 \(w(p_k)\) 就是我们假设的选择模型。

虽然从技术上说 \(w(p_k)\) 可以具有任何形状，但实际中人们常常把它实现为一个阶梯函数。如果 \(w(p_k)\) 是阶梯函数，意思是说：落在同一 \(p\) 值区间内的所有研究，具有相同的被选择概率。每个区间的选择概率记作 \(\omega_i\)，不同区间之间可以不同。也就是说，我们把 0 到 1 的可能 \(p\) 值范围切成若干段，并给每一段赋予自己的选择概率 \(\omega_1,\omega_2,\dots,\omega_c\)。

这些区间由若干切点（cut-points，记作 \(a_i\)）定义。切点的数量及其具体数值，都由研究者自行设定。假如一个阶梯函数有 4 个区间，那么它可以写成：

也就是说，对任意一个 \(p\) 值，该函数都根据它落入的区间返回相应的选择概率 \(\omega\)。为了更具体地理解这一点，下面设定一个带有实际数值的选择模型。

例如，我们可以假设某个 Meta 分析中的发表偏倚机制由 3 个切点描述。首先令 \(a_1 = 0.025\)。这对应于单侧 \(p = 0.025\)、双侧 \(p = 0.05\)，也就是最常见的显著性阈值。因此可以合理假定：所有小于 0.025 的结果，其被选中概率 \(\omega_1\) 为 100%。毕竟，如果研究结果是正向的，就没什么理由把它扔进文件抽屉。接着令 \(a_2 = 0.05\)，并假定此区间的选择概率为 80%；它仍然很高，但比明确显著的结果低一些。再往后设定一个较大的区间，从 0.05 到 0.5，假定选择概率为 60%；最后，对 \(p \ge 0.5\) 的研究，赋予更低的选择概率 \(\omega_4 = 35\%\)。

这就得到如下所示的阶梯型选择模型：

图 9.2：基于阶梯函数的选择模型。

在以阶梯函数定义选择模型时，我们通常只需设定切点 \(a_i\)。这些是模型中固定的部分，而各区间的选择概率 \(\omega = \omega_1,\omega_2,\dots,\omega_c\) 则由数据估计。随后，基于式 9.20 所定义的模型可以通过极大似然拟合到数据上，同时估计 \(\omega\)，以及考虑了不同选择概率之后修正的 \(\mu\) 与 \(\tau^2\)。其中修正后的 \(\mu\) 就代表在假设的发表偏倚机制下，真实平均效应量的估计。

需要说明的是，前文把 \(\omega\) 表示为 0% 到 100% 的绝对概率，只是为了帮助理解。实际拟合时，\(\omega_i\) 估计的不是绝对发表概率，而是相对选择可能性。这意味着：第一个区间会被固定为参考值 1，其他区间的 \(\omega\) 则表示相对于参考区间的选择可能性。例如，如果第二个区间估计得到 \(\omega_2 = 0.5\)，就表示该区间的研究被选中的可能性只有参考区间的一半。

当然，只有当我们设定的选择模型本身足够合适时，修正后的平均效应 \(\mu\) 才有意义。一个粗略的判断线索，是选择模型参数的似然比检验（LRT）是否显著。该检验的零假设是不存在选择性发表，即所有区间的相对选择可能性都相同：\(\omega_1 = \omega_2 = \dots = \omega_c\)。不过需要注意，这一显著性检验被发现往往偏保守性不足，因此解释时也要谨慎。

理论上，选择模型中的切点数量 \(a_i\) 可以任意设定；但每增加一个切点，就多需要估计一个 \(\omega_i\)。如果 Meta 分析规模本就不大，就很容易出现某些区间里研究数太少甚至为空的问题，导致各区间的 \(\omega_i\) 难以稳定估计。因此，切点很多的复杂选择模型通常只适用于研究数非常多的 Meta 分析（例如 \(K \ge 100\)）。

遗憾的是，大多数 Meta 分析只包含较少研究。这意味着，现实中往往只能使用切点很少的简单选择模型。其中一种特别简洁的变体，就是下面要介绍的三参数选择模型。它的优点是，即使纳入研究数量不多（例如 \(K = 15\) 到 20），也能使用。

9.2.3.1.1 三参数选择模型

三参数模型是一种只有一个切点的选择模型。之所以叫“三参数”，是因为只需估计 3 个参数：真实效应 \(\mu\)、研究间异质性方差 \(\tau^2\)，以及第二个区间的相对选择可能性 \(\omega_2\)。（不需要估计 \(\omega_1\)，因为第一个区间被固定为参考组，选择可能性恒定为 1。）

在三参数选择模型中，唯一切点 \(a_1\) 设为 0.025，也就是单侧 \(p = 0.05\)。这样就把全部 \(p\) 值分成两个区间：统计学显著与不显著。于是，\(\omega_2\) 表示“不显著结果被选中发表”的相对概率。实际上，p-curve 可以被看作三参数选择模型的一种特例：它也用 \(p = 0.05\) 作为切点，但只关注显著结果，并假定 \(\tau^2 = 0\)，因此真正由数据估计的只剩下真实效应量 \(\mu\)。

在 R 中，{metafor} 包的 selmodel 函数可以拟合多种选择模型，包括三参数模型。需要说明的是，在本书写作时，selmodel 仅存在于 {metafor} 的开发版本中；当时需要通过 remotes::install_github("wviechtb/metafor") 安装。不过如今你阅读本章时，它很可能已被整合进正式版本，无需额外安装开发版。

由于 selmodel 只接受由 {metafor} 的 rma 函数生成的 Meta 分析对象，因此我们首先需要创建这样一个对象。这里的设置尽量与之前用 {meta} 包 metagen 所做的分析保持一致：

library(metafor)

# 新对象命名为 m.rma
m.rma <- rma(yi = TE,
             sei = seTE,
             data = ThirdWave,
             slab = Author,
             method = "REML",
             test = "knha")

有了 m.rma 后，就可以通过 selmodel 拟合一个三参数选择模型。把 type 设为 "stepfun" 表示使用阶梯函数，在 steps 中填入切点 \(a_1 = 0.025\)：

selmodel(m.rma,
         type = "stepfun",
         steps = 0.025)

## [...]
##
## Model Results:
##
## estimate      se    zval    pval   ci.lb   ci.ub
##   0.5893  0.1274  4.6260  <.0001  0.3396  0.8390  ***
##
## Test for Selection Model Parameters:
## LRT(df = 1) = 0.0337, p-val = 0.8544
##
## Selection Model Results:
##
##                      k  estimate      se    pval   ci.lb   ci.ub
## 0     < p <= 0.025  11    1.0000     ---     ---     ---     ---
## 0.025 < p <= 1       7    1.1500  0.8755  0.8639  0.0000  2.8660
##
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

在 Model Results 中可以看到，选择模型给出的真实平均效应估计为 \(g = 0.59\)（95% CI: 0.34–0.84）。有意思的是，这几乎与我们之前算出的综合效应量 \(g = 0.58\) 完全一样。

这意味着：在这种设定下，我们的 Meta 分析并没有显示出因为不显著结果较少被选中而受到明显扭曲。这一点也得到 Test for Selection Model Parameters 的支持：检验不显著（\(\chi^2_1 = 0.034\), \(p = 0.85\)），说明 \(\omega_1\) 与 \(\omega_2\) 并没有显著差异。

在 Selection Model Results 中，可以看到两个区间的相对选择可能性估计。这里 \(\omega_2 = 1.15\)，意味着第二个区间的研究相对选择可能性甚至略高于第一个区间。如果发表偏倚非常严重，我们本应看到相反的结果：不显著结果的相对选择可能性应明显低于显著结果。

作为敏感性分析，我们可以把切点 \(a_1\) 从 0.025 改成 0.05 再重新运行。设为 0.05 的含义是：我们假定双侧 \(p\) 值处于 0.05 到 0.10 之间的结果，与 \(p < 0.05\) 的结果一样“可发表”。这可能对应于这样一种情况：原始研究中“趋势显著”的结果仍足以发表，或者某些研究使用了单侧检验。下面看看切点改变后结果是否不同：

selmodel(m.rma,
         type = "stepfun",
         steps = 0.05)

## [...]
##
## Model Results:
##
## estimate      se    zval    pval   ci.lb   ci.ub
##   0.3661  0.1755  2.0863  0.0370  0.0222  0.7100  *
##
## Test for Selection Model Parameters:
## LRT(df = 1) = 3.9970, p-val = 0.0456
##
## Selection Model Results:
##
##                    k  estimate     se   pval   ci.lb   ci.ub
## 0    < p <= 0.05  15    1.0000    ---    ---     ---     ---
## 0.05 < p <= 1      3    0.1786 0.1665 <.0001  0.0000  0.5050  ***
##
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

这一次就出现了不同的模式。新的平均效应估计为 \(g = 0.37\)，比前面明显更小；同时，似然比检验显著，说明两个区间之间存在差异。可以看到，\(\omega_2 = 0.18\)，意味着 \(p > 0.1\)（双侧）结果的选择可能性远低于那些（边缘）显著结果。

这提示我们的综合效应可能确实受到了选择性报告的轻度扭曲，尤其是那些明显不显著的研究，很可能更容易被塞进文件抽屉。

9.2.3.1.2 固定权重选择模型

在前面的三参数模型中，只设置了一个切点 \(a_1\)，而选择可能性由模型自由估计。正如前面所说，三参数模型由于只用一个切点，因此即便在研究数不多的 Meta 分析中也能使用，因为它较少发生某个区间“没有研究可估”的问题。

不过，如果不再让模型从数据中估计各区间的选择概率 \(\omega_i\)，而是直接由我们指定固定值，那么就能在不增加样本需求的前提下，使用更复杂、切点更多的选择模型。换言之，我们可以直接为每个区间设定一个固定的 \(\omega_i\)，然后检验：在这个假定的选择机制下，\(\mu\) 的估计会变成什么。

当然，这种固定权重选择模型也有代价。因为一旦我们事先设定的 \(\omega_i\) 并不符合真实情况，得到的平均效应估计就不再可信。因此，固定权重模型更适合作为一种敏感性分析：看看如果发表偏倚过程确实像我们假设的那样，那么真实效应量会是什么样子。

Vevea 与 Woods 提供了一些多切点固定权重选择模型的典型示例。下图展示了两种代表“中等选择”与“严重选择”的阶梯函数模型：

作为敏感性分析，我们可以检验：如果把这些选择模型视为适用于当前 Meta 分析，\(\mu\) 的估计会如何变化。首先要把图中的所有切点，以及每个区间的选择可能性 \(\omega_i\) 写出来：

# 定义切点
a <- c(0.005, 0.01, 0.05, 0.10, 0.25, 0.35, 0.50,
       0.65, 0.75, 0.90, 0.95, 0.99, 0.995)

# 定义各区间选择可能性
# （中等选择 / 严重选择）
w.moderate <- c(1, 0.99, 0.95, 0.80, 0.75, 0.65, 0.60,
                0.55, 0.50, 0.50, 0.50, 0.50, 0.50, 0.50)
w.severe <- c(1, 0.99, 0.90, 0.75, 0.60, 0.50, 0.40, 0.35,
              0.30, 0.25, 0.10, 0.10, 0.10, 0.10)

随后，把这些参数传给 selmodel。新的切点放在 steps，固定的选择可能性放在 delta：

# 假定中等选择
selmodel(m.rma, type = "stepfun", steps = a, delta = w.moderate)

## [...]
##
## Model Results:
##
## estimate      se    zval    pval   ci.lb   ci.ub
##   0.5212  0.0935  5.5741  <.0001  0.3380  0.7045  ***
## [...]

# 假定严重选择
selmodel(m.rma, type = "stepfun", steps = a, delta = w.severe)

## [...]
## Model Results:
##
## estimate      se    zval    pval   ci.lb   ci.ub
##   0.4601  0.1211  3.8009  0.0001  0.2229  0.6974  ***
## [...]

可以看到，在假定“中等选择”的模型下，综合效应估计为 \(g = 0.52\)；在假定“严重选择”的模型下，估计略降为 \(g = 0.46\)。

这两组结果都表明：即便考虑选择性发表，我们观察到的效应总体仍然比较稳健。不过务必牢记，这些估计只有在——并且只在——我们所指定的选择模型确实近似真实发表过程时才成立。

其他选择模型函数

上面只讨论了以阶梯函数为基础的选择模型。但要注意，这并不是唯一可用的选择函数形式。selmodel 还包含若干基于连续分布的选择模型，例如 half-normal、logistic 或 negative-exponential 选择模型。只需改变 type 参数即可调用。

这些模型超出了本书范围，不过 selmodel 的帮助文档本身就是一个很好的起点。加载 {metafor} 后，可直接在 R 中运行 ?selmodel 查看。

9.3 我该使用哪种方法？

至此，我们对统计学发表偏倚方法的讨论就告一段落了。本章内容很长，读者很自然会问：为什么要讲这么多不同的方法？难道不能只选一种方法，评估一下发表偏倚风险，然后继续往下做吗？

简短的回答是：不能。

发表偏倚方法至今仍是一个非常活跃的研究领域，过去很多研究都评估过不同方法的表现。然而，迄今并没有一个明确无争议的“赢家”。相反，越来越多的证据表明：并不存在一种方法能在所有情境下持续优于其他方法。

不同的发表偏倚方法得出截然不同的结果，不仅是可能的，而且相当常见。本章中的示例就是一个很好的说明。尽管我们每次都用的是同一个数据集，但对“真实偏倚校正后效应”的估计却从几乎为零一直到 \(g = 0.59\) 不等。这说明：方法的选择本身就会深刻影响结果，从而影响我们最终的结论。有些方法提示，一旦控制了小样本研究效应，综合效应几乎完全消失；而另一些方法则基本支持我们最初的发现。

因此，我们建议：在评估发表偏倚时，最好始终使用多种方法。通常很难甚至根本无法提前知道哪种方法最适合当前数据，也无法完全确定它给出的结果是否可信。正如前面反复强调的，选择性报告究竟在多大程度上影响了结果，这一点永远不可能被完全知晓。但通过并行使用多种发表偏倚方法，我们至少可以得到一个类似于“可信真实效应范围”的东西。

这个范围的宽窄本身就可以帮助解释结果。比如，如果 PET-PEESE、p-curve 和三参数选择模型都给出与原始综合效应相近的估计，我们对结果稳健性的信心就会增强；反之，如果这些方法彼此明显冲突，就说明发表偏倚和小样本研究效应对结果的影响远不确定，原始综合效应的可信度也随之下降。

无论如何，发表偏倚方法的结果都应当谨慎解释。现实中已经有不少例子表明，发表偏倚分析的结果可能引发激烈争论，例如有关 “ego depletion” 文献的争议。必须记住，控制发表偏倚的最好办法始终是：认真检索未发表证据，并从根本上改善科研发表实践。作为 Meta 分析者，我们所能提出的各种统计学“证据”，充其量也只能算有限而间接的线索。

为了帮助读者更快判断不同方法适用于什么场景，下面给出一个简要的优缺点概览。

这张表同时包含统计和实践层面的考虑，但既不能视为穷尽，也不能视为定论。发表偏倚方法仍是持续发展的研究领域，随着证据积累，未来的判断很可能会发生变化。

最后还要指出：目前还没有哪一种方法在研究间异质性非常高时（例如 \(I^2 \approx 75\%\)）能够始终提供令人满意的结果。这意味着，对于异质性极高的 Meta 分析，最好尽量避免直接开展发表偏倚分析。去除离群值后重做分析，或在更同质的亚组内分别分析，都是现实中常见的折中办法，但它们并不能从根本上解决问题。

9.4 问答

检验一下你的理解

1. “发表偏倚”可以如何定义？为什么它在 Meta 分析中是个问题？
2. 还有哪些其他报告偏倚？请至少列举并解释三种。
3. 请说出两种可疑研究实践（QRP），并解释它们如何威胁 Meta 分析的有效性。
4. 请解释小样本研究效应方法背后的核心假设。
5. 如果发现数据中存在小样本研究效应，这是否自动意味着存在发表偏倚？
6. P-curve 估计的是 Meta 分析中所有研究的真实效应，还是仅估计所有显著效应的真实效应？
7. 哪一种发表偏倚方法表现最好？

这些问题的答案见本书末尾附录 A： https://bookdown.org/MathiasHarrer/Doing_Meta_Analysis_in_R/qanda.html#qanda9

9.5 小结

• 发表偏倚是指：某些研究系统性地缺失于已发表文献中，因此也缺失于我们的 Meta 分析中。严格地说，只有当研究被发表的概率取决于其结果时，才称为发表偏倚。不过，除此之外还存在一系列其他报告偏倚，它们同样会影响某项研究最终是否进入我们的 Meta 分析，例如引文偏倚、语言偏倚和结局报告偏倚。
• 即便研究已经发表，证据仍可能因为可疑研究实践（QRP）而带有偏倚。两种常见 QRP 是 p-hacking 与 HARKing，它们都会增加在 Meta 分析中高估效应的风险。
• 许多发表偏倚方法都建立在小样本研究效应这一思想之上。这类方法假定：只有那些效应量异常高的小研究才能达到显著，从而被发表。因此它们会在漏斗图中造成不对称。不过，这种不对称并不一定由发表偏倚引起，也可能来自其他相对“无害”的原因。
• 相对较新的 p-curve 方法则认为，只看数据中显著（\(p < 0.05\)）效应的分布模式，就足以评估证据价值。它既可以检验真实效应是否存在，也可以检验其是否缺失，还能估计真实效应量的大小。
• 选择模型是一类高度灵活的方法，可以模拟不同形式的发表偏倚过程。但它们只有在所假定的模型足够合理时才会给出有效结果，而且通常需要相当多的研究数。较简单的三参数模型则也可用于较小的数据集。
• 目前没有任何一种发表偏倚方法能在所有情境下持续优于其他方法。因此，更稳妥的做法是始终并行使用多种技术，并谨慎解释校正后的效应量。相较于现有统计方法，认真检索未发表证据往往更能有效降低发表偏倚风险。