I'm trying to do meta-regression with a lot of trials (>40 trials with >100 'arms') investigating the efficacy of a procedure (abl) and any 'addon' procedure. Each trial will have 2 or more arms. In some trials, the 'control' arm is no procedure; in others is a 'vanilla' procedure, with the active arm adding in various 'extras'.

I want to see not only if the procedure improves outcomes, but also if various 'extras' to the procedure also influence outcome.

My columns which are predictors in my model start with 'pred_'. If they are to do with the procedure used (abl), they start 'pred_abl_'

> head(af_dat)
        id arm  xi  ni pred_age pred_gender pred_antiarrhythmics pred_abl_pvi_done pred_abl_pvi_type pred_abl_cryo pred_abl_cfea pred_abl_gp pred_abl_lines pred_abl_mapping                      pred_endpoint pred_at_included pred_blanking        pi
1 23094720   A 112 146     56.0        68.0                    0                 1                 1             0             0           0              1                1                             Holter                0             3 0.7671233
2 23094720   B 103 148     54.0        72.0                    1                 0                 0             0             0           0              0                0                             Holter                0             3 0.6959459
3 21539635   A  45  58     57.6        66.0                    0                 1                 1             0             1           0              0                1                            Holter                 1             3 0.7758621
4 21539635   B  14  24     56.4        67.0                    0                 0                 0             0             1           0              0                1                             Holter                1             3 0.5833333
5 21539635   C  27  35     52.2        71.0                    0                 1                 1             0             0           0              0                1                             Holter                1             3 0.7714286
6 24549549   A  50  66     56.3        77.3                    0                 1                 1             0             0           0              0                0 Mobile telemetry (transtelephonic)                1             3 0.7575758

The column which indicates if they had a procedure at all is pred_abl_pvi_done. The other columns beginning pred_abl_* (of which there are 6 shown above) describe the presence of any 'extras' to the procedure, and must be 0 if the pred_abl_pvi_done is 0, by definition.

Initially, I was going to do the meta-regression as follows:

rma(measure='PLO', xi=xi, ni=ni, data=af_dat, mods=~pred_age+pred_gender+pred_antiarrhythmics+pred_abl_pvi_done+pred_abl_pvi_type+pred_abl_cryo+pred_abl_cfea+pred_abl_gp+pred_abl_lines+pred_abl_mapping+pred_endpoint+pred_at_included+pred_blanking)

However, I'm not sure if this is valid; the pred_abl_* columns will of course be very strongly correlated with pred_abl_pvi_done as they will be 0 if the latter is 0. If the 'abl' procedure (pred_abl_pvi_done=1) is very efficacious, I'm worried the other pred_abl_* columns will also appear correlated even if they add nothing prognostically, as many of the '0's for these columns will merely be indicating no procedure has been done at all.

One way to get around this, I suppose, would be to only look at the 6 procedural columns if pred_abl_pvi_done is 1.

I could therefore do 2 analysies:

  1. Multivariate regression of all the studies, excluding the 6 procedural columns. I would then see if a procedure versus no procedure (pred_abl_pvi_done 1 or 0) is significant
  2. Select only cases where pvi_abl_pvi_done = 1 (ie the procedure is done (90% of the trials), and exclude medical treatment), and THEN do a further multivariate regression where I include all the procedural questions.

Do you think I need to do these 2 analyses separately due to the correlation?

  • 1
    $\begingroup$ Not sure about the answer to your questions, but do you not need to use rma.mv to take account of the multiple arms per study? You do refer to it as a multivariate analysis. $\endgroup$
    – mdewey
    May 22, 2016 at 20:42

1 Answer 1


First, let me answer the specific concern you have. It should be fine to use all of the data and include pred_abl_pvi_done and the variables that indicate the "extras" in the same model. Yes, pred_abl_pvi_done and the pred_abl_* variables will be highly correlated, but the coefficients you obtain indicate the "incremental" value of adding the extras on top of the procedure. If the extras don't add anything, the coefficients for the pred_abl_* variables will hover around 0 and any significant results will be due to chance.

Try out the following code:


k <- 100

vi <- runif(k, .05, 1)
trt <- c(rep(0,k/2), rep(1,k/2))
ext <- sample(c(0,1), k, replace=TRUE) * trt
yi <- 0.5 * trt + rnorm(k, 0, .25) + rnorm(k, 0, sqrt(vi))

rma(yi, vi, mods = ~ trt + ext)

Variable trt indicates whether the treatment was performed or not and ext whether some kind of extra was added. The two variables are highly correlated. But in the true model, only trt exerts an influence on the observed outcomes (plus there is some heterogeneity and of course sampling error). If you run the code many times, you will see that trt is usually correctly identified as a significant predictor while ext is not.

Of course, if ext adds something that is not already accounted for by trt, then this should also start to show up in the results.

So, that part is fine as far as I am concerned. However, I am a bit more concerned about the fact that you are not accounting for the fact that the arm-level outcomes have a multilevel structure (@mdewey also hinted at this in the comments). Right now, you are treating the data essentially cross-sectionally, when I think you should really be looking at the influence of your predictors within trials.

To do so, you can either add factor(id) as a predictor (so, you add fixed study effects) or you could use random study effects with rma.mv(..., random = ~ 1 | id/arm). There is a bit of a debate in the literature which of the two is more appropriate, but you could explore (and report) both approaches and hopefully the conclusions are consistent.

  • $\begingroup$ Any pointers to the literature to which you refer in last para? $\endgroup$
    – mdewey
    May 23, 2016 at 20:59
  • $\begingroup$ You could start here: dij.sagepub.com/content/34/2/535.abstract $\endgroup$
    – Wolfgang
    May 23, 2016 at 21:41
  • 1
    $\begingroup$ Many thanks! I've posted the results here: pastebin.com/8U67tKGU - I didn't seem to be able to get the second method working as you can see (even with "id/arm" as you wrote). Is it a fundamentally different function? It's not clear how to specify I want to use the proportional measure as I did with PLO in rma()... $\endgroup$
    – James
    May 24, 2016 at 15:14
  • 1
    $\begingroup$ You want to account for study- and arm-level heterogeneity. $\endgroup$
    – Wolfgang
    May 25, 2016 at 11:33
  • 1
    $\begingroup$ The results you posted say otherwise - the arm-level variance component is (at least) as large as the study-level one. $\endgroup$
    – Wolfgang
    May 25, 2016 at 12:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.