\pkgrobumeta: An \proglangR-package for robust variance estimation in meta-analysis

In traditional meta-analysis, inverse variance weights play two roles, First, they produce the most efficient estimate of $\bm{\beta}$. Second, the statistical estimator of the variance requires the weights to be inverse; when they are not, the estimate of the variance can be severely biased. In RVE, inverse variance weights are also used, but here they are only important for efficiency.

Weighting improves statistical efficiency by allocating more weight to studies that are more precise (i.e., have a smaller variance). Importantly, RVE provides asymptotically accurate estimates of the standard errors and valid inferential methods for any set of weights. This means the question of weights in RVE is only about efficiency. The most efficient weights would be $\mathbf{W}_{j}=\bm{\Sigma}_{j}^{-1}$. Clearly this is problematic since RVE is used when $\bm{\Sigma}_{j}^{-1}$ is unknown. For this reason, in RVE weights are developed based on a working model of the unknown covariance structure $\bm{\Sigma}_{aj}$ so that $\mathbf{W}_{j}\approx\bm{\Sigma}_{aj}^{-1}$. Two working models (hierarchical and correlated effects) have been proposed to model the types of dependency most commonly encountered in meta-analytic research. Both of these weighting methods are only approximately efficient and in practice, researchers are urged to choose the weighting method based on the most common source of dependence in their data.

The first adjustment is to the $\mathbf{V}^{R}$ estimator which can be written more generally as

where in general $\mathbf{A}_{j}\mathbf{e}_{j}=\mathbf{A}_{j}(\mathbf{I}-\mathbf{H})_{j}\bm{\epsilon}$, $\mathbf{H}=\mathbf{X}\mathbf{Q}\mathbf{X}^{\prime}\mathbf{W}$, $\mathbf{Q}=(\mathbf{X}^{\prime}\mathbf{WX})^{-1}$, and $\bm{\epsilon}$ is a $\Sigma\mathrm{k}_{j}\>\mathrm{x\>1}$ vector of true residuals. Hedges et al. (2010) proposed that this estimator could be improved by using

where $\mathbf{I}_{j}$ is a $k_{j}$ by $k_{j}$ identity matrix. Simulations indicate that this adjustment is inadequate except in very specific cases and Tipton (in press) develop a new solution based on extensions to the Bias Reduced Linearization (BRL) method proposed by McCaffrey et al. (2001). Here since $\mathbf{V}^{*}$ also depends on $\mathbf{W}$ and working covariance matrices $\mathbf{\Sigma}_{aj}$, the $\Sigma\mathrm{k}_{j}$ by $\Sigma\mathrm{k}_{j}$ adjustment matrix $\mathbf{A}_{j}$ is specified separately for the hierarchical

and correlated effects weighting model

where $\mathbf{H}_{jj}=\mathbf{X}_{j}\mathbf{Q}\mathbf{X}_{j}^{\prime}\mathbf{W}_{j}$. Additionally, if non-efficient weights are required (in \pkgrobumeta non-efficient weights can be specified using the \codeuserweights option) the adjustment matrix is specified as follows

where $(\mathbf{I}-\mathbf{H})_{j}$ includes the $k_{j}$ rows of $(\mathbf{I}-\mathbf{H})$ associated with study $j$ and the average variance of the study $j$’s effect sizes, $\nu_{\cdot j}=\left(1/k_{j}\right)\Sigma\nu_{ij}$.

In the above, for a general matrix $\mathbf{B}$ note that $\mathbf{I}=\mathbf{B}^{-1/2}\mathbf{BB}^{-1/2}$, and $\mathbf{B}^{-1/2}=\mathbf{P}\bm{\Lambda}^{-1/2}\mathbf{P}^{\prime}$ where $\mathbf{P}\bm{\Lambda}\mathbf{P}^{\prime}$ is the eigen decomposition. See the Appendix of Tipton (in press) for a detailed account of the development of these adjustments and their use with non-inverse variance weighting schemes.

The second adjustment proposed by Hedges et al. (2010) was to use the t-distribution with $df=m-p$ when making statistical inferences using $\mathbf{V}^{R}$. Building off of Bell and McCaffrey (2002), however, Tipton (in press) suggests using the Satterthwaite approximation (Satterthwaite, 1946) to estimate the distribution of $\mathrm{\mathbf{V}}^{*}$ with $\chi^{2}_{df_{Sk}}$, where the degrees of freedom $df_{Sk}$ is equal $2/c\nu^{2}$, and $c\nu$ is the coefficient of variation.

These degrees of freedom are calculated separately for each regression coefficient under studying using

where $\lambda_{jk}$ are the eigenvalues of $\mathbf{W}^{-1/2}(\sum\mathbf{g}_{jk}\mathbf{g}_{jk}^{\prime})\mathbf{W}^{-1/2}$ where $\mathbf{g}_{jk}=(\mathbf{I}-\mathbf{H})_{j}^{\prime}\mathbf{A}_{j}\mathbf{W}_{j}\mathbf{X}_{j}\mathbf{Q}l_{k}$. Here $l_{k}$ is a vector of $0$’s and $1$’s that indicates the regression coefficient, $\mathbf{A}_{j}$ and $\mathbf{W}_{j}$ are defined as before. Further details can be found in Tipton (in press). Simulations by Tipton (in press) show that the Satterthwaite approximation is valid so long as the $df>4$. This means that RVE can be effectively used to make inferences even with a very small number of studies.

In the first example we use data from a meta-analysis conducted by Oswald et al. (2013) examining the predictive validity of the Implicit Association Test (IAT) and various explicit measures of bias for a variety of criterion measures of discrimination. The \codeoswald2013 dataset provided in \pkgrobumeta contains $308$ effect sizes collected from $m=46$ individual studies.

For this example, we will look at a subset of the \codeoswald2013.ex1 dataset which contains only the studies where the criterion measure of discrimination used was either neurological activity or response latency. The (\codeoswald2013.ex1) dataset contains $32$ effect sizes collected from $m=9$ individual studies. Additionally, the reported correlations have been transformed to Fisher’s $Z_{r}$ (Gleser and Olkin, 2009). As the majority of effect-sizes in our dataset result from each study including multiple measures taken on the same individuals, we will use correlated effects weights to fit our model.

When performing a meta-regression it is standard practice to first calculate an average effect size without conditioning on the study covariates. In RVE this can be followed with a sensitivity analysis which will be discussed in more detail in the following section. The fitting function for the intercept-only model is specified with the \coderobu() function

oswald_intercept <- robu(formula = effect.size   1, data = oswald2013.ex1, studynum = Study, var.eff.size = var.eff.size, rho = .8, small = TRUE)

and the results displayed with \codeprint() as follows

print(oswald_intercept)

RVE: Correlated Effects Model with Small-Sample Corrections

Model: effect.size   1

Number of studies = 9 Number of outcomes = 32 (min = 1 , mean = 3.56 , median = 2 , max = 11 ) Rho = 0.8 I.sq = 82.05272 Tau.Sq = 0.1849812

Estimate StdErr t-value df P(|t|>) 951 intercept 0.277 0.181 1.53 7.84 0.164 -0.141 0.695 — Signif. codes: < .01 *** < .05 ** < .10 * — Note: If df < 4, do not trust the results

In addition to the coefficient table, the output from \codeprint() contains the formula call used to fit the model, summary information for the total number of studies and outcomes included in the analysis, and a coefficient table. In addition, output for the correlated effects model contains the user-specified $\rho$ value, the estimate of between-study variance in study-average effect sizes $\tau^{2}$, and a descriptive statistic for the ratio of true heterogeneity to total variance across the observed effect sizes, $I^{2}$ (Higgins et al., 2003). $I^{2}$ is calculated using

where $\mathrm{Q_{E}}$ is defined as the weighted residual sum of squares as in equation 8. Unlike $\tau^{2}$, $I^{2}$ is not on the same scale as the effect size index used in the analysis and is better suited to anwsering questions regarding the proportion of observed variance resulting from real effect size differences.

As mentioned previously, an estimate of $\tau^{2}$ is required to develop approximately efficient weights for the meta-regression model. The method of moments estimator for $\tau^{2}$ requires a value is specified for $\rho$. For this reason Hedges et al. (2010) suggests conducting a sensitivity analysis to determine the effect of $\rho$ on $\tau^{2}$. Generally, results from simulation studies suggest that $\tau^{2}$ and the meta-regression coefficients are relatively insensitive to changes in $\rho$ (Ishak et al., 2008; Williams, 2012; Tipton, 2013); however, if between-study covariances were to be considerably smaller than within-study covariances, these results might not hold. It should also be noted that a conservative approach and an external information approach have also been advocated for estimating $\rho$ in terms of developing efficient weights. The conservative approach involves setting $\rho$ to $1$, which ensures studies don’t receive additional weight due to having more effect sizes. Additional information on these strategies can be found in Hedges et al. (2010).

The \codesensitivty() function in \pkgrobumeta computes $\tau^{2}$, the average effect size and associated standard error for different values of $\rho$ in the interval $(0,1)$.

R > sensitivity(oswald_intercept)

Type Variable rho=0 rho=0.2 rho=0.4 rho=0.6 rho=0.8 rho=1 1 Estimate intercept 0.277 0.277 0.277 0.277 0.277 0.277 2 Std. Err. intercept 0.181 0.181 0.181 0.181 0.181 0.181 3 Tau.Sq 0.184 0.184 0.184 0.185 0.185 0.185

From the sensitivity analysis it seems $\tau^{2}$ and subsequently the average effect size are relatively robust to different values of $\rho$.

In meta-analysis, forest plots provide a graphical depiction of effect size estimates and their corresponding confidence intervals. The \codeforest.robu() function in \pkgrobumeta can be used to produce forest plots for RVE meta-analyses. The function requires the \pkggrid (R Core Team, 2013) package and is based on examples provided in Murrell (2011). As is the case in traditional forest plots, point estimates of individual effect sizes are plotted as boxes with areas proportional to the weight assigned to that effect size. Importantly, here the weight is not necessarily proportional to the effect size variance or confidence intervals, since the combined study weight is divided evenly across the study effect sizes. Two-sided 95% confidence intervals are calculated for each effect size using a standard normal distribution and plotted along with each block. The overall effect is included at the bottom of the plot as a diamond with width equivalent to the confidence interval for the estimated effect.

The RVE forest function is designed to provide users with forest plots which display each individual effect size used in the meta-analysis, while taking into account the study- or cluster-level properties inherent to the RVE analysis. As such, the user must specify columns from their original dataset that contain labels for the study or cluster and for the individual effect sizes. For instance, labels for the study column might be author names with corresponding publication years, and would be specified in the \codeforest.robu function using the \codestudy.lab argument. Labels for individual effect sizes might be “Math Score” or “Reading Score” for a meta-analysis that included such measures or as simple as “Effect Size 1” and “Effect Size 2,” and can be specified with the \codees.lab argument. Any number of additional columns can be specified to be plotted along side the confidence interval column and can be specified with the following syntax \code"arg1" = "arg2" where \code"arg1" is the title of the column on the forest plot, and \code"arg2" is the name of the column from the original \codedata.frame that contains the information to be displayed. An RVE forest plot with study labels corresponding to the author name and publication year, and effect size labels corresponding to the criterion for the Oswald data (see Figure 1 on page 1) is specified using the syntax below. In addition, a column of the RVE weights titled “Weight” and a column of effect sizes titled “Effect Sizes” is also included to illustrate the syntax for adding additional columns.

forest.robu(oswald_intercept, es.lab = "Crit.Cat", study.lab = "Study", "Effect Size" = effect.size, # optional column "Weight" = r.weights) # optional column

In the following section we walk through an example where \pkgrobumeta is used to fit simulated data described in Tanner-Smith and Tipton (2013). The dataset in question resulted from a fictional meta-analysis that examined the effectiveness of a brief alcohol intervention on alcohol intake. Intervention effectiveness was determined by differences in alcohol consumption between treatment and controls following the intervention. In this example standardized mean difference effect sizes are provided using Hedges’ g, however, any effect size can be used with the \coderobu() function. The \codehierdat dataset contains $m=15$ unique treatment centers with $68$ effect sizes. The covariates of interest in this example are \codefollowup, which indicates the length of the follow-up period before alcohol consumption was measured, in months and \codebinge which indicates whether or not the alcohol consumption measure targeted binge drinking specifically. As the majority of dependency in this meta-analysis resulted from studies nested within treatment centers, we use the hierarchical effects weights to fit our model.

In RVE dependencies result from hierarchically structured data. Borrowing terminology from the multilevel and hierarchical linear modeling literature, this hierarchical structure takes the form of lower-level observations nested within higher-level groupings. Within the RVE framework we are generally concerned with two-level structures where level-1 (the lowest level) represents effect sizes in the correlated effects case and studies in the hierarchical effects case, while level-2 represents the study-level in the correlated effects case and a higher-order grouping level such as research lab in the hierarchical effects case.

Since covariates can vary both between-level-2 and within-level-2 it is an important pre-processing step in any RVE analysis to identify the nature of this variation prior to model fitting. For example, at first glance the \codefollowup variable appears to vary both within- and between-level-2. This indicates the length of follow-up periods varies both across studies conducted by the same treatment center, and then between the treatment centers themselves. Failure to distinguish between these effects can complicate interpretation later on as the regression coefficient actually represents a weighted combination of both between- and within-level-2 effects, which can be problematic when the effects differ in direction or magnitude. For instance, if treatment centers which produce higher quality studies also tend to employ longer follow-up periods the effect of follow-up length on treatment effect sizes might differ when examined within a given treatment center’s studies, compared to the effect of follow-up length across treatment centers.

To partition these two types of effect we can create group mean $X_{\bullet j}$ and group centered $X^{C}_{ij}$ versions of \codefollowup using the convenience functions \codegroup.mean() and \codegroup.center() provided in \pkgrobumeta.

R > hierdat$followup_{m}<-group.mean(hierdat$followup, hierdat$studyid)R>hierdat$followup_c <- group.center(hierdat$followup,hierdat$studyid)

Subsequently, including both \codefollowup_m and \codefollowup_c in our meta-regression model as follows

allows us to decompose the effects of follow-up length into two independent estimates where $\beta_{1}$ and $\beta_{2}$ represent the effect of a 1-month increase in $X_{\bullet j}$ and a 1-month increase in $X^{C}_{ij}$ on $T_{ij}$, respectively. Another important feature of centering is that it often indicates that a variable only varies at the study level, not effect size level. We therefore encourage researchers to examine these variables before proceeding with their analysis.

Now that we have completed the necessary data pre-processing we can fit our meta-regression model using hierarchical effects weighting with the following command

R > model_hier <- robu(effectsize   followup_c + followup_m + binge, data = hierdat, modelweights = "HIER", studynum = studyid, var.eff.size = var, small = TRUE)

The \codeprint() method can be used to present the results

R > print(model_hier)

RVE: Hierarchical Effects Model with Small-Sample Corrections

Model: effectsize   followup_c + followup_m + binge

Number of clusters = 15 Number of outcomes = 68 (min = 1 , mean = 4.53 , median = 2 , max = 29 ) Omega.sq = 0.1650524 Tau.Sq = 0.02479249

Estimate StdErr t-value df P(|t|>) 951 intercept -0.154226 0.147671 -1.044 6.07 0.3361 -0.51461 0.20615 2 followup_c -0.000162 0.000695 -0.233 1.30 0.8470 -0.00534 0.00502 3 followup_m 0.003467 0.002320 1.495 3.28 0.2243 -0.00357 0.01050 4 binge 0.666645 0.115623 5.766 4.33 0.0035 0.35514 0.97815 *** — Signif. codes: < .01 *** < .05 ** < .10 * — Note: If df < 4, do not trust the results

The results of the meta-regression suggest the p-value for \codefollowup_c and \codefollowup_m are not reliable, as the Satterthwaite degrees of freedom are less than 4. This may result from the fact that center is contributing more than half of the effect sizes included in our meta-analysis did not actually vary the followup interval across studies, resulting in an unbalanced covariate. This is in line with results from simulation studies which suggest the Satterthwaite degrees of freedom are typically smaller for unbalanced covariates (Tipton, in press). The effect of whether or not the study used a binge drinking measure (\codebinge), however, was found to be positive and significant at $\alpha=0.05$. This suggests that studies which used binge drinking as their criterion measure were associated with increased treatment effect sizes compared to those studies that used some other measure of alcohol consumption.