metaplus: An R Package for the Analysis of Robust Meta-Analysis and Meta-Regression

2.1 Models

The random effect meta-analysis model assumes that the observed treatment effect \(Y_i\) for study \(i\) is \[\begin{equation*} Y_i=\mu+E_i+\epsilon_i, \end{equation*}\] where \(\mu\) is the overall mean for the studies, \(E_i\) is a random effect with mean zero, and \(\epsilon_i\) is a normally distributed error with variance \(\sigma^2_i\) for study \(i\), where the within study variance \(\sigma^2_i\) is assumed to be known.

An extension to the random effect meta-analysis model is to include covariates to explain (Berkey et al. 1995) the heterogeneity. Incorporating this into the meta-analysis model we obtain:

\[ \begin{equation*} Y_i=\mu+X_i^T\beta+E_i+\epsilon_i, \end{equation*} \] where \(X_i\) is a vector of covariate values for study \(i\), and \(\beta\) is a vector of the corresponding parameters.

In metaplus there are three available random effect distributions:

• Normal The probability density function for study \(i\) is

\[ \begin{equation*} f\left(Y_i|X_i;\mu, \tau\right)=\frac{1}{\sqrt{2\pi \left(\sigma_i^2+\tau^2\right)}} \exp\left( -\frac{\left(Y_i-\mu-X_i^T\beta\right)^2} {2\left(\sigma_i^2+\tau^2\right)} \right) . \end{equation*} \]

• Robust \(t\)-distribution Introduced as one of a number of distributions for robust meta-analysis by Lee and Thompson (2008) and Baker and Jackson (2008). This replaces the normal random effect distribution with a \(t\)-distribution. The degrees of freedom (\(\nu\)) of the \(t\)-distribution controls the heaviness of the tails, and is estimated from the data. The probability density function no longer has a closed-form expression, requiring integration over the \(t\)-distribution random effect as

\[ \begin{equation*} f\left(Y_i|X_i;\mu, \tau, \nu \right)= \frac{1}{\sqrt{2 \pi \sigma_i^2}} \int_{-\infty}^{\infty} \exp \left(\frac{-\left(Y_i-\mu-X_i^T\beta-\eta\right)^2}{2 \sigma_i^2}\right) g\left(\eta|\tau, \nu\right)d\eta , \end{equation*} \] where \(g\left(\eta | \tau, \nu \right)\) is the distribution function of a \(t\)-distribution with \(\nu\) degrees of freedom \[ \begin{equation*} g\left(\eta | \tau, \nu \right)=\frac{\Gamma\left(\left(\nu+1\right)/2\right)} {\tau\sqrt{\pi\nu}\Gamma\left(\nu/2\right)} \left(1+\frac{\eta^2}{\nu\tau^2}\right)^{-\left(\left(\nu+1\right)/2\right)}. \end{equation*} \]

• Robust mixture

This assumes that a study can belong to one of two classes, with each class a standard random effect model with the same mean but different random effect variance, which is higher for the outlier class (Beath 2014). The robust meta-analysis model takes the form

\[ \begin{equation*} Y_{i|k}=\mu+E_{i|k}+\epsilon_i, \end{equation*} \]

where \(\epsilon_i\) is as for the standard model, but \(E_{i|k}\) is now a random effect dependent on the class, where \(k=1, 2\) indexes the classes, with \(k=1\) corresponding to standard studies and \(k=2\) to outlier studies, with random effect variances \(\tau_1^2,\tau_2^2\) respectively, with the restriction that \(\tau_2^2>\tau_1^2\), and again zero mean. The probability density function becomes the weighted sum of the probability density function for each class, with weights the proportion of studies in each class \(\pi_1, \pi_2\) for the standard and outlier studies, respectively, as

\[ \begin{equation*} f\left(Y_i|X_i;\mu, \tau_1, \tau_2, \pi_1, \pi_2\right)=\sum_{k=1}^2 \pi_k \frac{1}{\sqrt{2\pi}} \left( \frac{1}{\sigma_i^2+\tau_k^2} \right) ^{1/2}\exp\left(-\frac{1}{2} \frac{\left(Y_i-\mu-X_i^T\beta\right)^2}{\sigma_i^2+\tau_k^2}\right). \end{equation*} \] with the constraint that \(\pi_1+\pi_2=1\)

2.2 Profile likelihood confidence intervals

A difficulty with the use of standard maximum likelihood techniques for random effect models is that it produces biased estimates for the variance of the random effect, which results in biased estimates of the standard errors for the parameters of interest, and therefore poor coverage. The solution for meta-analysis has been the use of REML, however this is difficult for the robust models. However, profile likelihood based confidence intervals [Pawitan (2001)][p 61] have been found to be superior (Hardy and Thompson 1996) to those obtained using REML, and they are used for all fitted models. The profile likelihood confidence intervals are obtained from routines based on the mle2 procedure in the package bbmle which provides an extended version of mle.

2.3 Parametric boostrap

Testing for the need for the robust distributions requires a test of \(\nu=0\) for the \(t\)-distribution and \(\pi_2=0\) for the robust mixture. Both tests involve a test of a parameter on the edge of the parameter space, so the usual asymptotic theory cannot be used. The solution is to use the parametric bootstrap, which involves simulating data sets under the null hypothesis and calculating the test statistic for each. The observed test statistic is then compared to the simulated test statistics to calculate a \(p\)-value. Application of the parametric bootstrap to mixture models is described in McLachlan (1987).

2.4 Other computational details

For both the robust models the starting values are important, as the optimisation may converge to a local minima. For the \(t\)-distribution, starting values are obtained by first fitting a standard normal random effect model, and the starting values from this model used with a range of starting values for the \(t\)-distribution degrees of freedom, and the maximum likelihood chosen. Numerical integration was used to obtain the likelihood. One difficulty is that the model is not identifiable when \(\tau^2=0\) as the \(t\)-distribution degrees of freedom doesn’t affect the likelihood, and this causes problems with the fitting optimiser. To avoid this a model was fitted with \(\nu=0\), to allow \(\tau^2=0\), and the likelihood from this model used if it was equal or larger than given by the optimisation with \(\nu\) unconstrained.

For the robust mixture model a generalized EM (GEM) algorithm was used. Starting values were obtained by partitioning the model into outliers and non-outliers, starting with no outliers. Then each study in turn was selected as an outlier, approximate starting values determined, and the model fitted. The model with the largest likelihood for all models with none or a single outlier was chosen. The process was then repeated adding each possible outlier in turn for two outlier models. This process is repeated until the likelihood cannot be increased, giving the final model. This approach was significantly faster than generating a random selection of outliers. It may be possible to improve the process further by only looking at studies with high residuals as possible outliers.

4.1 Intravenous magnesium in acute myocardial infarction

A number of studies have been performed to determine the effectiveness of intravenous magnesium in acute myocardial infarction, and the data obtained from Sterne et al. (2001). The studies have caused considerable controversy, as the results of a single large study ISIS-4 (ISIS-4: Collabarative Group 1995) contradicts the results of a meta-analysis. Higgins and Spiegelhalter (2002) discusses some of the history and some suggested methods from a Bayesian perspective. Woods (2002) comments on the variability between studies due to timing of infusion, and Downing (1999) on the higher level of dose used in ISIS-4, with a more recent meta-analysis by Li et al. (2009). Of interest is whether, given the heterogeneity between studies, the ISIS-4 study is unusual. The data have been obtained in the form of effect estimates and standard errors, but if raw data in the form of number of events per number of patients, then these can be converted using the escalc procedure in the metafor package, for example. We can perform the standard random effect meta-analysis, and obtain parameter estimates as follows:

data(mag)
mag.meta <- metaplus(mag$yi, mag$sei, slab=mag$study)
summary(mag.meta)
#>          Est. 95% ci.lb 95% ci.ub   pvalue
#> muhat -0.7463   -1.2583   -0.3428 0.000501
#> tau2   0.2540                             
#> 
#>      logLik      AIC      BIC
#>   -19.68459 43.36918 44.91436

Adding plotci=TRUE will add a plot giving details of the profile confidence intervals, and is shown in Figure 1. When the quadratic approximation to the likelihood holds the shape of the curve will be in the form of a “V”. In this case, the shape is not symmetric, so this does not hold. This is confirmed by the lack of symmetry of the confidence interval for muhat.

mag.meta <- metaplus(mag$yi, mag$sei, slab=mag$study, plotci=TRUE)

Figure 1: Profile plot.

The forest plot showing the studies and overall effect can be shown using plot(mag.meta), and is shown in Figure 2. The metaplus package makes extensive use of the metafor package which allows the parameters for the forest plot in metafor to be used when creating forest plots. As the results for the magnesium studies are actually log odds ratios it is more useful to produce results with units of odds ratios. This can be obtained by annotating the horizontal axis with odds ratios corresponding to the log odds, and requesting an exponential transformation for the coefficients, as shown in the following code. The plot is shown in Figure 3.

plot(mag.meta)

Figure 2: Forest plot for magnesium studies.

plot(mag.meta, atransf=exp, at=log(c(.01, .1, 1, 10, 100)), xlab="Odds Ratio")

Figure 3: Forest plot for magnesium studies (odds ratios).

Repeating the meta-analysis but now allowing for a \(t\)-distribution for the random effect by adding the random=“t-dist” option. From the summary the estimate of vinv, the inverse degrees of freedom, is zero indicating infinite degrees of freedom, or a normal distribution. The BIC is also a guide, with an increase for the \(t\)-distribution model indicating that it is likely that a standard normal is the correct model.

mag.tdist <- metaplus(mag$yi, mag$sei, slab=mag$study,
        random="t-dist")
summary(mag.tdist)
#>          Est. 95% ci.lb 95% ci.ub   pvalue
#> muhat -0.7463   -1.2583   -0.3430 0.000501
#> tau2   0.2540                             
#> vinv   0.0000                             
#> 
#>      logLik      AIC      BIC
#>   -19.68459 45.36918 47.68695

We can confirm this with the testOutliersMeta command, which performs a parametric bootstrap to obtain the null distribution of the likelihood ratio statistic for the test that \(\nu^{-1}=0\), required as the test is on the boundary of the parameter space. Note that this may take several minutes for the default of 999 simulations.

summary(testOutliersMeta(mag.tdist))
#> Observed LRT statistic 0.0 p value 1

We can repeat using the robust mixture distribution for the random effect. The variance of both the random effect for standard studies (tau2) and for outlier studies (tau2out) are very close indicating that there are no outlier studies and this is confirmed by the outlier test.

mag.mix <- metaplus(mag$yi, mag$sei, slab=mag$study,
      random="mixture")
summary(mag.mix)
#>                     Est.  95% ci.lb  95% ci.ub   pvalue
#> muhat         -7.463e-01 -1.258e+00 -3.427e-01 0.000501
#> tau2           2.540e-01                               
#> tau2out        2.540e-01                               
#> Outlier prob.  6.412e-05                               
#> 
#>      logLik      AIC      BIC
#>   -19.68459 47.36918 50.45954
summary(testOutliersMeta(mag.mix))
#> Observed LRT statistic 0.0 p value 1

4.2 CDP choline for cognitive and behavioral disturbances

This meta-analysis evaluates the effect of CDP choline for cognitive and behavioral disturbances associated with chronic cerebral disorders in the elderly (Fioravanti and Yanagi 2005). A study (Bonavita 1983, note that this was Senin 2003 in the original package due to an error in assigning the study labels) was previously determined to be an outlier by Gumedze and Jackson (2011). We can first fit a standard random effect meta-analysis, as previously.

data(cdp)
cdp.meta <- metaplus(cdp$yi, cdp$sei, slab=cdp$study)
summary(cdp.meta)
#>          Est. 95% ci.lb 95% ci.ub pvalue
#> muhat 0.38944   0.07269   0.76634 0.0218
#> tau2  0.14666                           
#> 
#>      logLik      AIC      BIC
#>   -8.198544 20.39709 21.00226

Note that Bonavita 1983 does have an unusually high value, as shown in Figure 5. We can then fit a robust model using the \(t\)-distribution.

cdp.tdist <- metaplus(cdp$yi, cdp$sei, slab=cdp$study,
        random="t-dist")
summary(cdp.tdist)
#>            Est. 95% ci.lb 95% ci.ub  pvalue
#> muhat 1.945e-01 5.294e-02 3.611e-01 0.00899
#> tau2  4.504e-05                            
#> vinv  2.023e+00                            
#> 
#>      logLik      AIC      BIC
#>   -4.058334 14.11667 15.02442
summary( testOutliersMeta(cdp.tdist))
#> Observed LRT statistic 8.3 p value 0.001

As a rough guide, the decrease in AIC and BIC demonstrate that the model is an improvement, and this is confirmed with the outlier test. We can also perform with the robust mixture.

cdp.mix <- metaplus(cdp$yi, cdp$sei, slab=cdp$study,
        random="mixture")
summary(cdp.mix)
#>                 Est. 95% ci.lb 95% ci.ub  pvalue
#> muhat         0.1910    0.0563    0.3479 0.00711
#> tau2          0.0000                            
#> tau2out       3.1558                            
#> Outlier prob. 0.1237                            
#> 
#>      logLik      AIC      BIC
#>   -3.007145 14.01429 15.22463
summary( testOutliersMeta(cdp.mix))
#> Observed LRT statistic 10.4 p value 0.001

The output from the robust mixture model has an interesting feature. For standard studies the estimated random effect variance is zero, indicating that only the outlier studies are contributing to the heterogeneity. The posterior probability of each study being an outlier can be obtained as

cdp.mix.outlierProbs <- outlierProbsMeta(cdp.mix)

and plotted using plot(cdp.mix.outlierProbs) in Figure 4. This shows clearly that Bonavita 1983 does have a posterior probability of nearly 1.0 of being an outlier. The other studies have a non-zero posterior probability of being outliers, as there is an overlap between the distribution of the standard and outlier studies.

plot(cdp.mix.outlierProbs)

Figure 4: Outlier probabilities for CDP studies.

Lastly, we can now produce a forest plot with the results of all three models, using the extrameta parameter to add the robust models, as plot(cdp.meta, extrameta=list(cdp.tdist, cdp.mix)), and these are shown in Figure 5. The effect of the robust models is to down-weight the Senin 2003 study, which has the consequence of both reducing the overall effect estimate and it’s standard error.

plot(cdp.meta, extrameta=list(cdp.tdist, cdp.mix))

Figure 5: Forest plot for CDP Studies.

4.3 Exercise for depression

This example is a meta-analysis of trials of exercise in the management of depression (Lawlor and Hopker 2001). Higgins and Thompson (2004) used the data as an example of meta-regression using a number of covariates, which will be limited here to a single covariate, the duration of trial. First perform the meta-analysis using standard normal random effect and the robust mixture model.

data(exercise)
exercise.meta <- metaplus(exercise$smd, sqrt(exercise$varsmd),
      mods=exercise[, c("duration"), drop=FALSE], slab=exercise$study)
summary(exercise.meta)
#>             Est. 95% ci.lb 95% ci.ub   pvalue
#> muhat    -2.8994   -4.3006   -1.5222 0.000884
#> tau2      0.1171                             
#> duration  0.2078    0.0584    0.3632 0.011570
#> 
#>      logLik      AIC      BIC
#>   -8.133435 22.26687 23.17462
exercise.mix <- metaplus(exercise$smd, sqrt(exercise$varsmd),
      mods=exercise[, c("duration"), drop=FALSE],
      slab=exercise$study, random="mixture")
summary(exercise.mix)
#>                   Est. 95% ci.lb 95% ci.ub   pvalue
#> muhat         -2.88472  -4.10690  -1.48262 0.000649
#> tau2           0.59398                             
#> tau2out        0.00000                             
#> Outlier prob.  0.74831                             
#> duration       0.21086   0.07824   0.34617 0.007123
#> 
#>     logLik      AIC     BIC
#>   -7.69139 25.38278 26.8957
exercise.testOutliers <- testOutliersMeta(exercise.mix)
summary(exercise.testOutliers)
#> Observed LRT statistic 0.9 p value 0.068
exercise.outlierProbs <- outlierProbsMeta(exercise.mix)

The test for outliers was close to significant (p=0.068), however a conservative approach seems appropriate, by using the robust model where the presence of outliers is not conclusive but there is a reasonable amount of evidence that there are outliers, as in this case. Note also that the \(p\)-value is different from that obtained in Beath (2014), due to the use of randomly generated data in the parametric bootstrap. Running the parametric bootstrap with a large number of simulations showed that the \(p\)-value was in fact near 0.04. Using plot(exercise.outlierProbs) the outlier probabilities are shown in Figure 6 where the study by Reuter is an obvious outlier with a posterior probability greater than 0.9. This study is a dissertation and was not published in a peer-reviewed journal, and was not included in a later meta-analysis by Krogh et al. (2011). There is also strong evidence of the effect of trial duration.

plot(exercise.outlierProbs)

Figure 6: Outcome probabilities for depression versus exercise.

To calculate the effect at each of Weeks 4, 8 and 12 it is easiest to centre the data at those times and fit a meta-regression for each. The intercept for each meta-regression will then be the estimated mean effect at that time. We also fit a model without including the covariate for study duration. The forest plot is shown in Figure 7. This shows that the effect of exercise decreases rapidly the longer the trial, possibly indicating a placebo effect that rapidly wears off. It would also be possible to plot the results from the standard random effect models as well.

exercise$duration4 <- exercise$duration-4
exercise$duration8 <- exercise$duration-8
exercise$duration12 <- exercise$duration-12

exercise.nodurn <- metaplus(exercise$smd, sqrt(exercise$varsmd),
    label="Random Mixture (No Duration)", slab=exercise$study,
    random="mixture")
exercise.wk4 <- metaplus(exercise$smd, sqrt(exercise$varsmd),
    mods=exercise[, c("duration4"), drop=FALSE],
    label="Random Mixture (Week 4)", slab=exercise$study,
    random="mixture")
exercise.wk8 <- metaplus(exercise$smd, sqrt(exercise$varsmd),
    mods=exercise[, c("duration8"), drop=FALSE],
    label="Random Mixture (Week 8)", slab=exercise$study,
    random="mixture")
exercise.wk12 <- metaplus(exercise$smd, sqrt(exercise$varsmd),
    mods=exercise[, c("duration12"), drop=FALSE],
    label="Random Mixture (Week 12)", slab=exercise$study,
    random="mixture")
plot(exercise.nodurn, extrameta=list(exercise.wk4, exercise.wk8,
      exercise.wk12))

Figure 7: Forest plot for exercise versus depression studies.