Meta-analysis of estimates based on correlated outcomes in brms

Interfaces brms
1

Hello,
This is related to the issue I raised which has been closed (Multiple response variables in meta-analysis. · Issue #1212 · paul-buerkner/brms · GitHub) which Paul has already kindly reviewed. I suggested a feature to allow users to meta-analyse results which were based on correlated outcome (response) variables.

Paul explained that the feature was not needed as the desired model could be obtained by arranging the data into long format and using fcor. However, I have been unable to do get this to work (I read the response to a previous question (Brms: Multivariate meta-analysis syntax) to help with the model specification on the RHS). Can someone please advise where I am going wrong in trying to specify the known standard errors and known correlation?

Commented reproducible example below.


## Dataset in long format. Two trials, two response variables. Effect estimates and standard errors for each combination (4 effect estimates)
mydf_lng <- structure(list(trial = c(1L, 1L, 2L, 2L), 
                           response = c("1", "2", "1", "2"), 
                           y = c(1.405, -1.652, 0.047, -0.169), 
                           se = c(0.642,0.697, 0.055, 0.114)),
                      row.names = c(NA, -4L), class = "data.frame")

## Correlation matrix indicating correlation between each effect estimate. Effect estimates within trials are correlated. Those between trials are not.
cors <- structure(c(1, -0.921, 0, 0,
                    -0.921, 1, 0, 0,
                    0, 0, 1, -0.48,
                    0, 0, -0.48, 1), .Dim = c(4L, 4L), .Dimnames = list(NULL, NULL))

## Works, specify se but not cors
mod1 <- brm(y|se(se) ~ 0 + response + (0 + response | trial) , data = mydf_lng, data2 = list(cors = cors))

## Works, specify cors but not se
mod2 <- brm(y ~ 0 + response + (0 + response | trial) + fcor(cors), data = mydf_lng, data2 = list(cors = cors))

## Error message "Error: Invalid addition arguments for this model" when try to specify both the se and correlation
mod3 <- brm(y|se(se) ~ 0 + response + (0 + response | trial) + fcor(cors), data = mydf_lng, data2 = list(cors = cors))

Many thanks!

2

You need to pass the covariance matrix to fcor() in this case. That is, the input to fcor needs to contain the SE information already. Then, no se() is needed anymore. If you then want to avoid estimating sigma, you additionally need to set it to for example via prior = prior(constant(1), class = "sigma")

1 Like
3

Thanks for the amazingly quick reply. I will do both.
That actually makes it more straightforward for me to implement as I have the coefficients and variance-covariance matrix from the original (IPD-derived) models.

4

PS presumably if I wanted to nest one outcome within a hierarchy but not others (eg response 1 but not response 2, I would be able to do so by setting-up indicator variables and writing

~ 0 + response1 + response2 + (0 + response1 | trial) + fcor(cors)
5

Yes. Indeed.

1 Like
7

Very enlightening, Paul. Of note, functions that construct variance-covariance matrices (eg, metafor::vcalc()) usually request the know sampling variances , not SE.

When you say

Does it matter if the covariance matrix was constructed based on the variance, not SE?

I guess not based on the Stan code that brms uses for meta-analyses (see below), but would like to confirm.

Thanks!

transformed data {
  vector<lower=0>[N] se2 = square(se);
}
2 Likes
8

Does anyone else know? :)

1 Like
9

The sampling variance is just the SE squared, so there shouldn’t be a difference

EDIT: As in, just pass the same matrix to fcor() that you would to `vcalc°

2 Likes
10

Sorry to engage in thread necromancy, but as a quick follow-up – I am trying to fit a meta-analysis with a known covariance structure between effect sizes, as described above. Based on what I see, I believe the code should be:

m1 = brm(SMD ~ 1 + IV + fcor(V) + (1 | Study_Name/id),
prior = c(
prior(normal(0, 1), class=‘b’),
prior(normal(0, 1), class=‘sd’),
prior(normal(0, 1), class=‘Intercept’),
prior(constant(0), class = “sigma”)
),
data2 = list(V = V),
data=dat, chains=4, cores = 4, iter = 1000)

This incorporates the SE associated with each estimate via fcor(V)… using a covariance matrix V calculated via impute_covariance_matrix from the clubSandwich package… but typically, if you fit a meta-analysis using se(), sigma is set to 0. I try to do so via the final line of the prior specification. However, if I set sigma to 0, the log probability evaluates to log(0) for each chain during initialization so sampling is not done. I haven’t been able to sort it out. Is there something I’m missing? Does the usual se(SE, sigma=FALSE) do something special rather than simply setting sigma to constant(0)? Or perhaps I am being dense.