Cluster-Bootstrapping a meta-analytic selection model

Preliminary Analysis

As a little warm-up exercise, here is a basic random effects meta-analysis of these data, fit via the metafor package (Viechtbauer, 2010). We use cluster-robust standard errors to account for effect size dependency.

# Estimate random effects model
RE_mod <- rma.uni(yi, vi = vi, data = lehmann_dat, method = "ML")

# Calculate cluster-robust standard errors
RE_robust <- robust(RE_mod, cluster = study, clubSandwich = TRUE)
RE_robust

## 
## Random-Effects Model (k = 81; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.1009 (SE = 0.0245)
## tau (square root of estimated tau^2 value):      0.3176
## I^2 (total heterogeneity / total variability):   81.54%
## H^2 (total variability / sampling variability):  5.42
## 
## Test for Heterogeneity:
## Q(df = 80) = 246.9683, p-val < .0001
## 
## Number of estimates:   81
## Number of clusters:    41
## Estimates per cluster: 1-6 (mean: 1.98, median: 1)
## 
## Model Results:
## 
## estimate      se¹    tval¹     df¹    pval¹   ci.lb¹   ci.ub¹     
##   0.2069  0.0569   3.6331   23.41   0.0014   0.0892   0.3245   ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 1) results based on cluster-robust inference (var-cov estimator: CR2,
##    approx t-test and confidence interval, df: Satterthwaite approx)

The random effects model indicates an average effect size of about 0.21 standard deviations and substantial heterogeneity, with \(\hat\tau = 0.32\). The cluster-robust standard error is about 27% bigger than the model-based standard error (not shown) because the latter does not account for dependent effect sizes.

Here is a contour-enhanced funnel plot of the data:

funnel(RE_mod, refline = 0, level=c(90, 95, 99), shade=c("white", "gray55", "gray75"))

The funnel plot shows some asymmetry, suggesting that there is reason to be concerned about selective reporting bias in these data.

A Selection Model

For starters, we will fit a very simple selection model, with a single step in the probability of reporting at \(\alpha = .025\). This is what’s come to be called the three-parameter selection model (McShane et al., 2016). The step is defined in terms of a one-sided p-value, so \(\alpha = .025\) corresponds to the point where an effect size estimate in the theoretically expected direction would have a regular, two-sided p-value of .05, right at the mystical threshold of statistical significance. We fit the model using metafor’s selmodel() function:

RE_sel <- selmodel(RE_mod, type = "stepfun", steps = .025)
RE_sel

## 
## Random-Effects Model (k = 81; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0811 (SE = 0.0261)
## tau (square root of estimated tau^2 value):      0.2848
## 
## Test for Heterogeneity:
## LRT(df = 1) = 37.3674, p-val < .0001
## 
## Model Results:
## 
## estimate      se    zval    pval   ci.lb   ci.ub    
##   0.1328  0.0655  2.0267  0.0427  0.0044  0.2612  * 
## 
## Test for Selection Model Parameters:
## LRT(df = 1) = 1.7646, p-val = 0.1840
## 
## Selection Model Results:
## 
##                      k  estimate      se     zval    pval   ci.lb   ci.ub    
## 0     < p <= 0.025  25    1.0000     ---      ---     ---     ---     ---    
## 0.025 < p <= 1      56    0.5485  0.2495  -1.8097  0.0703  0.0594  1.0375  . 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The selection parameter represents the probability that an effect size estimate that is not in the theoretically expected direction or not statistically significant at the conventional level would be included in the synthesis, relative to the probability that an affirmative, statistically significant effect size estimate would be included. In this example, the probability of selection is estimated as 0.55, meaning only about 55% of the non-significant results that we would expect were generated are actually reported. Adjusting for this selection bias, the model estimates an overall average effect size of 0.13 SD—smaller than the unadjusted random effects estimate—with heterogeneity of \(\hat\tau = 0.28\).

The problem with this analysis is that the selection model is set up under the assumption that the effect size estimates are all independent. As a result, the reported standard errors are probably smaller than they should be and the confidence intervals are narrower than they should be. We’ll use cluster bootstrapping to get standard errors and confidence intervals that should better account for effect size dependency.

Clustering and unclustering

The Lehmann dataset has one row per effect size, sometimes with multiple rows per study, so we need to modify the data to have one row per study. There are at least two ways to accomplish this. One option is to create a dataset consisting only of study-level IDs, then merge it back on to the full data to get the effect-size level data:

# Make a dataset of study IDs
cluster_IDs <- data.frame(study = unique(lehmann_dat$study))

glimpse(cluster_IDs)

## Rows: 41
## Columns: 1
## $ study <chr> "Banas, 2014 ", "Berthold, 2013 ", "Bigelow et al., 2013 ", "Ble…

# Merge with full data
full_dat <- merge(cluster_IDs, lehmann_dat, by = "study")

full_dat %>% select(study, yi, vi) %>% glimpse()

## Rows: 81
## Columns: 3
## $ study <chr> "Banas, 2014 ", "Berthold, 2013 ", "Bigelow et al., 2013 ", "Big…
## $ yi    <dbl> 0.05716727, 0.55411916, 0.31467980, -0.73259462, 0.07921700, -0.…
## $ vi    <dbl> 0.10267348, 0.06030579, 0.29520322, 0.53354343, 0.02692608, 0.05…

Another option is to use the dplyr::nest_by() function to nest the data by cluster (Wickham et al., 2019). Then, we can use tidyr::unnest() to recover the effect size level data (Wickham et al., 2019). Like so:

# Nest the data for each study
lehmann_nested <- nest_by(lehmann_dat, study, .key = "data")

glimpse(lehmann_nested)

## Rows: 41
## Columns: 2
## Rowwise: study
## $ study <chr> "Banas, 2014 ", "Berthold, 2013 ", "Bigelow et al., 2013 ", "Ble…
## $ data  <list<tibble[,49]>> [<tbl_df[1 x 49]>], [<tbl_df[1 x 49]>], [<tbl_df[2…

# Recover the full dataset
full_dat <-
  lehmann_nested %>%
  unnest(data)

full_dat %>% select(study, yi, vi) %>% glimpse()

## Rows: 81
## Columns: 3
## Groups: study [41]
## $ study <chr> "Banas, 2014 ", "Berthold, 2013 ", "Bigelow et al., 2013 ", "Big…
## $ yi    <dbl> 0.05716727, 0.55411916, 0.31467980, -0.73259462, 0.07921700, -0.…
## $ vi    <dbl> 0.10267348, 0.06030579, 0.29520322, 0.53354343, 0.02692608, 0.05…

We will follow the latter strategy for the remainder of our example.

With this nest-and-unnest approach, we can fill in a little bit more of our function skeleton:

fit_selmodel <- function(
    dat,   # dataset with one row per cluster
    index, # vector of indexes used to create the bootstrap sample
    ...    # any further arguments
) { 
  
  # take subset of data
  boot_dat_cluster <- dat[index, ]
  
  # expand to one row per effect size
  boot_dat <- tidyr::unnest(boot_dat_cluster, data)
  
  # fit selection model
  
  # compile parameter estimates into a vector
  
}

Selection model function

Next, we need to complete the function by writing code to fit the selection model. This is a little bit involved because of the way the metafor package implements the Vevea-Hedges selection model. We first need to fit a regular random effects model using metafor::rma.uni(), then pass the result to the metafor::selmodel() function to fit a selection model, and then pull out the parameter estimates as a vector. To make the code clearer, we will move this step out into its own function:

run_sel_model <- function(dat, steps = .025) {
  
  # initial random effects model
  RE_mod <- metafor::rma.uni(
      yi = yi, vi = vi, data = dat, method = "ML"
  )
  
  # fit selection model
  res <- metafor::selmodel(
    RE_mod, type = "stepfun", steps = steps,
    skiphes = TRUE, # turn off SE calculation
    skiphet = TRUE # turn off heterogeneity test
  )
  
  # compile parameter estimates into a vector
  c(beta = res$beta[,1], tau = sqrt(res$tau2), delta = res$delta[-1])
  
}

Note the use of skiphes and skiphet arguments in the selmodel() call, which skip the calculation of standard errors and skip the calculation of the test for heterogeneity. We don’t need the standard errors here because we’re going to use bootstrapping instead, and we’re not interested in the heterogeneity test. Turning off these calculations saves computational time.

A further complication with fitting a selection model is that the parameter estimates are obtained by maximum likelihood, using an iterative optimization algorithm that sometimes fails to converge. To handle non-convergence, we will pass our function through purrr::possibly() so that errors are suppressed, rather than causing everything to grind to a halt (Wickham et al., 2019). We set the otherwise argument so that non-convergent results are returned as NA values:

run_sel_model <- purrr::possibly(run_sel_model, otherwise = rep(NA_real_, 3))

A first bootstrap

Here is the completed fitting function called fit_selmodel():

fit_selmodel <- function(dat, 
                         index = 1:nrow(dat), 
                         steps = 0.025) {
  
  # take subset of data
  boot_dat_cluster <- dat[index, ]
  
  # expand to one row per effect size
  boot_dat <- tidyr::unnest(boot_dat_cluster, data)
  
  # fit selection model, return vector
  run_sel_model(boot_dat, steps = steps)
  
}

First, we take a subset of the data based on the index argument. This generates a bootstrap sample from the original data based on re-sampled clusters. We then use tidyr::unnest() to get the effect size level data for those re-sampled clusters. We then re-fit the model using our run_sel_model() function.

Now let’s apply our function to the Lehmann dataset. We will first need to create a nested version of the dataset, with one row per study:

lehmann_nested <- nest_by(lehmann_dat, study, .key = "data")

Now we can fit the selection model using fit_sel_model():

fit_selmodel(lehmann_nested)

## beta.intrcpt          tau        delta 
##    0.1327996    0.2848302    0.5484540

The results reproduce what we saw earlier when we estimated the three parameter selection model.

Now we can bootstrap using the boot::boot() function. The inputs to boot() are the nested dataset, the function to fit the selection model, and then any additional arguments passed to fit_selmodel()—here we include an argument for steps—and finally, the number of bootstrap replications:

# Generate bootstrap
set.seed(20230321)

tic()

boots <- boot(
  data = lehmann_nested,            # nested dataset
  statistic = fit_selmodel,         # function for fitting selection model
  steps = .025,                     # further arguments to the fitting function
  R = 1999                          # number of bootstraps
)

time_seq <- toc()

## 265.2 sec elapsed

This code takes a while to run, but we can speed it up with parallel processing.

Parallel processing

The boot package has some handy parallel processing features, but they can be a bit finicky to use here because of how R manages environments across multiple processes. With the above code, we can’t simply turn on parallel processing because the worker processes won’t know where to find the run_sel_model() function that gets called inside fit_selmodel(). To fix this, we include the run_sel_model() function inside our fitting function, as follows:

fit_selmodel <- function(dat, 
                         index = 1:nrow(dat), 
                         steps = 0.025) {
  
  # take subset of data
  boot_dat_cluster <- dat[index, ]
  
  # expand to one row per effect size
  boot_dat <- tidyr::unnest(boot_dat_cluster, data)
  
  # build run_selmodel
  run_sel_model <- function(dat, steps = .025) {
  
    # initial random effects model
    RE_mod <- metafor::rma.uni(
      yi = yi, vi = vi, data = dat, method = "ML"
    )
    
    # fit selection model
    res <- metafor::selmodel(
      RE_mod, type = "stepfun", steps = steps,
      skiphes = TRUE, # turn off SE calculation
      skiphet = TRUE # turn off heterogeneity test
    )
    
    # compile parameter estimates into a vector
    c(beta = res$beta[,1], tau = sqrt(res$tau2), delta = res$delta[-1])
    
  }
  
  p <- 2L + length(steps)  # calculate total number of model parameters
  
  # error handling for run_sel_model
  run_sel_model <- purrr::possibly(run_sel_model, otherwise = rep(NA_real_, p))
  
  # fit selection model, return vector of parameter estimates
  run_sel_model(boot_dat, steps = steps)
  
}

Now we can call boot() with options for parallel processing. The machine we used to compile this post has 12 cores. We will use half of the available cores for parallel processing. The configuration of parallel processing will depend on your operating system (different options are available for Mac), so you may need to adapt this code a bit.

ncpus <- parallel::detectCores() / 2

# Generate bootstrap
set.seed(20230321)

tic()

boots <- boot(
  data = lehmann_nested,
  statistic = fit_selmodel, steps = .025,
  R = 1999,
  parallel = "snow", ncpus = ncpus # parallel processing options
)

time_par <- toc()

## 63.08 sec elapsed

Parallel processing is really helpful here. We get 2000 bootstraps in 63 seconds, 4.2 times faster than sequential processing.

Standard Errors

The standard deviations of the bootstrapped parameter estimates can be interpreted as standard errors for the parameter estimates that take into account the dependence structure of the effect size estimates. Here is a table comparing the cluster-bootstrapped standard errors to the model-based standard errors generated by selmodel():

est <- boots$t0 # original parameter estimates

# calculate bootstrap SEs
boot_SE <- apply(boots$t, 2, sd, na.rm = TRUE)  

# calculate model-based SEs
model_SE <- with(RE_sel, c(se, se.tau2 / (2 * sqrt(tau2)), se.delta[-1]))

# make a table
res <- tibble(
  Parameter = names(est),
  Est = est,
  `SE(bootstrap)` = boot_SE,
  `SE(model)` = model_SE,
  `SE(bootstrap) / SE(model)` = boot_SE / model_SE
)

The standard errors based on cluster bootstrapping are all substantially larger than the model-based standard errors, which don’t account for dependence.

Confidence Intervals

For reporting results from this sort of analysis, it is useful to provide confidence intervals along with the model parameter estimates and standard errors. These can be calculated using the boot::boot_ci() function. This function provides several different types of confidence intervals; for illustration, we will stick with simple percentile confidence intervals, which are calculated by taking percentiles of the bootstrap distribution of each parameter estimate. To use the function, we’ll specify the type of confidence interval and the index of the parameter we want. An index of 1 is for the overall average effect size:

boot.ci(boots, type = "perc", index = 1) # For overall average ES

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1995 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boots, type = "perc", index = 1)
## 
## Intervals : 
## Level     Percentile     
## 95%   (-0.0038,  0.4171 )  
## Calculations and Intervals on Original Scale

Here is the confidence interval for between-study heterogeneity:

boot.ci(boots, type = "perc", index = 2) # For heterogeneity

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1995 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boots, type = "perc", index = 2)
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 0.001,  0.489 )  
## Calculations and Intervals on Original Scale

And for the selection weight:

boot.ci(boots, type = "perc", index = 3) # For selection weight

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1995 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = boots, type = "perc", index = 3)
## 
## Intervals : 
## Level     Percentile     
## 95%   ( 0.0574,  2.7536 )  
## Calculations and Intervals on Original Scale

Percentile confidence intervals can be asymmetric, and here the confidence interval for the selection weight parameter is notably asymmetric. The end-points of the confidence interval range from 0.057 (which represents very strong selective reporting, with only 6% of non-significant results reported) to 2.754 (which represents very strong selection against affirmative results).² So, overall we can’t really draw any strong conclusions about the strength of selective reporting.