[R-sig-ME] metafor: estimate correlation between response variables (meta-analysis)
Viechtbauer Wolfgang (SP)
wolfgang.viechtbauer at maastrichtuniversity.nl
Thu Dec 15 09:50:17 CET 2016
After you fit the "res.with.covar" model, you have:
print(res.no.covar, digits=3)
which should be:
print(res.with.covar, digits=3)
which shows a correlation of -0.967.
I haven't really delved into the rest of the code, so I cannot say much about that.
Best,
Wolfgang
> -----Original Message-----
> From: Gustaf Granath [mailto:gustaf.granath at gmail.com]
> Sent: Wednesday, December 14, 2016 15:56
> To: r-sig-mixed-models at r-project.org
> Cc: Viechtbauer Wolfgang (SP)
> Subject: Re: metafor: estimate correlation between response variables
> (meta-analysis)
>
> I tried to add between response covariances with each study, here using
> the 'cor*sqrt(v1*v2)' approach. I think I got it right but the estimated
> correlation coef (between responses across studies) dont change between
> a model without correlation within studies, compared to a model with
> this correlation. Not sure if it is just the way I set up my toy
> example, or if Im doing something wrong.
>
> require(metafor)
>
> # make data
> set.seed(1)
> dat <- data.frame(exp = gl(8,4), resp = gl(2,2, length=32), yi =
> rnorm(32,10,4),
> vi = rnorm(32,4,0.5), vi.cov = rnorm(32,0.5,0.25) )
> dat$vi.cov <- rep(c(0.1,0.15),8, each=2)
> dat$st.depend <- interaction(dat$exp, dat$resp)
>
> # known var-cov matrix (within-study dependence for each response)
> # load function 'v_func' first (see below)
> V = v_func(dat = dat, study.dep = "st.depend", var.y = "vi" , var.y.cov
> = "vi.cov")
>
> # run model
> res.no.covar <- rma.mv(yi, V, mods = ~ resp - 1,
> random = ~ resp | exp,
> struct="UN", data=dat, method="REML")
> print(res.no.covar, digits=3)
> # correlation = -0.613
>
> # test with correlation between responses within studies
> cor.r = 0.5
>
> covaris = list()
> for (i in 1:nlevels(dat$exp)) {
> studid <- dat$exp == levels(dat$exp)[i]
> vari <- dat[ studid, "vi"]
>
> # make cor matrix
> a <- matrix(nrow=length(vari), ncol=length(vari))
> diag(a) <- 1
> a[lower.tri(a)] <- cor.r
> a[upper.tri(a)] <- cor.r
>
> # calculate covariance for the given cor coef
> b <- sqrt(vari) %*% t(sqrt(vari))
> a_cov <- b * a # covariance matrix
> # cov2cor(a_cov) #test
>
> # store var-cov block for each study
> covaris[[i]] <- a_cov
> }
>
> # var-cov matrix for the between response correlation
> # within each study
> V2 <- as.matrix(bdiag(lapply(covaris, function (x) x) ))
> # cov2cor(V2) #test if cor is correct
>
> # add correlation between responses to
> # within-response matrix
> V = ifelse(V==0, V2, V)
>
> # run mode with complete var-covar matrix
> res.with.covar <- rma.mv(yi, V, mods = ~ resp - 1,
> random = ~ resp | exp,
> struct="UN", data=dat, method="REML")
> print(res.no.covar, digits=3)
> # correlation still -0.613....
>
> # function to create var-covariance matrix
>
> v_func <- function (dat, study.dep, var.y , var.y.cov) {
>
> dat.agg <- data.frame(st = unique(dat[[study.dep]]), len =
> rle(as.numeric(dat[[study.dep]]))$lengths,
> var = unique(dat[[var.y.cov]]))
> require(Matrix)
> ll = list()
> for (i in 1:NROW(dat.agg)) {
> ll[[i]] <- matrix(dat.agg$var[i], ncol = dat.agg$len[i], nrow =
> dat.agg$len[i])
> }
> mat <- as.matrix(bdiag(lapply(ll, function (x) x) ))
> diag(mat) <- dat[[var.y]]
> return(mat)
> }
>
> On 2016-12-14 11:14, Viechtbauer Wolfgang (SP) wrote:
> > So there really is a 4x4 var-cov matrix for each study then. So, like
> this:
> >
> > trt1-var1 trt1-var2 trt2-var1 trt2-var2
> > trt1-var1 v_11 cov_11,12 cov_11,21 cov_11,22
> > trt1-var2 v_12 cov_12,21 cov_12,22
> > trt2-var1 v_21 cov_21,22
> > trt2-var2 v_22
> >
> > (leaving out subscript i -- but there is such a 4x4 block for each
> study).
> >
> > And the covariance for different treatments (for the same variable?) is
> known. So that would be cov_11,21 and cov_12,22. That leaves 4 unknown
> covariances. So yes, one could impute those (by making some reasonable
> assumptions about the correlations and then back-computing the
> covariances using appropriate equations) followed by a sensitivity
> analysis.
> >
> > I agree that implementing something like this can be a bit tricky, but
> can be done.
> >
> > Best,
> > Wolfgang
> >
> >> -----Original Message-----
> >> From: Gustaf Granath [mailto:gustaf.granath at gmail.com]
> >> Sent: Tuesday, December 13, 2016 15:55
> >> To: r-sig-mixed-models at r-project.org
> >> Cc: Viechtbauer Wolfgang (SP)
> >> Subject: Re: metafor: estimate correlation between response variables
> >> (meta-analysis)
> >>
> >> Hi,
> >>
> >> Im sorry if my toy example wasnt clear. You almost got it right
> though.
> >> But yi is not independent for different treatments within a study. I
> >> have the estimated covariance though, so no problem to account for
> this.
> >> In my example I did include this covariance when I constructed V (I
> used
> >> the outcome covariance in the berkey1998 data to illustrate this). If
> I
> >> understand you correctly, the best way check the influence of
> >> within-study correlation between responses, is to manually add
> >> covariances (try a range of reasonable values) between responses in V.
> I
> >> will try that - although, Im not sure that it is easy to implement.
> >>
> >> Thanks,
> >>
> >> Gustaf
> >>
> >> On 2016-12-13 13:41, Viechtbauer Wolfgang (SP) wrote:
> >>> I have a hard time making sense of that toy example.
> >>>
> >>> But if I understand you correctly, you have data like this:
> >>>
> >>> study trt respvar yi vi
> >>> -----------------------
> >>> 1 1 1 . .
> >>> 1 1 2 . .
> >>> 1 2 1 . .
> >>> 1 2 2 . .
> >>> 2 1 1 . .
> >>> 2 1 2 . .
> >>> 2 2 1 . .
> >>> 2 2 2 . .
> >>> ...
> >>>
> >>> where 'yi' and 'vi' are the observed outcomes and corresponding
> >> variances.
> >>> Within studies, can we assume that 'yi' is independent for different
> >> treatments? Then V (the var-cov matrix of the 'yi' vector) will be
> block-
> >> diagonal, each block being a 2x2 matrix (since you only need to
> consider
> >> the covariance between 'yi' for respvar 1 and 'yi' for respvar 2). And
> >> the problem is that the covariances are unknown. Correct so far?
> >>> If so, the 'R' argument has nothing to do with this. If you want to
> >> approach this by means of a sensitivity analysis, then just impute the
> >> unknown covariances directly into 'V' and analyze by means of an
> >> appropriate multilevel/multivariate model. Something like this:
> >>> dat$study.trt <- interaction(dat$study, dat$trt)
> >>> rma.mv(yi, V, mods = ~ factor(trt) - 1, random = list(~ 1 | study, ~
> >> respvar | study.trt), struct="UN", data=dat)
> >>> To impute the covariances, you may be able to use cor*sqrt(v1*v2),
> >> where 'cor' is some kind of assumed correlation (maybe constant across
> >> studies, maybe not). However, whether this is appropriate depends on
> your
> >> outcome measure. For example, this would be fine for means or mean
> >> differences, but the covariance between standardized mean differences
> >> cannot be computed that way (see Gleser & Olkin, 2009, for the correct
> >> equation).
> >>> Gleser, L. J., & Olkin, I. (2009). Stochastically dependent effect
> >> sizes. In H. Cooper, L. V. Hedges, & J. C. Valentine (Eds.), The
> handbook
> >> of research synthesis and meta-analysis (2nd ed., pp. 357-376). New
> York:
> >> Russell Sage Foundation.
> >>> Best,
> >>> Wolfgang
> >>>
> >>>> -----Original Message-----
> >>>> From: R-sig-mixed-models [mailto:r-sig-mixed-models-bounces at r-
> >>>> project.org] On Behalf Of Gustaf Granath
> >>>> Sent: Tuesday, December 13, 2016 11:42
> >>>> To: r-sig-mixed-models at r-project.org
> >>>> Subject: [R-sig-ME] metafor: estimate correlation between response
> >>>> variables (meta-analysis)
> >>>>
> >>>> Hi,
> >>>> Im trying the estimate the correlation between two response variable
> >> in
> >>>> a meta-analysis. Basically, are effects on y1 associated with
> effects
> >> on
> >>>> y2 across studies. Unfortunately, I dont have information on y1-y2
> >>>> correlations within each study. In addition, each study contain
> >> multiple
> >>>> treatments, adding within-study dependence for each response
> variable.
> >>>>
> >>>> Because I dont have the y1-y2 covariance in each study, my idea is
> to
> >>>> run analyses with different covariance/correlation values to explore
> >> how
> >>>> this covariance affect the result. Reading the rma.mv()
> documentation,
> >>>> this seems possible through the R= argument, but I cant figure out
> >> how.
> >>>> Alternatively, I can add covariances in the matrix describing the
> >> known
> >>>> var-covariance matrix of the within-study dependence, but it seemed
> >> like
> >>>> the R argument is a easier solution (I may be wrong though). Below
> >>>> follow code illustrating my problem using the dat.berkey1998
> example.
> >>>>
> >>>> Cheers
> >>>>
> >>>> Gustaf
> >>>>
> >>>> # Estimate correlation between two response variables
> >>>> # with within-study dependence
> >>>> library(metafor)
> >>>>
> >>>> # Make up an example based on the berkey1998 data.
> >>>> # Multiple outcomes (multiple treatments)
> >>>> # within each study, for each response variable, are created.
> >>>> # Covariances between response variables (within studies) are
> unknown.
> >>>> dat <- get(data(dat.berkey1998))
> >>>> st.vcov <- dat[, c("v1i", "v2i")] # save vcov matrices to add later
> >>>> dat <- rbind(dat, dat)
> >>>> dat <- dat[order(dat$trial, dat$outcome),] # fix order
> >>>> dat[,c("v1i", "v2i")] <- rbind(st.vcov,st.vcov) # put back matrices
> >>>> dat$trial.treat <-paste(rep(1:2, nrow(dat)/2), dat$trial, sep="_") #
> >> id
> >>>> for treatments within trials
> >>>>
> >>>> # Covariances within studies, for each response variable, is known.
> >>>> Hence, a
> >>>> # varcov-matrix for each study and response, can be made and added
> as
> >>>> blocks
> >>>> # in a large varcov-matrix for the data set.
> >>>> # First a within-study dependence dummy must be added
> >>>> dat$stud.unit <- interaction(dat$trial, dat$outcome) # within-trial
> >>>> dependence dummy
> >>>> # Put together known varcov matrix
> >>>> V <- bldiag(lapply(split(dat[,c("v1i", "v2i")], dat$stud.unit),
> >>>> as.matrix))
> >>>>
> >>>> # plot relationship between response variables
> >>>> dat.wide <- reshape(dat, direction="wide", v.names = "yi", timevar
> =
> >>>> "outcome", idvar = "trial.treat")
> >>>> plot(yi.AL ~yi.PD, dat.wide)
> >>>> cor(dat.wide$yi.AL, dat.wide$yi.PD)
> >>>> # r = 0.39, if using the study outcomes ignoring all
> >>>> dependence/uncertainty.
> >>>>
> >>>> # Run analysis
> >>>> res <- rma.mv(yi, V, mods = ~ outcome - 1,
> >>>> random = ~ outcome | trial,
> >>>> struct="UN", data=dat, method="REML")
> >>>> print(res, digits=3)
> >>>> # correlation = 0.51, when accounting for dependence etc. But, y1
> and
> >> y2
> >>>> are assumed to be independent
> >>>> # within each study. How to perform sensitivity analyses by adding
> >>>> different values
> >>>> # on this within-study correlation between the two response
> variables?
> >>>>
> >>>> --
> >>>> Gustaf Granath (PhD)
> >>>> Post doc
> >>>> Swedish University of Agricultural Sciences
> >>>> Department of Ecology
More information about the R-sig-mixed-models
mailing list