r语言如何计算t分布临界值_如何利用R语言进行meta分析—详细教程（2）

如何利用R语言进行meta分析—详细教程（2）

--（如何在meta分析中进行离群值检测与meta回归）

• 小伙伴们大家好，在我上一篇的推文里详细讲解了如何利用R语言的meta包并以OR值和 95%CI作为效应尺度进行合并作meta分析的文章，大家有需要的可以去看看。由于最近我做了不少关于meta分析的研究，也有了更多的经验，现在计划一点点向大家分享出来。其实meta分析并不只是用R才能做，其他的还有SAS,STATA,RevMan,MIX,Meta-DiSc等等许多，但是可以说R语言是个meta分析的全才，不仅可以完成经典的meta分析功能，一些新近出现的meta分析方法如网状meta分析等，也都可以用它完美实现，而且我也发现用R做的话，它的代码其实并不多，就那么几行核心代码，但是能做的事情却很多，许多的数据一丢进来，跑个代码就OK。所以我还是推荐大家用R做meta分析，它并不像大家想象的那么难。不过当然，这里面不同的分析方法和程序却有很多，在这篇文章例我主要向大家介绍metaplus包的使用方法，请小伙伴们结合自己的实际需要进行使用。
• 对于meta分析的理论部分还有R软件包以及RStudio的获取大家可以看我公众号（全哥的学习生涯）的上一篇推文《如何利用R语言进行meta分析—详细教程（以OR值和 95%CI作为效应尺度进行合并）》，在这里就不赘述了。metaplus 程序包是 R 专用于实现 robust meta 分析和 Meta 回归的程序包metaplus 程序包，总体上可以分为 Meta 分析和Meta回归两大功能。其主要执行函数命令有metaplus 函数、testOutliers 函数和 outlierProbs 函数。其中，metaplus 函数用于模型拟合、选择合适模型及绘制森林图等；testOutliers 及 outlierProbs函数用于检测数据中的离群值。我们在这里以图1的数据为例（数据来源：The effectiveness of exercise as an intervention in the management of depression: systematic review and meta-regression analysis of randomised controlled trials, BMJ,2001, 322(7289): 763-767.）,以data <- read.csv()命令录入到R中。[图1中的数据说明--study：纳入研究的作者名；smd（study effect estimate）：研究效应量；varsmd（study effect variance）：研究效应量的方差；sesmd（study effect standard error）：研究效应量的标准误；duration：干预时间，单位为周；duration4：干预时间与 4 周差值；duration8：干预时间与8周差值；duration12：干预时间与 12 周差值] 。
• 若有读者想与我取得联系，可以在关注我的公众号：全哥的学习生涯，公众号内有我的个人联系方式~

图1

1 模型的拟合与选择

安装包与录入数据：（整理后的数据见图2）

install.packages(“metaplus”)

library(“metaplus”)

对”duration”列按升序排列：data<-data[order(data$duration), ]

图2

1.1 标准正态随机效应模型

图3

data.meta<-metaplus（smd, sesmd, mods=duration, slab=study,

random="normal", data=data）

summary（data.meta）

命令中 smd 为观察的效应量，sesmd 为效应量的标准误，mods 为各研究相对应的协变量数据框架，slab 为与每个研究相对应的字符串向量（Vector of character strings corresponding to each study），random 为随机效应的类型（normal 表示正态分布，t-dist 表示 t 分布，mixture 表示混合分布）,data为数据。上述数据的标准正态随机效应模型 Meta 分析结果见表2和表3。结果如图3。

1.2 t 分布随机效应模型

图4

metaplus 软件包进行 t分布随机效应模型 Meta 分析与标准正态随机效应模型类似，只需要将参数 random 赋值为“t-dist”即可，具体命令如下：

data.meta<-metaplus（smd,sesmd ,mods= duration, slab= study, random= "t-dist",data=data）

summary（data.meta）

上述数据的标准正态随机效应模型 Meta 分析结果见表4和表5。结果如图4。

1.3 混合随机效应模型

图5

metaplus 进行混合随机效应模型合并只需要参数 random 赋值为

“mixture”即可，具体命令如下：

data.meta<-metaplus（smd,sesmd ,mods= duration, slab= study, random= "mixture",data=data）

summary（data.meta）

上述数据的标准正态随机效应模型 Meta 分析结果见表6和表7。结果如图5。

2 离群值检验

对于 t 分布随机效应模型，离群值的检验需要自由度趋近为无穷大；对于混合随机效应模型，离群值的检验需要离群研究概率密度函数的权重趋近为零。由于这两个检验都涉及参数空间边界上的参数，渐近理论不适用，所以在零假设下运用参数自举的方法来检验统计量的经验分布。然后将观察到的统计量似然比与此分布进行比较，以确定P 值。由于参数自举中使用大量随机生成的数据，P 值是不同的。metaplus 提供了 testOutliers 函数进行离群值检验，命令如下:

data.testOutliers<- testOutliers(data.meta)

summary(data.testOutliers)

图6

离群值检验的结果见图6。

此外，metaplus 还提供 outliersProbs 用于计算离群概率，对于混合随机效应模型和 t 分布随机效应模型可以计算出每个研究的离群概率，并使用plot 命令绘制离群概率图像，结果见图7。命令如下：

data.outlierProbs<- outlierProbs(data.meta)

plot(data.outlierProbs)

图7

离群概率图结果见图7。

3 模型选择

对于相关模型的拟合与选择，总体上依据AIC（赤池信息量准则）和 BIC（贝叶斯信息规则）的数值。AIC 和 BIC 都是对模型的拟合效果进行评价的指标，若 AIC 和 BIC 值越小，则模型对数据的拟合越好。当 t 分布随机效应模型和混合随机效应模型比标准正态随机效应模型的 AIC 和 BIC 的值低，那么倾向于选择 t 分布随机效应模型和混合随机效应模型。结合 testOutliers 和 outlierProbs 函数，可以对研究数据中的离群值进行检测，若数据中存在明显的离群值，则倾向于选择混合随机效应模型。因本例中数据存在明显的离群值，故选择应用混合随机效应模型。

4 输出meta分析结果

经过前期的模型拟合与选择，并完成离群值检验后，metaplus 程序包可以调用 metafor 程序包的绘制森林图功能，完成森林图的绘制。命令如下：

plot（data.meta，extrameta=NULL，xlab=“Forset Plot-全哥制作”）

森林图结果见图9。

图9

5 meta回归

目前，metaplus 程序包尚未有预测方法，计算4、8 和 12 周效应的替代方法是在这些时间将数据置于中心，并对每个数据进行 Meta 回归。命令如下：

data$duration4<-data$duration4

data$duration8<-data$duration8

data$duration12<-data$duration12

data.nodurn<-metaplus（smd, sesmd, label ="Random Mixture-全哥制作 （No Duration）", slab=study, random="mixture", data=data）

data.week4<-metaplus(smd,sesmd, mods=duration4,label="Random Mixture-全哥制作 (Week4)", slab=study, random="mixture", data=data)

data.week8<-metaplus（smd, sesmd, mods=duration8,label ="Random Mixture-全哥制作 （Week8）", slab=study, random="mixture", data=data）

data.week12<-metaplus（smd, sesmd, mods=duration8,label ="Random Mixture-全哥制作 （Week12）", slab=study, random="mixture", data=data）

plot(data.nodurn,extrameta=list（data.week4,data.week8,data.week12)，xlab="Standardised mean difference-全哥制作"）

Meta 回归森林图见图10，Meta 回归的截距可以估计对应的平均效应。图10结果表明，研究效应在持续时间延长时会下降，可能与安慰剂效应减低或消失有关。

图10

6 模型的比较

在 Meta 分析中，最常见的两种模型是固定效应模型和随机效应模型。固定效应模型假设所有纳入的研究拥有共同的真实效应量，而随机效应模型中的真实效应量随研究的不同而改变。基于不同模型的运算，所得到的合并后的效应量的均数值也不相同。尽管不一定适用于实践，但为了数据处理的方便，传统上假定随机效应模型的真实效应量呈正态分布。而 t 分布可以更灵活地来替代这个假设，通过马尔科夫链蒙特卡罗方法实现正态分布和 t 分布的偏态扩张，比常规正态分布更加灵活。混合随机效应模型可以通过比较具有和不具有离群值的模型来确定是否存在异常值，并且使用后验概率来确定，其优点在于 Meta 分析时考虑了含离群值的研究，且具有适当的权重。

7 总结

metaplus 运用的是标准 Meta分析法和基于轮廓似然置信区间的 robust 方法。metaplus 程序包的目的不是为了取代一般的 Meta分析程序包，而是为了提供额外的专业分析，如离群值的检测和离群值的概率计算等。metaplus 程序包可以绘制离群概率图，并调用 Metafor 程序包绘制森林图等结果图形。然而，值得注意的是，metaplus 程序包的不足在于进行大量数据分析计算时会延长计算时间，特别是进行testOutliers() 命令时。另一个不足之处是 metaplus程序包只适用于标准正态随机效应、t 分布随机效应和混合正态随机效应 3种模型，没有扩展到效应量的其他分布类型。总之，metaplus 程序包可以在实现基本 Meta分析功能的基础上进行离群值的检测和概率计算。

最后，如果屏幕前的你对R语言学习还有什么问题或者看法，可以在我的公众号（全哥的学习生涯）给我留言，公众号里也有我的个人联系方式，我也希望可以结合更多志同道合的伙伴。

感谢你的阅读。

参考文献

[1]Beath KJ. Metaplus: an R package for the analysis of robust metaanalysis and meta-regression. R Journal, 2016, 8(1): 5-16.

[2] Lawlor DA, Hopker SW. The effectiveness of exercise as an intervention in the management of depression: systematic review and meta-regression analysis of randomised controlled trials. BMJ,2001, 322(7289): 763-767.(本文数据来源)