[R-sig-ME] metafor: estimate correlation between response variables (meta-analysis)
Gustaf Granath
gustaf.granath at gmail.com
Wed Dec 14 15:56:21 CET 2016
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