[R-meta] Between study correlation between two different outcomes

[James Pustejovsky]

The correlation between random effects is a notoriously difficult thing to
estimate well. Here are two citations that discuss the point:

* Riley, R. D., Abrams, K. R., Sutton, A. J., Lambert, P. C., & Thompson,
J. R. (2007). Bivariate random-effects meta-analysis and the estimation of
between-study correlation. BMC Medical Research Methodology, 7(3), 1–15.
https://doi.org/10.1186/1471-2288-7-3

* Chen, Y., Hong, C., & Riley, R. D. (2014). An alternative
pseudolikelihood method for multivariate random-effects meta-analysis.
Statistics in Medicine, (October), n/a-n/a. https://doi.org/10.1002/sim.6350

Two things you could do to get a better sense of how well rho can be
estimated:

1. Plot the profile log-likelihood for rho (example code below). The REML
estimate of rho is based on maximizing the restricted log likelihood. If
the profile likelihood is quite flat, then there's just not much
information available in the data to estimate this correlation.

library(metafor)
> dat <- get(data(dat.berkey1998))
> cor <- 0.7
> V <- bldiag(lapply(split(dat[,c("v1i", "v2i")], dat$trial), function(vi)
> {vi[1,2] <- vi[2,1] <- sqrt(vi[1,1] * vi[2,2]) * cor; as.matrix(vi)}))
> ### multiple outcomes random-effects model
> res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial,
> struct="UN", data=dat)
> profile(res, rho = 1)


2. Calculate the number of studies in your meta-analysis that include
estimates of both outcomes. If it's small, then rho will not be well
estimated.

There are a few potential remedies:
- Get more data that includes effect size estimates for both outcomes.
- Sophia Rabe-Hasketh and some of her colleagues have proposed using
penalized log likelihood to estimate variance component in multi-level
models, using priors that keep the random effect variance structure from
degenerating. The blme R package will fit those models easily, although I'm
not sure if it can do meta-analysis models.
- If you can develop a reasonable prior for rho you could go full Bayes.
The brms package is one good tool for fitting such models.


On Wed, May 16, 2018 at 4:29 AM, Aulbach, Matthias B <
matthias.aulbach at helsinki.fi> wrote:

> Hello once again,
>
> thanks to Wolfgang's most valuable replies, I have run the analyses as
> suggested. The problem occurring now is that the estimates for rho don't
> seem to make a lot of sense (point estimates of 1 and/or huge confidence
> intervals). I found this post from James E. Pustejovsky about it in the
> archives of this mailing list [http://jepusto.github.io/
> imputing-covariance-matrices-for-multi-variate-meta-analysis] where I
> read this:
>
> "The metafor fit is also a bit goofy because the correlation between the
> random effects for math and verbal scores is very close to -1, although
> evidently it is not uncommon to obtain such degenerate estimates of the
> random effects structure."
>
> How is it that the estimation of rho is "goofy" and could someone point to
> any other examples on the literature to underline that "is not uncommon"?
> Does that imply that the method doesn't work properly here? Or did I
> probably do something wrong?
>
> To clarify: What I want to know is to what degree studies that find larger
> effects in outcome 1 also find larger effects in outcome 2.  In my initial
> solution to this problem, I had conducted univariate analyses and then ran
> a meta-regression with the outcome 1 as dependent variable and outcome 2 as
> a moderator. This was criticized by reviewers which is why I now conducted
> the multivariate analysis with both outcomes, trying to determine rho.
>
> Thanks a lot in advance (once more)!
>
> Best,
> Matthias
>
>
> -----Original Message-----
> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces at r-project.org> On
> Behalf Of Aulbach, Matthias B
> Sent: torstai 10. toukokuuta 2018 9.16
> To: Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer@
> maastrichtuniversity.nl>
> Cc: r-sig-meta-analysis at r-project.org
> Subject: Re: [R-meta] Between study correlation between two different
> outcomes
>
> Dear Wolfgang,
>
> thank you again for your really detailed answer. This is helping me a lot!
> I will definitely check the mailing list archives.
>
> All the best,
> Matthias
>
> -----Original Message-----
> From: Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer@
> maastrichtuniversity.nl>
> Sent: keskiviikko 9. toukokuuta 2018 11.46
> To: Aulbach, Matthias B <matthias.aulbach at helsinki.fi>
> Cc: r-sig-meta-analysis at r-project.org
> Subject: RE: Between study correlation between two different outcomes
>
> Dear Matthias,
>
> Please keep the list in CC when replying.
>
> V should have the sampling *variances* along the diagonal. And yes, the
> off-diagonal values will be the covariances. Since you have two values per
> study, V will therefore be block-diagonal with 2x2 blocks along the
> diagonal. So your case is basically the same as this one:
>
> http://www.metafor-project.org/doku.php/analyses:berkey1998
>
> The issue of missing information about the covariances has been discussed
> on this mailing list quite extensively, so it would be good to browse
> through the archives.
>
> One solution you will find mentioned there is ignoring the covariances
> (i.e., assuming that they are 0) and then using cluster-robust inference
> methods. This approach should be ok when interest is only on the fixed
> effects of the model. In your case, this approach is *NOT* appropriate,
> since you are specifically interested in the correlation between the random
> effects. If you assume that the covariances between the sampling errors is
> 0, then this will (usually) drive up the covariance between the random
> effects, which is going to lead to an overestimate of the correlation
> between the random effects. The cluster-robust approach fixes up the
> standard errors of the fixed effects, but won't do anything to correct for
> the bias in the estimated correlation between the random effects.
>
> So, if you do not have any information to compute the covariances, you
> will have to make an educated guess and then do a sensitivity analysis.
>
> Using the Berkey et al. (1998) data as an example, here is an illustration
> of what this might look like:
>
> library(metafor)
>
> ### load data
> dat <- get(data(dat.berkey1998))
>
> cors <- seq(0, .99, length=100)
> rhos <- rep(NA, length(cors))
>
> for (j in 1:length(cors)) {
>
>    print(j)
>
>    ### construct V matrix with guestimated correlation
>    V <- bldiag(lapply(split(dat[,c("v1i", "v2i")], dat$trial),
> function(vi) {vi[1,2] <- vi[2,1] <- sqrt(vi[1,1] * vi[2,2]) * cors[j];
> as.matrix(vi)}))
>
>    ### multiple outcomes random-effects model
>    res <- rma.mv(yi, V, mods = ~ outcome - 1, random = ~ outcome | trial,
> struct="UN", data=dat)
>    rhos[j] <- res$rho
>
> }
>
> plot(cors, rhos, type="o", pch=19)
>
> So what I am doing here is constructing the V matrix based on the sampling
> variances and different assumed values for the correlation between the
> sampling errors (cors). Then I fit the model and save the correlation
> between the random effects (rhos). As you can see in the plot, as I assume
> a higher value for cor, the value for rho goes down. Assuming a correlation
> of 0 leads to an estimate of rho equal to 0.78. On the other hand, assuming
> a correlation of 0.9 leads to an estiamte of rho equal to 0.36. That's
> quite different.
>
> I don't know to what extent rho will depend on the assumed correlation for
> your data. Also, I used a very wide range for 'cors' just for illustration
> purposes. In practice, you should be able to narrow down the range a bit
> more.
>
> Best,
> Wolfgang
>
> -----Original Message-----
> From: Aulbach, Matthias B [mailto:matthias.aulbach at helsinki.fi]
> Sent: Wednesday, 09 May, 2018 10:07
> To: Viechtbauer, Wolfgang (SP)
> Subject: RE: Between study correlation between two different outcomes
>
> Dear Wolfgang,
>
> thank you very much for your fast and very helpful answer!
>
> About the data structure: basically everything you assumed is correct. I
> have two rows for each study ("ID") ,first row outcome 1, second row
> outcome 2. "type_outcome" is a binary variable defining which kind of
> outcome is entered in which line.
>
> And yes, the measurements are from the same subjects. The problem is that
> I don't have data about within-study correlations (or covariance), so I
> will have to make a guess about how large that is (or try out different
> values?). Is that right?
> So as I understand, I need a 32x32-matrix with the diagonal being the
> standard errors from the studies and the off-diagonal values being my
> educated guess about the within-study covariance between the outcomes. Is
> that right?
>
> Thank you once more for the help!
>
> Best,
> Matthias
>
> -----Original Message-----
> From: Viechtbauer, Wolfgang (SP) <wolfgang.viechtbauer@
> maastrichtuniversity.nl>
> Sent: torstai 3. toukokuuta 2018 20.45
> To: Aulbach, Matthias B <matthias.aulbach at helsinki.fi>;
> r-sig-meta-analysis at r-project.org
> Subject: RE: Between study correlation between two different outcomes
>
> Dear Matthias,
>
> Can you explain the data structure a bit more? I assume you have two rows
> for each level of 'ID', the first row for outcome 1 and the second row for
> outcome 2 and that 'type_outcome' is a dummy variable to indicate the
> outcome. Is that correct?
>
> Then the output should include the correlation between the underlying true
> effects (rho). You can get a CI for this with the confint() function. To
> test its significance, you can conduct a likelihood ratio test. This should
> do it:
>
> rp <- rma.mv(yi, vi, data=pd, mod = ~ factor(type_outcome) -1, random= ~
> type_outcome | ID)
> rp0 <- rma.mv(yi, vi, data=pd, mod = ~ factor(type_outcome) -1, random= ~
> type_outcome | ID, rho=0) anova(rp, rp0) confint(rp, rho=1)
>
> Two notes:
>
> The model you are fitting uses struct="CS" by default. This assumes that
> the amount of heterogeneity is the same for the two outcomes, which may not
> be appropriate. So you might want to use struct="UN".
>
> Also, I assume the two outcomes are measured in the same subjects. In that
> case, the sampling errors of the two outcomes are correlated. So the V
> matrix (the second argument for the rma.mv()) function is not diagonal,
> but also should include the covariances. If you do not account for this,
> the correlation between the underlying true effects is very likely an
> overestimate.
>
> Best,
> Wolfgang
>
> -----Original Message-----
> From: R-sig-meta-analysis [mailto:r-sig-meta-analysis-
> bounces at r-project.org] On Behalf Of Aulbach, Matthias B
> Sent: Thursday, 03 May, 2018 16:52
> To: r-sig-meta-analysis at r-project.org
> Subject: [R-meta] Between study correlation between two different outcomes
>
> Hi,
>
> I am conducting a meta-analysis, using the great metafor package. I have
> run into the problem of dependency when using more than one outcome from
> the same set of studies. In an earlier attempt, I had run two univariate
> meta-analyses and then used the outcome of one outcome as a predictor in a
> meta-regression with the other outcome as the dependent variable. But that
> ignores the within-study correlation between the two, so I'd like to
> improve that. So I'd like to handle that using a multivariate
> meta-analysis, using rma.mv and this line of code (with yi denoting the
> effect sizes, vi their standard errors, type_outcome the kind of outcome
> that was measured, and ID as the study identifier):
>
> rp <- rma.mv(yi, vi, data=pd, mod = ~ factor(type_outcome) -1, random= ~
> type_outcome | ID)
>
> This nicely gives me the different effects for the two kinds of outcomes.
> However, what I am so desperately interested in is the between-study
> correlation between the two outcomes, i.e. if there's a strong effect on
> one outcome, is there also a strong effect on the other (or not)? Is there
> a way to get that information, including confidence intervals and a
> significance test for the correlation coefficient?
>
> Any kind of advice is deeply appreciated!
>
> Best,
>
> Matthias
>
> _______________________________________________
> R-sig-meta-analysis mailing list
> R-sig-meta-analysis at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>
> _______________________________________________
> R-sig-meta-analysis mailing list
> R-sig-meta-analysis at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-meta-analysis
>

	[[alternative HTML version deleted]]