[R-meta] Network meta-analysis with randomized and non-randomized studies in metafor

[Chris Rose]

Hi Michael, all

I think I had actually tried your suggestion of constructing a matrix but did not include it my attempt in my original email. After exhausting a few alternatives that seemed sensible based on my reading of the docs, I started to try some “random stuff" in the hope that I would learn something purely by experimentation — unfortunately not.

I’ve included below an example from the beginning, which includes what I think you mean towards the end. It’d be great if you could have a look and tell me where I’m going wrong.

Best,

— Chris


### load the library
library(metafor)

### copy data into 'dat'
dat <- dat.hasselblad1998

### calculate log odds for each study arm
dat <- escalc(measure="PLO", xi=xi, ni=ni, add=1/2, to="all", data=dat)

### convert trt variable to factor with desired ordering of levels
dat$trt <- factor(dat$trt, levels=c("no_contact", "self_help", "ind_counseling", "grp_counseling"))

### add a space before each level (this makes the output a bit more legible)
levels(dat$trt) <- paste0(" ", levels(dat$trt))

### simulate a RS/NRS moderator (type)
set.seed(1234)
dat$type <- rep(sample(c("RS", "NRS"), length(unique(dat$study)), replace=TRUE), times=tapply(dat$study, dat$study, length))

### fit a meta-regression model
res <- rma.mv(yi, vi, mods = ~ trt + type, random = list(~ type | study, ~ type | id), struct=c("DIAG","DIAG"), data=dat)

### call predict without additional arguments
predict(res) # This works!

### make a matrix of moderator values at which to predict, and try to predict
newvals <- as.matrix(c(trt = " self_help", type = "RS"))
predict(res, newmods = newvals) # fails with...
# Error in predict.rma(res, newmods = newvals) : 
#   Dimensions of 'newmods' do not match dimensions of the model.

### try to specify the tau and gamma levels
predict(res, newmods = newvals, tau2.levels = "RS", gamma2.levels = "RS") # fails with...
# Error in predict.rma(res, newmods = newvals, tau2.levels = "RS", gamma2.levels = "RS") : 
#   Dimensions of 'newmods' do not match dimensions of the model.

### is the matrix moderators incorrectly transposed?
predict(res, newmods = t(newvals), tau2.levels = "RS", gamma2.levels = "RS") # fails with...
# Error in dimnames(x) <- dn : 
#   length of 'dimnames' [2] not equal to array extent


> On 14 Feb 2019, at 18:05, Michael Dewey <lists using dewey.myzen.co.uk> wrote:
> 
> Dear Chris
> 
> I suspect the problem may be that when it converts the character vector which you are supplying as newmods it makes it a 1-column matrix rather than a 1-row matrix. What happens if you explicitly create a matrix?
> 
> Michael
> 
> On 13/02/2019 10:28, Chris Rose wrote:
>> Hi Wolfgang
>> A couple of months ago you provided some help with a network meta-analysis that adjusts for possible differences between randomised and non-randomised studies (for which, thanks again). I’m now at the stage of reporting results, and would like to provide prediction intervals. I’m aware of the predict.rma function, but I’m struggling to get it to work!
>> Taking an example from your original reply (based on the Hasselblad 1998 data, to which you added a study type factor [called type] to identify randomised [“RS”] vs nonrandomised [“NRS”] studies), one of the models you suggested was:
>> res <- rma.mv(yi, vi, mods = ~ trt + type, random = list(~ type | study, ~ type | id), struct=c("DIAG","DIAG"), data=dat)
>> I can run predict(res) and see predictions for each arm of each study (i.e., each row of the original data frame, dat). This works fine.
>> However, I want to make predictions for specific values of the trt and type variables. I have tried:
>> predict(res, newmods = c(trt = "trt self_help", type = "RS”))
>> and
>> predict(res, newmods = c(trt = "trt self_help", type = "RS"), tau2.levels = "RS", gamma2.levels = "RS”)
>> and
>> predict(res, newmods = c(trt = " self_help", type = "RS"), tau2.levels = "RS", gamma2.levels = "RS”)
>> # Note the space at the start of " self_help”.
>> and
>> predict(res, newmods = res$X) # Here I’m trying to use the same dummy coding as rma.mv.
>> But these all give the error:
>> Error in dimnames(x) <- dn :
>>   length of 'dimnames' [2] not equal to array extent
>> I have read the documentation for the predict.rma function, but I must be misunderstanding a key piece of information (or missing something very obvious and stupid!). I am aware that newmods should be a matrix, and have tried various versions of the above by constructing matrices. I have also searched the web and this list, but have not found something that has helped.
>> Thanks very much in advance,
>> — Chris
>>> On 19 Dec 2018, at 10:42, Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer using maastrichtuniversity.nl> wrote:
>>> 
>>> Hi Chris,
>>> 
>>> Why don't you just include a moderator that distinguishes RS from NRS? The model should already include fixed/random effects for studies and random effects for estimates, so adding a moderator is sufficient. Here is an example, using the data from help(dat.hasselblad1998):
>>> 
>>> ### copy data into 'dat'
>>> dat <- dat.hasselblad1998
>>> 
>>> ### calculate log odds for each study arm
>>> dat <- escalc(measure="PLO", xi=xi, ni=ni, add=1/2, to="all", data=dat)
>>> dat
>>> 
>>> ### convert trt variable to factor with desired ordering of levels
>>> dat$trt <- factor(dat$trt, levels=c("no_contact", "self_help", "ind_counseling", "grp_counseling"))
>>> 
>>> ### add a space before each level (this makes the output a bit more legible)
>>> levels(dat$trt) <- paste0(" ", levels(dat$trt))
>>> 
>>> ### network meta-analysis using an arm-based model
>>> res <- rma.mv(yi, vi, mods = ~ trt, random = ~ 1 | study/id, data=dat)
>>> res
>>> 
>>> ### simulate a RS/NRS moderator (type)
>>> set.seed(1234)
>>> dat$type <- rep(sample(c("RS", "NRS"), length(unique(dat$study)), replace=TRUE), times=tapply(dat$study, dat$study, length))
>>> 
>>> ### adding type (RS/NRS) to the model
>>> res <- rma.mv(yi, vi, mods = ~ trt + type, random = ~ 1 | study/id, data=dat)
>>> res
>>> 
>>> This last model assumes that the amount of heterogeneity at the study and estimate level is the same for RS and NRS. This might not be true; we may suspect that NRS are more heterogeneous. We can allow for this:
>>> 
>>> ### allowing for different amounts of heterogeneity for RS vs NRS
>>> res <- rma.mv(yi, vi, mods = ~ trt + type, random = list(~ type | study, ~ type | id), struct=c("DIAG","DIAG"), data=dat)
>>> res
>>> 
>>> This is already a pretty complex model, so let's check on the identifiability of the variance components:
>>> 
>>> profile(res)
>>> 
>>> Looks good -- all profile plots have a clear peak (one being at 0, but that is fine).
>>> 
>>> Finally, one may not just want to distinguish between RS/NRS, but actually allow the effect of the treatment to differ depending on the type of study. Then we allow type and trt to interact:
>>> 
>>> res <- rma.mv(yi, vi, mods = ~ type * trt, random = list(~ type | study, ~ type | id), struct=c("DIAG","DIAG"), data=dat)
>>> res
>>> 
>>> I hope this gives you some ideas on how to model your data.
>>> 
>>> Best,
>>> Wolfgang
>>> 
>>> -----Original Message-----
>>> From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-bounces using r-project.org] On Behalf Of Chris Rose
>>> Sent: Wednesday, 21 November, 2018 10:25
>>> To: r-sig-meta-analysis using r-project.org
>>> Subject: [R-meta] Network meta-analysis with randomized and non-randomized studies in metafor
>>> 
>>> Hi
>>> 
>>> I'm working on a network meta-analysis (NMA) that will include randomized studies (RS) and non-randomized studies (NRS). I anticipate that the NRS may be biased, and that magnitude and direction of bias may vary unpredictably across the NRS. I am less concerned that the NRS may be overly precise (though I may be wrong!). It is also possible that I will have to include fixed effects for moderators (i.e., do a network meta-regression). I'm considering using the metafor package (see the end of this email for my reasoning).
>>> 
>>> How should I model the distinction between the RS and NRS in metafor? I’m considering using a random effect (intercept) for each NRS. I.e., each effect estimate would be supported by some mixture of direct and indirect evidence that is (hopefully) unbiased in the case of the RS, and potentially biased in the case of the NRS (with that bias explained by the random effects).
>>> 
>>> How would I specify this model in metafor? I.e., how do I specify that the random effects are conditional on a study being a NRS? Could I use a factor variable that has the same value for all the randomized studies, and distinct values for the non-randomized studies? Is there a better way to do this?
>>> 
>>> Thanks in advance,
>>> 
>>> Chris
>>> 
>>> Why metafor?
>>> 
>>> I’m aware of various NMA packages available for R, and I'm able to implement bespoke Bayesian models in JAGS and Stan. I’m considering metafor because:
>>> 
>>> 1. I need to use R.
>>> 2. I have multiple outcomes to analyze and need to do various sensitivity analyses, so run-time is a consideration and MCMC is less attractive for this reason.
>>> 3. I'd prefer to use a well-tested model implementation (I'm experienced enough with JAGS and Stan to know it's possible to make programming errors).
>>> 4. I think the existing NMA-specific packages do not support modeling the distinction between RS and NRS (but I could be wrong).
>>> 5. I think the existing NMA packages either do not support regression, or only support a single covariate or continuous covariates. I think this means they do not allow me to adjust for the randomized/non-randomized designs as well as other covariates if that was necessary. Again, I could be wrong.
>> 	[[alternative HTML version deleted]]
>> _______________________________________________
>> R-sig-meta-analysis mailing list
>> R-sig-meta-analysis using r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
> 
> -- 
> Michael
> http://www.dewey.myzen.co.uk/home.html


	[[alternative HTML version deleted]]