[R-sig-ME] metafor: estimate correlation between response variables (meta-analysis)
Gustaf Granath
gustaf.granath at gmail.com
Thu Dec 15 10:03:01 CET 2016
Oh, that is embarrassing. Thanks for spotting my mistake. It also
explains why I got more sensible results when I tried it on my real data.
Short additional question. Is the plan to implement BLUPs for rma.mv?
Cheers
Gustaf
On 2016-12-15 09:50, Viechtbauer Wolfgang (SP) wrote:
> 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