dmetar.Rmd
dmetar
The dmetar
package serves as the companion R package for the online guide Doing MetaAnalysis in R  A Handson Guide written by Mathias Harrer, Pim Cuijpers, Toshi Furukawa and David Ebert. This freely available guide shows how to perform metaanalyses in R from scratch with no prior R knowledge required. The guide, as well as the dmetar
package, have a focus on biomedical and psychological research synthesis, but methods are applicable to other research fields too. The guide primarily focuses on two widely used packages for metaanalysis, meta
(Schwarzer, 2007) and metafor
(Viechtbauer, 2010), and how they can be applied in realworld use cases. The dmetar
package thus aims to provide additional tools and functionalities for researchers conducting metaanalyses using these packages and the Doing MetaAnalysis in R guide.
In this vignette, we provide a rough overview of the core functionalities of the package. An indepth introduction into the package and how its functions can be applied to “realworld” metaanalyses can be found in the online version of the guide. To get detailed documentation of specific functions, you can consult the dmetar
reference page.
Currently, the dmetar
package is still under development (version 0.0.9000). This means that, despite intense testing, we cannot guarantee that functions will work as intended under all circumstances and for all environments used. To report a bug, or ask a question, please contact Mathias (mathias.harrer@fau.de) or David (d.d.ebert@vu.nl)
Given that dmetar
is currently under development, the package is only available from GitHub right now. To install the development version, you can use the install_github
function from the devtools
package. Given that the package already passes the R CMD Check
, we aim to submit the package to CRAN in the near future after the development process has been completed.
Use the code below to install the package:
if (!require("devtools")) { install.packages("devtools") } devtools::install_github("MathiasHarrer/dmetar")
The package can then be loaded as usual using the library()
function.
library(dmetar)
The dmetar
package requires that R Version 3.6.3 or greater is installed and used in RStudio. You can check your current R version by running:
R.Version()$version.string
If you have an R version below 3.6.3 installed, installing dmetar
will likely cause an error (e.g, because the dependency metafor
was not found). This means that you have to update R. A tutorial on how to update R on your system can be found here.
The dmetar
package provides tools for different stages of the metaanalysis process. Functions cover topics such as power analysis, effect size calculation, smallstudy effects, publication bias, metaregression, subgroup analysis, risk of bias assessment, and network metaanalyses. Many dmetar
functions heavily interact with functions from the meta
and metafor
package to improve the work flow when conducting metaanalysis. Therefore, the meta
and metafor
package should be loaded from the library first.
To show some of the core functionality of the dmetar
package, we will use three datasets which come shipped with the package itself: ThirdWave
, MVRegressionData
and NetDataNetmeta
.
The dmetar
package contains two functions for a priori power analyses of a metaanalysis: power.analysis
and power.analysis.subgroups
. Let us assume that researchers expect to have approximately 18 studies in their metaanalysis, with moderate betweengroup heterogeneity and about 50 participants per arm and study. Will there be sufficient power to detect an assumed minimally important difference of \(d=0.18\)? The power.analysis
function can be used to answer this question.
power.analysis(d = 0.18, k = 18, n1 = 50, n2 = 50, heterogeneity = "moderate")
## Randomeffects model used (moderate heterogeneity assumed).
## Power: 83.86%
The dmetar
package includes several functions to calculate effect sizes needed for metaanalyses: NNT
, se.from.p
and pool.groups
.
Using the NNT
function, we can calculate the number needed to treat \(NNT\) for the first effect size in ThirdWave
(\(g\)=0.71). In this example, we use Furukawa’s method (Furukawa & Leucht, 2011), assuming a control group event ratio (CER
) of 0.2
NNT(0.71, CER = 0.2)
## Furukawa's method used.
## [1] 4.038088
We can also pool together two groups of a study into one group using the pool.groups
function once we have obtained the \(n\), mean and SD of each arm. This can be helpful if we want to avoid a unitofanalysis error. Here is an example:
pool.groups(n1 = 50, n2 = 65, m1 = 12.3, m2 = 14.8, sd1 = 2.45, sd2 = 2.89)
## Mpooled SDpooled Npooled
## 1 13.71304 2.969564 115
When extracting effect size data, studies sometimes only report an effect size of interest, and its \(p\)value. To pool effect sizes using functions such as the metagen
function, however, we need some dispersion measure (e.g., \(SE\), \(SD\) or the variance). The se.from.p
function can be used to calculate the standard error (\(SE\)), which can the be used directly for pooling using, for example, the metagen
function. Here is an example assuming an effect of \(d=0.38\), a \(p\)value of \(0.0456\) and a total \(N\) of \(83\):
se.from.p(effect.size = 0.38, p = 0.0456, N = 83)
## EffectSize StandardError StandardDeviation LLCI ULCI
## 1 0.38 0.1904138 1.734752 0.006795843 0.7532042
In biomedical literature, it is common to assess the Risk of Bias of included studies using the Cochrane Risk of Bias Tool. Such Risk of Bias assessments can be directly performed in RevMan, but this comes with certain drawbacks: RevMan graphics are usually of lower quality, and journals often require highresolution charts and plots; using RevMan along with R to perform a metaanalysis also means that two programs have to be used, which may consume unnecessary extra time; lastly, using RevMan to generate Risk of Bias summary plots also reduces the reproducibility of your metaanalysis if you decide to make all your other R analyses fully reproducible using tools such as RMarkdown.
The rob.summary
function allows you to generate RevManstyle Risk of Bias charts based on ggplot2
graphics directly in R from a R data frame. Here is an example:
rob.summary(data, studies = studies, table = TRUE)
The dmetar
package contains two functions related to the topic of moderator variables of metaanalysis results: subgroup.analysis.mixed.effects
and multimodel.inference
.
MixedEffects Subgroup Analysis
The first function, subgroup.analysis.mixed.effects
performs a subgroup analysis using a mixedeffects model (fixedeffects plural model; Borenstein & Higgins, 2013), in which subgroup effect sizes are pooled using a randomeffects model, and subgroup differences are assessed using a fixedeffect model. The function was built as an additional tool for metaanalyses generated by meta
functions. In this example, we therefore perform a metaanalysis using the metagen
function first.
meta < metagen(TE, seTE, data = ThirdWave, studlab = ThirdWave$Author, comb.fixed = FALSE, method.tau = "PM")
We can then use this meta
object called meta
as input for the function, and only have to specify the subgroups coded in the original data set we want to consider.
subgroup.analysis.mixed.effects(x = meta, subgroups = ThirdWave$TypeControlGroup)
## Subgroup Results:
## 
## k TE seTE LLCI ULCI p Q I2
## information only 3 0.4015895 0.1003796 0.205 0.598 6.315313e05 1.144426 0.00
## no intervention 8 0.5030588 0.1185107 0.271 0.735 2.187522e05 16.704467 0.58
## WLC 7 0.7822643 0.2096502 0.371 1.193 1.905075e04 22.167163 0.73
## I2.lower I2.upper
## information only 0.00 0.82
## no intervention 0.08 0.81
## WLC 0.42 0.87
##
## Test for subgroup differences (mixed/fixedeffects (plural) model):
## 
## Q df p
## Between groups 2.723884 2 0.2561628
##
##  Total number of studies included in subgroup analysis: 18
##  Tau estimator used for withingroup pooling: PM
Multimodel Inference
The multimodel.inference
function, on the other hand, can be used to perform Multimodel Inference for a metaregression model.
Here is an example using the MVRegressionData
dataset, using pubyear
, quality
, reputation
and continent
as predictors:
library(metafor) multimodel.inference(TE = 'yi', seTE = 'sei', data = MVRegressionData, predictors = c('pubyear', 'quality', 'reputation', 'continent'))
##

  0%

====  6%

=========  12%

=============  19%

==================  25%

======================  31%

==========================  38%

===============================  44%

===================================  50%

=======================================  56%

============================================  62%

================================================  69%

====================================================  75%

=========================================================  81%

=============================================================  88%

==================================================================  94%
##
##
## Multimodel Inference: Final Results
## 
##
##  Number of fitted models: 16
##  Full formula: ~ pubyear + quality + reputation + continent
##  Coefficient significance test: knha
##  Interactions modeled: no
##  Evaluation criterion: AICc
##
##
## Best 5 Models
## 
##
##
## Global model call: metafor::rma(yi = TE, sei = seTE, mods = form, data = glm.data,
## method = method, test = test)
## 
## Model selection table
## (Intrc) cntnn pubyr qulty rpttn df logLik AICc delta weight
## 12 + + 0.3533 0.02160 5 2.981 6.0 0.00 0.536
## 16 + + 0.4028 0.02210 0.01754 6 4.071 6.8 0.72 0.375
## 8 + + 0.4948 0.03574 5 0.646 10.7 4.67 0.052
## 11 + 0.2957 0.02725 4 1.750 12.8 6.75 0.018
## 15 + 0.3547 0.02666 0.02296 5 0.395 12.8 6.75 0.018
## Models ranked by AICc(x)
##
##
## Multimodel Inference Coefficients
## 
##
##
## Estimate Std. Error z value Pr(>z)
## intrcpt 0.38614661 0.106983583 3.6094006 0.0003069
## continent1 0.24743836 0.083113174 2.9771256 0.0029096
## pubyear 0.37816796 0.083045572 4.5537402 0.0000053
## reputation 0.01899347 0.007420427 2.5596198 0.0104787
## quality 0.01060060 0.014321158 0.7402055 0.4591753
##
##
## Predictor Importance
## 
##
##
## model importance
## 1 pubyear 0.9988339
## 2 continent 0.9621839
## 3 reputation 0.9428750
## 4 quality 0.4432826
To obtain a more robust estimate of the pooled effect size, especially when the betweenstudy heterogeneity of a metaanalysis is high, it can be helpful to search for outliers and recalculate the effects when excluding them.
The find.outliers
function automatically searches for outliers (defined as studies for which the 95%CI is outside the 95%CI of the pooled effect) in your metaanalysis and recalculates the results without these outliers. The function works for metaanalysis objects created with functions of the meta
package as well as the rma.uni
function in metafor
.
meta < metagen(TE, seTE, data = ThirdWave, studlab = ThirdWave$Author, method.tau = "SJ", comb.fixed = FALSE) find.outliers(meta)
## Identified outliers (randomeffects model)
## 
## "DanitzOrsillo", "Shapiro et al."
##
## Results with outliers removed
## 
## 95%CI %W(random) exclude
## Call et al. 0.7091 [ 0.1979; 1.2203] 5.0
## Cavanagh et al. 0.3549 [0.0300; 0.7397] 6.9
## DanitzOrsillo 1.7912 [ 1.1139; 2.4685] 0.0 *
## de Vibe et al. 0.1825 [0.0484; 0.4133] 10.4
## Frazier et al. 0.4219 [ 0.1380; 0.7057] 9.1
## Frogeli et al. 0.6300 [ 0.2458; 1.0142] 7.0
## Gallego et al. 0.7249 [ 0.2846; 1.1652] 6.0
## HazlettStevens & Oren 0.5287 [ 0.1162; 0.9412] 6.4
## Hintz et al. 0.2840 [0.0453; 0.6133] 8.1
## Kang et al. 1.2751 [ 0.6142; 1.9360] 3.5
## Kuhlmann et al. 0.1036 [0.2781; 0.4853] 7.0
## Lever Taylor et al. 0.3884 [0.0639; 0.8407] 5.8
## Phang et al. 0.5407 [ 0.0619; 1.0196] 5.4
## Rasanen et al. 0.4262 [0.0794; 0.9317] 5.1
## Ratanasiripong 0.5154 [0.1731; 1.2039] 3.3
## Shapiro et al. 1.4797 [ 0.8618; 2.0977] 0.0 *
## Song & Lindquist 0.6126 [ 0.1683; 1.0569] 5.9
## Warnecke et al. 0.6000 [ 0.1120; 1.0880] 5.3
##
## Number of studies combined: k = 16
##
## 95%CI z pvalue
## Random effects model 0.4708 [0.3297; 0.6119] 6.54 < 0.0001
##
## Quantifying heterogeneity:
## tau^2 = 0.0361 [0.0000; 0.1032]; tau = 0.1900 [0.0000; 0.3213];
## I^2 = 24.8% [0.0%; 58.7%]; H = 1.15 [1.00; 1.56]
##
## Test of heterogeneity:
## Q d.f. pvalue
## 19.95 15 0.1739
##
## Details on metaanalytical method:
##  Inverse variance method
##  SidikJonkman estimator for tau^2
##  Qprofile method for confidence interval of tau^2 and tau
forest(find.outliers(meta), col.study = "blue")
Influence Analysis can be helpful to detect studies which:
The InfluenceAnalysis
is a wrapper around several influence analysis function included in the meta
and metafor
package. It provides four types of influence diagnostics in one single plot. The function works for metaanalysis objects created by meta
functions, which can then be directly used as input for the function:
infan = InfluenceAnalysis(meta)
## [===========================================================================] DONE
plot(infan)
Additionally, the gosh.diagnostics
can be used to analyze influence patterns using objects generated by metafor
s gosh
function. The gosh.diagnostics
function uses unsupervised learning algorithms to determine effect sizeheterogeneity patterns in the metaanalysis data. We can use dmetars
inbuilt m.gosh
data set, which has been generated using metafor
s gosh
function as an example:
data("m.gosh") res < gosh.diagnostics(m.gosh)
##
## Perform Clustering...
## ========================================================================================== DONE
summary(res)
## GOSH Diagnostics
## ================================
##
##  Number of Kmeans clusters detected: 3
##  Number of DBSCAN clusters detected: 4
##  Number of GMM clusters detected: 4
##
## Identification of potential outliers
## 
##
##  Kmeans: Study 3, Study 4
##  DBSCAN: Study 3, Study 4
##  Gaussian Mixture Model: Study 3, Study 4
plot(res)
The package contains two functions to assess the potential presence of publication bias in a metaanalysis: eggers.test
and pcurve
. Both methods are optimized for conducting metaanalyses using the meta
package, and only have to be provided with a meta
metaanalysis results object.
Here is an example output for eggers.test
:
eggers.test(meta)
## Eggers' test of the intercept
## =============================
##
## intercept 95% CI t p
## 4.111 2.39  5.83 4.677 0.0002524556
##
## Eggers' test indicates the presence of funnel plot asymmetry.
And here for pcurve
:
pcurve(meta)
## Pcurve analysis
## 
##  Total number of provided studies: k = 18
##  Total number of p<0.05 studies included into the analysis: k = 11 (61.11%)
##  Total number of studies with p<0.025: k = 10 (55.56%)
##
## Results
## 
## pBinomial zFull pFull zHalf pHalf
## Rightskewness test 0.006 5.943 0.000 4.982 0
## Flatness test 0.975 3.260 0.999 5.158 1
## Note: pvalues of 0 or 1 correspond to p<0.001 and p>0.999, respectively.
## Power Estimate: 84% (62.7%94.6%)
##
## Evidential value
## 
##  Evidential value present: yes
##  Evidential value absent/inadequate: no
The package contains two utility functions for network metaanalysis using the gemtc
(Van Valkenhoef & Kuiper, 2016) and netmeta
(Rücker, Krahn, König, Efthimiou & Schwarzer, 2019) packages.
The first one, direct.evidence.plot
creates a plot for the direct evidence proportion of comparisons included in a network metaanalysis model and displays diagnostics proposed by König, Krahn and Binder (2013). It only requires a network metaanalysis object created by the netmeta
function as input. We will use dmetar
s inbuilt NetDataNetmeta
dataset for this example:
library(netmeta) data("NetDataNetmeta") nmeta = netmeta(TE, seTE, treat1, treat2, data=NetDataNetmeta, studlab = NetDataNetmeta$studlab) direct.evidence.plot(nmeta)
The second, sucra
, calculates the \(SUCRA\) for each treatment when provided with a mtc.rank.probability
network metaanalysis results object, or a matrix containing rank probabilities. In this example, we will use dmetar
s inbuilt NetDataGemtc
data set.
library(gemtc) data("NetDataGemtc") # Create Network MetaAnalysis Model network = mtc.network(data.re = NetDataGemtc) model = mtc.model(network, linearModel = "fixed", n.chain = 4, likelihood = "normal", link = "identity") mcmc = mtc.run(model, n.adapt = 5000, n.iter = 100000, thin = 10) rp = rank.probability(mcmc) # Create sucra sucra(rp, lower.is.better = TRUE)
Furukawa, T. A., & Leucht, S. (2011). How to obtain NNT from Cohen’s d: comparison of two methods. PloS one, 6(4), e19070.
Harrer, M., Cuijpers, P., Furukawa, T.A, & Ebert, D. D. (2019). Doing MetaAnalysis in R: A Handson Guide. DOI: 10.5281/zenodo.2551803.
König J., Krahn U., Binder H. (2013): Visualizing the flow of evidence in network metaanalysis and characterizing mixed treatment comparisons. Statistics in Medicine, 32, 5414–29
Rücker, G., Krahn, U., König, J., Efthimiou, O. & Schwarzer, G. (2019). netmeta: Network MetaAnalysis using Frequentist Methods. R package version 1.01. https://CRAN.Rproject.org/package=netmeta
Schwarzer, G. (2007), meta: An R package for metaanalysis, R News, 7(3), 4045.
Van Valkenhoef, G. & Kuiper, J. (2016). gemtc: Network MetaAnalysis Using Bayesian Methods. R package version 0.82. https://CRAN.Rproject.org/package=gemtc
Viechtbauer, W. (2010). Conducting metaanalyses in R with the metafor package. Journal of Statistical Software, 36(3), 148. URL: https://www.jstatsoft.org/v36/i03/