R语言元分析专题第六章：异质性分析
在元分析中，我们的目的是将多个研究中的效应整合为一个效应。例如，尽管元分析的整体效应可能很小，但可能有几个研究的效应量很大。这种信息在混合效应中就会丢失，有个别研究效应量很大、有的很小，知道这些情况非常重要。因为元分析中有一些极端的效应值，也称为outliers，这些outliers会扭曲整体效应，我们需要知道没有这些outliers时整体效应如何。在一个元分析内，效应量的变化程度被称为异质性（heterogeneity）。评估异质性非常重要，因为如果研究中有两个或更多亚团体（subgroups)，它们有不同的真实效应（true effect），那么这些研究的异质性会很高。异质性信息对于研究来说非常重要，因为它可以让我们找到特定的干扰、效应量较低或较高的群体。从统计的角度来说，异质性过高也是存在问题的，非常高的异质性意味着这些研究没有共同之处，数据没有**真实的**正确效应，报告混合效应也就没有任何意义（Borenstein etal., 2011)。
6.1 异质性思想
Rücker和其同事命名了元分析中的三种异质性：
第1和第3点可以通过将研究范围限定在特定的干预类型、人群和结果上而得到一定程度的控制。第2点，必须是我们在混合（pooling）研究时进行评估，这是本章的重点内容。
6.2 异质性指标
6.2.1 选择异质性指标
一般来说，当我们在元分析中评估并报告异质性时，需要一个稳健的、不易受统计功效影响的值。当研究的数量（k）提高时、或者当精度（precision，即一个研究的样本量N）提高时，Cochran’s Q都会提高。因此，Q值以及它是否显著很大程度上依赖于元分析的大小以及统计功效。因此，当评估异质性时我们不应该仅仅依赖于Q值。
I2值对元分析中研究数量（K）的变化不敏感。因此广泛地应用于医学和心理学研究，特别是存在“rule of thumb”可以对它进行解释（Higgins et al., 2003）：【rule of thumb,大拇指原则，又称“经验法则”，是一种可用于许多情况的简单的、经验性的、探索性的但不是很准确的原则】
尽管广泛地应用于文献当中，它也不是评估异质性非常充分的指标，因为它仍然依赖于元分析中研究的统计功效（Rücker et al. 2008; Borenstein et al. 2017）。如前所述，是由非采样误差引起的变异量。如果我们的研究样本越来越大，采样误差趋近于0，然而同时会趋近于100%，这仅仅是因为单个研究的样本量（N）更大（不是衡量采样误差引起的变异，但样本量增大却变大）因此，仅仅依赖于也不是合适的做法。
τ2对研究的数量和精度(统计功效)均不敏感。然而有时很难解释与实际的观点（practical standpoint）有多大的相关。
**预测区间**（Prediction intervals）是一个很好的可以克服该缺陷（IntHout et al., 2016）的方法(IntHout et al., 2016)，因为预测区间考虑了研究间的变异。预测区间给出了一个范围，我们可以基于**现在元分析的证据**预测**未来研究的效应落在**这个范围内。例如，如果预测区间完全落在支持干预效果的“阳性”那边，我们可以非常有信心地说，尽管效应会有变化，但在以后的各种研究情景中，干预措施至少在某些程度上是有效的。如果预测区间包括0，我们就不那么确定了，尽管在医学和心理学研究中预测区间很宽这种现象非常普遍。
6.3评估混合效应量的异质性
在元分析中使用metagen(), metabin(), 或者metacont()计算好混合效应量后，可以很容易地提取**三个最常见的异质性指标**。我们以章节4.2.2中的混合效应模型元分析为例，结果存储在结果m.hksj中。
得到元分析异质性指标的方法之一是输入元分析结果：
print(m.hksj)结果如下：
我们可以看到输出中包含**3个异质性指标**（H也是，但在这里不作陈述）：
如何解释这些指标呢？上述三个指标说明，该元分析有**中等到很强的异质性**。**预测区间非常宽**，下限小于0，我们不能非常自信地说干预在每种情景中都有阳性效果。有可能干预措施在以后的一些情况下不会产生阳性效果，甚至有可能产生较小的负性效果，即使效应量很高（0.5935），也存在这种可能。
6.3.1当输出中没有呈现异质性指标
在函数metagen, metabin,或者metacont中的参数设置可以决定是否输入某些指标，但所有的指标都会存储在元分析客体中，只是有些没有呈现出来而已。我们可以使用$操作输出想要的指标，非常容易。例如下表：
6.4 检测异常值&有影响力的个案
如前所述，**研究间的异质性**（between-study heterogeneity ）也可能由另一项或多项研究造成的，这些研究不是存在**极端效应值**（extreme effect sizes ）的情况。特别是，当这些研究的**质量很低**或者**研究的数量很小**的时候，这可能会扭曲我们对混合效应的估计。如果在元分析中把这些异常值排除，再求合并效应值是比较合理的做法。
另一方面，我们也想知道计算的合并效应值是否稳健，也就是说该效应并不严重依赖于某一项单独的研究。因此，我们也想知道是否有研究把我们的分析结果明显地推向一个方向。这样的研究被称为**有影响力的案例**（influential cases ）。
需要注意的是，有许多方法可以检测异常值和有影响力的案例，这里不会对这些方法进行全面地描述。如果你想了解更多有关这个主题的基础知识，我们推荐阅读Wolfgang Viechtbauer和MikeCheung 的论文(Viechtbauer and Cheung, 2010 )。
6.4.1 检测极端效应值（异常值）
直接探测异常值的一种常用方法是，如果某个研究的置信区间与混合效应的置信区间没有重叠则被认为是异常值。这意味着，当某一项研究的估计效应量过于极端，使我们很有把握确定该研究不能成为所汇总的结果中的效应值“总体”的一部分时，我们将其定义为一个异常值（I.e.,个别研究与整体效应存在显著差异）。为了探测这样的异常值，我们使用dplyr包中的`filter()`函数，使用该函数可以搜索所有的研究：
95%置信区间的上限低于混合效应值置信区间的下限(I.e., 极端小的效应)
95%置信区间的下限高于混合效应值置信区间的上限(I.e., 极端大的效应)
以混合效应模型元分析的结果m.hksj为例 。输出元分析混合效应量置信区间的上限和下限，如下代码及结果。如果是固定效应，则输出这些值使用代码：$lower.fixed 和$upper.fixed。
m.hksj$lower.randomm.hksj$upper.random置信区间为[0.389147, 0.7979231]。使用dplyr包中的函数**spot.outliers.random**或** spot.outliers.random**，以及置信区间的上限、下限值我们可以过滤掉异常的研究。先加载该包：Library(dplyr)
以混合效应模型元分析m.hksj为例，使用函数**spot.outliers.random()**。目前，R还不知道这个函数，我们可以将代码全部复制并粘贴到RStudio左下窗格的控制台console中，然后点击回车。
spot.outliers.random<-function(data){data<-dataAuthor<-data$studlablowerci<-data$lowerupperci<-data$upper m.outliers<-data.frame(Author,lowerci,upperci) te.lower<-data$lower.random te.upper<-data$upper.random dplyr::filter(m.outliers,upperci < te.lower) dplyr::filter(m.outliers,lowerci > te.upper)}现在，我们唯一要告诉`spot.outliers.random()` 的是——我们要检查异常值的元分析输出，这里我们检查混合效应元分析m.hksj的异常值。
spot.outliers.random(data = m.hksj)结果：
我们可以看到，该函数已经检测到了两个异常值。从较低的ci值（lowerci），即两个研究置信区间的下限，我们可以看到这两个研究都有极高的正向效应，因为它们的下限比混合效应值的置信区间的上限都要高得多，即g = 0.798。
因此，我们可以进行一个**敏感性分析**（sensitivity analysi），排除这两个异常值。这可以通过元分析中的`update.meta()`函数来实现，更新m.hksj中未包含异常值的元分析输出结果。我们可以看出异质性小了很多。
m.hksj.outliers <-update.meta(m.hksj,subset=Author !=c("DanitzOrsillo","Shapiro et al."))m.hksj.outliers如果你进行的是一个**固定效应模型**的元分析，那么整个过程都是类似的，仅仅需要将random改为fixed。spot.outliers.fixed函数的代码，具体如下：
spot.outliers.fixed<-function(data){ data<-data Author<-data$studlab lowerci<-data$lower upperci<-data$upper m.outliers<-data.frame(Author,lowerci,upperci) te.lower<-data$lower.fixed te.upper<-data$upper.fixed dplyr::filter(m.outliers,upperci < te.lower) dplyr::filter(m.outliers,lowerci > te.upper)}6.5 影响力分析
我们已经说明了如何在元分析中探测和剔除极端效应值（异常值）。然而，正如我们之前所提到的，不仅是统计上的异常值会影响混合效应值的稳健性。元分析中的有些研究也可能对我们的整体结果产生非常大的影响。例如，原本我们可能会发现整体效应（overall effect ）并不显著，但如果移除某一个特定的研究，分析就显著了。当我们报告元分析的结果时，这些信息也非常重要。
比排除异常研究更深入的探测数据的方法，如**Leave-One-Out-method**，该方法每次排除一项研究，重新计算元分析K-1次。通过这种方式，我们可以更容易检测到那些对元分析的结果影响最大的研究，以及评估这种影响是否会扭曲混合效应量（Viechtbauer and Cheung, 2010 ）。因此，这种分析被称为影响力分析（Influence Analyses）。可以使用**dmetar**包中的`InfluenceAnalysis()`函数进行并使结果可视化。
Library(dmetar)如果您不想使用dmetar包，可以使用函数的源代码（书中有链接）。该函数依赖于ggplot2、ggrepel、forcats、dplyr、grid、gridExtra、metafor和meta包。
我们需要定义InfluenceAnalysis 函数中的几个参数。
影响力分析代码：
Library(dmetar)InA <- InfluenceAnalysis(x = m.hksj, random = TRUE)plot(InA)plot(InA, “influence”) #输出influence 图plot(InA, "Baujat") #Baujat图plot(InA, "I2") #输出按I^2排序的图plot(inf.analysis, "es") #输出按Effect size排序的图输出结果:
6.5.1 影响力分析（Influence Analyses）
在第一个分析结果中，我们可以看到不同的影响指标，这些图中包含元分析中的每个研究。影响力分析的类型由Viechtbauer和Cheung提出(Viechtbauer和Cheung，2010)。接下来，让我们讨论一下最重要的子图（subplots ）：
然而，通常情况下，您没有必要深究单个指标的计算。根据经验法则（arule of thumb），有影响力的研究是一个在图中具有非常极端值的研究。Viechtbauer和Cheung还提出了临界值（cut-offs），根据此临界值判断一个研究是否是有影响的研究（influential case）（其中p是模型系数的数量，k是研究的数量）：
如果采用这些临界值将一个研究确定为有影响力的研究，它的值将显示为红色（在例子中“Dan”或Danitz-Orsillo的研究标记为红色）。
请注意，正如Viechtbauer和Cheung所强调的那样，这些临界值在某种程度上是一种随意设定的阈值。因此，我们不能只看研究的颜色，也要看图表的结果，综合起来解释结果。在例子中，折线图中只有Danitz-Orsillo的研究被定义为有影响力的研究，但实际上在大多数图形中是有两个峰值的（另一个是Shapira等人的研究）。因此，我们还可以将Shapiro等人的研究定义为一个有影响力的个案。
前述分析中，我们发现“Danitz-Orsillo”和“Shapiro et al.的研究可能是有影响力的。这是一个有意思的发现，因为先前我们在观察统计异常值（outliers）的情况下也检测到了同样的研究。这进一步证实了这两项研究可能扭曲了对混合效应值估计，并可能引起元分析中的部分组间异质性。
6.5.2 Baujat图
Baujat图(Baujat etal., 2002 )是一种一个诊断性图，用于检测影响一个元分析的异质性的研究。该图横轴上的Cochran’s Q值显示了每项研究对整体异质性的贡献，纵轴表示它对混合效应值的影响。由于我们想要评估的是异质性及对异质性有影响的研究，所以我们需要重点查看位于图右侧的研究，主要右侧的研究导致了异质性。如果一项研究对总体异质性有很大影响，而同时对综合效应值却影响不大（e.g., 因为该研究的样本量非常小），这一点就更加需要重要。因此，**Baujat 图右侧下面的研究，对我们来说十分重要**。你可能已经发现，在这一区域发现的仅有的两项研究正是我们之前已经检测到的两项研究（Danitz & Orsillo, Shapiro et al ）。这些研究对总体效应的影响不大(大概是因为它们的样本量非常小)，但它们大大增加了元分析的异质性。
6.5.3“去一”分析（Leave-One-Out Analyses）
在森林图中，每次去掉一个研究，重新计算混合效应，产生两个图，它们提供了相同的数据，但按不同的值排序。第一个图按异质性排序（从低到高）。我们在图中可以看到，排除Danitz & Orsillo和Shapiro等人的研究，异质性最低，再次证实两项研究是引起异质性的罪魁祸首。第二个图按效应值大小排序（从低到高）。在这里，我们可以看到排除一项研究后，混合效应值的估计是如何变化的。同样，由于这两个异常研究的效应值非常高，我们发现当它们被移除时，混合效应值变得最小。总而言之，例子的异常值探测和影响分析的结果均指向相同的两个研究，它们很可能是会扭曲效应值估计及精确度。因此，我们应该报告排除这些研究的敏感性分析（sensitivity analysis）。
6.6 GOSH Plot Analysis
上一节，我们使用影响分析（influence analyses）和“去一法”（leave-one-out method）探究了元分析结果的稳健性。还有一种更为复杂的方法探究效应量和异质性的模式，即GOSH 图（GraphicDisplay of Heterogeneity plots）（Olkin et al., 2012）。即拟合所有研究子集的元分析模型。相较于“去一法”，我们不仅要拟合K-1个模型，还要拟合所有2k-1个可能的研究组合。因此，当参与元分析的研究数量比较多的时候，跑出GOSH图会让计算量变得异常庞大。在R中，可以最多拟合一百万个随机选择模型。
一旦模型被计算，就可以绘制它们，该图以混合效果量作为x轴，研究间的异质性作为y轴。这可以让我们找到一些特定的模式，例如不同效果量的子集群。这将说明实际上在我们的数据中，并不是只有一个效果量的总体，提醒我们要做一个亚群体分析。如果在我们的样本中的效果量是同质的，那异质性图将显示为有一个峰值的对称分布。Gosh函数在数据包**metaphor**中。
我们使用m.hksj数据绘制异质性图。在绘制前，我们需要将meta包创造的对象（object）转化成可以使用gosh函数的metafor元分析对象。在metafor中进行元分析的函数是`rma()`。只需要提供三个数据：效应量（TE）、标准误（seTE）和研究间异质性估计量（method.tau）。因为我们想完全重复结果，因此在`rma()`中设置test=“knha”校正。其语句如下：
m.rma <- rma(yi = m.hksj$TE, sei = m.hksj$seTE, method = m.hksj$method.tau, test = "knha")如果使用的是固定效应模型，需要设置method=“FE”。
使用m.rma对象来产生异质性图，完成时间取决于参与分析的研究数量多少。
dat.gosh <- gosh(m.rma)plot(dat.gosh, alpha=0.1, col="blue")我们可以看到数据存在两种模式：低效果量和低异质性，高效果量和相当大的异质性：似乎数据存在两个子集。
现在我们已经清楚了数据的效果量-异质性的模式，那么真正重要的问题是：哪些研究导致了这些模式，哪些研究属于一个子集群？回答这一问题可以使用`gosh.diagnostics()`函数。这个函数使用三个集群算法（类似于监督的机器学习）来探测在GOSH图数据中的集群，并且确定是哪些研究的贡献最大。该函数在dmetar包中。这个函数也需要安装dplyr,cluster, mvtnorm , factoextra, fpc, cowplot, reshape2 和 flexmix包。代码如下：
gosh.diagnostics(dat.gosh)结果：
我们看到这里有三种算法：k-means, DBSCAN 和Gaussian Mixture Model，已经检测到有三个研究（研究3、10、16）可能导致了集群的不平衡。这些研究看起来几乎完全可以解释第二个高效果量和高异质性集群。
如果我们移除这三个研究后再进行元分析的话，将会发生什么呢？
metagen(TE, seTE, data=madata, studlab=paste(Author), comb.fixed = FALSE, comb.random = TRUE, method.tau = "SJ", hakn = TRUE, prediction=TRUE, sm="SMD", exclude = c(3, 10, 16))结果：
我们可以看到研究间的异质性已经降到了0%，这意味着现在参与元分析的这些研究源于一个同质性的总体或亚总体。
R元分析专题介绍：本活动是由荷兰心理统计联盟组织的第三次线上同辈互助学习活动，小组成员每周学习一章，并整理笔记录屏内容共享，从R基础到元分析效应量计算，作图，元分析高阶，如多水平meta, MetaSEM, network Meta-Analysis 。学习资料为《Doing Meta-Analsis in R 》这本书。笔记为小组成员整理（付费与否全由小组成员决定）