[R-meta] rma.mv metafor models for genome-wide meta-analyses are (too) conservative

[Viechtbauer, Wolfgang (NP)]

Dear Beate,

Just so I understand this correctly:

1) You have 130 estimates (of something) within each of 23 cohorts, so the input to yi in the model below is of length 130*23=2990, correct?

2) And Vsampoverlap is therefore a 2990x2990 variance-covariance matrix, which is block-diagonal with 23 130*130 blocks, correct?

3) Then you use the results from this model to compute 8 million predicted values (for various combinations of COV1, COV2, and COV3) and compute 'z = pred / se', correct?

4) Then you examine the distribution of these 8 million z-values and find that they are too conservative (i.e., they are not quite standard normal, but have an SD that is somewhat larger than 1 / more than 5% of the z-values are >= +-1.96), correct?

I hope I got all of this correct. I am not familiar with the specifics of this application, but here are some thoughts.

First of all, it is not entirely clear to me why you would expect the predicted values to be standard normal in the first place. Maybe this is similar to GWAS (in a single sample) where one would expect the vast majority of SNPs not to be associated with some phenotype of interest so the p-values should form a roughly uniform distribution, but this may not be true when there is population stratification (or other things) going on and one can use the distribution of p-values to check for this.

You also wrote that the correlation matrix is 'scaled to the SE of each BETA'. So does the diagonal of Vsampoverlap correspond to the standard errors of the estimates or their sampling variances (squared standard errors)? It should be the latter.

The random effects structure of the model is quite simple, assuming a constant correlation for every pair of estimates and a constant amount of heterogeneity for the various estimates. These may be fairly sensible or rather strict assumptions, depending on the application. Relaxing these assumptions is tricky though, because with 130 different estimates, even allowing the amount of heterogeneity to differ would introduce 129 additional parameters into the model and let's not even think about a fully unstructured var-cov matrix for the random effects.

An interesting attempt might be to try:

res <- robust(model_cov, cluster=COHORT)

(or possibly also adding clubSandwich=TRUE) and then using 'res' for obtaining the predicted values, but I am not sure how well cluster-robust inferences methods will handle 130 estimates with 'only' 23 cohorts. The idea is that even if Vsampoverlap and/or the random effects structure is misspecified, cluster-robust inferences methods should give you asymptotically correct standard errors / tests. But I somewhat doubt that this works with so many estimates relative to the number of cohorts. Still worth a try though.

Maybe some of these comments will be helpful or at least spark a bit more discussion.

Best,
Wolfgang

> -----Original Message-----
> From: R-sig-meta-analysis <r-sig-meta-analysis-bounces using r-project.org> On Behalf
> Of St Pourcain, Beate via R-sig-meta-analysis
> Sent: Friday, January 26, 2024 19:14
> To: r-sig-meta-analysis using r-project.org
> Cc: St Pourcain, Beate <Beate.StPourcain using mpi.nl>
> Subject: [R-meta] rma.mv metafor models for genome-wide meta-analyses are (too)
> conservative
>
> Hi,
>
> We are a group of geneticists using meta-regression for genome-wide meta-
> analysis and encountered a hidden thorny issue.
> We use metafor rma.mv models to meta-analyse:
>
> * 130 correlated input statistics
> * each statistics has BETAs for 8 million variants (reflecting genetic
> association effects)
>
> capturing:
> * 70000 individuals
> * 400000 repeat observations
>
> from
> 23 cohorts.
>
> Across the 8 million variants, the following model fits quite well for each
> variant, based on a (fairly) well-known phenotypic correlation matrix for sample
> overlap (Vsampoverlap) scaled to the SE of each BETA.
>
> model_cov <- rma.mv(yi = BETA, V = Vsampoverlap, mods = COV1+ COV2+COV3, random=
> list(~ 1|COHORT/VAR), data = df)}, silent = F)
>
> where VAR specifies the input statistics, COHORT represents cohorts, and COV
> represents fixed effects.
>
> However, when predicting BETAs for each of the 8 million variants across a grid
> of fixed effect predictors (COV1, COV2, COV3) we note that, on average, derived
> genome-wide Z-scores based on predicted BETAs and SEs are too conservative and
> deviate from the expected null distribution in a quantile-quantile plot,
> affecting all predictions. We can also quantify this deviation from a null
> distribution across all 8 million variants using the LDSC intercept (i.e. the
> Linkage disequilibrium Score regression intercept for a variant-based
> heritability estimation; Bullik-Sullivan 2015:
> https://pubmed.ncbi.nlm.nih.gov/25642630/), which should be one if unbiased, but
> we observe ~0.9(SE=0.005). Note that the intercept would be above one, if test
> statistics were inflated.
>
> Thus, we think we subtly over-correct for relatedness and quench the power of
> our analysis. Would there be any thoughts within the group as to how to relax
> the adjustment for relatedness in metafor?
>
> Thanks so much for any input,
> Beate
>
> PS: Note to random geneticists in the group: If we could, we would replace the
> phenotypic correlation matrix with the non-genetic part of the phenotypic
> correlation matrix, but such an estimation is unreliable given the small sample
> size of most cohorts.
>
> R version 4.1.3 (2022-03-10)
> Platform: x86_64-w64-mingw32/x64 (64-bit)
> Running under: Windows 10 x64 (build 19045)
>
> Matrix products: default
>
> locale:
> [1]
> LC_COLLATE=English_Europe.1252  LC_CTYPE=English_Europe.1252    LC_MONETARY=Engl
> ish_Europe.1252 LC_NUMERIC=C
> [5] LC_TIME=English_Europe.1252
>
> attached base packages:
> [1] parallel  stats     graphics  grDevices
> utils     datasets  methods   base
>
> other attached packages:
> [1] data.table_1.14.2   metafor_4.4-0       numDeriv_2016.8-1.1 metadat_1.2-
> 0       Matrix_1.4-0
>
> Beate St Pourcain, PhD
> Senior Investigator & Group Leader
> Room A207
> Max Planck Institute for Psycholinguistics | Wundtlaan 1 | 6525 XD Nijmegen |
> The Netherlands
>
> @bstpourcain
> Tel: tel:+31%2024%20352%201964
> Fax: tel:+31%2024%20352%201213
> ORCID: https://orcid.org/0000-0002-4680-3517
> Web: https://www.mpi.nl/departments/language-and-genetics/projects/population-
> variation-and-human-communication/
> Further affiliations with:
> MRC Integrative Epidemiology Unit | University of Bristol | UK
> Donders Institute for Brain, Cognition and Behaviour | Radboud University | The
> Netherlands